Nanoscale soil-water retention curve of unsaturated clay via molecular dynamics

. This paper characterizes nanoscale soil-water retention mechanism of unsaturated clay through molecular dynamics simulation. Series of molecular dynamics simulations of clay at low degrees of saturation were conducted. Soil water was represented by a point cloud through the centre-of-mass method. Water-air interface area was measured numerically by the alpha shape method. Spatial variation of water number density is characterized and used to determine the adsorbed water layer. The soil-water retention mechanism at the nanoscale was analysed by distinguishing adsorptive pressure and capillary pressure at different mass water contents and considering apparent interface area (water-air interface area per unit water volume).


Introduction
The physics and mechanics of unsaturated soils are significant in geotechnical and geo-environmental engineering. Soil water retention curve is fundamental for understanding the relationship between water content and matric suction [1]. For instance, soil water retention curve is required in modelling multiphase fluid flow, shear deformation, and stress-strain relationships of unsaturated soils. The adsorptive pressure might be ignored at a high degree of saturation. However, it has been found that adsorption plays an important role at low degrees of saturation [2]. Furthermore, both experimental and theoretical studies have suggested that wetting-nonwetting (e.g., waterair) interface should be considered to better describe soil water retention curves of unsaturated soils. The existence of wetting-nonwetting interface could cause changes in soil volume and the stress states [3].
Due to the strong justification for including a fourth independent phase namely the water-air interface in unsaturated soil, several techniques have been developed to make water-air interfacial area measurement experimentally accessible. These technologies are mainly divided into two main categories: interfacial active tracer method and imaging techniques such as synchrotron X-ray microtomography and environmental scanning electron microscope. In addition, computational modelling through physics-based numerical methods has gained success in resolving and quantifying water-air interface in porous media. One common method is the porenetwork modelling technology [3]. However, it remains challenging to quantify the impact of adsorption on the soil water retention curve and explain the mechanism of soil-water adsorption at the nanoscale. This opening question matters because of the multiscale and multiphysics nature of clays. At the nanoscale adsorptive forces in fine-grained clay become pronounced and could modify the water structure, e.g., adsorptive water film tightly attached to clay surface [5]. The mechanical and physical behaviour of clay at the continuum scale are related to its physical behaviour at the nanoscale. For example, the interatomic and intermolecular bonding forces at the nanoscale hold clay particles together. These unbalanced forces that exist at phase boundaries could influence the formation of clay minerals, the structure, size, and shape of soil particles, and further influence the macroscopic properties [2]. The laboratory measurement techniques suffice to account for capillary effects in the water retention mechanism, however, indeed no viable experimental technique exists to quantify adsorption and measure water-air interfacial area in unsaturated soils at the nanoscale. For this reason, molecular dynamics (MD) simulations are widely used to gain detailed insights into the physics and mechanics of clay and related phases at the atomistic scale [6][7][8][9][10][11][12][13]. For example, Song & Zhang [6] measured contact angle between clay and water at elevated temperatures from molecular dynamics simulations. The nanoscale hydrodynamics of claywater systems have been also investigated through nonequilibrium MD simulations [8].
In this study, we conduct full scale MD simulations of unsaturated clay at various water contents. By using the alpha shape method, we measure the water-air interfacial area of unsaturated clay at the nanoscale and incorporate it in soil water retention curves. Adsorption and capillarity are considered explicitly and temperature effect on pore water pressure is investigated.

Materials and methods
The MD model of unsaturated clay consists of two parallel equal size pyrophyllite platelets and confined water between them. The chemical formula of pyrophyllite is 2 [ 4 10 ]( ) 2 . The dimensions of pyrophyllite unit cell are 5.28 Å × 9.14 Å × 6.56 Å in the x-y-z Cartesian coordinate system. Each clay platelet consists of 36 × 21 × 1 unit cells, which corresponds to 190.08 Å × 191.94 Å in the x-y plane. The interlayer spacing between two clay platelets is 5 nm.
We use CHARMM force field to describe intermolecular and intramolecular interactions. Rigid water molecules are described using the TIP3P model. During simulation, the − bond length and − − bend angle are kept constant as 0.9572 Å and 104.52° using the SHAKE algorithm. The total potential energy is the combination of bonded energy and nonbonded energy. In this study, clay platelets are immobilized, and water molecules are free to move. Therefore, only non-bonded interactions between clay and water are considered. The non-bonded interaction is expressed as where is the well-depth of Lennard-Jones (LJ) potential, the distance at the minimum LJ interactions energy, and the charge of atoms and , the distance between atom and .
The cut-off radii for van der Waals and Coulombic interactions are 10 Å. Table 1 lists the nonbonded parameters for clay-water system. ℎ and represent water hydrogen and water oxygen, respectively. and ℎ are hydrogen groups in the octahedral layer in clay. and are silicon and oxygen atoms in the tetrahedral layer.
is the bridging oxygen connecting octahedral layer and tetrahedral layer, and is the aluminium atoms in the octahedral layer. All MD simulations were performed on LAMMPS, a large scale atomic/molecular massively parallel simulator [14] using NVT ensemble with a time step of 0.5 fs. Langevin dynamics was used to maintain the temperature. The Verlet algorithm was employed to integrate the equations of motion of water molecules. Periodic boundary conditions were applied to all directions. The simulation was first run for 1 ns to bring the system into equilibrium. Then the production simulation was run for 0.2 ns to output the averaged water trajectories. In this study, mass water content ( ) varies with different numbers of water molecules in the clay nanopore. Mass water content is defined as the ratio between the weight of water and the weight of dry clay = , (2) where and is the number of water molecules and clay unit cells, respectively. = 18 g/mol is the molar mass of water, and = 360 g/mol is the molar mass of pyrophyllite unit cell. Fig. 1 shows the clay-water system at equilibrium with mass water content 31.3%.

