Molecular Dynamics Investigation of Pressure-Driven Water Transport in Kaolinite Nanopores

1 Introduction

Clay is one of the most widespread sediments on Earth. It is extensively distributed in soils, rocks, and sedimentary layers, and plays a crucial role in environmental, geological, and petroleum engineering [1,2]. As a natural layered silicate mineral composed of alternating stacks of silica tetrahedral sheets and alumina octahedral sheets [3], the microscopic pore structure and surface properties of clay exert a decisive influence on the occurrence state and flow behavior of water molecules within it. These characteristics not only govern the migration and distribution of fluids in pores, but also directly affect the permeability and mechanical response of geomaterials. For example, in hydraulic fracturing operations for shale gas development, the behavior of pore water in clay layers influences fracture propagation and reservoir seepage efficiency [4]; in subsurface contaminant barriers and nuclear waste disposal projects, the adsorption and seepage properties of clay materials determine the long-term reliability of barrier systems [5]. Therefore, revealing the microscopic occurrence characteristics and flow mechanisms of water molecules in clay nanopores is of great scientific and practical significance for understanding seepage laws in low-permeability media and guiding the design of related engineering projects [6,7].

Traditional geotechnical seepage theories are primarily based on macroscopic experiments and the continuum hypothesis. Since it was proposed in 1856, Darcy’s law has become the fundamental theory for describing fluid flow in porous media [8], which is widely applied in macroscopic engineering fields such as enhanced oil recovery [9], heterogeneous porous media flow [10], and geothermal energy exploitation [11]. It assumes steady, incompressible flow in a homogeneous and isotropic medium. However, in nanoscale pores, enhanced interfacial effects and slip behavior induce flow characteristics that deviate significantly from those predicted by classical continuum hydrodynamics [12,13,14,15]. This implies that seepage is not only influenced by porosity and permeability but is also significantly governed by solid-liquid interfacial forces, surface adsorption, and variations in fluid properties [16,17]. The microstructure of clay minerals significantly influences their macroscopic geomechanical properties. The conventional Darcy model neglects these microscopic effects and therefore fails to accurately capture the complex non-Darcy flow characteristics in clays [18], shale oil reservoirs [19], and other complex thermohydraulic coupled processes [20].

To address this limitation, various modified models have been proposed [21,22,23]. However, experimental validation at the nanometer scale remains challenging. The requirements in measurements of ultra-low flow velocities, precise control of hydraulic gradients, and mitigation of sample disturbance, end effects, and interfacial leakage severely limit their accuracy and number [24]. Given these constraints, molecular dynamics (MD) simulation has attracted increasing attention in recent years. As a microscopic simulation method based on Newtonian mechanics and statistical physics, MD can resolve fluid-solid interfacial interactions and structural evolution at the atomic scale and has been widely applied to study the structure and flow characteristics of fluids under nanoscale confinement [25,26,27].

With respect to interfacial wettability and adsorption structure, Šolc et al. [28] employed MD simulations to reveal pronounced wettability differences between the two kaolinite (001) basal surfaces. Their results showed that the aluminum-oxygen octahedral surface is strongly hydrophilic, whereas the silicon-oxygen tetrahedral surface exhibits hydrophobic characteristics. Building on these findings, Liu et al. [29] further combined MD simulations with the Young-Dupré equation to quantitatively characterize the wettability of kaolinite surfaces. Zhang et al. [30] combined MD and grand canonical Monte Carlo simulations to systematically analyze carbon dioxide adsorption in kaolinite nanopores under varying pore sizes and water contents, revealing the influence of confinement and hydration on adsorption capacity and spatial distribution. Subsequently, Zhang and Tang [31] investigated the density distribution and multilayer adsorption mechanisms of methane in silica nanopores. In addition, Lu et al. [32] explored the mechanical softening and structural degradation of kaolinite under hydrated conditions from a molecular perspective.

Based on the heterogeneity of interfacial wettability and adsorption structure, the molecular layering of water near pore walls plays a critical role in momentum transfer and boundary shear, and consequently exert a strong influence on flow mechanisms in nanoporous media. In this context, Boţan et al. [33] provided the first systematic characterization of water structure, diffusion, and transport in clay nanopores, and evaluated the applicability of continuum hydrodynamics at the nanoscale. Subsequently, Wu et al. [34] incorporated wettability effects by combining MD simulations with experimental observations and proposed a nonlinear correction model based on an effective slip length to improve macroscopic flow predictions. Zhan et al. [35] further compared flow boundary conditions across several representative clay nanopore systems, significantly enhancing the accuracy of microscale water flow predictions. Wei et al. [36] employed nonequilibrium MD simulations to analyze the coupling between boundary slip, viscosity evolution, and water-clay interactions from a molecular mechanics perspective, providing deeper mechanistic insight into changes in apparent viscosity and boundary behavior under nanoconfinement. As research increasingly approaches realistic geological conditions, growing attention has been paid to the effects of pore geometry complexity and microstructural evolution on flow behavior. Kang et al. [37] constructed flocculated kaolinite microstructures that more closely resemble natural systems and demonstrated that water adsorption, diffusion, and pore-scale flow can drive transitions from flocculated to dispersed states, suggesting a dynamic feedback between seepage processes and pore-structure evolution. Bi et al. [38] further showed that pore-wall roughness and geometric configuration exert a significant influence on nanoscale two-phase water-gas flow behavior, indicating that boundary conditions and effective transport parameters may vary spatially in complex pore networks.

