Thermo-Osmosis in Charged Nanochannels: Effects of Surface Charge and Ionic Strength

Thermo-osmosis refers to fluid migration due to the temperature gradient. The mechanistic understanding of thermo-osmosis in charged nano-porous media is still incomplete, while it is important for several environmental and energy applications, such as low-grade waste heat recovery, wastewater recovery, fuel cells, and nuclear waste storage. This paper presents results from a series of molecular dynamics simulations of thermo-osmosis in charged silica nanochannels that advance the understanding of the phenomenon. Simulations with pure water and water with dissolved NaCl are considered. First, the effect of surface charge on the sign and magnitude of the thermo-osmotic coefficient is quantified. This effect was found to be mainly linked to the structural modifications of an aqueous electrical double layer (EDL) caused by the nanoconfinement and surface charges. In addition, the results illustrate that the surface charges reduce the self-diffusivity and thermo-osmosis of interfacial liquid. The thermo-osmosis was found to change direction when the surface charge density exceeds −0.03C · m–2. It was found that the thermo-osmotic flow and self-diffusivity increase with the concentration of NaCl. The fluxes of solvent and solute are decoupled by considering the Ludwig–Soret effect of NaCl ions to identify the main mechanisms controlling the behavior. In addition to the advance in microscopic quantification and mechanistic understanding of thermo-osmosis, the work provides approaches to investigate a broader category of coupled heat and mass transfer problems in nanoscale space.


INTRODUCTION
Thermo-osmosis refers to the motion of a fluid due to a temperature gradient. The interest in this phenomenon in micro-and nano-scale spaces has recently increased due to its potential application in a wide range of fields. For example, in the area of low carbon energy conversion, nano-porous membranes are being explored for converting low-grade waste heat into mechanical energy, which in turn can be converted into electrical energy 1,2 or in geo-materials where the transport of fluids in tight porous media under thermal effects is of importance in the energy geo-structures (e.g., energy foundations, geological disposal of radioactive waste). 3 In the emerging field of micro/nano-fluidics, thermo-osmosis is important to describe the fluid transport and mixing, with many promising applications such as the chip-level cooling without extra pumps 4 and accurate control of (micro)nanoconfined cells, flexible fibers, colloids, etc., via shear flows. 5 The understanding of thermo-osmosis in nano-porous media, including the underlying microscopic mechanisms, has been substantially improved as a result of recent theoretical, numerical, and experimental studies, e.g., refs 2 and 4−8. The thermo-osmotic component of fluid flow and transport in many natural and technical systems occurs under the effects of surface charges and saline environment (e.g., charged nano-porous membranes for ionic selectivity, 9 osmotic energy recovery, 10 electric energy storage technologies, 11 and clay nanopores for saline geofluids are mostly negatively charged 12 ). Thermo-osmosis in such systems is predominantly controlled by the pore size (nano-confinement effect), the charge of the solid surface (surface charge effect), and the salinity of pore fluid (ionic strength effect). Quantification of the nano-confinement effects on thermo-osmosis have progressed in literature. 6,7 However, understanding the surface charge effects and ionic strength effects on thermo-osmosis is still limited, in contrast to the understanding of these two effects on other transport phenomena such as thermal diffusivity and Fickian diffusivity. 13−15 This paper aims to provide new insights into the surface charge effects and ionic strength effects on thermo-osmosis and the mechanisms controlling the phenomenon by isolating the coupled processes of thermo-osmosis and thermal diffusion (or Soret−Ludwig effect).
The focus of this work is on (1) thermo-osmosis of pure water confined by negatively charged silica nanochannels, where the pore fluid systems contain pure water and the surface charges are compensated by a small number of monovalent sodium counterions, and (2) thermo-osmosis of aqueous NaCl solutions confined by uncharged silica nanochannels, where the pore fluid systems contain NaCl with several concentrations. The choice of these systems was based on the fact that silica phases have become crucial constituents of many nanodevices for drug delivery, 16 water desalination, 17 and biomolecule detection. 18 The work addresses questions of long-standing scientific interest and current technological relevance, namely, quantification of the interfacial structure, e.g., aqueous electrical double layer (EDL) on the charged surface, and the dynamical properties of interfacial liquid. 19,20 Using nonequilibrium molecular dynamics (NEMD) simulations, it is shown that the surface charges and ionic strength strongly affect the thermo-osmotic response of the interfacial liquid, which can be attributed to the alteration of aqueous EDL.

