Thermo-osmosis in hydrophilic nanochannels: mechanism and size effect

Understanding thermo-osmosis in nanoscale channels and pores is essential for both theoretical advances of thermally induced mass ﬂ ow and a wide range of emerging industrial applications. We present a new mechanistic understanding and quanti ﬁ cation of thermo-osmosis at nanometric/sub-nanometric length scales and link the outcomes with the non-equilibrium thermodynamics of the phenomenon. The work is focused on thermo-osmosis of water in quartz slit nanochannels, which is analysed by molecular dynamics (MD) simulations of mechano-caloric and thermo-osmotic systems. We investigate the applicability of Onsager reciprocal relation, irreversible thermodynamics, and continuum ﬂ uid mechanics at the nanoscale. Further, we analyse the e ﬀ ects of channel size on the thermo-osmosis coe ﬃ cient, and show, for the ﬁ rst time, that these arise from speci ﬁ c liquid structures dictated by the channel size. The mechanical conditions of the interfacial water under di ﬀ erent temperatures are quanti ﬁ ed using a continuum approach (pressure tensor distribution) and a discrete approach (body force per molecule) to elucidate the underlying mechanism of thermo-osmosis. The results show that the ﬂ uid molecules located in the boundary layers adjacent to the solid surfaces experience a driving force which generates the thermo-osmotic ﬂ ow. While the ﬁ ndings provide a fundamental understanding of thermo-osmosis, the methods developed provide a route for analysis of the entire class of coupled heat and mass transport phenomena in nanoscale structures.


Introduction
Thermo-osmosis is a non-equilibrium heat and mass transfer cross-phenomenon. It is one of the main physical mechanisms of thermally induced particle motion. While the other observed mechanisms, generically referred to as the thermophoresis (also thermal diffusion or the Ludwig-Soret effect), involve thermally induced migration of particle mixtures, e.g., thermally induced flux of aqueous species relative to the fluid, thermo-osmosis describes the migration of a single-component fluid due to thermal gradients. This mechanism is encountered in many areas of biology, 1 chemistry, 2 energy, 3,4 separation and membrane technologies 5 and energy geo-structures. 6,7 Although its discovery has been reported as early as 1907, 8 the need for an in-depth understanding of thermo-osmosis has re-emerged recently due to the important role of thermo-osmosis in a wide range of engineering processes and technologies (e.g. membrane systems, thermo-osmotic energy conversion heat recovery; geo-energy systems). In the fields of energy and environment, the thermo-osmotic response of nanoporous membranes has been used for low-grade waste heat harvesting 3,4 and recovery of water from organic wastewater solution. 5 Thermo-osmosis has also been reported to be an important process for fuel cell operation. 2 In addition, quantifying the role of thermo-osmosis in tight rocks (e.g. shale) and clay soils is essential for underpinning the design and safety assessments of many geotechnical infrastructures such as nuclear waste repositories and groundsource heat or energy foundations. 6,7 Compared to the theoretical studies of thermophoresis, in particular thermal diffusion, [9][10][11][12][13] the microscopic mechanistic understanding of thermo-osmosis is rather incomplete, although its macroscopic thermodynamical description has been widely discussed and used both in modelling multiscalemultiphase couplings problems 14,15 and in experimental works. 16,17 In fact, there is a debate between two approaches for describing the underlying mechanism for thermo-osmosis, namely the interfacial approach and the energetic approach. 18 Further understanding of thermo-osmosis at the nanoscale is of particular interest because the very limited studies of thermo-osmosis at solid-fluid interfaces [19][20][21][22][23] do not allow for linking the two approaches and for a complete understanding of the underlying mechanisms.
The parameter controlling thermo-osmosis is the thermoosmotic coefficient. It describes the magnitude of fluid flux induced by a unit temperature gradient. Experimental and theoretical works have shown that the thermo-osmotic coefficient can be described as a function of the interfacial excess enthalpy distribution of the fluid and its hydrodynamic properties at the interface. 24 However, the influence of various factors on the mechanism of thermo-osmosis is still unclear. For example, the excess enthalpy distribution is affected by temperature, pressure, liquid polarity, 25 and liquid structure in the boundary layer. 19 This makes it difficult to be experimentally measured. 21,26 The experimental results from previous studies 16,17 also show that the magnitude and the direction of thermo-osmotic flow can vary for different types of chemicals. Furthermore, the existing formulas for the thermo-osmotic coefficient are derived on the basis of a local continuum theory, 21 e.g. Navier-Stokes equation, where assumptions are made on the fluid properties near the solid surface including the excess enthalpy distribution and fluid viscosity. Such continuum approximations may be sufficient for solid-fluid interfaces at larger length scales. However, for the solid-fluid interfaces at the nanometric/sub-nanometric scales, the validity of the continuum description for the interfacial liquid structure or hydrodynamics is not proven. Understanding thermo-osmosis at such length scales requires a mechanistic approach.
Analysis of thermo-osmosis at solid-fluid interfaces by molecular dynamics (MD) is the natural approach to address the knowledge gap, as it provides direct, comprehensive and detailed information on the microscopic thermodynamic states and hydrodynamic characteristics at the interfaces. Several groups have started utilising the approach. For instance, Ganti et al. 21 reported that different schemes of MD simulations to measure the thermo-osmotic coefficient would yield consistent results. Fu et al. 20 investigated the effect of interfacial wettability and interfacial hydrodynamic characteristics on the thermoosmotic coefficient. However, little progress has been made on the understanding and quantification of the relation between the thermo-osmotic response and the interfacial fluid structure.
The aims of this work are to develop a new mechanistic understanding and theoretical description of thermo-osmosis at nanometric/sub-nanometric length scales and to link with the non-equilibrium thermodynamic formulation of the coupled phenomenon. The system selected for investigation is quartz and water. This choice is dictated by two constraints. Firstly, the molecular dynamics limits the spatial dimension of the system that is analysed to nanochannels of several nanometres width. Secondly, the inclusion of charged minerals (e.g., kaolinite) and aqueous solutes will not produce pure thermo-osmotic flow, as surface charges and concentration gradient will create unwanted coupling effects. Quartz is electrically neutral and accounts for approximately 20% of the Earth's exposed crust, 27 making it a major constituent of soils, clays, and rocks. It is acknowledged that such systems contain a multitude of minerals, but it is reasonable to consider that nanoscale pores are formed in single minerals. Moreover, understanding water interaction in quartz or other silica phases is of interest to many fields such as geology, 28 biology, 29 physics 30,31 and chemistry. 32,33 This specific system is also relevant to membrane science as there are membranes made up quartz and other silica phases such as hydrophilic quartz fibre membranes 34 and silica membranes. 35 While a number of studies have been dedicated to the interactions of water with quartz and other silica phases, research on temperature effects, such as thermo-osmosis, in quartz-water systems is very limited. This is in contrast with the need to understand such effects, since natural or artificial temperature gradients are ubiquitous in the real quartz-water system, and they may generate substantial, otherwise unaccounted for mass transport. Besides, quartz surfaces can be fully hydroxylated within the water environment, exhibiting strong hydrophilicity. 36 Therefore, investigation on the quartz-water interfaces can naturally provide fundamental insights to the studies on fluid properties in the hydrophilic environment.
The aims of the work are achieved by molecular dynamics simulations of thermo-osmosis in quartz slit nanochannels. Microscopic thermo-osmotic coefficients are measured by mechano-caloric and thermo-osmotic systems, introduced in sections 2.2. Results presented in sections 3.1 demonstrates the validity of Onsager reciprocal relation at nanoscale and consistency between the microscopic thermo-osmotic coefficient and the macroscopic thermo-osmotic coefficient based on irreversible thermodynamics. The dependence of the thermo-osmotic coefficient on the channel size is analysed in section 3.2. It is quantified and found to be a consequence of the changing liquid structure in boundary layers adjacent to the solid surfaces. The changing liquid structure alters either the excess enthalpy distribution or the hydrodynamic properties of the interfacial water. The mechanism of thermo-osmosis in quartz nanochannels is revealed in section 3.3 by quantifying the mechanical interactions of the interfacial water under different temperatures using a continuum approach (pressure tensor distribution) and a discrete approach (body force per molecule). It is found that the fluid molecules experience a driving force for thermo-osmotic flow only in the boundary layers adjacent to the solid surfaces. The systems analysed, the analysis methods, and the insights of this work provide a solid foundation for future research in the field of thermophoretic phenomena.