Overall, existing studies have systematically investigated the structure and seepage behavior of water in clay nanopores from the perspectives of interfacial wettability contrast, nanoconfined fluid occurrence, molecular layering, and boundary-condition modification. However, most previous work has focused on macroscopic or effective flow responses, while the molecular-scale origins of slip flow and non-Darcy behavior remain insufficiently elucidated. Current models still struggle to explicitly capture how dynamic processes such as hydrogen bond network restructuring and molecular orientational evolution regulate boundary slip and seepage behavior, and how these microscopic mechanisms are intrinsically linked to macroscopic seepage parameters such as effective permeability and threshold hydraulic gradients.

To fill these gaps, this study focuses on kaolinite and constructs three parallel-plate kaolinite nanochannel models. Section 2 describes the model construction procedure and the parameter settings of the MD simulations. Section 3 systematically analyzes the density distribution, hydrogen bond structure, and molecular orientation of water under different surface-wettability conditions. Section 4 investigates the flow behavior of water in kaolinite pores and elucidates the influence mechanisms of surface wettability and interfacial effects on flow regimes. The results provide microscopic theoretical support for understanding the seepage mechanisms of water in clay-like media.

2 Models and Simulation Methods

2.1 Model Construction

Kaolinite is a typical 1:1 layered silicate mineral composed of one alumina octahedral sheet and one silica tetrahedral sheet interconnected through bridging oxygen atoms [Fig. 1a]. Its chemical formula is Al2Si2O5(OH)4, and the crystallographic parameters determined from experiments are: a = 5.15 Å, b = 8.93 Å, c = 7.38 Å; α = 91.9°, β = 105.0°, and γ = 89.8° [39]. Based on these structural parameters, a kaolinite supercell was constructed by replicating the unit cell 12, 16, and 2 times along the x, y, and z directions, forming a layered kaolinite slab [Fig. 1b]. Two basal surfaces of the kaolinite slab are exposed along the z direction:

• •(001) surface: the alumina octahedral surface, which is strongly hydrophilic.
• •(001¯) surface: the silica tetrahedral surface, which is hydrophobic.

Water molecules were represented using the Simple Point Charge (SPC) model [Fig. 1c], a rigid three-site configuration with fixed bond lengths and angles that balances computational efficiency with acceptable accuracy for clay-water systems [40].

Figure 1: Molecular structures and seepage channel models: (a) Kaolinite crystal structure; (b) Surface structure of kaolinite; (c) SPC water molecule; (d) Seepage channel structure of the (001)-(001) model.

To investigate the effects of wettability on water occurrence and flow, three slit-pore models were constructed:

• 1.(001)-(001) model: hydrophilic pore with both walls exhibiting (001) surfaces.
• 2.(001¯)-(001¯) model: hydrophobic pore with both walls exhibiting (001¯) surfaces.
• 3.(001)-(001¯) asymmetric model: one hydrophilic and one hydrophobic wall.

A vacuum region was inserted between kaolinite layers to form the slit pore, which was filled with a water layer of uniform size and 60 Å thickness, with an average density of 1 g/cm3, corresponding to 17,736 water molecules [Fig. 1d]. The channel dimensions of all three models were kept identical for comparison. It is assumed that the upper and lower surfaces remain rigid in this study.

2.2 Potential Energy Parameters

All MD simulations were performed using the LAMMPS software package [41]. Interactions between kaolinite and water were described using the CLAYFF force field [42], which was specifically developed for clay minerals and accurately captures both bonded and nonbonded interactions in layered silicates. The force field primarily consists of electrostatic and Van der Waals components.

In the CLAYFF force field, the total potential energy of the system can be classified as: Etotal=Evdw+Ecoul+Ebond+Eangle(1) where Evdw is the Van der Waals potential energy, Ecoul is the Coulombic electrostatic potential energy, Ebond is the bond stretching potential energy, and Eangle is the bond angle bending potential energy.

Non-bonded interactions include Van der Waals and electrostatic forces. The Van der Waals interactions were represented using the Lennard-Jones (12-6) potential [43]. EvdW=4εijσijrij12−σijrij6(2) where εij and σij are the Lennard-Jones energy and size parameters for the atom pair i-j, and rij is the interatomic distance.

The cross-term parameters for unlike atom pairs are obtained using the Lorentz-Berthelot mixing rules [44].

The electrostatic interactions between charged particles are described by the Coulomb potential. ECoul=keqiqjrij(5) where ke is the electrostatic constant and qi and qj are the charges of particles i and j, respectively. The Van der Waals and charge parameters for all atomic species in the system are listed in Table 1.

The intramolecular bond-stretching and angle-bending potential energies are both described using harmonic oscillator models. Ebond=k1rij−r02(6) Eangle=k2θijk−θ02(7) where, k1 and k2 are the force constants for bond stretching and angle bending, respectively; r0 and θ0 denote the equilibrium bond length and equilibrium bond angle; and θijk represents the angle formed by atoms i, j, and k. The corresponding potential parameters are summarized in Table 2.

