Thermal diffusion of ionic species in charged nanochannels

Diffusion of ions due to temperature gradients (known as thermal diffusion) in charged nanochannels is of interest in several engineering fields, including energy recovery and environmental protection. This paper presents a fundamental investigation of the thermal diffusion of sodium chloride in charged silica nanochannels performed by molecular dynamics (MD). The results reveal the effects of nanoconfinement and surface charges on the sign and magnitude of the Soret coefficient. It is shown that the sign and magnitude of the Soret coefficient are controlled by the structural modifications of the interfacial solutions. These modifications include the ionic solvation and hydrogen bond structure induced by the nanoconfinement and surface charges. The results show that both nanoconfinement and surface charges can make the solutions more thermophilic. Furthermore, the thermal diffusion of solutions in boundary layers is significantly different from that of solutions in bulk fluid, contributing to the overall difference between the thermal diffusivity of pore fluid and that associated with bulk fluid. The findings provide further understanding of thermal diffusion in nano-porous systems. The proposed MD simulation methodology is applicable to a wider category of coupled heat and mass transfer problems in nanoscale spaces.


Introduction
Thermal diffusion in fluid mixtures and thermophoresis in suspensions, which are also referred to as the Ludwig-Soret effect, describe the migration of solutes due to temperature gradients. Thermal diffusion participates in various natural processes and is of interest for a range of engineering applications such as thermal field-flow fractionation of colloids and synthetic polymers 1, 2 , spatial composition analysis of hydrocarbon reservoirs 3 , combustion processes 4 , and behaviour of solutions encompassing ions and polar solutes. It has been shown that the thermal diffusion of biomolecules is highly sensitive to their structural modifications in aqueous solutions. This has been used to monitor the protein-ligand binding reactions under temperature gradients 5 . Thermal diffusion can also be used to microscopically manipulate the motions of biomolecules 6 and colloidal particles, which can work as nanocarrier agents for drug delivery 7 . In low carbon energy conversion, understanding the thermodiffusive response of electrolyte solutions and other charged fluids is particularly useful in quantifying the thermoelectric Seebeck effect. This is being extensively studied to develop devices capable of recovering waste heat or solar energy and converting it to electricity [8][9][10][11][12] . Further, evaluating the thermal diffusion of geofluids is vital for the long-term safety of many geotechnical infrastructures such as nuclear waste repositories, landfill containers, and geothermal energy facilities 13,14 .
Despite significant developments in the understanding of thermal diffusion and the successful utilization of this physical phenomenon in aqueous systems in many fields, a full microscopic explanation of the phenomenon is missing. As a result, there is no general theory that can successfully evaluate the magnitude and direction of diffusive flux due to thermal diffusion in aqueous solutions in a wide range of situations. A number of experimental observations have revealed peculiar phenomena, 3 such as a change of the sign of the Soret coefficient at a specific temperature (a negative coefficient means that migration occurs from colder to hotter regions) and the existence of a minimum Soret coefficient 15 . Application of different membranes, microfluidic and nanofluidic devices by utilizing thermal diffusion of aqueous solutions, e.g., thermophoretic molecule trap for DNA replication and accumulation 16 , and temperature-controlled nano-porous membrane for molecular transport and characterization 17 have emerged. Such technological developments highlight the need for further understanding of thermal diffusion in microscale and nanoscale spaces and channels. It has been suggested in the literature that in addition to the nano-confinement effect on thermal diffusion, the effect of charged solid surfaces should also be considered, as it exists in various industrial applications and natural systems, e.g., negatively/positively charged nano-porous membranes for ion transport and ionic selectivity 18 , osmotic energy harvesting based on pressure-retarded osmosis and reversed electrodialysis methods 19 , electric energy storage technologies such as alkaline zinc-based aqueous flow batteries 20 , clay minerals forming the nanochannels for geofluids are mostly negatively charged 21 . While the understanding of the thermo-diffusive response of bulk aqueous solutions has been improved recently by experimental, theoretical, and numerical studies, e.g., 6,13,14,[22][23][24] , the research on thermal diffusion in nano-confined aqueous solutions has been very limited, particularly in charged media. It is known that nanoconfinement and surface charges will significantly change the transport properties of the interfacial liquid, e.g., diffusive transport [25][26][27][28][29] . Therefore, the existing knowledge on the thermo-diffusive properties of bulk liquid cannot be directly applied.
The aim of this work is to advance the mechanistic understanding of thermal diffusion in charged nanoscale structures. To the best of our knowledge, the combined effects of nanoconfinement and surface charges on thermal diffusion have not been studied from a molecular-level perspective before, whilst such understanding is critical to facilitate the technological and scientific advances in the areas mentioned above. To this end, the thermal diffusion of NaCl aqueous solutions confined in charged quartz nanochannels is investigated by equilibrium molecular dynamics (EMD) and non-equilibrium molecular dynamics (NEMD) simulations. The choice of the solid phase is dictated by practical considerations. Silicon-based nanomaterials are appropriate for device applications due to their superb 4 biocompatibility, biodegradability, and unique optical, electronic, and mechanical properties 30 . These make quartz and other silica phases vital constituents of micro/nanodevices for water filtration [31][32][33][34] , drug delivery 30 , biomolecule detection 35,36 , etc. The paper is organized as follows. The theory and simulation approaches are discussed in Section 2. The results of the simulations are presented in Section 3. They show how the thermo-diffusive response of the solutions is affected by the solution-quartz interactions and the nanoconfinement. The findings are explained by the alteration of the interfacial liquid structure, which is also discussed in Section 3. Finally, conclusions are drawn in Section 4.