MOLECULAR DYNAMICS MODEL
Two MD systems have been studied in this work. The first system (system A) was designed to study the nanofluidic systems containing pure water and confined by charged silica substrates. Figure 1a shows the first system that includes a nano-channel of silica with a channel size of 2.4 nm. The model contains around 1400 water molecules, 4400 silica surface atoms, and between 0 and 32 counterions (Na + ) depending on the surface charge density of the charged silica substrate. The system consists of pure water and a small number of counterions confined in a charged silica slit nanochannel. Periodic boundary conditions are imposed in all directions. It has been previously shown that the Onsager reciprocal relations are valid at the nanoscale, 2,6,7,21 i.e., the thermo-osmotic coefficient relating the mass flux to the temperature gradient equals the mechano-caloric coefficient relating the heat flux to the pressure gradient. Because the application of the pressure gradient is computationally simpler, and the steady-state flow can be achieved faster, the thermoosmotic coefficient is determined here by calculating the mechano-caloric coefficient. A constant pressure gradient, is created by applying an external particle force, f i , in the x direction on every water molecule in the system with N molecules (see Figure 1a).
The second system (system B) was designed to study the effects of ionic solutes (NaCl) on thermo-osmosis. The system considered is neutral silica substrates containing aqueous NaCl solutions (see Figure 1b). The concentrations studied are 0.00, 1.00, 1.23, 1.60, 2.29, and 4.04 mol/kg. Upon applying the constant pressure gradient, the system reaches a steady state with a constant planar Poiseuille flow after a period of equilibration time. The velocity profiles v x (z) at the steady state are recorded and fitted to the Stokes where d = L z − 2z s is the actual channel width after deducting the thickness of immobile layers z s . η and b are the fitted liquid viscosity and interlayer slippage length, respectively. The thermo-osmotic coefficient m 21 is computed by 6  where J h is the heat flux. It is calculated from the excess specific enthalpy δh(z) and the velocity profile v x (z) by where the specific excess enthalpy is calculated by where h b is the specific enthalpy of the bulk liquid and h(z) is the local specific enthalpy. The latter is calculated by 23 where p xx (z) is the calculated pressure component of interfacial liquid along the x direction using the atom-based virial equation 24 and u(z) is the specific internal energy (per unit volume) of interfacial liquid. For sufficiently wide channels, h b = h(z = 0). 2,23,25 Considering the hydrodynamic interlayer slippage length, b, and the thickness of the immobile layer, z s , the theoretical value of thermo-osmotic coefficient, M 21 , can be obtained by 2 The difference between m 21 and M 21 is that the former uses the measured velocity profile of interfacial liquid while the latter assumes a linear velocity profile in the boundary layer with the same hydrodynamic boundary condition as the former.
The EDL is the origin of many electrokinetic transport phenomena at charged interfaces such as electro-osmosis, electrophoresis, streaming potential, and streaming current. Like the streaming current experiments, the induced fluid flow of interfacial liquid allows for quantifying the widely used zeta potential, ζ, which is defined as the electrostatic potential at the position of the plane of shear. The zeta potential is an important quantity in electrokinetic transport because it quantifies the coupling between electrostatic properties of the ion cloud and hydrodynamics of the solvent within the EDLs. Following previous studies 26,27 where is the streaming electric current density, ρ q (z) is the charge density, − ∇ p is the applied pressure gradient, and η b and ε b are the bulk viscosity and dielectric permittivity of the liquid, respectively. Similarly, the theoretical prediction of the zeta potential, ζ*, can be obtained by 28,29

MOLECULAR DYNAMICS SIMULATIONS
All NEMD simulations were implemented with LAMMPS. 30 The Velocity-Verlet integration is used to solve the Newton's equations of motion with a time step of 3.0 fs.