Periodic boundary conditions were applied in all three directions, and a vacuum region of 100 Å was introduced along the z-axis to eliminate spurious interactions between periodic images. A cutoff distance of 10 Å was used for Van der Waals interactions, and the integration time step was set to 1 fs. Long-range Coulomb interactions were evaluated using the Particle-Particle Particle-Mesh method (PPPM) with a target accuracy of 99.99% [45]. Intramolecular geometry in water molecules was constrained using the SHAKE algorithm [46]. All simulations were performed under the NVT ensemble with temperature controlled by the Nosé-Hoover thermostat [47,48]. In the Non-equilibrium molecular dynamics (NEMD) simulations, to avoid the artificial contribution of the imposed flow velocity to the calculated temperature, the center-of-mass velocity of water molecules in the flow direction was removed during temperature control [33].

To verify the reliability of the adopted force field parameters and the constructed models, contact angle simulations of water on kaolinite surfaces were carried out using the aforementioned force field and models. The model temperature is set at 300 K, and it contains 500 water molecules. The results show that the hydrophilic surface exhibits a contact angle close to 0° [Fig. 2a], while the hydrophobic surface exhibits a contact angle of approximately 104° [Fig. 2b]. These values are consistent with the MD simulations reported by Šolc et al. [28], demonstrating that the models constructed in this study possess high accuracy and physical validity.

Figure 2: Wettability of kaolinite surface: (a) Kaolinite hydrophilic surface contact angle; (b) Kaolinite hydrophobic surface contact angle.

2.3 Simulation Methods

Before the production simulations, energy minimization was performed using the Conjugate Gradient (CG) algorithm to eliminate any unreasonable contacts or stress concentrations in the initial configuration. After minimization, the normal pressure P t o p , which acted as the formation overburden pressure, was applied by exerting an equivalent external force on the atoms of the top surface. and the force F t o p is expressed as [38]: Ftop=−Ptop·AtopNtop(8) where A t o p is the area of the top surface, and N t o p is the number of atoms involved in the loading on the top surface. After the 1 ns loading simulation, the system reached equilibrium, and the stable pore size of the water molecules at the current temperature and pressure was determined.

Subsequently, the applied normal pressure was removed, and the kaolinite layer was fixed in place. Equilibrium molecular dynamics (EMD) simulations were then carried out to analyze the structural distribution and diffusion characteristics of the confined water molecules. The total duration of the equilibrium simulation was 10 ns, with the first 5 ns used to allow the system temperature, energy distribution, and molecular configurations to stabilize, ensuring thermal equilibrium. The trajectory data from the remaining 5 ns were used for statistical analysis.

In addition, to elucidate flow behavior, separate NEMD simulations were conducted under an applied pressure gradient. Acceleration in the y direction is applied to each atom of the water to simulate the flow of water within the kaolinite pores. The acceleration is 0.03 Å/ps2, which means the driving pressure is 42.9 MPa. Such a high-driving-force setup is commonly adopted in NEMD simulations to overcome thermal fluctuations and obtain a statistically stable directional flow within the limited simulation time, and the resulting flow velocities should be interpreted as nanoscale simulation responses rather than direct counterparts of macroscopic flow velocities. The relationship between acceleration (ay) and the driving pressure Δ P across the flow domain is given as follows [49]: ΔP=NmayAxz(9) where N is the number of particles within the pore, m is the particle mass, Axz is the cross-sectional area perpendicular to the flow direction.

The initial simulation temperature was set to 313 K (40°C), and the applied normal pressure during the loading stage was 23.5 MPa, corresponding to the high-temperature and high-pressure environment at an approximate burial depth of 1000 m in shale-gas reservoirs. These conditions represent typical subsurface environments and allow examination of the effects of temperature and pressure on the occurrence and flow characteristics of water in nanopores. Specifying non-zero initial temperature and pressure ensures that the system begins from a stable thermodynamic state, providing appropriate initial conditions for the subsequent analyses.

3 Microscopic Mechanisms of Water Occurrence and Diffusion in Kaolinite Nanopores

3.1 Density Analysis

In the MD simulations, the slit-like kaolinite nanopores were divided along the z direction into sublayers with a thickness of 0.1 Å. After the system reached equilibrium, the number density of atoms within each sublayer was calculated and time-averaged to obtain the density distribution of water molecules along the z-axis. Because the characteristic size of a water molecule is significantly larger than 0.1 Å, the hydrogen and oxygen atoms of an individual water molecule may fall into different sublayers. The density values are converted to units of g/cm3 (where 1 g/cm3 corresponds to a particle number density of 0.1002 particles/Å3). To ensure consistency in comparison, the z-coordinate origin for the (001) surface is defined at the position of the hydroxyl oxygen atoms, whereas the origin for the (00 1 ¯ ) surface is defined at the position of the silicon-oxygen atoms.

The z direction density profiles of water molecules exhibit a clear interfacial layering characteristic of confined fluids (Fig. 3). Strong density oscillations appear near the pore walls, whereas the density becomes uniform toward the center of the pore. This behavior indicates that, under the influence of solid-liquid interfacial interactions, water molecules form a multilayered structure composed of adsorbed layer, diffuse layer, and bulk region. Within approximately 1.8 Å of the pore wall, virtually no oxygen atoms are present, while hydrogen atoms can approach the solid surface more closely. A distinct interfacial density peak appears at the pore wall, after which the density gradually decays with increasing distance until reaching the stable value of the bulk region (approximately 1 g/cm3). This characteristic “oscillatory-to-flat” distribution reflects the structural ordering of water molecules induced by the electrostatic field and hydrogen bond interactions at the kaolinite surface.