Theory and method
Thermal diffusion and the associated MD systems to measure the Soret effect of sodium chloride (NaCl) in water are presented in this section. Aqueous NaCl solutions are binary mixtures consisting of salt (component 1, solute) and water (component 2, solvent).

Thermodynamic formulation and molecular dynamics setup
Based on non-equilibrium thermodynamics 37 , the diffusive mass flux of solute (component 1), 1 , in a non-reacting binary mixture under a thermal gradient and in mechanical equilibrium is given by 38 where and 12 are the thermal and the effective mutual mass (Fickian) diffusion coefficients of the mixtures, respectively, 1 is the weight fraction of component 1, is the mass density of the mixture, and is the temperature. Under a given temperature gradient, the system tends towards an equilibrium state of zero mass flux. The Soret coefficient, generally used to quantify thermal diffusion, characterises the equilibrium state. In the case of dilute solutions, it is evaluated by: where 1 is the molar fraction of component 1 (solute). This equation is an approximation, as it is 5 derived with the assumption that the quantity of solvent significantly exceeds the quantity of the solute.
Thus, a positive implies that the solute tends to move with respect to the solvent from the hotter to the colder region. Notably, the existing experimental approaches do not explicitly distinguish the Soret coefficient for the cations and anions and provide Soret coefficients of individual ions. In this work, it is assumed that , + = , − = , which is also assumed in previous MD simulations 15 .
Previous studies have found that the following empirical relation describes the temperature dependence of the Soret coefficient of alkali halide solutions 15, 22, 39 and aqueous colloidal suspensions where ∞ is the asymptotic limit of the Soret coefficient at high temperatures, 0 is the inversion temperature where switches sign, and the rate determines the strength of the temperature effect on . Combining Eq. (2) and Eq. (3) provides the following relation where 1 ( ) is the molar fraction of component 1 at temperature , and is a fitting parameter, in addition to the previously defined ∞ , , 0 .
The thermodynamic formulation shows that the computation of the Soret coefficient ( ) can be achieved by imposing a known temperature gradient and measuring the induced concentration gradient. This is used in the MD simulation approaches described below and employed for the analysis of thermal diffusion. It is noted that the use of solid mineral surfaces and alkali halide solutions in the MD setup of this study mimics the environments of pores in hydrothermal vents/rocks.