Theory and method
In this section, we describe the macroscopic/thermodynamic and the microscopic/mechanistic perspectives to thermoosmosis in nano-sized systems. We are interested in the coefficients of proportionality between different gradients ( pressure and temperature) and fluxes (fluid and heat). In order to distinguish between the two perspectives, we will denote the coefficients of the thermodynamic formulation (macroscopic coefficients) with uppercase letters, and the coefficients of the mechanistic formulation (microscopic coefficients) with lowercase letters. The thermodynamic formulation, introduced in section 2.1, shows the dependence of the macroscopic coefficients on the system size, interfacial hydrodynamics and the excess enthalpy density. Determining excess enthalpy for nano-sized systems requires molecular-level analysis, which leads to the mechanistic formulation by MD models intro-duced in section 2.2. The MD simulation details are provided in section 2.3, while the calculation of various parameters from these simulations is presented in section 2.4.

Irreversible thermodynamics of thermo-osmosis in nanochannels
Based on irreversible thermodynamics, coupled fluid flow and heat transfer can be described by the following system of equations 22,26,37,38 where, J f is the fluid flux, J h is the heat flux, M ij are macroscopic ( phenomenological) coefficients, p is fluid pressure, and T is temperature. The off-diagonal parameter M 12 in eqn (1) is the macroscopic thermo-osmotic coefficient, which controls the fluid flux induced by a temperature gradient. The off-diagonal parameter M 21 in eqn (2) is the macroscopic mechano-caloric coefficient, which controls the heat flux induced by a hydraulic pressure gradient. The Onsager reciprocal relation for this coupling states that M 12 = M 21 . 39, 40 We consider a quartz slit nanochannel and introduce x-axis parallel to and z-axis normal to the quartz surfaces (z = 0 at a surface). Further, we consider coupled fluid flow and heat transfer described by eqn (1) and (2), resulting from a non-zero applied temperature gradient along the x-axis and a zero applied hydraulic pressure gradient. In this situation, the fluid flux is induced by the thermal gradient alone and governed by the thermo-osmotic coefficient.
The changes of temperature, dT, and pressure, dp, are related to the changes of the chemical potential of a single component system, dμ, via the Gibbs-Duhem equation where, N is the number of moles of the component, S is the entropy and V is the volume of the system. We apply the Gibbs-Duhem equation to the fluid in the nanochannel. Dividing eqn (3) by the fluid volume, V, differentiating with respect to x, and rearranging terms, provides the derivative of the pressure in the form where, ρ = N/V is the number density and s = S/N is the specific entropy of the interfacial fluid molecules. The variation of the chemical potential with temperature in eqn (4) can be simplified under the condition of hydrostatic equilibrium in the bulk fluid, 21,41 i.e. for z → +∞. In this case p equals the bulk fluid pressure p b , and where, s b is the specific entropy of bulk fluid. Further, it can be assumed that the chemical potential and temperature are constant through the thickness of the nanochannel (z-direction), while pressure, entropy and density may vary. This leads to the following expression for the pressure derivative: The normal variation of the specific entropy of the interfacial fluid molecules can be described by 21 where, h(z) is the specific enthalpy of the interfacial fluid molecules. Similarly, the specific entropy of the bulk fluid can be represented by the specific enthalpy of the bulk fluid, h b . This allows for rewriting eqn (6) in the form 20,21,42 where, δh(z) is the excess enthalpy density at a distance z from the surface. Then we can then solve the linearized (Navier-)Stokes equation as where, η is the interfacial fluid viscosity (assumed to be homogeneous) and v x (z) is the fluid velocity along the direction of the temperature gradient. For the calculation of the fluid flux, we consider the general case of partial slip boundary condition at the solid-fluid interface. This is illustrated in Fig. 1, where the slip length, b, depends on the nature of the intermolecular interactions between solid and fluid. 43 The thermo-osmotic velocity in the bulk fluid is then obtained as 20 Eqn (8) shows that a temperature gradient applied to the interfacial fluid results in a pressure gradient parallel to the surface, which drives the thermo-osmosis with fluid flux given by eqn (11). Comparing eqn (1) and (11) and recalling the condition of zero applied hydraulic pressure, leads to the following expression for the macroscopic thermo-osmotic coefficient, and via Onsager reciprocal relations the macroscopic mechano-caloric coefficient When no-slip (stick) boundary condition is applied, 44 i.e. when b = 0 (or when the fluid/solid friction coefficient κ → ∞), eqn (12) reduces to which is equivalent to that of Derjaguin et al., 37 based on Onsager's linear non-equilibrium thermodynamics. 45,46 We note eqn (13) indicates that the phenomenological coefficients are inversely proportional to the fluid viscosity, which is also found in recent rigorous theoretical derivations based on mixture coupling theory and irreversible thermodynamics. 14,15 The calculation of the excess enthalpy density in the nano-channel requires atomistic analysis to which we turn in the next sub-section.