Figure 3: Density distributions along the z direction for water molecules and their atomic species in the three models: (a) (001)-(001) model; (b) (00 1 ¯ )-(00 1 ¯ ) model; (c) (001)-(00 1 ¯ ) model; (d) Density distributions of water for the three models.

The adsorbed layer generally begins at the clay mineral surface and extends to the position of the first major peak in the water density distribution, with a thickness approximately equal to the effective collision diameter of a water molecule. The region beyond the adsorbed layer and up to the fourth density trough is defined as the diffuse layer, where the density oscillations gradually decay. Outside the diffuse layer, the density becomes stable and invariant with distance, forming the bulk region, which may be regarded as the domain of free or unconfined water.

For the (001)-(001) model, the first small density peak corresponds to the orientational arrangement of water molecules near the surface, in which one hydrogen atom points toward the surface while the other points away (Fig. 4). The second peak represents the interfacial density peak of the adsorbed layer, located approximately 3.0 Å from the surface, with a maximum density of 1.78 g/cm3. Within the adsorbed layer, hydrogen and oxygen atoms exhibit pronounced stratification, indicating that the hydrophilic surface exerts strong adsorption and orientational ordering effects on water molecules. In contrast, for the (00 1 ¯ )-(00 1 ¯ ) model, the interfacial density peak is lower (1.61 g/cm3), and the stratification of hydrogen and oxygen atoms within the adsorbed layer is not pronounced, indicating that the hydrophobic surface imposes relatively weak adsorption and orientational constraints on water molecules. This comparison highlights the critical role of wettability in controlling interfacial water structure and ordering.

Figure 4: (001)-(001) model surface and water molecules.

The (001)-(00 1 ¯ ) asymmetric model exhibits both hydrophilic and hydrophobic surface characteristics. On the hydrophilic (001) side, the interfacial density peak reaches 1.86 g/cm3, whereas on the hydrophobic (00 1 ¯ ) side, it is lower at 1.43 g/cm3. The bulk region stabilizes at a density of 0.985 g/cm3, slightly lower than in the other two models. Compared with the single-surface models, the mixed-wettability configuration amplifies the contrast between the two surfaces. Adsorption is enhanced on the hydrophilic side and weakened on the hydrophobic side. This coupling between the surfaces also modifies the bulk region, reflecting the asymmetric influence of surface wettability on the overall water structure.

3.2 Orientation Analysis

To quantitatively characterize the orientational behavior of water molecules near the kaolinite surfaces, the orientation order parameter SZ is employed to describe the degree of molecular alignment [50]: Sz=1.5×〈cos2α〉−0.5(10) where α is the angle between the considered directional vector and the z-axis, and the bracketed term represents the ensemble average of the MD simulations. When SZ ≈ 0, the water molecules exhibit a random orientational distribution; when SZ > 0, the orientation of the molecules tends to align parallel to the z-axis (with SZ = 1 representing perfect alignment). When SZ < 0, the molecular orientation tends to be perpendicular to the z-axis (with SZ = −0.5 corresponding to perfectly perpendicular alignment).

In the orientation analysis, two types of directional vectors are defined (Fig. 5):

• •OM vector: the vector from the oxygen atom (O) to the midpoint (M) between the two hydrogen atoms, representing the direction of the molecular dipole moment.
• •HH vector: the vector connecting one hydrogen atom to the other, representing the orientation of the molecular plane.

Figure 5: OM and HH direction diagram.

The system was analyzed using a spatial layering approach, similar to the density calculations, with each water molecule assigned to a layer according to the z-coordinate of its oxygen atom. The variations of SZ(OM) and SZ(HH) along the z direction for the three nanopore models are shown in Fig. 6. All three models exhibit a similar trend: water molecules show pronounced orientational ordering within the adsorbed layer, whereas in the diffuse layer, SZ gradually oscillates and decays toward zero, indicating a progressive loss of orientational alignment; in the bulk region, water molecules exhibit a random orientation distribution. In general, water molecules near hydrophilic surfaces show stronger orientational ordering than those near hydrophobic surfaces.

Figure 6: Variations of the two S_Z parameters along the z direction in the three models: (a) (001)-(001) model; (b) (00 1 ¯ )-(00 1 ¯ ) model; (c) (001)-(00 1 ¯ ) model.

For the (001)-(001) model, SZ(OM) is strongly negative within the adsorbed layer, reaching peak values of approximately −0.32, indicating that the dipole moment of the water molecule tends to orient perpendicular to the z-axis [Fig. 6a]. Meanwhile, SZ(HH) is strongly positive, surpassing 0.72 near the surface, showing that the H-H vector tends to align parallel to the z-axis. This “OM-perpendicular and HH-parallel” configuration corresponds to the hydrogen-oxygen stratification observed in the density profiles and suggests that the strong interfacial electric field of the hydrophilic surface induces a stable orientational arrangement of interfacial water molecules. At the boundary between the adsorption and diffuse layers, both SZ(OM) and SZ(HH) exhibit noticeable turning points, after which they oscillate and decay toward zero. The system exhibits pronounced symmetry about the two opposing surfaces, indicating that the dual hydrophilic surfaces form a mirror-symmetric polarized interfacial water structure. This structure reflects the strong adsorption effect of hydroxyl groups on hydrophilic surfaces.