Silica Substrates.
The silica substrates with surface charge densities of 0.00, −0.03, −0.06, and −0.12 C · m −2 were built following our previous work. 15 The partial charges on atoms were determined based on the study of Kroutil et al. 31 The interatomic interactions were computed by the ClayFF force field. 32 3.2. Nanoconfined Liquid. Sufficient water molecules were packed into the silica nanochannels, and 0−32 sodium counterions with positive charges were dissolved in the interfacial water to compensate for the negatively charged silica surfaces. For system B, NaCl ions were further dissolved into the interfacial liquid to achieve the specified ion concentrations. The formed interfacial liquids were described with the rigid SPC/E water model 33 and the force field by Smith and Dang 34 for the ion−ion(water) interactions. The Coulombic and Lennard-Jones 12−6 potentials were adopted to describe the silica−liquid interactions, and standard Lorentz−Berthelot mixing rules were used to determine the cross-interaction parameters between various atomic types. A cutoff radius of r c = 1.5 nm was used for short-range interactions, and the long-range electrostatic interactions were determined by the particle−particle−particle−mesh (PPPM) method. 35 The bonds involving hydrogen atoms in the system were constrained during simulation to enable a longer timestep, i.e., 3.0 fs. The same large time step of 3 fs, together with an increased interatomic force cutoff, has been used by some previous MD studies, e.g., ref 36, to perform NEMD simulations of rigid SPC/E water. In our work, the bulk atoms of silica slabs are tethered to their respective initial positions by independent springs during simulation and the surface atoms are flexible with constrained O−H bonds in silanol groups to enable a longer timestep. We repeated the simulations for several cases with a smaller time step of 1 fs. The results showed a negligible difference; therefore, a time step of 3 fs was used.

NEMD Simulations.
Energy minimization was first performed on the established system using the steepest descent algorithm. A relaxation stage composed of a 0.1 ns run in an NVT ensemble and a 1 ns run in the NPT ensemble, using a Nose−Hoover thermostat and a barostat, was then carried out to achieve a system pressure of P = 600 bar with a coupling time of 3 ps and a system temperature of T = 300 K with a coupling time of 0.3 ps. The adopted thermodynamic conditions are consistent with our previous study 15 on the thermal diffusion of aqueous NaCl solutions occurring at the same nanofluidic system. This enabled to directly use the thermal diffusion coefficient obtained in that study and further decouple the fluxes of solvent and solute (see details in Discussion). Following the previous studies and simulation procedures of Kroutil et al. 31 and Quezada et al. 37 on the same systems where silica surfaces are immersed in aqueous NaCl solutions, the semi-isotropic pressure coupling condition is used in the NPT simulation run, where the Nose−Hoover barostat set to 600 bar was used with scaling only in the z direction of the simulation box, which is due to the strong structural heterogeneity of the simulation box in the z direction and x/y direction.
The simulation was followed by applying the external particle force and creating a pressure gradient equaling to − ∇ p = 0.055 GPa/nm inside the nanochannel, and an equilibrium run was performed. The steady state was achieved after approximately 9 ns. Finally, a production run was performed for 18 ns. The system configuration was recorded and analyzed every 0.3 ps to produce thermodynamic data, density, charge, excess enthalpy, velocity profiles, etc., by dividing the simulation cell into typically 600 sampling layers along the z direction. During the equilibrium and production runs, the system temperature was modulated by a canonical sampling thermostat adopting global velocity rescaling with Hamiltonian dynamics. 38 The surface atoms (see Figure 1c) were thoroughly flexible throughout the simulation. In the NVT simulations, bulk atoms were constrained in all directions (see Figure 1c), while in the NPT simulations, they were permitted to move only in the z direction to enable the pressure equilibration. These fixations were achieved by applying a spring with a force constant of 1000 kJ mol −1 nm −2 to tethered atoms in the specified directions. For better comparison between different cases, the zero line was specified as the averaged position of all surface silicon atoms (see Figure 1c).