Molecular dynamics setup for nanochannel simulations
Two molecular dynamics systems are used to calculate the phenomenological coefficients controlling thermo-osmosis in quartz slit nanochannels. These include a thermo-osmotic system, illustrated in Fig. 2(a), which follows a similar carbon nanotube membrane and water system devised by Fu et al., 42 and a mechano-caloric system, illustrated in Fig. 2(b), which is developed based on a similar system by Fu et al. 20 primarily for Lennard-Jones liquid and Einstein solids, as well as graphene-water systems. Both systems contain a slit nanochannel of width L z between two quartz slabs of 1.5 nm thickness each. The size of the simulation domains in y direction is 5.0 nm (L y = 5.0 nm).
In the thermo-osmotic system, Fig. 2(a), the nanochannel is connected to two water reservoirs with pressure applied by two rigid graphene pistons. MD simulations with this system are performed with two sets of boundary conditions. The first set prescribes equal pressure at the pistons, i.e. zero applied pressure gradient, and different temperatures in the two reservoirs, i.e. non-zero applied thermal gradient along the quartz nanochannel. This allows for obtaining the microscopic thermo-osmotic coefficient m 12 , as well as m 22 if required, from calculated fluid and heat fluxes, respectively. The second set of boundary conditions prescribes different pressures at the pistons, i.e. applied non-zero pressure gradient, and equal temperatures in the two reservoirs, i.e. zero applied thermal gradient. This allows for obtaining the microscopic mechano- caloric coefficient m 21 , as well as m 11 , from calculated fluid and heat fluxes at the two ends of the nanochannel, respectively. In both sets, periodic boundary conditions are considered in y and z-direction. The MD simulations provide further the excess enthalpy density, required for the calculation of the macroscopic coefficients by eqn (12) or (13). The thermo-osmotic system is adopted only to demonstrate the thermo-osmosis phenomenon and validate the Onsager reciprocity relations, while the majority of the results are obtained with the mechano-caloric system. Details of the use of the thermo-osmotic system are given in Appendix A.
In the mechano-caloric system, Fig. 2(b), a graphene piston, which is composed of a two-dimensional monolayer of sp 2bonded carbon (C) atoms in a hexagonal honeycomb lattice, is set at the top of the upper quartz slab with a pressure of 1 atmosphere in the negative z-direction. In the simulations, we apply periodic boundary conditions in x and y directions. A body force per particle f i is exerted to every water molecule to generate an external pressure gradient ∇p parallel to the x-axisan approach commonly employed in previous studies. [47][48][49] Notably, this approach is different from the thermo-osmotic system, where the external pressure was applied by two rigid graphene pistons. Previous study 50 suggested that these two different approaches should lead to highly consistent outcomes. The generated Poiseuille flow induces an excess heat flux, which allows for independent from the thermo-osmotic system calculation of the microscopic mechano-caloric coefficient m 21 . The reason for considering the mechano-caloric system is that the simulations with the thermo-osmotic system require a relatively high computational cost. Demonstrating that the two systems lead to the same coefficient and that the Onsager reciprocal relations are fulfilled will allow for subsequent parametric studies with the computationally simpler mechano-caloric system. Importantly, the mechano-caloric system is used for explaining the mechanism of thermo-osmosis in section 3.3.
To observe the structural and dynamical behaviours of this quartz-water system, the simulation domains were segmented into equal bins of size 0.004 nm. The selection of this size of bins is justified in Appendix B. The bins have been used to calculate average temperature, mass density, number density, charge, excess enthalpy, velocity, and pressure profiles.

Molecular dynamics simulation procedure
All MD simulations have been performed with the LAMMPS 51 and visualized with the VMD 52 package. The crystal structure of α-quartz (α-SiO 2 ) 27,53 is adopted to build the quartz slit nanochannels, and sufficient water molecules are added to fill the nanochannels and reservoirs. The typical surface d = (010) 54,55 of bulk α-quartz structure is cleaved as the quartzwater interface. This is illustrated in Fig. 3 together with the distribution of partial charges at the quartz surface sites. Due to the discontinuity of the bulk α-quartz structure at the quartz-water interface, the under-coordinated surface oxygen (O) and silicon (Si) atoms are capped with hydrogen (H) or hydroxyl (OH) groups as in previous studies. 30,56,57 Simulations are conducted in the NVT ensemble with a 1.0 fs time step, and the long-range coulombic (electrostatic) interactions are resolved by the particle-particle and particle-mesh (PPPM) method 58 and temperature is controlled by a canonical sampling thermostat that uses global velocity rescaling with Hamiltonian dynamics. 59 Regarding the force field, ClayFF, parameterized by Cygan et al. 60 for describing the interaction of clays with water, is adopted to describe the quartz-water system in our simulations and the flexible single point charge (SPC) 61,62 water model is used in this forcefield. 60 Note that in the ClayFF force field, the bond interactions only define the hydroxyl (-OH) on the quartz surfaces and O-H bonds in water molecules. In our simulations, some silicon (Si) atoms in the upper and lower quartz slabs are fully or partially constrained to achieve the stability of the nanochannels during the simulations. The carbon (C) atoms in the graphene pistons are only allowed to move in one direction to apply constant pressure to the fluid. Fig. 4 illustrates the degrees of freedom of atoms in mechano-caloric system and thermoosmotic systems. Thermostats are applied to all atoms except for the constrained silicon (Si) atoms.
The graphene pistons are maintained rigidly in order to restrict potential distortions and to exert a constant and uniform normal pressure on the water. A Lennard-Jones 12-6 potential and Lorentz-Berthelot combining rules 60,63 are used for the description of the intermolecular interactions between graphene pistons and water molecules. Table 1 lists the adopted potential parameters for graphene pistons, which are similar to those used by Wang et al. 64 A cutoff radius of 1.5 nm is used for both Lennard-Jones interactions and coulombic interactions in the system.