Computational details
LAMMPS 43 was adopted for all non-equilibrium molecular dynamics (NEMD) and equilibrium molecular dynamics (EMD) simulations, and the VMD 44 package was used for the trajectory 8 visualization. The Velocity-Verlet algorithms were adopted to integrate Newton's equations of motion with a time step of 3.0 fs.
Charged quartz substrates. The charged quartz substrates with different surface charge densities were built by following the previous study of Kroutil et al. 45 and Quezada et al. 46 . The surface = (101) of α-quartz (α-SiO2) crystal is cut as the quartz-solution interface, and the under-coordinated superficial O and Si atoms are capped with H or OH groups to form outer and inner silanol groups (SiOH) (see Fig. 1b) as in previous studies 47,48 . As shown in Fig. 1c The solution-quartz interactions were depicted utilizing the Coulombic and Lennard-Jones potentials, and the cross-interaction parameters between various atom types were obtained by adopting Lorentz-Berthelot rules. The cut-off distance of = 1.5nm was adopted for the short-range interactions, and 9 the long-range electrostatic forces were calculated utilizing the Particle-Particle-Particle-Mesh (PPPM) method 57 . RATTLE algorithm 58 was adopted to constrain the bonds involving hydrogen atoms.
NEMD Simulation. For NEMD simulation, first, systems were energetically minimized by adopting the steepest descent integrator. Second, a pre-equilibrium stage containing 0.1 ns NVT ensemble and 1 ns NPT ensemble was completed to reach a system pressure of = 600 bar, and a system temperature of = ( H + C )/2, where H and C were the thermostat temperatures in the subsequent NEMD simulations (see Fig. 1a). The system after this stage is referred to as the "pre-equilibrium system".
Third, by switching on the two thermostats and creating a thermal gradient inside the nanochannel, the system was driven out of equilibrium, where water and ions migrate and gradually reach the steady state after several nanoseconds (9 ns in the case studies here). Finally, the production stage lasting 18 ns was run, and the configurations were collected every 100 steps (~0.3 ps) to extract the temperature, density, concentration profiles, and thermodynamic averages by separating the simulation cell into 120 statistical bins in the direction of the temperature gradient, x. The temperature distribution of the solution was computed based on the equipartition principle by collecting the velocities of the ions and the water molecules. Note that in the NEMD simulations, the NVE ensemble is applied to the system, and the temperatures of the water molecules within the two thermostatting layers were further rescaled by the local thermostats every timestep.
By applying a harmonic potential with a spring constant equal to 1000 kJ mol −1 nm −2 , the positions of water oxygen atoms in the hot and cold thermostatted regions were restrained in the direction, but the water molecules were allowed to move freely in the plane. The restrained water molecules could rotate freely and exchange momentum with the unrestrained water molecules and ions. The temperatures of the hot and the cold regions were rescaled to the specific values every timestep by a canonical sampling thermostat implementing global velocity rescaling with Hamiltonian dynamics 59 , and the linear momentum of the system was reset after velocity rescaling. The streaming velocity of the system was monitored, and a statistically averaged zero value was observed.
The silanol atoms, as well as the surface silicon and oxygen atoms (see Fig. 1b), were kept fully flexible. During the NPT simulations, the bulk quartz atoms were only permitted to move in the zdirection (fixed in the xand y-directions) to equilibrate the system pressure and maintain the surface geometry, while in NVT and NVE simulations, these atoms were fixed in all directions (see Fig. 1b).
These fixations were also implemented by utilizing a harmonic potential with spring constant of 1000 kJ mol −1 nm −2 to the fixed atoms in the corresponding directions.

EMD Simulation.
To obtain the structural modifications of confined NaCl aqueous solutions, additional EMD simulations were performed on these systems without the thermal gradient applied.
The production runs in the NVT ensemble lasting 20 ns were performed on the pre-equilibrium system as mentioned above to analyse the equilibration phase of the systems.

Bulk Solutions.
To better demonstrate the effect of nanoconfinement and the surface charge on the thermo-diffusive response of the NaCl aqueous solutions, the same NEMD and EMD simulations were subsequently performed on identical NaCl aqueous solutions in bulk conditions.