Measurement of water-air interfacial area using alpha shape method
In this study, we first converted water molecules to a 3-D point cloud using the centre-of-mass method. Fig. 2 shows this process. Having obtained the water point cloud, we reconstruct the surface of point cloud using alpha shape method [15] to measure the surface area. The general idea is finding piece-wise triangles (the so-called alpha shapes) to represent the surface of water point cloud. Fig. 3 shows the schematic procedure of interfacial area calculation using the alpha shape method [12]. We could define a point set where each point represents the center-of-mass of one water molecule. Let be an open ball with radius . We restrict to be empty such that it can occupy its space without enclosing any of the points of , i.e., ∩ = ∅. Let = { , , } be a subset of . Given , we could define a 2-simplex Δ (i.e., a triangle) as the convex hull of . Here, the convex hull is the smallest convex set that contains all points in . The 2 simplex is said to be -exposed if = ∩ . Here is the boundary of open ball . For implementation, we first perform the Delaunay triangulation of the surface of the water point cloud and then define the alpha complex that associates Delaunay triangulation with the alpha shape. The algorithm adopted is summarized as follows: a) Compute the Delaunay triangulation of , knowing that the boundary of -shape is contained in Delaunay triangulation. b) Determine -complex by inspecting all simplices Δ in Delaunay triangulation. If the circumsphere of Δ is empty and the radius of circumsphere is smaller than , we accept Δ as a member of -complex. c) All simplices on the boundary of -complex form the -shape. In Fig. 3 (2), the grey object is an empty open ball , the red sphere is the boundary , and = ∩ = { 1 , 2 , 3 }. Meanwhile, we must have ∩ = ∅ and is exterior to . The spherical cap is straightened by a 2-simplex (i.e., triangle) connected by points 1 , 2 , and 3 . Thus, the alpha shape of is the polytope whose boundary consists of all the 2simplices/triangles. The interfacial area is the sum of the area of each boundary triangle.

Numerical results
In this section, we present the numerical results of water-air interfacial area, pore water pressure, capillarity and adsorption at various mass water contents. The simulated mass water content ranges from 6.2% to 31.3%. Fig. 4 presents four unsaturated clay MD models with different mass water contents.

Pore water pressure
In MD simulation, the pressure tensor (or stress tensor) can be expressed by the virial form where is atom index, is the number of atoms in the system, is the system volume. , , , and represent the mass, velocity, position, and force of atom , respectively. The pore water pressure can be determined from = 1 3 ( 11 + 22 + 33 ),

Fig. 5.
Variation of pore-water pressure with respect to mass water content. Fig. 5 plots the pore-water pressure with respect to mass water content. It shows that as mass water content increases from 6.2% to 31.3%, the negative pore-water pressure decreases from -208.6 atm to -256 atm (1 atm = 101.325 kPa).
Here, the negative water consists of two parts -adsorptive water pressure and capillary water pressure. Thus, the decrease of negative water pressure with increasing mass water content may be due to the increase of absolute adsorptive water pressure, which deserves further study. The effect of temperature on pore water pressure is briefly investigated. Fig. 6 plots the variation of the negative pore water pressure at elevated temperatures from 298 K to 348 K. The results demonstrate that the magnitude of pore water pressures at constant mass water content decreases as temperature increases. For example, the magnitude of pore water pressure at T = 298 K is 256 atm. As temperature increases to 348 K, the magnitude of pore water pressure decreases to 245.8 atm. This observed phenomenon may be due to the decrease of the adhesion between the meniscus water and the clay platelet and the decrease in surface tension of the meniscus water by the temperature increase. The general temperature dependence of pore water pressure at the nanoscale from MD modelling is consistent with the trend observed at the continuum scale [16].