Calculations of fields and coefficients
The pressure tensor of interfacial water is calculated by using the atom-based virial equation. 65 The pressure variation of interfacial water along x-direction (or y, z-direction if required) is calculated by 21 where, a = x, y, z, and 〈·〉 is the time average of the variable. The first term on the right-hand side of eqn (14) describes the kinetic energy contribution, where T (x) is the water temperature in the bin at position x, ρ(x) is the number density of water in the bin at position x, and k B is the Boltzmann constant. The second term is the virial contribution due to intramolecular and intermolecular interactions, where V(x) is the bin volume at position x, and W aa is given by 65 In eqn (15), N(x) is the number of atoms in the bin at position x, N p is the number of neighbours of atom i, N b is the number of bonds containing atom i, N a is the number of angles containing atom i. r 1 and r 2 represent the positions of the pairwise interacting two atoms. F 1 and F 2 stand for the forces acting on the two atoms as the result of the pairwise interaction. The first term on the right-hand side of eqn (15) is a pairwise energy contribution from all N p neighbours. The second term is a bond contribution from all N b bonds. The third term is an angle contribution from all N a angles. The fourth term is the KSpace contribution which is as the result of the long-range coulombic interactions.
Based on eqn (2), the mechano-caloric system, Fig. 2(b), subjected to pressure gradient is used to calculate the microscopic mechano-caloric coefficient m 21 by where, the heat flux is calculated from the velocity profile and the excess specific enthalpy by The local specific enthalpy is found by 21 where, u(z) is the specific internal energy ( per unit volume) of interfacial water in the bin at position z, and p xx (z) is the pressure of interfacial water along x-direction in the bin at position z, which are calculated by eqn (14). The excess specific enthalpy δh(z) required in eqn (17) is obtained by where, h b is the specific enthalpy of the bulk water. When the channel width is large enough, h b can be obtained by 20,21,42 Graphene pistons Carbon 12.01115 3.39848 0.29288 0.00000 The fluid velocity profile of the planar Poiseuille flow, which is generated by a pressure gradient ∇p in the nanochannel can be theoretically computed by solving the Stokes where, L z is the channel width, η is the fluid viscosity, and b is the slip length (see Fig. 1). Fitting eqn (21) to the numerically obtained velocity profiles provides the hydrodynamic parameters η and b. These values, together with the calculated excess specific enthalpy are used in eqn (12) to obtain the macroscopic value of the thermo-osmotic and the mechanocaloric coefficients M 12 = M 21 . The pressure gradient −∇p in eqn (16) is determined from the applied body forces on the water molecules f i by the expression where, N is the total number of the interfacial water molecules, and V is the nanochannel volume.

Results and discussion
Due to the high computational cost of the thermo-osmotic system, Fig. 2(a), we use mechano-caloric system, Fig. 2(b), in our studies. Therefore, in section 3.1 we first justify the equivalence of the two systems in the sense that they yield equal reciprocal coefficients obeying Onsager relations. Then we use the mechano-caloric system to study the nanochannel size effects on thermo-osmotic coefficient and the mechanism of thermoosmosis in section 3.2 and section 3.3, respectively.

Equivalence of mechano-caloric and thermo-osmotic systems
Details of calculations with the thermo-osmotic system are given in Appendix A. In summary, it has been found by molecular dynamics simulations that m 11  To demonstrate the equivalence of the two systems, the mechano-caloric system shown in Fig. 2(b) with channel width L z = 2.0 nm is used here to calculate the mechano-caloric coefficient m 21 . The system is subjected to the following simulation steps. The water and quartz slabs are maintained at 300 K. The system is firstly equilibrated for 1.2 ns. External particle forces f i are then exerted to water molecules in the x-direction to produce a Poiseuille flow across the channel length. A constant acceleration a 0 is applied to all oxygen (O) and hydrogen (H) atoms via f i . 67 The pressure gradient is calculated by eqn (22). The hydrodynamic resistance may increase as the channel size decreases according to the previous study. 43 Therefore, the pressure gradient was tested for different channel widths to ensure that a steady flow was generated. After 0.4 ns, the fluid flow becomes stable, and the numerical velocity profile is recorded every 0.001 ns for 0.8 ns afterwards. The calculations use a similar approach to Appendix A to obtain the heat flux from the excess enthalpy density and velocity profiles. The profiles are not shown here as they will be presented in the next sub-section. The calculated microscopic mechano-caloric coefficient is m 21 = (1.2 ± 0.06) × 10 −8 m 2 s −1 , which is very close to the microscopic coefficients obtained with the thermo-osmotic system. This provides further support for the validity of the Onsager reciprocal relations at molecular level, and more importantly allows us to use the computationally simpler mechano-caloric system for studying size effects on the microscopic coefficients controlling thermo-osmosis.

Nanochannel size effect on the thermo-osmotic coefficient
Previous studies suggest that the channel size can change the liquid structure of the interfacial fluid, causing a potential significant alteration to the transport properties of the channel. 31,43 This section will present microscopic evidence for the effects of channel size L z on thermo-osmotic coefficient. The mechano-caloric system, Fig. 2(b), with several channel sizes (widths) ranging from 0.8 to 3 nm (L z = 0.8, 1, 1.5, 2, 2.5, 3 nm) is analysed to calculate the corresponding microscopic mechano-caloric coefficients, which are equal to the microscopic thermo-osmotic coefficients by the established validity of the Onsager reciprocal relation at the molecular level.
The liquid structure of the interfacial water can be effectively analysed by the spatial arrangement and the solvation (hydration) mode of the interfacial water molecules. 30 We follow a previous investigation on silicon-water interfaces 68 to study the spatial arrangement of the interfacial water molecules by looking at the mass density and the charge distribution profiles of the interfacial water for different channel sizes.
Firstly, the calculated mass density distribution profiles, ρ m (z), are shown in Fig. 5. The results highlight significant mass density fluctuations/variations near the quartz-water interfaces and approximately constant density in the midsection of the nanochannel. The constant mid-section density is found to be equal to the density of bulk water: 1000 kg m −3 . As discussed in previous studies, 41,68,69 these variations are related to interactions between the interfacial water molecules and the quartz surfaces. Analogous interfacial fluid density variations at other solid-fluid interfaces have been reported previously. 33,41,63,70 These variations indicate that in the vicinity of the quartz-water interface, the interfacial water molecules tend to occupy some fixed layers which are close and parallel to the quartz surfaces. 41,68 The peaks (maxima) in the mass density curves show clearly the locations of these layers.