For the (00 1 ¯ )-(00 1 ¯ ) model, SZ(OM) becomes negative within the adsorbed layer, reaching only about −0.28, which are significantly closer to 0 than those near hydrophilic surfaces [Fig. 6b]. Simultaneously, SZ(HH) shows only a weak negative deviation, with a peak of approximately −0.1. These values provide a clear quantitative confirmation that the hydrophobic surface exerts minimal orientational constraints on water molecules. At the boundary between the adsorbed and diffuse layers, the SZ(HH) curve still shows a slight turning pause, whereas no clear turning behavior is observed for SZ(OM). Both parameters gradually decay within the diffuse layer, and their values approach zero in the bulk region, demonstrating that the water molecules become randomly oriented and lose directional preference. These findings are consistent with the density distribution curves and confirm that hydrophobic surfaces impose only weak structural ordering, with no significant hydrogen-oxygen separation within the adsorbed layer.

For the asymmetric (001)-(00 1 ¯ ) model, the upper and lower surfaces exhibit a distinct asymmetry. On the hydrophilic (001) side, SZ(OM) approaches its theoretical limit for perpendicular alignment, dropping close to −0.42, indicating a stronger tendency for the water dipole moment to orient perpendicular to the z-axis [Fig. 6c]. Meanwhile, SZ(HH) is positive, signifying that the H-H vector tends to align parallel to the z-axis. On the hydrophobic (00 1 ¯ ) side, SZ(OM) remains negative and shows a small “turning pause” at the adsorbed-diffuse boundary, while SZ(HH) transitions from a slight negative deviation to a slight positive deviation. This behavior suggests that the presence of the opposite hydrophilic surface may contribute to partial orientational ordering on the hydrophobic side [51,52,53], making the overall orientational behavior tend to resemble that associated with the hydrophilic surface. As a result, the nanopore develops an asymmetric polarized water layer, which is consistent with the trends observed in the density distributions.

3.3 Hydrogen Bond Analysis

The adsorption of water molecules on the alumina (Al-O) surface is significantly stronger than on the silica (Si-O) surface, which also explains why the interfacial water density near the alumina layer is higher than that near the silica layer. To elucidate the mechanisms of water-water and clay-water interfacial interactions, a quantitative analysis of the hydrogen bond distribution was conducted. The hydrogen-bond criterion is defined as follows: a hydrogen bond is considered to exist when the distance between the donor oxygen atom (D) and acceptor oxygen atom (A), RDA, is less than 3.5 Å and the angle ∠HDA is less than 30° [54]. Based on the donor-acceptor origin, hydrogen bonds in the nanochannels can be classified into three categories: Owater-H…Owater (water-water hydrogen bonds), Owater-H…Oclay (hydrogen bonds donated by water to the clay surface), and Oclay-H…Owater (hydrogen bonds donated by clay hydroxyl groups to water molecules), which occur only on the hydrophilic (001) surface.

The system was analyzed using a spatial layering approach, similar to the density calculations, with each water molecule assigned to a layer according to the z-coordinate of its oxygen atom. The hydrogen-bond statistics were calculated by averaging 50 configurations sampled every 0.1 ns from the equilibrated trajectory of the last 5 ns. The hydrogen bond distributions for the three pore models are shown in Fig. 7, and the corresponding total hydrogen bond counts are summarized in Table 3. Regarding spatial distribution, the hydrogen bond density of water molecules increases significantly near the solid surface, exhibiting a characteristic “interfacial peak-decay-stabilization” pattern that is consistent with the layering structure observed in the density profiles. Therefore, the increase in the number of hydrogen bonds near the interface is partly related to the enhanced local water density, but is more importantly attributed to the additional hydrogen-bonding sites provided by the kaolinite surface, especially the surface hydroxyl groups on the hydrophilic (001) surface.

Figure 7: Variation of the number of hydrogen bonds along the z direction in the three models: (a) (001)-(001) model; (b) (00 1 ¯ )-(00 1 ¯ ) model; (c) (001)-(00 1 ¯ ) model.

On the (001) hydrophilic surface, out of the 1152 available surface hydroxyl groups, approximately 60% participate in constructing the interfacial hydrogen bond network, forming approximately 688 Oclay-H…Owater hydrogen bonds between surface hydroxyl groups and water molecules. In addition, a noticeable number of 312 Owater-H…Oclay hydrogen bonds appear within the adsorbed layer, indicating the presence of bidirectional donor-acceptor interactions at the interface. These interactions give rise to a highly ordered interfacial water structure, allowing hydrogen atoms to approach the solid surface more closely. The combined effects of interfacial electric fields and hydrogen bonding produce the interfacial density peak in the density profile.

On the hydrophobic (00 1 ¯ ) surface, the absence of surface hydroxyl groups results in 0 Oclay-H…Owater hydrogen bonds. Furthermore, the number of Owater-H…Oclay hydrogen bonds within the adsorbed layer sharply drops to only 111, which represents a 64.4% reduction compared to the hydrophilic channel. This suggests that the hydrophobic surface imposes weaker constraints on water molecules. Overall, the hydrogen bond peak at hydrophilic surfaces is substantially higher than at hydrophobic surfaces, reflecting stronger interfacial polarization effects.

In the asymmetric (001)-(00 1 ¯ ) model, the numbers of Owater-H…Oclay hydrogen bonds on the hydrophilic and hydrophobic surfaces become more comparable: the hydrophilic side exhibits a reduction to 241, whereas the hydrophobic side shows a corresponding increase to 148. This trend suggests a tendency toward interfacial polarity equalization within the confined slit pore. The average number of hydrogen bonds per water molecule (Hbondwater) is approximately 3.2 in all three models, which is close to the typical value for liquid water at 313 K [55]. In the bulk region, we calculated the average hydrogen bond density (Hbondbulk) per 0.1 Å height. The results showed that the Hbondbulk values of the three models were similar, with the (001)-(00 1 ¯ ) model being slightly lower, which is consistent with the feature that the density in the bulk region is slightly smaller.