Water-air interfacial area
Water-air interfacial area is measured using alphashape method. Fig. 7 plots the variation of water-air interfacial area as a function of mass water content. The interfacial area increases from 5445 to 12768 Å 2 as mass water content varies from 6.2% to 31.3%. Fig. 8 plots the apparent interfacial area as a function of mass water content. The apparent water-air interfacial area is defined as the water-air interfacial area per unit water volume. The apparent interfacial area increases with the decrease of mass water content. As mass water content increases from 6.2% to 31.3%, the apparent interfacial area decreases from 0.15 A −1 to 0.05 A −1 . It is noted that the apparent interfacial area-water content relationship obtained from MD simulation is consistent with the laboratory testing results [17], which show that the air-water interfacial areas are inversely proportional to water content. The nanoscale soil-water retention surface incorporating water-air interfacial area is finally presented in Fig. 9.

Adsorptive and capillary pressure at the nanoscale
Under the assumption of a local thermodynamic energy equilibrium, soil water is one of two types (capillary and adsorptive). To distinguish capillarity and adsorption, we characterize spatial variation of water number density from MD simulations [8]. To do this, the nanopore space is divided into several layers of equal thickness along the z direction. The number density of each water layer is computed as the ratio between the number of water molecules in the layer and the layer volume. Fig. 10. Number density profiles of soil water at different mass water contents. Fig. 10 shows the water number density profiles at different mass water contents. The first and second density peaks are found 3.75 Å and 6.75 Å from clay surface, respectively. Further away from clay surface, water number density gradually decreases. At the centre of the clay nanopore (i.e., z = 0), water number density reaches the local minimum. Large number density fluctuations in the vicinity of clay-water interface could indicate the strong effect of soil water adsorption. It also results in the layered water structure. The thickness of adsorbed water layer is a variable ranging from only a fraction of the surface being covered by one molecular layer up to many tens of molecular layers. It is a function of the electrochemical environment and interparticle stress. Based on numerical results of water number density, we assume the thickness of adsorbed water layer in this work is approximately the distance between the outermost water layer and the local minimum after the second density peak. The adsorptive water layer is bounded by two dash horizontal lines in Fig. 8. The average thickness of adsorptive water layer is 7 Å under different mass water contents.
Based on this criterion, the soil water could be partitioned into two regions, e.g., adsorptive region and capillary region. Fig. 11 plots the variation of percentages of adsorptive and capillary water pressures versus mass water content. At low mass water content, Fig. 11. Percentage of adsorptive and capillary water pressure as a function of mass water content. e.g., =6.2%, the adsorptive water pressure occupies 62.9% of the total pore water pressure. With the increase of mass water content, the effect of adsorption is gradually weakened. When the mass water content is around 15%, the effect of adsorption and capillarity are similar. The percentage of adsorptive pressure fluctuates around 45% when mass water content exceeds 20% and capillarity becomes a dominant factor in the total negative porewater pressure (i.e., matric suction). During wetting process, one sees the significant increase in capillary pressure as mass water content increases to 20%. This observation is consistent with the rate of change in water-air interfacial area before and after 20% water content as shown in Fig. 7. The increasing rates of interfacial area before and after 20% water content are 338 Å 2 and 237 Å 2 per 1% mass water content. It may be concluded from the MD results that the adsorption plays a significant role in the soil water retention mechanism at low degree of saturation in clay. An empirical formula of adsorptive pressure induced by long-range molecular (van der Waals) force [5] is expressed as where is the Hamaker constant and = 7 Å is the thickness of adsorptive water layer. For a soil-water system, is in the order of -10 20 to -10 19 Joules. In this work we take = -6×10 20 Joules suggested by [5]. Thus, the calculated adsorptive pressure is about -143.1 atm, which is in the same order with the numerical results from MD simulation (as shown in Fig. 12). Though the effect of water content is not considered in the above equation, the general consistence could imply that the adsorption potential is dominated by van der Waals force at the nanoscale.

Summary
MD simulations are conducted to investigate the soil water adsorptive and capillary mechanism of clay. The MD model of the unsaturated clay consists of two parallel clay plates and confined soil water. For postprocessing, soil water was represented by a point cloud using the centre-of-mass method. Water-air interfacial area was measured using the alpha shape method. We characterize water number density profile and calculate the thickness of adsorbed water layer. Based on this, we distinguish adsorption from capillarity. The nanoscale soil water retention curves were plotted in terms of mass water content, apparent interfacial area, and pore water pressure. The MD simulation results have demonstrated that at low mass water contents, adsorption is the dominant mechanism on soil water retention. For instance, the adsorptive water pressure accounts for more than 60% of the total pore water pressure at low mass water content around 6%. The apparent interfacial area (water-air interfacial area per unit water volume) decreases during the imbibition of water. The numerical results also showed that the magnitude of pore water pressures at a constant mass water content decreases as temperature increases which agrees with the laboratory testing result.