Surface Charge Effect.
For every simulation case of system A, the mass/number and charge density profiles of water and Na + as a function of distance from the silica surface were calculated. The results are shown in Figure 2. The zero plane of the silica surface defined in Figure 1c is fixed at the abscissa x = 0 for a clearer comparison. The molecular organization of the interfacial water into layers near the silica surface can be clearly identified by the peaks and valleys shown in Figure 2a,b. This layering becomes increasingly pronounced with decreasing surface charge density. The cation number density near the silica surface in Figure 2a shows that Na + form three categories of adsorbate species: (i) inner-sphere surface complexes (ISSC), (ii) outer-sphere surface complexes (OSSC), and (iii) diffuse swarm (DS) species. 39 The corresponding distinct counterion adsorption planes, namely, 0-plane, β-plane, and d-plane of the widely used triple-layer model (TLM 40  shows that with decreasing surface charge density, more water molecules and sodium counterions are attracted to the charged surfaces, forming ISSC and OSSC with different ionic solvation structures and less mobility compared with those in the bulk liquid, where ISSC is partially hydrated and adsorbed directly on the silica surfaces and OSSC is fully hydrated. The accumulation of sodium counterions on the surface increases with the surface charge and is limited mainly to ISSC and OSSC regions, where the local viscosity will increase accordingly. An enhanced local fluid viscosity in the EDLs near the solid surface compared with that in the bulk fluid away from the solid surface has been reported by previous studies, e.g., refs 41 and 42. This results in an apparently enhanced effective viscosity of the interfacial fluid. It was also previously found that increasing the surface charge density can further increase this apparent viscosity, due to the increasing adsorption and entrapment of counterions on charged surfaces, and resultant lateral hindrance for the fluid flow of surrounding solvent molecules. 41 A similar observation is made in Section 4.3. Figure 2b shows the number density of water hydrogen atoms (H w ) and water oxygen atoms (O w ), which indicates that the water mass density in Figure 2a is dictated by the number density of water oxygen atoms due to its much larger mass relative to the water hydrogen atoms. Figure 2c shows the charge density distribution of interfacial liquid, which can be used to evaluate the electro-kinetic responses at the interface in the subsequent analysis. The results indicate that the charge density distribution is concurrently determined by the distribution of sodium counterions, water hydrogen atoms, and water oxygen atoms, and its peaks and valleys are consistent with those in Figure 2b. The charge density is zero in the bulk liquid due to the random arrangement and distributions of water molecules and sodium counterions.
The thickness of the interfacial EDLs (including Stern layer and diffuse ion swarm in Figure 2a) can be characterized by the Debye length, λ D . Within the classical mean-field theories and for dilute electrolyte (low ion density) regime, the Debye length can be calculated by 43 where ε r is the relative dielectric permittivity of interfacial fluid, ε 0 is the vacuum permittivity, k B is the Boltzmann's constant, T is the temperature, ρ ∞i is the number density of ion type i in the bulk region in Figure 2a Substituting ε r = 70.8 44 into eq 9 gives λ D in the range of 0.49−0.67 nm for charged silica surfaces in Figure 2a. This theoretical screening length is in good agreement with the one obtained by MD simulations.

Concentration Effect.
The EDL structure of the interfacial liquid for system B is presented in Figure 3. Figure  3a−c shows that the shape of the number/mass density profile of water molecules and NaCl ions and the position of TLM planes are hardly modified by the addition of NaCl ions. This is confirmed by the unchanged positions of peaks and valleys in the profiles and the vertical blue dashed lines. Accordingly, the charge density distribution in Figure 3d is also kept constant with the addition of NaCl ions. In addition, Figure 3c shows an obvious cation condensation and anion exclusion (repulsion) at the silica surface, where the number density distribution profile of chloride ions roughly starts beyond the 0 plane. The distribution of sodium and chloride ions becomes equal in the bulk liquid. The increasing concentration also leads to enhanced ion pairing between sodium and chloride ions, which might concurrently decrease the mobility of each ion. Note that in system B, where concentrated aqueous NaCl solutions (beyond the Debye−Huckel region) are confined by two planar charged surfaces, a previous experimental study 45 with a similar setup has pointed out that in such concentrated electrolytes, the theoretical Debye screening length of eq 9 will become inapplicable and substantially underestimate the real thickness of EDL, i.e., the electrostatic screening length, which was also found to increase with concentration. This is also consistent with our MD results in Figure 3, where λ D is estimated to be only within 0.13−0.26 nm according to eq 9.  Figure 4a and b, respectively, where the latter is calculated using eq 4. The regions of the boundary layer and the bulk liquid are consistent with those in Figure 2a. The excess enthalpy profiles in Figure  4b oscillate notably in the boundary layers and diminish to zero in the bulk liquid. It is found that the excess specific enthalpy tends to be more negative with decreasing surface charge density, due to the enhanced attractive interactions relative to the weakened repulsive interactions between interfacial liquid particles and the negatively charged surfaces, 46 e.g., the electrostatic attractions between sodium counterions and surfaces. Following previous studies, 6,23,47 thermo-osmotic forces acting on interfacial liquid under a unit temperature gradient, i.e., − ∇ T = 1, are calculated and shown in Figure 4c According to the TLM defined in the previous section, it is found that the interfacial liquid located in the β-plane and dplane is thermophilic, i.e., the liquid will migrate from colder to hotter areas under a temperature gradient, while the liquid in the 0 plane is thermophobic, i.e., the liquid will migrate from hotter to colder areas. In addition, Figure 4c,d shows that the thermo-osmotic force acting on the interfacial liquid tends to be more negative with decreasing surface charge density, which indicates that thermo-osmotic flow will tend toward the −x direction when the surface charge density is applied. This tendency can be presented more clearly by comparing the integrals ∫ −Ld z /2 L z /2 f x (z)dz and ∫ −Ld z /2 L z /2 f x P (z)dz in Figure 4c and d, respectively, which decrease with increasing surface charge density. In fact, the effective excess enthalpy, ∫ −Ld z /2 L z /2 δh(z)dz, has been used in many previous studies, e.g., ref 5, to determine the direction of thermo-osmosis, as is shown in Figure 4b. This is the same as our integrals of thermo-osmotic forces. Therefore, it can be expected that the thermo-osmotic  coefficient may become negative with decreasing surface charge density. This expectation will be proved shortly.

Concentration Effect.
For system B, the thermodynamic properties of interfacial liquid are computed and presented in Figure 5. It is found from Figure 5a,b that p xx (z) generally keeps constant with the addition of NaCl ions, while δh(z) varies significantly, indicating that the internal energy of the interfacial liquid changes a lot with the addition of NaCl ions according to eq 5. Figure 5c,d shows that the addition of NaCl ions increases the magnitude of the thermoosmotic force on the interfacial liquid in either the +x or −x direction. The position of these enhancements is consistent with that of TLM planes, indicating that the addition of NaCl ions facilitates the thermo-osmosis in either the +x or −x direction depending on their occupied TLM planes. Furthermore, the integral of the thermo-osmotic force over the channel section monotonically increases with the addition of NaCl ions, indicating an overall flowing trend toward the +x direction.

Hydrodynamic Properties of Nanoconfined Liquid. 4.3.1. Surface Charge Effect.
For system A and by applying an external particle force on each water molecules, the generated Poiseuille-type parabolic velocity profiles of interfacial liquid confined by silica surfaces are plotted in Figure 6b−e. These can be described well by the continuum hydrodynamics in eq 1, even in the range of EDLs. Figure 6a shows the definitions of the fitting parameters in eq 1, where the shear plane position, z s , in which the slip boundary condition is imposed, is determined as the positions of first peaks of the water density profiles in Figure 2a. The derived fitting parameters in Figure 6f show that with decreasing surface charge density, the viscosity of confined liquid increases and the thickness of (stagnant) immobile layers decreases. These observations are consistent with the results in the previous subsection that more counterion complexes are formed near the charged silica surfaces with decreasing surface charges. The vertical dashed lines in Figure 6b−e indicating the TLM planes show that the interfacial liquid with less mobility, i.e., being nearly stagnant, is located between the silica surfaces (zero line in Figure 1c) and 0 planes. The description of the immobile layer is from classical continuum theories, which define a fluid adsorption layer near the solid surface with infinite fluid viscosity and zero fluid mobility. 42 Notably, the MD data in Figure 6 suggests that the fluid molecules located within the defined immobile layer is nearly immobile, but not completely immobile. Due to the strong solid−fluid friction and higher local fluid viscosity, they move relatively slowly and nearly synchronically as is observed by the relative flattening of velocity profiles. These findings are consistent with a previous MD study of Moh et al. 48 on the decane transport through slit calcite pores. In addition, due to the hydrophilic nature of the silica surface, 6,7 the interlayer slippage length, i.e., the velocity jump, is small, and its magnitude decreases with the surface charge density becoming more negative. This is consistent with a previous study on the interface between aqueous sodium chloride and charged graphene. 26 The decreasing interlayer slippage length can be attributed to the enhanced binding (trapping) effect of Na + counterions on the negatively charged sites of the silica surfaces with decreasing surface charge. The trapped counterions stick out from the silica surface inducing a viscous Stokes drag.

Concentration Effect.
For system B, by applying an external particle force on every water molecule, the generated velocity profiles of interfacial liquid for different cases are presented in Figure 7a−f, which can be well described by the continuum hydrodynamics. The fitting parameters in Figure 7g show that increasing NaCl concentration leads to monotonically increasing liquid viscosity η and interlayer slippage length b, with the thickness of immobile layer z s being almost constant. Overall, a general decrease of liquid mobility can be observed when the concentration increases. Note that even the channel size is fixed at around 2.4 nm to exclude the channel size effects, whereas the addition of NaCl ions may slightly enlarge the channel size, e.g., see Figure 7a−f.

Self-Diffusivity of the Interfacial Liquid.
The selfdiffusion coefficient is an important quantity, which can be used to estimate the transport distance of mass diffusion within a given time, e.g., ref 49. For systems A and B, the onedimensional self-diffusion coefficient of the interfacial liquid, including water molecules (at the center of mass) and NaCl ions (if exist), in the x direction (D x ) were calculated for the entire interlayer space by using the Einstein relation: 50 ) where d = 1 is the diffusion dimensionality, τ is the simulation time, and MSD x (τ) = ⟨(x(τ) − x(0)) 2 ⟩ is the mean-square displacement of interfacial liquid along the x direction. The self-diffusivity of the interfacial liquid under bulk conditions is determined for comparison by averaging the mean-square displacements over three dimensions and setting d = 3. The dependance of D x on the surface charge density and ionic strength is presented in Figure 8a and b, respectively. It can be observed from Figure 8a that both the nanoconfinement and the surface charges can reduce the self-diffusivity of the interfacial liquid. The obtained bulk value, i.e., (2.91 ± 0.01) × 10 −9 m 2 /s, is consistent with a previous MD study for the same rigid SPC/E water model under a similar thermodynamic condition, i.e., (2.97 ± 0.05) × 10 −9 m 2 /s. 51 This provides further confidence in the accuracy of the simulations presented here. Furthermore, Figure 8b shows that the increasing ionic strength can also reduce the self-diffusivity of the interfacial liquid. The same varying trends of interlayer diffusion with increasing surface charge density and ion concentration have also been reported by Greathouse et al. 49 for MD systems consisting of montmorillonites and interlayer aqueous NaCl/CaCl 2 solutions. They explained that a higher surface charge density makes the interlayer environment more hydrophilic with slower diffusivity of interfacial liquid, which well agrees with our results in Figures 8a and 6f. Since the diffusion coefficient is inversely proportional to the viscosity according to the Stokes−Einstein relation, our results for interlayer self-diffusivity in Figure 8 agree well with the calculated viscosity in Figures 6f and 7g. 4.4. Thermo-Osmotic Properties of Nanoconfined Liquid. 4.4.1. Surface Charge Effect. For system A, the thermo-osmotic coefficient m 21 is calculated using eq 2. Figure  9a shows the variation of the thermo-osmotic coefficient with surface charge density. Three independent simulations with  various initial states of interfacial liquid were performed to determine the variability of the computed m 21 . The variabilities are shown by the error bars in the figure. It was found that decreasing the surface charge can significantly reduce the thermo-osmotic coefficient. The coefficient changes sign to negative when the surface charge density is below −0.03C · m −2 , indicating that the thermo-osmotic flow would change its direction and liquid would migrate toward the hot region. At this threshold, there is no thermo-osmotic flow. This conclusion is consistent with the analysis of thermo-osmotic force in Figure 4. Note that the calculated thermo-osmotic coefficients are comparable to those from previous experiments 52,53 on silica-based materials, ±(10 −10 ∼10 −9 )m 2 /s, and that both positive and negative signs of the thermo-osmotic coefficient have been observed in those experiments. Bregulla et al. 52 have attributed the negative thermo-osmotic coefficient to the electrostatic contribution of EDL caused by surface charges. In addition, it was found that the thermo-osmotic coefficient (10 −10 − 10 −8 m 2 /s) is of the same order of magnitude or larger than the self-diffusion coefficient (10 −10 − 10 −9 m 2 /s). This shows that the thermo-osmotic flow at the nanoscale ought not be overlooked when considering coupled transport. Substituting the corresponding values in Figures 4b  and 6f into eq 6, the computed M 21 is presented in Figure 9a. This is found to agree well with the calculated m 21 . The discrepancies come from the linear assumption for the velocity profile in the boundary layers when deriving eq 6. 2 A high consistency between m 21 and M 21 is found at low surface charge densities, i.e., below −0.03C · m −2 , while their discrepancy increases with the increasing surface charge density.

