Thermal diffusion of ionic species in charged nanochannels

Abstract

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.

---

1. 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,³ combustion processes,⁴ 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.⁵ Thermal diffusion can also be used to microscopically manipulate the motions of biomolecules⁶ and colloidal particles, which can work as nanocarrier agents for drug delivery.⁷ In low carbon energy conversion, understanding the thermo-diffusive 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–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, 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.¹⁵ 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,¹⁶ and temperature-controlled nano-porous membrane for molecular transport and characterization¹⁷ 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,¹⁸ osmotic energy harvesting based on pressure-retarded osmosis and reversed electrodialysis methods,¹⁹ electric energy storage technologies such as alkaline zinc-based aqueous flow batteries,²⁰ clay minerals forming the nanochannels for geofluids are mostly negatively charged.²¹ 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–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–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 biocompatibility, biodegradability, and unique optical, electronic, and mechanical properties.³⁰ These make quartz and other silica phases vital constituents of micro/nanodevices for water filtration,^31–34 drug delivery,³⁰ biomolecule detection,^35,36etc. 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.

2. 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).

2.1 Thermodynamic formulation and molecular dynamics setup

Based on non-equilibrium thermodynamics,³⁷ the diffusive mass flux of solute (component 1), J₁, in a non-reacting binary mixture under a thermal gradient and in mechanical equilibrium is given by³⁸

J₁ = −ρD₁₂∇w₁ − ρw₁(1 − w₁)D_T∇T, (1)

where D_T and D₁₂ are the thermal and the effective mutual mass (Fickian) diffusion coefficients of the mixtures, respectively, w₁ is the weight fraction of component 1, ρ is the mass density of the mixture, and T 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 x₁ is the molar fraction of component 1 (solute). This equation is an approximation, as it is derived with the assumption that the quantity of solvent significantly exceeds the quantity of the solute. Thus, a positive S_T 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 S_T,Na⁺ = S_T,Cl⁻ = S_T, which is also assumed in previous MD simulations.¹⁵

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 S^∞_T is the asymptotic limit of the Soret coefficient at high temperatures, T₀ is the inversion temperature where S_T switches sign, and the rate τ determines the strength of the temperature effect on S_T. Combining eqn (2) and (3) provides the following relationwhere x₁(T) is the molar fraction of component 1 at temperature T, and k is a fitting parameter, in addition to the previously defined S^∞_T, τ, T₀.

The thermodynamic formulation shows that the computation of the Soret coefficient (S_T) 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.

Fig. 1a shows the MD simulation cell adopted in the work. It contains ∼4800 water molecules, ∼450 NaCl ion pairs, and ∼8900 quartz surface atoms. The system consists of NaCl aqueous solutions confined in a charged quartz slit nanochannel. Periodic boundary conditions are applied at all boundaries of the cell. Thermodynamic ensembles including the microcanonical (NVE), canonical (NVT), and isothermal-isobaric (NPT) ensemble are used in different stages of our modellings. The methodology for calculating the Soret coefficient of confined NaCl aqueous solutions is like the one presented in,^15,41,42 where boundary-driven molecular dynamics simulations are used. The system first undergoes a relaxation stage in NPT ensemble to achieve the desired pressure by dynamically adjusting the cell's dimensions. After that, two regions of the interfacial solutions with widths ∼0.15 nm are thermostatted, i.e., the temperature of the water molecules in these regions is regulated to prescribed higher and lower values relative to the rest of the system. These are the hot region and the cold region in Fig. 1a. The water molecules in the rest of the system, the ions, and the quartz surface atoms are not thermostatted. These particles regulate their temperature via interactions with the thermostatted water molecules by thermal conduction. A constant thermal gradient in the x direction can be created by directly conditioning the water molecules in the hot region at high temperatures (T_H), and thermostatting the ones in the cold region at low temperature (T_C). It has been found that a constant heat flux is established in a few hundred picoseconds. After equilibration, a steady-state condition of zero mass flux is reached, and a constant concentration gradient along the concomitant temperature gradient is established. The obtained composition and temperature profiles are then used to compute the Soret coefficient based on the theory introduced above. For comparison of different profiles, the zero planes of the quartz-water interface are specified as shown in Fig. 1b, where the average position of silicon atoms in the superficial silicon atoms (O–Si–O) and the ones in the silanol groups (Si–OH) has been adopted.