Paper Nanoscale
The interval of adjacent peaks stems from the strong Lennard-Jones repulsion of oxygen (O) atoms along z-direction in adjacent water layers. 68 In contrast, the troughs (minima) of the mass density curves show that the interfacial water molecules have the lowest distribution probability at these interlayer regions. Fig. 6 presents the charge distribution profiles, Q(z). The results show that the spatial/molecular orientation of the interfacial water molecules is simultaneously under the strong influence in the vicinity of the quartz-water interfaces, as also shown by Markesteijn et al. 68 for silicon-water interfaces. The first peak (maximum) in charge distribution profiles shows that larger number of hydrogen (H) atoms instead of oxygen (O) atoms are located in the vicinity of the quartzwater interface because the hydrogen (H) atom carries a charge of +0.41e while the oxygen (O) atom carries a charge of −0.82e in the ClayFF force field. 60 Similarly, the first trough (minimum) in Fig. 6 describes a layer consisting mainly of oxygen (O) atoms in place of hydrogen (H) atoms. The position z of the first trough (minimum) in the charge profiles (Fig. 6) is consistent with that of the first peak in the mass density curve (Fig. 5). In addition, the second minimum and maximum indicate a similar trend for the peak and trough in the mass density profile. The charge distribution profiles are approximately zero in the middle of the channel if channel size is large enough (i.e. L z ≥ 1.5 nm). This shows a random spatial distribution and orientation (disordered configuration) of the interfacial water molecules, 68 which behave as bulk (free) water. 30 The results presented in Fig. 5 and 6 show the formation of a thin region in the vicinity of the quartz-water interface where the liquid structure of the interfacial water is significantly different from the bulk (free) water in terms of distribution and orientation. This layer can be called a boundary or an interaction layer. 22 In the boundary layer, the interfacial water molecules form a layer-by-layer, highly ordered liquid structure near the solid surface, which is reflected by the alternating peaks and troughs in the mass density and charge distribution profiles. Contrary to the boundary layer, the bulk liquid in the middle of a sufficiently wide nanochannel is composed of randomly distributed liquid molecules, which is reflected by the absence of large fluctuations ( peaks and troughs) in the corresponding distribution profiles for channels wider than 1.5 nm. Therefore, the boundary layer can be estimated to extend to the position where the profile fluctuations become negligible compared to the first peaks and troughs. Based on the results presented in Fig. 5 and 6, the thickness of the boundary layer is estimated to be ∼0.75 nm for all simulated cases. Notably, when the channel width is less than two times the boundary layer thickness, i.e. L z < 1.5 nm, the boundary layers at the two quartz-water interfaces overlap and the middle of the channel does not contain bulk water.
One potential difference between water behaviour in channels of different widths could arise from differences in water solvation mode. The variation of the solvation mode with channel size was obtained by calculating the spatial distribution function (SDF) of the interfacial water molecules. Details of SDF can be found in Chen at al. 30 The SDF provides an alternative method for analysis of the liquid structure, to the frequently-used radial distribution function (RDF), 71 specifically suitable for 3-D representations. The SDF describes statistically the relative positions between any central water molecule and neighbouring water molecules. The positions of neighbouring water molecules are determined by the molecular mass centres. 30,63 The spatial probability is calculated by averaging over the solvation modes of all water molecules within a searching radius for neighbouring water molecules. Fig. 7 provides a visualisation of SDF with the searching radius r = 0.5 nm.
In more detail, Fig. 7(a) shows a 3-D plot of the probable positions, calculated by SDF, of the water molecules surrounding a central water molecule. Three most favourable positions of the adjacent water molecules are depicted as I-III. Molecules in positions II and III belong to the first solvation shell and are most strongly bonded. 72 These are less affected by the channel size. Molecules in the position I belong to the second solvation shell. 72 The projections on 2-D planes, shown in Fig. 7(b) and (c), indicate that the molecules in the channel with L z = 1.0 nm have different arrangement from the larger channels. Firstly, molecules in position I are missing, and the  characteristics of positions IV and V are different from the other cases. This is an indication that a unique liquid structure is formed for this specific channel width. It is anticipated that a distinctive value of the thermo-osmotic coefficient for the channel width L z = 1.0 nm will be obtained due to possibly distinct thermodynamic properties and hydrodynamic response of interfacial water. Note that a persistent solvation mode is created for interfacial water molecules when L z ≥ 1.5 nm, indicating that the effect of channel size on liquid structure diminishes for larger channels.
The excess enthalpy distribution profiles of interfacial water for different channel sizes are calculated by eqn (19) and shown in Fig. 8. The presence and sizes of the boundary layers and the bulk water section can be clearly identified and coincide with those in Fig. 5 and 6. The excess enthalpy profiles vary significantly in the boundary layers and are zero in the bulk water, when present. In smaller channel sizes, e.g. L z < 1.5 nm, the interference of the overlapping boundary layers results in a complex variation of the excess enthalpy. Interestingly, the first trough is found before the first peak, which differs from some previous studies 20, 42 reporting excess enthalpy distribution profiles. The difference arises from different substrates used. In our case, there are dangling hydroxyl groups on the quartz surfaces, seen in Fig. 3, which create spaces for some liquid water molecules to penetrate into the quartz. 69 The penetrating water molecules are strongly attracted and trapped within these spaces, which can be seen in Fig. 8 for the example of L z = 2.0 nm, where the vertical dotted lines representing H position are inside the physical quartz boundary. This results in negative potential energy and according to eqn (18) in negative excess enthalpy density, which explains the appearance of trough before the first peak.
We note that the positions of the first troughs in mass density profiles are different from those of the first peaks in charge density and excess enthalpy distribution profiles. This can be explained by the following observation. The mass density distribution profiles, ρ m (z), shown in  distribution profiles in Fig. 25 of Appendix C, which is also found to have the same coincidence.
In addition to the excess enthalpy, the existence of interfacial regions may significantly affect the hydrodynamic properties of the interfacial water. 44 Fig. 9 shows the velocity profiles of the interfacial water molecules along the x-direction, v x (z), for different channel sizes. It is seen that the velocity distribution profiles vanish at the quartz-water interfaces and have maxima at the middle of the channel. The fitting curves shown in the plots are based on eqn (21) and are used to determine the hydrodynamics parameters η and b. The curves are fitting the simulation results with high correlation, indicating that the induced pressure-driven flow is described well by planar Poiseuille flow, and the continuum hydrodynamics remains valid in nanometric/sub-nanometric channels. The fitting results for b indicate that there is basically no slippage at the quartz-water interfaces (b ≈ 0 nm), due to the large friction coefficient κ between interfacial water and quartz surfaces (κ ≈ 10 11 N s m −3 from the fitting results). The fitting results for η show that for L z ≥ 1.5 nm, η gradually reduces to a constant value (∼0.60 mPa s) which is expected to further converge to the shear viscosity of SPC bulk water at 300 K and 1 atm (∼0.475 mPa s) 73 when the channel width continues to increase, while for L z < 1.5 nm, the overlap and interference of the two boundary layers leads to a higher value and a complex variation with channel size of water viscosity.
The microscopic mechano-caloric (or thermo-osmotic) coefficient m 21 for different channel sizes L z is calculated by eqn (16). The corresponding macroscopic coefficient M 21 is also calculated by eqn (12), using the calculated excess enthalpy density and fitted viscosity and slip length. For each nanochannel size, we performed three independent simulations with different initial configurations of interfacial water molecules, from which we computed the error bars on m 21 and M 21 . Fig. 10 shows the variation of the thermo-osmotic coefficient with channel size L z .
The results present new findings and observations: (1) The microscopic thermo-osmotic coefficients measured by MD simulations (m 21 ) and the macroscopic coefficients calculated by thermodynamics (M 21 ) are approximately the same for all channel sizes, with the exception of the smallest.