By combining the density distribution and orientation analyses, it is evident that the formation of the adsorbed layer involves not only density enrichment but also strong orientational polarization and enhanced hydrogen bonding among water molecules. The hydrophilic surface provides stable hydrogen bond donor sites through its surface hydroxyl groups, promoting the formation of an ordered interfacial hydrogen bond network. As a result, the interfacial water exhibits pronounced orientational alignment and a highly compact structure. In contrast, the hydrophobic surface lacks functional groups capable of forming hydrogen bonds, and the interfacial water structure relies primarily on water-water hydrogen bonding and surface charge effects. Consequently, the orientational ordering is weaker and the interfacial density is slightly lower.

3.4 Other Factors Affecting the Density Distribution of Water Molecules

To further investigate the influence of factors other than interfacial wettability on the structural characteristics of confined water, the density distribution of water molecules in the (001)-(001) kaolinite model was analyzed under different pore widths, temperatures, pressures, and externally applied driving pressure.

As shown in Fig. 8a, when the pore height increases from 20 Å to 100 Å, the overall shape of the water density distribution remains the same, exhibiting a characteristic three-layer structure consisting of an adsorbed layer, diffuse layer, and bulk region. As the pore height increases, the thicknesses of the interfacial adsorbed and diffuse layers remain nearly unchanged, whereas the width of the bulk region expands significantly. This suggests that the layering structure is governed primarily by surface interactions rather than by pore size. In the 20 Å channel, the interfacial water layer from the two surfaces almost overlaps, and the density of water molecules inside the pores shows a trough in the middle of the pores, indicating that the water molecules within the channel are all in the interfacial water layer. When the pore height increases to 40–100 Å, the density in the central region approaches 1.0 g/cm3, which is close to that of bulk water. Therefore, the main influence of pore size on confined water lies in changing the volume fraction of the bulk region, whereas the structural characteristics of the interfacial water layer exhibit a degree of scale invariance.

Figure 8: Effects of different factors on the water density distribution in the (001)-(001) model: (a) channel heights; (b) temperatures; (c) normal pressures; (d) driving pressures.

As shown in Fig. 8b, when the temperature increases from 298 K to 328 K, the peak densities in both the adsorbed and diffuse layers decrease to varying extents. The elevated temperature enhances thermal motion, enabling water molecules to overcome the interfacial adsorbed potential wells more easily, thereby weakening the density peaks of the interfacial water layer. With increasing temperature, the equilibrium pore height increases slightly, accompanied by a slight increase in the bulk-region thickness and a slight decrease in the average density of the bulk water region. Therefore, when the profiles are aligned near the lower surface for comparison, the adsorbed and diffuse layers near the upper surface do not completely overlap among different temperature conditions.

As shown in Fig. 8c, the density profiles obtained at three normal pressures (0.1 MPa, 23.5 MPa, and 40 MPa) exhibit an overall upward shift with increasing pressure, indicating an increase in the average system density. High pressure strengthens intermolecular interactions and local packing among water molecules, leading to a denser arrangement and a slight elevation of the density peaks.

As shown in Fig. 8d, when different magnitudes of driving pressure Δ P (0–79.3 MPa) are applied along the y direction, the density curves nearly overlap, indicating that the driving pressure has a negligible effect on the density distribution of water molecules. Both the peak positions and amplitudes remain essentially unchanged in the adsorbed, diffuse, and bulk regions. This suggests that the driving pressure primarily affects the velocity distribution of water flow, while its influence on the static density distribution and structural layering is negligible. Therefore, the subsequent analyses of velocity fields and slip phenomena can be carried out under the assumption of a stable density distribution.

4 Flow Mechanisms of Water Molecules in Kaolinite Nanopores

4.1 Analysis of Slip-Flow Velocity

To investigate the flow mechanisms of water molecules in kaolinite nanopores, the slip-flow behavior under different surface-wettability conditions was analyzed using a spatial layering approach. The simulation domain was divided into layer units of 1 Å thickness. The average displacement of the oxygen atoms in water molecules within each layer was computed over time intervals of 100 ps, and the time-averaged values during the steady-state stage after 5 ns were used as the fluid velocity under the driving pressure. To ensure statistical reliability, only layers containing more than 100 water molecules (approximately one-third of the average density) were included in the calculation.

Fig. 9 shows the velocity distributions of the three pore models under conditions of 313 K and 23.5 MPa with the driving pressure of 42.9 MPa. All velocity curves exhibit a characteristic Poiseuille-type parabolic shape superimposed with slip effects [56]. The velocity near the pore walls does not drop to zero, indicating the presence of significant boundary slip.

Figure 9: Velocity distributions of the three models under the driving pressure of 42.9 MPa.

In the (001)-(001) model, the velocity profile is relatively smooth, with a maximum velocity of approximately 24.43 m/s and a boundary slip velocity of only 3.53 m/s. This suggests that the hydrophilic surface exerts strong adsorption and hindering effects on water molecules, thereby yielding flow behavior that closely resembles classical Poiseuille flow.