Fig. 1 Schematic of the MD system setup: (a) simulation box showing dimensions of constituents and thermostatted regions; (b) magnification focused on the quartz surface showing the inner and outer silanol groups. The solid black line denotes the zero plane, and the dashed black line separates the surface atoms (flexible region) from bulk atoms (fixed region); (c) locations of the deprotonated outer silanol groups for the four investigated surface charge densities (note that only half of the surface along the x direction is shown); and (d) definition of ϕ, which quantifies the spatial orientation of the water dipole with respect to the quartz surface.

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.

2.2 Computational details

LAMMPS⁴³ was adopted for all non-equilibrium molecular dynamics (NEMD) and equilibrium molecular dynamics (EMD) simulations, and the VMD⁴⁴ package was used for the trajectory visualization. The Velocity-Verlet algorithms were adopted to integrate Newton's equations of motion with a time step of 3.0 fs.

2.2.1 Charged quartz substrates. The charged quartz substrates with different surface charge densities were built by following the previous study of Kroutil et al.⁴⁵ and Quezada et al.⁴⁶ The surface d = (101) of α-quartz (α-SiO₂) 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, 0, 4, 8, and 16 outer silanol groups on each quartz surface were deprotonated by directly deleting the dangling hydrogen atoms to produce surface charge densities of approximately 0.00, −0.03, −0.06, and −0.12 C m⁻², matching along with pH values of around 2.0–4.5, 7.5, 9.5, and 11, respectively.⁴⁵ The partial charge values of the quartz surface atoms are reassigned according to the quantum mechanical calculations of Kroutil et al.⁴⁵ The quartz interactions are described by the ClayFF,⁴⁹ where the same nonbonded and bonded parameters were adopted for both uncharged and charged quartz substrates.

2.2.2 NaCl aqueous solutions. A certain amount of water molecules and sodium and chloride ions were packed into the charged quartz nanochannels to achieve an average salt concentration of 4.0 mol kg⁻¹ (an intermediate concentration where non-ideal effects become important³⁹). This concentration enables direct comparison of the results in this work for nano-confined solutions with previous studies for bulk ones. The formed NaCl aqueous solutions were simulated utilizing the rigid SPC/E model⁵⁰ for water and the force field of Dang et al.^51–54 for the ion-ion/water interactions, whose accuracy in simulations of bulk and interfaces has been demonstrated.⁵⁵ Previous works^15,22,56 have shown that these force–field parameterizations accurately reproduce the experimental thermo-diffusive response of NaCl solutions, including the concentration and temperature effects on the Soret coefficient, particularly the change in sign and minimum.

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 r_c = 1.5 nm was adopted for the short-range interactions, and the long-range electrostatic forces were calculated utilizing the Particle–Particle-Particle-Mesh (PPPM) method.⁵⁷ RATTLE algorithm⁵⁸ was adopted to constrain the bonds involving hydrogen atoms.