Structural modifications of nano confined NaCl aqueous solutions
At the high concentration of the solution (4.0 mol kg-1) considered, the Coulombic interactions are effectively screened 22 . Therefore, the ion-ion/water dispersive interactions become increasingly Water molecules interact with the charged quartz surface establishing hydrogen bonds. The comparisons between the number density profiles of Fig. 2a-d and those of Fig. 2e-f show a considerable degree of overlap between solutions and the quartz surfaces, indicating penetrations of both water molecules and ions into the quartz surface. A noticeable layer-by-layer arrangement of water molecules next to the surface can be observed from the peaks and dips in Fig. 2a-b. This result, together with the charge density profiles of the interfacial water molecules presented in Fig. 3b, shows that the water molecules reorient as they meet the quartz surface, with the hydrogen atoms positioned closer to the surface than the oxygen atoms. In the middle of the nanochannel, the number density of O w , H w , Na + , and Cl − converges to their bulk value at 300/350 K and 600 bars in about 1.2 nm from the quartz surface (see Fig. 2a-d).

Fig. 2g
shows the number fraction of NaCl ions NaCl = ( Na + + Cl − )/( Na + + Cl − + water ) at the charged quartz surface is significantly different from that in the interior of the nanochannel, underlining the impact of the quartz-ion interactions. This leads to the generation of an electrical double layer (EDL) in which the Na + ions absorb on the quartz surface first, resulting in a number density peak as shown in Fig. 2c and a number fraction peak as shown in Fig. 2h, while the Cl − ions absorb immediately above Na + (see Fig. 2d and i) with a concurrent depletion of Na + ions (see Fig. 2c and h). In addition, with increasing surface charge density, more cations are absorbed on the quartz surface to compensate for the negative surface charges (see Fig. 2c and h), while the anions are less affected by the surface charge density (see Fig. 2d and i). In short, the results highlight a stronger affinity of cations towards the negatively charged quartz surfaces.
The effects of nanoconfinement and surface charge on the spatial orientation of the interfacial water are quantified. The spatial orientation is characterized by cos ( ), in which is the angle between the normal to the quartz surface and a vector opposite to a water dipole (see Fig. 1d): cos( ) = −1 represents a water dipole pointing away from the quartz; cos( ) = 1 represents a water dipole pointing towards the quartz; −1 < cos( ) < 1 represents partial orientation of the water dipole. Fig. 2j shows the average of cos ( ) for the interfacial water molecules at different distances from the quartz surfaces for all simulated cases. The results in Fig. 2j show that without electrolytes, the quartz surface induces the spatial orientation of water dipoles at the quartz-solution interface, resulting in two water layers.

12
The spatial orientation of water dipoles in the first layer, near the surface, is mostly pointing away from the surface (cos( ) > 0); in the second layer, the orientation is mostly pointing toward the surface (cos( ) < 0). Apart from the effect of nanoconfinement, Fig. 2j shows that the surface charge can further enhance this spatial orientation. Generally, the magnitude of cos( ) increases with the surface charge increasing. According to these results, the spatial orientation of interfacial water molecules, which is vital for the mobility of particles and the effective anchoring of ions, is greatly affected by the Temperature, nanoconfinement, and surface charge modify the first solvation shell of the interfacial water molecules in the NaCl solutions. These structural modifications can be quantified by the radial distribution functions (RDF), ( ) (Fig. 4a-c and i-j), and the coordination numbers, (Fig. 5a).  Fig. 5c when surface charge density changes. In addition, an enhanced hydration shell of the ions Na + and Cl − can be observed when the system temperature is decreased from 350 K to 300 K (see Fig. 5b-c). In summary, a less tight solvation shell of NaCl ions can be obtained under lower temperature, nanoconfinement and 17 surface charge conditions. This can also be seen from the different heights of the main peaks in the RDF profiles of Fig. 4. The position of these peaks is not significantly affected by the confinement and the surface charges. Observed also is an enhancement of ion pairing in the solution under nanoconfinement conditions or under increased temperature (see Fig. 4h and Fig. 5d). However, the variation of the surface charge appears to have little effect on ion pairing.   in Fig. 6. The concentration profiles show that the distribution of the cations is strongly affected by the surface charge density, featuring a gathering around the deprotonated silanol groups (see Fig. 6c), while the one of the anions is basically independent of the surface charge density (see Fig. 6d). This shows the screening effect of the cations on the charged quartz surfaces. that changing the surface charge density can also achieve the same effect.