In the (00 1 ¯ )-(00 1 ¯ ) model, the velocity peak is the highest (36.39 m/s), and the boundary slip velocity reaches 16.35 m/s. The slip ratio is significantly larger, indicating that the hydrophobic surfaces markedly reduce interfacial friction, allowing water molecules to “slide” along the pore walls with much greater ease.

In the asymmetric (001)-(00 1 ¯ ) model, the velocity distribution lies between the two symmetric cases but clearly leans toward hydrophobic characteristics. The slip velocity on the hydrophilic side is relatively small (5.04 m/s), whereas the slip velocity on the hydrophobic side is much larger (18.68 m/s), causing the velocity peak to shift away from the channel center. The maximum velocity reaches 32.9 m/s. These results demonstrate that the velocity field within the pore exhibits pronounced asymmetry, with water molecules attaining higher velocities near the hydrophobic surface. This reflects the stronger adsorption on the hydrophilic side and weaker confinement on the hydrophobic side.

According to the modified Poiseuille flow model for laminar flow between parallel plates, the velocity distribution can be expressed as follows [36]: vyz=−12μΔPLyz2−h∗22−ls·h∗(11) where μ is the dynamic viscosity of water, ls is the slip length, h* is the effective pore-channel height, Ly is the length of the flow direction, and z is the distance from the center of the flow channel (where the velocity is at its maximum) to the calculation point. The calculated results obtained by fitting the parabolic portion of the velocity profiles are summarized in Table 4.

The viscosity of the simulated fluid was calculated to be approximately 0.65 cP, which is consistent with the macroscopic viscosity of water at 313 K. In hydrophilic channels, stable hydrogen bond networks form between surface hydroxyl groups and water molecules, leading to strong orientational ordering and restricted mobility of interfacial water. As a result, the slip length is relatively small (around 2.5 Å), corresponding to an “adhesive flow” regime. In contrast, in hydrophobic channels, the absence of polar functional groups results in weak interactions between water molecules and the channel walls, causing the slip length to increase to more than 12 Å, characteristic of a “weakly constrained flow” regime.

In the asymmetric (001)-(00 1 ¯ ) channel, the slip lengths on both surfaces increase, indicating mutual influence between the two interfacial properties. This may be associated with interfacial coupling across the hydrophilic-hydrophobic interface, where the hydrophilic side becomes comparatively more prone to slip, while the hydrophobic side also exhibits enhanced slip behavior. As a result, the overall flow characteristics of the channel tend to shift toward those observed in the hydrophobic system.

4.2 Influence of Channel Heights and Driving Pressure on Slip-Flow Velocity

As indicated by the density profiles in Fig. 8a, the width of the stable bulk region is directly determined by the height of the flow channel, and the proportion of the bulk region within channels of different sizes is a key factor governing the differences between microscopic flow behavior and macroscopic flow characteristics. The velocity distributions of the (001)-(001) model under the same driving pressure of 42.9 MPa but with different channel heights are illustrated in Fig. 10.

Figure 10: Velocity distributions in the (001)-(001) model with different channel heights under the driving pressure of 42.9 MPa.

As the channel height increases, the proportion of the bulk region increases, the overall velocity profile becomes smoother, and the peak velocity rises. This indicates that wider channels reduce the fraction of water molecules strongly influenced by wall-water interactions, thereby decreasing the effective viscosity of the system. In nanoscale channels, variations in pore size imply that different fractions of water molecules are subject to the potential field of the kaolinite surface. Water molecules located close to the surface reside in the adsorbed and diffuse layers, where their motion is constrained by both the surface potential and the hydrogen bond network, resulting in higher viscosity. In contrast, water molecules in the bulk region farther from the surface experience weaker confinement, leading to a pronounced increase in flow velocity.

Consistent with these observations, the calculated results in Table 5 show that when the channel height increases from 20 Å to 100 Å, the effective viscosity of water decreases from 1.60 cP to 0.55 cP, while the slip length decreases slightly. This behavior suggests that the flow system gradually approaches a free-flow regime with minimal wall confinement as the channel widens. At the channel height of 60 Å, the viscosity enhancement caused by interfacial water layers may be partly compensated by the intrinsic underestimation of water viscosity by the SPC model, resulting in an effective viscosity close to the macroscopic value.

To further verify the stability of fluid viscosity and slip length under different driving pressure conditions, simulations were performed on the three kaolinite nanopore models by applying driving pressure of varying magnitudes (3.6–114.4 MPa). The results indicate that within a certain driving pressure range, the velocity profiles consistently follow the Poiseuille-flow pattern superimposed with slip flow (Fig. 11). The velocity distributions maintain a parabolic shape, the boundary velocity remains nonzero, and the profiles become smoother with more pronounced peaks at higher driving pressure. This demonstrates that stronger driving pressures can effectively overcome disturbances caused by thermal fluctuations, causing the flow to resemble a stable laminar regime more closely.

Figure 11: Velocity distributions of the three models under different driving pressures: (a) (001)-(001) model; (b) (00 1 ¯ )-(00 1 ¯ ) model; (c) (001)-(00 1 ¯ ) model.

The corresponding calculated parameters for each model are summarized in Table 6. It is evident that within the tested driving pressure range, the variations in viscosity and slip length are relatively small for all models, remaining within a stable interval. This suggests that the intrinsic fluid properties of water within the confined nanochannels do not change significantly with the driving pressure, confirming the robust dynamical stability of the system.