2.2.3 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 P = 600 bar, and a system temperature of T = (T_H + T_C)/2, where T_H and T_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⁻¹ nm⁻², the positions of water oxygen atoms in the hot and cold thermostatted regions were restrained in the x direction, but the water molecules were allowed to move freely in the yz 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,⁵⁹ 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 z-direction (fixed in the x- and 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⁻¹ nm⁻² to the fixed atoms in the corresponding directions.

2.2.4 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.

2.2.5 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.

3. Results and discussion

3.1 Structural modifications of nano confined NaCl aqueous solutions

At the high concentration of the solution (4.0 mol kg⁻¹) considered, the coulombic interactions are effectively screened.²² Therefore, the ion-ion/water dispersive interactions become increasingly prominent. The structural properties of confined solutions were derived from the 20 ns of the trajectory at the production stage in the EMD simulations, using the charged quartz surfaces with four different surface charge densities, i.e., 0.00 C m⁻², −0.03 C m⁻², −0.06 C m⁻², and −0.12 C m⁻² (see Fig. 1c) at the temperature of 300/350 K and the pressure of 600 bar, which is consistent with the previous experimental and numerical study²² and facilitates the subsequent comparative analysis. At the same time, the bulk solutions under the same thermodynamic conditions were investigated for comparison.

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 and 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 and 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. 2 Calculated variation of quantities in a nano-channel of width L_z = 4.0 nm at temperature T = 300/350 K and pressure P = 600 bar: (a–f) number density profiles (ρ_N) of water oxygen atoms (O_w), water hydrogen atoms (H_w), sodium ions (Na⁺), chloride ions (Cl⁻), quartz surface oxygen atoms (O_s), and quartz surface hydrogen atoms (H_s) along the z direction for all surface charge densities; (g–i) number fraction of Nacl ions along the z direction; and (j) water orientation on the quartz surface. Vertical solid or dashed lines indicate the positions of two water layers (peaks in the water density profiles).

Fig. 3 (a–d) Charge density profiles (ρ_q) of the whole system, water molecules, sodium ions, and chloride ions along the z direction for all surface charge densities. All the data correspond to the quartz nanochannel of width L_z = 4.0 nm, temperature T = 350 K, and pressure P = 600 bar.

Fig. 2g shows the number fraction of NaCl ions x_NaCl = (N_Na⁺ + N_Cl⁻)/(N_Na⁺ + N_Cl⁻ + N_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. 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 nanoconfinement and surface charge. Previous studies^42,60 have found correlations between the thermomolecular orientation effects and the Soret effect. The results in Fig. 2j demonstrate that the nanoconfinement and surface charges further modify the spatial orientation of the interfacial water molecules, inducing the corresponding variations of the thermo-diffusive responses of the interfacial solutions.

The formations of electrical double layers on the charged quartz surfaces can be evaluated by the charge density distribution profiles (see Fig. 3). The charge separation is observed at the quartz-solution interface in Fig. 3a, which shows that the system reaches electroneutrality at ∼1.2 nm from the quartz surface. This is consistent with the number density profiles in Fig. 2. The fluctuating charge density profiles (see Fig. 3a–d) highlight the effects of nanoconfinement and surface charge on the interfacial liquid structure.

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), g(r) (Fig. 4a–c, i and j), and the coordination numbers, N_c (Fig. 5a). These show that inside the nanochannels, part of the water–water hydrogen bonds (O_w–H_w) are replaced by water–quartz ones (O_w–H_s and O_s–H_w), and that the replacement is enhanced with increasing surface charge, leading to the increase of the total coordination number N_c,total = N_c,Ow–HW + N_c,Os–HW + N_c,Ow–Hs in Fig. 5a. Specifically, when the system temperature (T) is equal to 350 K, in a system with zero surface charge, the total coordination number changes from 1.29 under the bulk conditions to 1.52 under the confined conditions, showing enhanced hydration of water molecules under confinement conditions. With increasing surface charge density, the total coordination number increases from 1.52 to 1.63 for the confined case with a surface density of −0.12 C m⁻², highlighting that the hydration of water molecules is further enhanced by the surface charge. The hydration of water molecules can be enhanced further, albeit not strongly, by decreasing the system temperature from 350 K to 300 K.

Fig. 4 (a–j) Normalized radial distribution functions g(r) for NaCl solutions confined in charged quartz nanochannels of the width L_z = 4.0 nm with different surface charge densities. The data correspond to 300/350 K, and 600 bar and the radial distribution functions g(r) are normalized by g(r = 1.5 nm).

Fig. 5 (a–d) Coordination numbers N_c for NaCl solutions confined in charged quartz nanochannels of the width L_z = 4.0 nm with different surface charge densities. The data correspond to 300/350 K, and 600 bar and the radial distribution functions g(r) are normalized by g(r = 1.5 nm).

Moreover, the (water solvation) hydration shell of the ions Na⁺ and Cl⁻ is modified by temperature, the nanoconfinement conditions, and the surface charges of the quartz substrates. When the system temperature (T) is equal to 350 K, on the one hand, the coordination number of Na⁺–water (the number of water molecules around the cation Na⁺) under the bulk conditions is 4.84 (quantified by Na⁺–O_w and see Fig. 4e and 5b), while under confinement conditions and zero surface charge density, the solvation shell of Na⁺ is slightly loosed, and the coordination number decreases to 4.77. This further decreases to 4.69 when the quartz surface charge density is −0.12 C m⁻². On the other hand, the partly dehydrated Na⁺ ions coordinate with the negatively charged chemical groups on the quartz surfaces, specifically, the siloxane bridges ≡Si–O–Si≡ and the dangling oxygen ≡SiO⁻ (quantified by Na⁺–O_s and see Fig. 4d and 5b). With surface charge density increasing from 0.00 C m⁻² to −0.12 C m⁻², more surface oxygen atoms O_s coordinate to Na⁺ ions, leading to the increase of the total coordinate number in Fig. 5b. The Cl⁻–water coordination number also decreases under the nanoconfinement conditions, from 6.1 in bulk to 5.9 for the cases of zero surface charge (quantified by Cl⁻–H_w and see Fig. 4g and 5c). Cl⁻ ions coordinate with the hydrogen atoms H_s in the silanol groups of the quartz surface instead. Unlike Na⁺ ions, the coordination number of Cl⁻–H_s decreases with increasing surface charge density because of the repulsive interactions between negatively charged surfaces and anions and due to less hydrogen atoms for the deprotonated silanol groups. However, no significant modifications can be observed for the total coordination number, the Cl⁻–H_w, and Cl⁻–H_s coordination shells in 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 and c). In summary, a less tight solvation shell of NaCl ions can be obtained under lower temperature, nanoconfinement and 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 5d). However, the variation of the surface charge appears to have little effect on ion pairing.