Paper Nanoscale
(2) For channels with L z ≥ 1.5 nm, m 21 is approximately constant, while for L z ≥ 1.5 nm, m 21 varies with the channel size. This is due to the overlap and interference of the boundary layers which affect the excess enthalpy distribution and hydrodynamic properties of interfacial water, causing variations of the thermo-osmotic coefficient.
(3) A unique value for thermo-osmotic coefficient was found for the case of L z = 1.0 nm, and shown to be associated with the unique liquid structure (i.e. solvation mode) observed in Fig. 7.
(4) The MD simulations with the mechano-caloric system have revealed the dependence of the thermo-osmotic coefficient on the channel size. It has been found that this dependence emerges from the changing liquid structure in boundary layers adjacent to the solid surfaces, which alters either the excess enthalpy distribution or the hydrodynamic properties of the interfacial water. The comparisons between number/mass density, charge, excess enthalpy, and pressure distribution profiles, as well as the coincidence of the corresponding peak/ trough values (Fig. 5,6,8,24 and 25) indicate that the excess enthalpy distribution profile is dictated by the layered liquid structure of interfacial water. The velocity distribution profiles and the fitting results of fluid viscosity in Fig. 9 indicate that the dynamic properties of the interfacial fluid are also strongly affected by the interfacial liquid structure. Specifically, with channel size increasing, the fluid viscosity firstly increases and then decreases to a steady value. This is in agreements with the observed variation of the thermo-osmotic coefficient, since according to eqn (12) and (13) the thermo-osmotic coefficient is inversely proportional to the fluid viscosity.