Thermo-diffusive response of confined NaCl aqueous solutions
The role of hydrogen bond structure in thermal diffusion has been previously discussed by Niether and Wiegand 62 . Their conclusions show that water-water hydrogen bonds are easier to form at low temperatures, and the solution is more thermophilic, while the water-water hydrogen bonds are destroyed at high temperatures, and the solution is more thermophobic. These conclusions agree well with the results for bulk NaCl solutions (see Fig. 7f and h), namely, under high-temperature conditions, the quantity of hydrogen bonds is lower in the bulk solution, which therefore is more thermophobic.
For solutions under nano-confinement and surface charge, the total quantity of hydrogen bonds is higher than that in bulk, which agrees well with the observation in Fig. 7f that under nanoconfinement and surface charge, the solution exhibits a more intense thermophilic behaviour, with the ions tend to accumulate in hot areas. It has been shown previously that the hydration structure and entropy of the ions have a considerable influence on the magnitude of the Soret coefficient 39 . Their conclusion is that a less tight water solvation shell of NaCl ions makes the solution more thermophilic and results in a more negative Soret coefficient.
This conclusion is in full agreement with the results in Fig. 5 and Fig. 7f and the associated discussions. while the negative one means that the NaCl is migrating from the cold toward the hot side (thermophilic, < 0). Therefore, it is found that with the surface charge increasing, the thermo-diffusive behaviour of some regions in the boundary layers is thermophilic instead of thermophobic, as shown in Fig. 8d.
where 1 and 2 are the particle numbers of component 1 (NaCl) and component 2 (water), 1 = 12, ( ) = 2 ∑ , ( ) is the relative velocity of concentration current in the direction, and , ( ) and , ( ) are the particle velocities. The calculated results are presented in Fig. 9a. They show that both nanoconfinement conditions and surface charges can lower the mutual diffusivity. Previous studies 26-29, 65, 66 have attributed this phenomenon to structural modifications, e.g., the extension of the size of electrical double layers (EDL) ion distribution. In addition, the calculated 12 enables evaluation of the thermal diffusivity of NaCl ions using Eq. (2), i.e., = 12 , where, is presented in Fig. 7f. The calculated is presented in Fig. 9b. It shows that the thermal diffusivity of NaCl ions is slightly lower than that of concentration-driven ionic diffusivity, i.e., 10 −10~1 0 −9 m 2 /s. This suggests that thermal diffusion should not be neglected when considering coupled nanoscale transport under temperature gradients. Moreover, it is found that, like ionic diffusivity, the thermal diffusivity of NaCl ions is lowered by nano-confinements and surface charges. In summary, based on the existing studies on concentration-driven diffusive transport in nanofluidic systems, this work shows that thermallyinduced diffusive transport is affected by nano-confinement and surface charges, which facilitate the structural modifications of EDL.

Conclusions
Molecular dynamics (MD) simulations were performed to investigate the thermal diffusion of NaCl solutions in quartz nanochannels with surface charge densities ranging from 0.00 C • m −2 to −0.12 C • m −2 . The results obtained with the non-equilibrium molecular dynamics (NEMD) analyses showed that the nano-confinement and the surface charges considerably reduced the thermo-phobic response of 4.0 mol/kg NaCl solutions. It was found that by increasing the surface charge density, the ionic solution became thermophilic, reflected in a negative Soret coefficient ( < 0 ) , i.e., the salt moved from cold to hot regions. The equilibrium molecular dynamics (EMD) simulations revealed that the structural modifications induced by the nano-confinement and the surface charge were strongly related to the observed thermo-diffusive response of NaCl solution in the NEMD simulations, and the thermal diffusion behaviour was different in the bulk liquid and the boundary layers.
The study revealed for the first time that the thermal diffusion of alkali solutions in nanopores can be controlled not only by varying the nanopore size, but also by changing the surface charge density. This outcome is important for the design of many devices/facilities, such as nanofluidic separation devices that utilise thermal fields to separate and enrich salts, nano-porous membranes for water desalination and low-grade waste heat recovery, engineering barriers for hazardous contaminants and high-level nuclear waste etc. The modelling methods presented in this paper, together with the insights derived from our results, show a path for further investigation of coupled transport phenomena in nanoscale structures.

Author contributions
W. Q. C., M. S., and A. J. conceived the idea. W. Q. C. performed the numerical simulations. All authors discussed the results and wrote and reviewed the manuscript. M. S. and A. J. supervised the entire project.