3.2 Thermo-diffusive response of confined NaCl aqueous solutions

Section 3.1 showed how the confinement conditions and surface charges modified the liquid structure of the solutions. This is expected to impact the thermo-diffusive response of the confined solutions. To corroborate the expectation, a temperature gradient along the x direction was used in the nanochannel of width L_z ≈ 4.0 nm. The stationary state of the confined solutions was reached in ∼9 ns. At this state where well-defined temperature and concentration profiles were developed. The profiles along the x direction in the stationary state for different confined cases and the reference bulk condition are shown 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.

Fig. 6 Temperature and concentration profiles along the x direction of the NaCl aqueous solutions inside the charged quartz nanochannels with different surface charge densities (a–d) and under the bulk conditions (e–h).

By eliminating the marginal regions where the impact of the thermostats is perceived (see the dashed thermostatted regions in Fig. 6), the temperature dependences of the solute molar fraction in the NEMD simulations are obtained from Fig. 6. These are presented in Fig. 7a–e by scattered symbols and fitted using eqn (4). The fitting curves represent well the simulation results, indicated by the high R-square values for all cases. The obtained fitting parameters, i.e., S^∞_T, T₀, and τ are further used in eqn (3) to calculate the Soret coefficients of the confined NaCl aqueous solutions at different temperatures. These are presented in Fig. 7f. Romer et al.²² have reported experimental and simulated Soret coefficients of bulk NaCl aqueous solutions at similar thermodynamic states and concentrations, which are also given in Fig. 7f. The Soret coefficients calculated here agree with the simulated values in ref. 22, which verifies the accuracy of the results in this work. The calculated Soret coefficients show the same temperature dependence as the ones obtained by experiments but with lower magnitudes. Currently, the existing MD simulation models for NaCl aqueous solutions can only reproduce Soret coefficients of the right order of magnitude and predict the correct dependence of the Soret coefficient with temperature and concentration, possibly due to the contradiction between the complexity of this non-isothermal transport phenomenon and the comparative simplicity of the interatomic interaction model. This discrepancy between MD and experimental results within ∼10⁻³ K has been thus deemed as acceptable in previous studies.^15,22 Indeed, the Soret coefficients measured by different experimental approaches under the same thermodynamic conditions also disagree with each other (e.g., the infrared thermal diffusion forced Rayleigh scattering experiment²² and thermo-gravitational column experiments⁶¹). Compared with the results in the bulk solutions, the measured Soret coefficient of the confined solutions in charged nanochannel is generally lower, indicating that the solution becomes more thermophilic. It is found that after the equilibrium stage, the steady states of the thermal gradient and concentration gradient have been reached, with respective characteristic times of 1 and 9 ns. The position of the maximum in the molar fraction profile (see Fig. 7a–e) indicates the sign-reversal temperature, T₀, where the Soret coefficient is zero. It is found that this sign-reversal temperature increases with the surface charge density. All cases studied here show a temperature inversion effect. The sign of the Soret coefficient changes at the temperature of T₀. At T₀ the thermo-diffusive response of the solution shifts from thermophilic (the salt accumulates in the cold region) at low temperature to thermophobic (the salt accumulates in the hot region) at high temperature. In addition, it is also found in Fig. 7f that the solution tends to be more thermophilic with increasing negative surface charges. The previous work by Di Lecce et al.⁴¹ has shown that by varying the nanopore size, the thermal diffusion of alkali halide solutions in silica nanopores can be effectively regulated. The results of the present work show, for the first time, that changing the surface charge density can also achieve the same effect.