Concentration Effect.
The thermo-osmotic response for different cases in system B is quantified. The thermoosmotic coefficient, the measured m 21 , and the predicted M 21 are presented in Figure 9b. It is found that the measured thermo-osmotic coefficient m 21 generally increases with the addition of NaCl ions, although it experiences a slight drop at the concentration of around 1.6 mol/kg. The predicted one, M 21 , has a trend consistent with m 21 but with a lower magnitude. It is found that the increasing NaCl concentration will enlarge the difference between m 21 and M 21 . This is consistent with the effect of surface charge density on the difference between m 21 and M 21 , as is presented in Figure 9a. Similarly, the difference is due to the linear assumption for the velocity profile in the boundary layers when M 21 is derived. According to the formulas for m 21 and M 21 , the thermoosmotic response is the combined action of thermodynamic properties in Figure 5b and hydrodynamic properties in Figure   7g, where the magnitude of the former one increases with increasing NaCl concentration, in either +x or −x direction, while the latter one decreases with increasing the NaCl concentration. Their competition finally determines the variation of the thermo-osmotic coefficient seen in Figure 9b 51 where ε 0 is the vacuum dielectric permittivity, the computed zeta potential from eq 7 is presented as solid symbols in Figure 10a, which shows a monotonic decrease of ζ potential with decreasing surface charge density from 0 to −0.12 C · m −2 . Since the surface charge density of 0.00, −0.03, −0.06, and −0.12 C · m −2 corresponds to the pH values of around 2.0−4.5, 7.5, 9.5, and 11, respectively, 31 the predicted decrease of zeta potential with increasing pH values is consistent with previous experimental and numerical measurements for silica/aqueous electrolyte solution interfaces by Brkljacǎ et al., 54 including the sign inversion at a specific pH value. The calculated values of the zeta potential are also in the same order of magnitude with the experimental observations. 54 A nonzero zeta potential for the uncharged surface is observed in Figure 10a, which has also been reported by numerical studies on the electro-osmotic flow in hydrophobic nanochannels. 28,29 The calculated ζ* from eq 8, shown in Figure 10a with hollow symbols, is in the same order of magnitude as ζ. However, it shows a negligible variation with the surface charge density. The discrepancies between ζ* and ζ also come from the linear assumption for the velocity profile in the boundary layers when deriving eq 8.