4.3 Microscopic Verification of Darcy’s Law and Equivalent Permeability Analysis

To establish a connection between molecular scale simulation results and macroscopic seepage behavior, this study attempts to verify the applicability of Darcy’s law within the molecular dynamics framework. For the present microscopic model, the classical Darcy’s law can be expressed as: v¯=kγwμNmayγwV(12) where v ¯ is the average fluid velocity, k is the permeability, γw is the unit weight of water, and V is the simulated fluid volume.

This expression represents an “equivalent Darcy form” under molecular dynamics conditions, in which the term Nmay/γwV may be regarded as the equivalent hydraulic gradient i. By fitting the linear relationship between v ¯ and i, the equivalent permeability coefficient K of the system is obtained to determine whether the microscale flow process obeys Darcy’s law.

The Darcy fitting curves for the three pore models under driving pressures ≥3.6 MPa are presented in Fig. 12. Within the investigated driving-force range, the average velocity of water molecules shows an approximately linear relationship with the equivalent hydraulic gradient. This result indicates that, under sufficiently large external driving forces where stable directional flow can be obtained, the water flow in the kaolinite nanopores exhibits a Darcy-type linear response. The fitted parameters are listed in Table 7.

Figure 12: The Darcy curves of the simulation results of the three models.

The permeability coefficients obtained from the slopes of the fitted curves follow the order: (001)-(001) model < (001)-(00 1 ¯ ) model < (00 1 ¯ )-(00 1 ¯ ) model. This trend corresponds to the increasing hydrophobicity of the channel surfaces. In the hydrophilic channel, surface hydroxyl groups form stable hydrogen-bond networks with interfacial water molecules, leading to stronger adsorption, more pronounced orientational ordering, and greater resistance to flow. In contrast, the hydrophobic channel exhibits weaker water-clay interactions and lower interfacial friction, which promotes boundary slip and results in a larger apparent permeability. Therefore, the fitted permeability coefficients reflect the important role of surface wettability in regulating nanoscale seepage behavior through interfacial water structure and slip effects.

It should be noted that the threshold hydraulic gradient was obtained by linear extrapolation of the fitted Darcy-type curve, rather than by direct simulation in the extremely low-driving-force regime. In NEMD simulations, molecular thermal motion strongly interferes with flow responses under very low driving forces, making it difficult to obtain a stable directional flow within accessible simulation times. Therefore, relatively high driving forces are commonly applied in NEMD simulations to overcome thermal fluctuations and improve statistical efficiency. As a result, the flow behavior near the true threshold-gradient was not directly quantified in this study. The apparent threshold hydraulic gradient obtained here should be regarded as a fitted nanoscale indicator of interfacial flow resistance.

In addition, the quantitative values obtained in this study are affected by the limitations of the molecular model. The SPC water model may not perfectly reproduce all transport properties of real water, and the assumption of rigid kaolinite walls neglects possible surface deformation and mineral-water coupling under flow. The present Darcy-type analysis is intended to provide a qualitative molecular-scale interpretation of how surface wettability, interfacial water structure, and boundary slip influence the seepage capacity of kaolinite nanopores, while direct quantitative extrapolation to natural shale or macroscopic seepage conditions should be made with caution.

5 Conclusions

This study systematically investigates the adsorption structures and flow behaviors of water in three types of kaolinite nanopores with different wettability using MD simulations. The main conclusions are as follows:

• (1)The wettability of kaolinite pore surfaces significantly influences the nanoscale structure of confined water by regulating the number and spatial distribution of interfacial hydrogen-bond interactions. Hydrophilic surfaces form a large number of clay-water hydrogen bonds between interfacial water molecules and surface hydroxyl groups, resulting in pronounced orientational ordering and more distinct layering structures in the interfacial region. In contrast, hydrophobic surfaces form only a limited number of clay-water hydrogen bonds, leading to reduced orientational order and lower density peak values. These results suggest that pore-surface properties are key microscopic factors governing the occurrence state and migration behavior of water in nanopores, providing a basis for understanding its subsequent transport behavior.
• (2)The distinct structural properties of the interfacial water layer directly dictate anomalous flow dynamics that deviate from classical continuum predictions. As the pore size decreases, the increasing volumetric fraction of interfacial water layer within the pore space elevates the effective viscosity and enhances solid-liquid coupling, driving pronounced deviations in viscosity, density and velocity. Consequently, a larger external driving pressure is required to achieve the same volume-averaged flow velocity, and the flow exhibits a Poiseuille-type velocity profile with pronounced velocity slip. Moreover, differences in interfacial interaction strength under varying surface-wettability conditions give rise to significant variations in both slip velocity and mean flow velocity.
• (3)The average flow velocity of water in nanopores exhibits an approximately linear relationship with the hydraulic gradient within the investigated range, which can be regarded as an effective nanoscale manifestation of Darcy’s law. Fitting results obtained under different wettability conditions indicate that the effective permeability of the system increases with increasing surface hydrophobicity, demonstrating that interfacial wettability plays a crucial role in controlling seepage capacity through its regulation of slip effects. Furthermore, the threshold hydraulic gradient for the three models was derived from the linear extrapolation of the Darcy fitting curves. This extrapolated result is qualitatively consistent with the nonzero threshold gradients commonly observed in macroscopic clay materials, suggesting that nanoscale pore channels, strong interfacial interactions, and size effects may contribute to the microscopic origin of threshold hydraulic gradients.