Fig. 7 (a–e) Temperature dependence of the molar fraction of NaCl ions at the stationary state in the bulk and confined in the charge quartz nanochannels, where full lines are using eqn (4). (f) Temperature dependence of the Soret coefficient for the NaCl solutions in the bulk and confined in the charge quartz nanochannels. Our MD simulations used a thermal gradient of ∇T = 34 K nm⁻¹ and the average temperature, pressure, and salt molality were set to 350 K, 600 bar, and 4.0 mol kg⁻¹, respectively. (g) The adopted criterion to define the hydrogen bonds between water molecules and silanol groups, as well as water molecules and water molecules. (h) Average hydrogen bond number per water molecule in NaCl solutions confined by charged quartz nanochannels and in bulk.

The role of hydrogen bond structure in thermal diffusion has been previously discussed by Niether and Wiegand.⁶² 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.³⁹ 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 7f and the associated discussions.

Fig. 8 shows the molar fraction distributions of the NaCl along the z direction in the hot and cold side (see Fig. 1a) for four kinds of charged quartz slit pores. It is found that the thermo-diffusive responses of the NaCl aqueous solutions are different in the bulk liquid and in the boundary layers, characterized by the differences between the hot and the cold sides (solid grey lines). The difference is constant in the bulk liquid and fluctuating in the boundary layers. Secondly, the positive value of the difference indicated that the NaCl is driven from the hot toward the cold side (thermophobic, S_T > 0) while the negative one means that the NaCl is migrating from the cold toward the hot side (thermophilic, S_T < 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.

Fig. 8 (a–d) Evolution of the molar fraction distribution of the NaCl along the z direction in the hot and cold side for four kinds of surface charges.

Using the atomic velocities from EMD simulations, the one-dimensional effective mutual diffusion (interdiffusion) coefficient in the x direction, i.e., D₁₂ in eqn (2), was quantified by the time integral of the autocorrelation function of the concentration current,^63,64

where N₁ and N₂ are the particle numbers of component 1 (NaCl) and component 2 (water), x₁ = N₁/(N₁ + N₂) and x₂ = N₂/(N₁ + N₂) are their molar fractions,is the relative velocity of concentration current in the x direction, and v_i,x(t) and v_j,x(t) 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 D₁₂ enables evaluation of the thermal diffusivity of NaCl ions using eqn (2), i.e., D_TT = S_TD₁₂T, where, S_T is presented in Fig. 7f. The calculated D_TT 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⁻⁹ m² 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 thermally-induced diffusive transport is affected by nano-confinement and surface charges, which facilitate the structural modifications of EDL.

Fig. 9 The (a) mutual diffusion coefficients (D₁₂) and (b) thermal diffusivity (D_TT) of NaCl ions under different surface charge conditions and at temperatures of 300 K and 350 K.

4. 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⁻² to −0.12 C m⁻². 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 (S_T < 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.

Conflicts of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

Chen acknowledges the President Doctoral Scholarship Award (PDS Award 2019) from The University of Manchester, UK. Jivkov acknowledges gratefully the financial support from the Engineering and Physical Sciences Research Council (EPSRC), UK, via Grant EP/N026136/1. The authors acknowledge the assistance provided by the Research IT team for the use of the Computational Shared Facility at The University of Manchester.

References

1. G. H. Thompson, M. N. Myers and J. C. Giddings, Sep. Sci., 1967, 2, 797–800.
2. J. C. Giddings, V. Kumar, P. S. Williams and M. N. Myers, Polymer separation by thermal field-flow fractionation: high-speed power programming, 1990, pp. 3–21.
3. M. F. Riley and A. Firoozabadi, AIChE J., 1998, 44, 452–464.
4. D. Rosner, R. Israel and B. La Mantia, Combust. Flame, 2000, 123, 547–560.
5. P. Baaske, C. J. Wienken, P. Reineck, S. Duhr and D. Braun, Angew. Chem., Int. Ed., 2010, 49, 2238–2241.
6. S. Duhr and D. Braun, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 19678–19682.
7. A. L. Sehnem, A. M. F. Neto, D. Niether and S. Wiegand, Phys. Rev. E, 2018, 98, 062615.
8. S. Di Lecce and F. Bresme, J. Phys. Chem. B, 2018, 122, 1662–1668.
9. R. F. Stout and A. S. Khair, Phys. Rev. E, 2017, 96, 022604.
10. Y. Liu, M. Cui, W. Ling, L. Cheng, H. Lei, W. Li and Y. Huang, Energy Environ. Sci., 2022, 15, 3670.
11. S. Zhou, Z. Qiu, M. Strømme and C. Xu, Energy Environ. Sci., 2021, 14, 900–905.
12. J. Bai, Y. Hu, T. Guang, K. Zhu, H. Wang, H. Cheng, F. Liu and L. Qu, Energy Environ. Sci., 2022, 15, 3086.
13. H. Yan, M. Sedighi, H. Ding, Z. Sun and H. Xie, J. Hydrol., 2021, 603, 126848.
14. M. Sedighi, Investigation of hydro-geochemical processes in coupled thermal, hydraulic, chemical and mechanical behaviour of unsaturated soils, Cardiff University, United Kingdom, 2010.
15. S. Di Lecce, T. Albrecht and F. Bresme, Sci. Rep., 2017, 7, 1–10.
16. C. B. Mast and D. Braun, Phys. Rev. Lett., 2010, 104, 188102.
17. M. Belkin, S.-H. Chao, G. Giannetti and A. Aksimentiev, J. Comput. Electron., 2014, 13, 826–838.
18. X. Zhu, Y. Zhou, J. Hao, B. Bao, X. Bian, X. Jiang, J. Pang, H. Zhang, Z. Jiang and L. Jiang, ACS Nano, 2017, 11, 10816–10824.
19. T. Ma, E. Balanzat, J.-M. Janot and S. Balme, ACS Appl. Mater. Interfaces, 2019, 11, 12578–12585.
20. Z. Yuan, X. Liu, W. Xu, Y. Duan, H. Zhang and X. Li, Nat. Commun., 2018, 9, 1–11.
21. Y. Yang, R. Qiao, Y. Wang and S. Sun, Appl. Clay Sci., 2021, 201, 105924.
22. F. Römer, Z. Wang, S. Wiegand and F. Bresme, J. Phys. Chem. B, 2013, 117, 8209–8222.
23. B. Arlt, S. Datta, T. Sottmann and S. Wiegand, J. Phys. Chem. B, 2010, 114, 2118–2123.
24. H. Ning, S. Datta, T. Sottmann and S. Wiegand, J. Phys. Chem. B, 2008, 112, 10927–10934.
25. Y.-C. Zou, L. Mogg, N. Clark, C. Bacaksiz, S. Milovanovic, V. Sreepal, G.-P. Hao, Y.-C. Wang, D. G. Hopkinson and R. Gorbachev, Nat. Mater., 2021, 1–6.
26. N. Di Trani, A. Pimpinelli and A. Grattoni, ACS Appl. Mater. Interfaces, 2020, 12, 12246–12255.
27. G. Bruno, N. Di Trani, R. L. Hood, E. Zabre, C. S. Filgueira, G. Canavese, P. Jain, Z. Smith, D. Demarchi, S. Hosali, A. Pimpinelli, M. Ferrari and A. Grattoni, Nat. Commun., 2018, 9, 1682.
28. N. F. Durand, C. Dellagiacoma, R. Goetschmann, A. Bertsch, I. Märki, T. Lasser and P. Renaud, Anal. Chem., 2009, 81, 5407–5412.
29. M. Taghipoor, A. Bertsch and P. Renaud, Nanoscale, 2015, 7, 18799–18804.
30. F. Peng, Y. Su, Y. Zhong, C. Fan, S.-T. Lee and Y. He, Acc. Chem. Res., 2014, 47, 612–623.
31. K. Leung, S. B. Rempe and C. D. Lorenz, Phys. Rev. Lett., 2006, 96, 095504.
32. Z. Wang, M. Sedighi and A. Lea-Langton, Water Res., 2020, 184, 116165.
33. Z. Wang, M. Sedighi, A. R. Lea-Langton and M. Babaei, J. Hydrol., 2022, 612, 128220.
34. A. Patti, A. D. Mackie, V. Zelenak and F. R. Siperstein, J. Mater. Chem., 2009, 19, 724–732.
35. C.-C. Chien, S. Shekar, D. J. Niedzwiecki, K. L. Shepard and M. Drndić, ACS Nano, 2019, 13, 10545–10554.
36. P. Karau and V. Tabard-Cossa, ACS Sens., 2018, 3, 1308–1315.
37. S. R. de Groot and P. Mazur, Non-equilibrium thermodynamics, Dover, New York, 1984.
38. J. Colombani, G. Galliéro, B. Duguay, J.-P. Caltagirone, F. Montel and P. A. Bopp, Phys. Chem. Chem. Phys., 2002, 4, 313–321.
39. S. Di Lecce, T. Albrecht and F. Bresme, Phys. Chem. Chem. Phys., 2017, 19, 9575–9583.
40. S. Iacopini, R. Rusconi and R. Piazza, Eur. Phys. J. E, 2006, 19, 59–67.
41. S. Di Lecce, T. Albrecht and F. Bresme, Nanoscale, 2020, 12, 23626–23635.
42. F. Römer, F. Bresme, J. Muscatello, D. Bedeaux and J. M. Rubí, Phys. Rev. Lett., 2012, 108, 105901.
43. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
44. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics Modell., 1996, 14, 33–38.
45. O. Kroutil, Z. Chval, A. Skelton and M. Predota, J. Phys. Chem. C, 2015, 119, 9274–9286.
46. G. R. Quezada, R. E. Rozas and P. G. Toledo, J. Phys. Chem. C, 2017, 121, 25271–25282.
47. S. J. Chen, W. Q. Chen, Y. B. Ouyang, T. Matthai and L. H. Zhang, Nanoscale, 2019, 11, 22954–22963.
48. W. Q. Chen, M. Sedighi and A. P. Jivkov, Nanoscale, 2021, 13, 1696–1716.
49. R. T. Cygan, J. J. Liang and A. G. Kalinichev, J. Phys. Chem. B, 2004, 108, 1255–1266.
50. H. J. C. Berendsen, J. R. Grigera and T. P. Straatsma, J. Phys. Chem., 1987, 91, 6269–6271.
51. L. X. Dang, J. Chem. Phys., 1992, 96, 6970–6977.
52. L. X. Dang and B. C. Garrett, J. Chem. Phys., 1993, 99, 2972–2977.
53. D. E. Smith and L. X. Dang, J. Chem. Phys., 1994, 100, 3757–3766.
54. L. X. Dang, J. Am. Chem. Soc., 1995, 117, 6954–6960.
55. A. Wynveen and F. Bresme, J. Chem. Phys., 2010, 133, 144706.
56. S. Di Lecce and F. Bresme, Mol. Simul., 2019, 45, 351–357.
57. R. W. Hockney and J. W. Eastwood, Computer simulation using particles, 2nd ed., IOP, Bristol, 1988.
58. H. C. Andersen, J. Comput. Phys., 1983, 52, 24–34.
59. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101.
60. F. Bresme, A. Lervik, D. Bedeaux and S. Kjelstrup, Phys. Rev. Lett., 2008, 101, 020602.
61. F. Gaeta, G. Perna, G. Scala and F. Bellucci, J. Phys. Chem., 1982, 86, 2967–2974.
62. D. Niether and S. Wiegand, J. Phys.: Condens. Matter, 2019, 31, 503003.
63. J. Hansen, F. Joly and I. McDonald, Phys. A, 1985, 132, 472–488.
64. A. P. Lyubartsev and A. Laaksonen, J. Phys. Chem., 1996, 100, 16410–16418.
65. S. Pennathur, F. Baldessari, J. G. Santiago, M. G. Kattah, J. B. Steinman and P. J. Utz, Anal. Chem., 2007, 79, 8316–8322.
66. M. A. Alkhadra, X. Su, M. E. Suss, H. Tian, E. N. Guyes, A. N. Shocron, K. M. Conforti, J. P. de Souza, N. Kim, M. Tedesco, K. Khoiruddin, I. G. Wenten, J. G. Santiago, T. A. Hatton and M. Z. Bazant, Chem. Rev., 2022, 122, 13547–13635.

---