Mechanism of thermo-osmosis
This section provides a mechanistic demonstration that the thermo-osmotic flow of interfacial water is driven by a pressure gradient induced by a temperature gradient. Though Ganti et al. 21,74 have revealed and confirmed the similar mechanisms for conceptual nonpolar Lennard-Jones system, this mechanism has not been found in a more complex and realistic system (e.g., quartz-water interfaces). Therefore, we followed the same method and used a more advanced MD configuration with a dynamically self-adjusted piston to achieve a con-stant bulk fluid pressure and revisit this mechanism in a realistic quartz-water system.
We use the mechano-caloric system, Fig. 2(b), with the channel width L z = 2.0 nm. This is subjected to zero external pressure gradient ∇p = 0 along the channel (x direction). A series of MD simulations with a range of temperatures from 270 K to 360 K are performed. The corresponding pressure distributions across the channel section are obtained by eqn (14). For each temperature, the simulation system is relaxed/equilibrated first for over 1.6 ns. The gap between the quartz surfaces (channel size) is adjusted dynamically by the graphene piston to ensure the bulk pressure of interfacial water remains at the designated value (P piston = 1 atm ≈ 0.0001 GPa). The measurement of pressure is then followed for 1.2 ns. Fig. 11 presents the number density profiles ρ(z) in the channel obtained for all temperatures studied. In the figure, z = 0 nm represents the middle point of the channel, and z ≈ 1.0 nm represents the quartz surface. It is found that the increasing temperature leads to increasing nanochannel gap due to the thermal expansion of the interfacial water, as the enclosed pressure applied by graphene piston is kept constant at 1 atm. Fig. 12 presents the temperature distributions of the interfacial water along z-direction T (z) at different temperatures, indicating that no temperature gradient exists in this direction. Fig. 13 presents the pressure tensor components of the interfacial water in the channel for the case of T = 300 K (as an example). It can be observed that in the middle of the nanochannel, the pressure components are equal, i.e. p xx = p yy = p zz = p bulk ≈ 0 GPa. However, near the quartz surfaces, the pressure components vary several orders of magnitude from the bulk. The component p xx is approximately equal to p yy over the entire channel section, indicating insignificant anisotropy on the plane parallel to the quartz surfaces. However, p zz deviates considerably from p xx and p yy near the quartz surfaces. The boundary layer discussed in section 3.2 can be identified in Fig. 13 as the region where the influences due to the quartz-water interactions are substantial on the properties of interfacial water.  all temperatures, p xx in the bulk is within 1% of the designated value of 1 atm, and deviates significantly from this value as the quartz surfaces are approached. The pressure starts oscillating with a wavelength corresponding roughly to the size of a single water molecule (∼0.28 nm) 41,63 and its amplitude increases with decreasing temperature. The latter observation may be related to a higher bulk density of fluid at a lower temperature, Fig. 11, and subsequently, a stronger effect of solid on the thermal mobility according to a previous study. 41 The pressure profiles in Fig. 14 suggest that a temperature gradient will not produce a pressure gradient in the x-direction in the bulk water, since the pressure is nearly independent of temperature away from the quartz surfaces (z = 0). The bulk water is in mechanical equilibrium. In the vicinity of the quartz-water interface, however, a finite pressure gradient along the channel will be generated by changing the temperature, seen by different pressure magnitudes at different temperatures in Fig. 14. Such pressure differential can be seen as the main driving force for the flow in the x-direction.
The pressure gradient can be obtained by eqn (8). Fig. 15 presents the pressure gradient of the interfacial water @p xx ðzÞ @x across the channel width under different temperatures when a unit temperature gradient is applied in the x-direction   Paper Nanoscale (∇T = 1 K m −1 . Under the same temperature gradient ∇T, the magnitude of the pressure gradient will decrease as the water temperature increases. The considerations above assumed that the interfacial water was a continuum, i.e. thermally induced pressure on a fluid element. The body force on each water molecules can be calculated microscopically. The body force on a fluid element at a distance z from the quartz surface is given by and the thermo-osmotic force per particle is given by Fig. 16 presents the distribution of particle forces across the channel width under different temperatures when a unit temperature gradient is applied in the x-direction (∇T = 1 K m −1 . It is seen that only in the boundary layer, the water molecules will experience a body force to drive the flow along the x-direction. In the bulk region, no particle force will act on the water molecules. In addition, under the same temperature gradient ∇T, the magnitude of the body force will decrease as the water temperature increases, which is in accord with the result in Fig. 14 and 15. By investigating the mechanical interactions of the interfacial water under different temperatures using a continuum approach (pressure tensor distribution) and a discrete approach (body force per molecule), we have revealed and confirmed the microscopic mechanism of thermo-osmosis of water in quartz nanochannels. It has been found that only in the boundary layers adjacent to the solid surfaces, the fluid molecules experience a driving force that produces a thermo-osmotic flow.

Conclusions
The thermo-osmotic flow of water in quartz slit nanochannels has been investigated by MD simulations. The microscopic thermo-osmotic coefficient has been measured by analysis of mechano-caloric and thermo-osmotic systems. The following conclusions are drawn from the results and discussion of the work: • The validity of Onsager reciprocal relation and the continuum hydrodynamics description of flow at the nanoscale has been confirmed, and the consistency between the microscopic and the macroscopic/thermodynamic thermo-osmotic coefficient has been demonstrated; • The channel size effect on the thermo-osmotic coefficient has been quantified by analysis of systems with different channel widths, ranging from 0.8 to 3.0 nm. It has been found that the liquid structure of the interfacial water changes with Fig. 15 The pressure gradient of the interfacial water across the channel width for a unit temperature gradient under different temperatures. decreasing channel size, due to the overlap and interference of liquid boundary layers adjacent to the solid surfaces. This leads to changes in the excess enthalpy distribution and in the hydrodynamic properties of interfacial water, causing variations of the thermo-osmotic coefficient; • The pressure tensor distributions of the interfacial water under different temperatures have been quantified. It has been found that the temperature gradient inside the nanochannel induces a pressure gradient parallel to the quartz surface. This pressure gradient is identified as the source of thermo-osmotic motion of the interfacial water. It has been found that, although the bulk water is at mechanical equilibrium under a thermal gradient, such an equilibrium state cannot be achieved in the vicinity of the quartz-water interfaces. Layers of water at these surfaces experience a pressure gradient parallel to the surfaces, which drives the thermo-osmotic flow.
The modelling and simulation approaches used in this work, as well as the insights provided by the results, form a good foundation for further investigation of coupled transport phenomena in nanoscale structures.

Author contributions
The manuscript was written through the contributions of all authors. All authors have given approval to the final version of the manuscript.

Conflicts of interest
There are no conflicts of interest to declare.

Appendix A. The thermo-osmotic system
For the studied thermo-osmotic system, the water flux is calculated by where, ΔN is the number of water molecules transported through the nanochannel within a time interval Δt, ρ is the number density of interfacial water molecules, and A is the cross-section area of the nanochannel (A = L y L z ). This can be further simplified as where, N = ρL x L y L z is the total water molecule number in the quartz nanochannel.
The thermo-osmotic system subjected to temperature gradient is used to calculate the microscopic thermo-osmotic coefficient by where, ∇T is the temperature gradient across the domain, T is the average temperature in the nanochannel, and J f is the water flux calculated by eqn (26). The thermo-osmotic system subjected to pressure gradient is used to calculate the microscopic coefficient m 11 by: where, ∇p is the pressure gradient along the nanochannel, and J f is the water flux calculated by eqn (26). The pressure gradient is obtained by calculating the pressure tensor of interfacial water using the atom-based virial equation, i.e. eqn (14) in the main text.
A.1. Determination of m 12 using the thermo-osmotic system with a temperature gradient The system shown in Fig. 2(a) of the main text with the channel width L z = 2.0 nm is subjected to the following simulation steps. The two water reservoirs are kept at the same pressure ( p 1 = p 2 = 1.0 atm) applied by two graphene pistons.
The system is first relaxed/equilibrated for 1.2 ns by maintaining both the water and quartz at an initial temperature (T 1 = T 2 = T 0 = 300 K). The temperature gradient in the nanochannel is then applied by increasing the temperature of reservoir 1 to 330 K (T 1 = 330 K) and decreasing the temperature of reservoir 2 to 270 K (T 2 = 270 K). The system is then stabilized for another 3.0 ns until the temperature profile reached steadystate and a relatively stable flow is achieved under the constant temperature gradient created by the two reservoirs. The calculation of the microscopic thermo-osmotic coefficient, eqn (27), requires the temperature gradient and the water flux. Fig. 17 shows the temperature profile of both water and quartz along the channel from the MD simulation. The results show that a constant temperature gradient ∇T is created along the quartz slit nanochannel. The quartz temperature profile is consistent with the interfacial water temperature profile. Linear regression of the water temperature profile inside the nanochannel gives −∇T = 4.815 K nm −1 . Note that the range of water temperature in the system is of ∼260-340 K, previous studies 75, 76 have indicated that SPC water is held in the stable Fig. 17 The temperature profile in the thermo-osmotic system.

Paper
Nanoscale liquid phase in this temperature range and atmospheric pressure. The calculation of water flux by eqn (26) requires the water molecule numbers transported through the nanochannel in a time interval. Fig. 18(a) presents the changes of water molecule numbers and temperatures with time in the reservoir 1 (ΔN 1 ), reservoir 2 (ΔN 2 ) and in the channel (ΔN 0 ). It can be observed that the temperature gradient between the two reservoirs drives the water molecules from the hot reservoir to the cold reservoir at a constant rate, demonstrating the establishment of a thermo-osmosis process in the system. The change of water molecule numbers, ΔN, required for eqn (26), is obtained as an average of the changes in the two reservoirs, i.e. ΔN = (|ΔN 1 | + |ΔN 2 |)/2, and the water flux is calculated to be J f = 0.191 m s −1 .
The microscopic thermo-osmotic coefficient is calculated by eqn (27) to be m 12 = 1.19 × 10 −8 m 2 s −1 . The (nearly) linear flow profile in Fig. 18(a) indicates that the thermo-osmotic coefficient m 12 is a constant, and the inequality m 12 > 0 represents a flow from the hot reservoir to the cold reservoir.
We extended our simulations to more scenarios of temperature gradients and calculated the corresponding microscopic thermo-osmotic coefficients. The results are shown in Table 2.
A.2. Determination of m 21 using the thermo-osmotic system with a pressure gradient The system shown in Fig. 2(a) of the main text with the channel width L z = 2.0 nm is subjected to the following simu-lation steps. The temperature in the two reservoirs and the nanochannel is set at T 1 = T 2 = T 0 = 300 K and the pressure on the two graphene pistons is set at p 1 = p 2 = 1.0 atm ≈ 1.01 × 10 −4 GPa. The system is initially relaxed/equilibrated for 1.2 ns, after that the temperature is kept constant everywhere (∇T = 0) and a pressure gradient is generated by increasing p 1 = 3001 atm ≈ 0.30 GPa, and maintaining p 2 = 1 atm. Fig. 19(a) shows the calculated changes of water molecule numbers in reservoir 1 (ΔN 1 ), reservoir 2 (ΔN 2 ) and in the channel (ΔN 0 ), as well as the pressure acting on the two pistons with time. As in sub-section A.1, the water flux is calculated from the average change of water molecule numbers and eqn (26).
The calculation of the microscopic coefficient m 11 by eqn (28) requires the pressure gradient along the nanochannel. The pressure tensor is calculated by eqn (14), and the pressure profile along the x-direction is shown in Fig. 20. The profile indicates a pressure drop induced by the viscous entrance effects. 71 Therefore, the generated −∇p = 0.042 GPa nm −1 is   The calculation of the microscopic mechano-caloric coefficient m 21 by eqn (16) requires the heat flux in addition to the pressure gradient. This is calculated by eqn (17) using the excess enthalpy density and the velocity profiles. Fig. 21 presents the excess enthalpy density profile over the channel section calculated by eqn (19). Fig. 22 presents the velocity profile of the interfacial water molecules along the x-direction v x (z). The hydrodynamics fitting curve based on eqn (21) is also shown. The fitting yields b ≈ 0 and η ≈ 0.61 mPa s. Using the calculated heat flux and pressure gradient in eqn (16) provides the microscopic mechano-caloric coefficient m 21 = 1.22 × 10 −8 m 2 s −1 .
We extended our simulations to more scenarios of pressure gradients and calculated the corresponding microscopic mechano-caloric coefficients. The results are shown in Table 3. They can be summarised by the following expression for the mechano-caloric coefficient m 21 = (1.21 ± 0.04) × 10 −8 m 2 s −1 . We note that the calculated microscopic thermo-osmotic coefficient in section A.1. is m 12 = (1.19 ± 0.03) × 10 −8 m 2 s −1 , which shows that the Onsager reciprocal relation (m 12 = m 21 ) is obeyed microscopically. Furthermore, using the excess enthalpy profile and the calculated slip length and viscosity in eqn (12) the value of the macroscopic thermo-osmotic (and mechano-caloric) coefficient is obtained as M 12 = M 21 ≈ (1.21 ± 0.02) × 10 −8 m 2 s −1 , which is consistent with the microscopic values.
By conducting the MD simulations of a thermo-osmotic system, we have shown how the thermo-osmosis emerges from the collective behaviour of microscopic constituents, allowing for the calculation of the microscopic coefficients m 11 , m 12 , m 21 . Further, we have shown that the calculated mixed coefficients obey the Onsager reciprocal relation (m 12 = m 21 ) and that they can be considered as emergent properties of the system since they are very close to the value based on irreversible thermodynamics (m 12 = m 21 = M 12 = M 21 ).

Appendix B. Effect of bin size on density profiles and their peak/trough values
The choice of bin size will influence the accuracy of the derived first peak value (z 0 , ρ 0 ): the smaller the bin size, the more accurate profiles and peak values will be obtained, as the error range can be narrowed down by decreasing the bin size. For example, for the position of the peak z 0 , a bin size of 0.04 nm causes an error range −0.04 nm-0.04 nm. However, a very small bin size can cause significant fluctuations of profiles as there is not enough data in each bin to be averaged for producing smooth profiles. To assess whether the peak value will change as the bin size changes, we take the mass density profiles for the case L z = 1.0 nm as an example. By varying the bin size, the corresponding profiles and peak values are   The pressure profile along x-direction ( p xx (x)). Table 3 Sets of pressures adopted to create pressure gradient inside the quartz nanochannels of the thermo-osmotic system, as well as the calculated microscopic and macroscopic mechano-caloric coefficient   derived and shown in Fig. 23. When the bin sizes are too large, e.g. 0.040-0.100 nm, the produced profiles are quite rough and the peak values cannot be accurately determined. When the bin size is too small, e.g. 0.002 nm, the lack of enough data in each bin results in fluctuating profiles. With intermediate bin sizes, e.g. 0.004-0.010 nm, the profile fluctuations are avoided with acceptable accuracy of peak values. Therefore, in our study, a bin size of 0.004 nm is adopted.
Appendix C. Number density and pressure distribution profiles for different channel size The mass density distribution profiles, ρ m (z), presented in Fig. 5 of the paper are disproportionally affected by the distribution of oxygen (O) atoms, because of their greater molar mass (∼16 g mol −1 ) compared with hydrogen (H) atoms (∼1 g mol −1 ). We present the number density distribution profiles, ρ(z), of O and H atoms in Fig. 24, with the pressure tensor distribution profiles, ρ xx (z), in Fig. 25, and mass density, charge density, excess enthalpy distribution profiles in Fig. 5, 6 and 8, respectively to clarify the coincidence of peak/trough values in different distribution profiles. The peaks and troughs in number density profiles show good agreements with other distribution profiles. It can be seen that: (1) the position of the first peak in the number density profile of O-atoms is the same as the position of the first peak in the mass density profile, the position of the first trough in the charge profile, and the position of the first peak in the pressure profile; (2) the position of the first peak in the number density profile of H-atoms is the same as the position of the first peak in the charge profile, and the position of the first trough in the pressure profile; (3) the position of the first peak in the excess enthalpy profile is between the first peak position of the H-atom number density profile and the first peak position of the O-atom number density profile, approximately coinciding with the first peak in O number density distributions. Other peaks/troughs on these profiles coincide with each other in the same way.