Concentration Effect.
For system B, the zeta potential, measured ζ, and predicted ζ* are quantified and shown in Figure 10b. It is found that ζ linearly and slightly decreases with increasing NaCl concentration, while ζ* shows a different variation and a larger magnitude. Similarly, the discrepancy is due to the linear assumption for the velocity profiles in the boundary layers when ζ* is derived.

DISCUSSION
The aim of this section is to decouple the fluxes of solvent and solute based on the results of system B. The movement of nanoconfined ionic solutions due to temperature gradients is not a pure thermo-osmosis problem. It is the result of thermal osmosis combined with thermally induced chemical transport processes (thermal diffusion). Therefore, the thermo-osmotic coefficient quantified in Figure 9b intrinsically includes the Ludwig−Soret effect of ionic species. Based on the thermoosmotic coefficient measured in this study and the thermal Figure 10. The variations of zeta potential with (a) surface charge densities and (b) NaCl concentrations, where ζ is computed using eq 7 and ζ* is calculated using eq 8. diffusion coefficient obtained from a previous one, 15 the fluxes of solvent and solute are de-coupled to assess the extent to which the thermo-osmosis and thermal diffusion co-exist in the nanoconfined fluid mixture when a temperature gradient is imposed. Figure 11 shows that thermo-osmosis is the creep fluid flow of the interfacial liquid propelled by the pressure gradient induced by the temperature gradient in the boundary layers and this phenomenon will not occur under bulk conditions because thermal gradients do not lead to pressure gradients in the bulk liquid. The thermo-diffusion (i.e., Ludwig−Soret effect) is the separation tendency of different components in the interfacial liquid under a temperature gradient, which occurs in both nanoconfined and bulk liquids. In the case of aqueous salt solution, the thermal diffusion is the flux of salt (solute) relative to water molecules (solvent). Figure 11b shows that the thermally induced solution flux in a nanochannel is the combination of thermo-osmotic flow and thermo-diffusive flow. Therefore, the thermo-osmotic coefficient quantified in our study includes the Ludwig−Soret effect of NaCl ions.
In the nonequilibrium thermodynamics and in the linear regimes, the flux of nanoconfined solution, J Solution (volumetric flow rate, unit: m 3 /s), under a temperature gradient and without the pressure and chemical concentration gradients can be described by where V̅ Solvent and V̅ Solute are the molar volumes of solvent and solute, J Solvent and J Solute the fluxes of solvent and solute, m 21 the thermo-osmotic coefficient of the solution as obtained in this study, and A the cross section for the solution flux. The excess solute flux (unit: mol/s, the solute flux relative to the solvent flux, computed as the difference between the observed solute flux J Solute and the solute flux predicted according to the total solute concentration, c Solute , and solution flux, J Solution ) is induced by the thermal diffusion: where ρ is the molar density of the mixture, x Solute the molar fraction of the solute, and D T T the thermal diffusivity of the solute.
Equations 11 and 12 show that the solute flux and solvent flux can be determined if all other data are known. Taking our case with the NaCl concentration of 4.04 mol/kg, the surface charge density of 0.00C · m −2 , and a unit temperature gradient as an example, the corresponding data from this study and the previous one on the thermal diffusion 15 in the same nanofluidic system is presented in Table 1. The obtained J Solvent and J Solute indicate that the thermally induced migration of water molecules and dissolved NaCl ions are in opposite directions due to the thermophilic thermal diffusion of NaCl ions, i.e., D T T < 0.

CONCLUSIONS
MD simulations were used to investigate the thermo-osmotic response of pure water confined in charged silica slit nanochannels with surface charge density varying from 0.00 to −0.12 C · m −2 . The NEMD results showed that decreasing surface charges reduces progressively the self-diffusivity and thermo-osmotic mobility of the nanoconfined liquid. The thermo-osmotic flow was found to change direction for surface charge densities less than −0.03C · m −2 . It was shown that the structural modifications caused by the nanoconfinement, and the surface charge were significantly correlated with the quantified thermo-osmotic properties of the confined liquid. In addition, by further dissolving NaCl ions to the interfacial water confined by uncharged silica nanochannels, the ion concentration effect was found to facilitate the interfacial thermo-osmotic flow. The fluxes of solvent and solute in the thermo-osmotic flow were decoupled by considering the Ludwig−Soret effect of NaCl ions. It was shown that these two fluxes were in opposite directions due to the thermophilic thermal diffusion of NaCl ions.
The presented study was carried out on realistic surfaces. It is the first one to reveal that the thermo-osmotic response of nanoconfined liquid can be manipulated by changing the surface charge density and ionic strength. The study provides molecular-level quantifications and explanations for the existing experimental observations and macroscopic analysis. The simulation methods adopted in this work and the insights