Pore network modeling approach for simulating soil water retention curve under different stress conditions

. The field of unsaturated soil mechanics has recently seen the introduction of pore network models, which attempt to replicate the void structure of porous materials. They are robust physically-based simulation tools and have been used to simulate constitutive relationships like soil water retention curves (SWRC) and unsaturated hydraulic conductivity functions. This work aims to present a pore network modeling approach for predicting hysteretic SWRC at various stress states using only grain size distribution and porosity data of the granular soils. The soil sample subjected to given stress conditions is simulated using the Discrete Element Method to obtain a stable packing of spherical particles representing the soil structure. From this packing, an algorithm based on the medial axis is availed to extract a network of pores and throats that describes the geometry and topology of the void structure of the granular soils. Various pore-scale mechanisms like the piston-like advance, corner flow, pore body filling, and snap-off are then used to model fluid displacements at the pore scale to simulate SWRC along the drying and wetting paths. The modeled SWRC under various stress conditions is compared with the measured curves obtained for granular soils from the literature, and the predictions are in good agreement with the experimental results.


Introduction
The soil water retention curve (SWRC) is an important constitutive relation between soil suction and its water content. The hydromechanical properties of unsaturated soils are related to SWRCs, which are hysteretic and consist of a drying branch where air invades the pore spaces of soil to replace water and a wetting branch where water replaces air in the soil pores [1]. SWRC serves as a foundation for modeling several unsaturated soil properties, including unsaturated hydraulic conductivity and shear-strength functions. Hence it is crucial to estimate SWRC correctly under varied in situ conditions.
The stress state of soil significantly affects the SWRC. Several experimental investigations have been conducted to study the influence of the stress state on SWRC [2,3]. Malaya and Sreedeep [4] objectively analyzed and reviewed the findings reported in the literature to reveal the impact of stress history on the SWRC. Most soils show a decrease in water content, an increase in air entry value, and a steeper SWRC at highstress levels. The differences in SWRC originate from variations in the void structure caused due to alteration of the stress state of soils. Therefore, the stress state impacts major parameters in the SWCR, and failing to take this into account may lead to incorrect conclusions about the soil's permeability or shear strength in the field.
Experimental procedures to determine the SWRC under various stress conditions are very timeconsuming. The researchers have proposed several constitutive models to predict SWRC at different net stresses [5,6]. The calibration of model parameters of these constitutive models also needs experimental measurements of SWRC for at least two stress state conditions for predictions at higher stress states. Although the constitutive models preclude conducting the experiments at all the stress states, the experiments required for calibration for stress dependency are also time-consuming due to the axis translation technique.
Recent advances in computational techniques have led to the development of pore network modeling (PNM) for simulating two-phase flow in porous media, such as soil. In this simulation approach, porous soil structure is modeled as a network of pores and the throats where pores occupy a sizeable volume of the network. Khaksar et al. [7] used PNM to study the effect of stress levels on SWRC. Experimental SWRCs data at varying stress levels are used to determine the optimal structural parameters for building a pore network that best matches the data. SWRCs at other stress levels are then generated using the calibrated pore network model. Rostami et al. [8] used PNM for simulating SWRC under various net stresses. Experimental SWRC at zero net confining pressure is used to construct a representative pore network of granular soil. The pore network model at different net confining pressure is obtained by applying the volume change to the pores and throats of the representative pore network at zero net confining pressure in accordance with the change in porosity of the soil when subjected to higher stresses.
However, the pore network modeling techniques also require experimental SWRCs for calibration. Integrating the Discrete Element Method (DEM) with PNM is another potential simulation approach that needs only rudimentary experimental data (i.e., grain size distribution and porosity) [9,10]. To simulate the void structure of the porous media, this method extracts a pore network from a stable packing of spheres modeled using DEM. A wide range of geo-mechanical processes has been simulated using DEM, which is based on modeling the interaction between individual particles of a granular system. A particle packing with the same particle size distribution and packing density (porosity) as the porous media to be modeled is generated using DEM. Flow simulations are carried out in the pore networks extracted from the packing using various network extraction techniques like triangulation, medial axis-based algorithms, and maximal ball-based network extraction algorithms. Yuan et al. [11] generated packing of spheres representing a system of glass beads using DEM to simulate drainage capillary pressure-saturation curves of glass beads. SWRC and unsaturated hydraulic conductivity of granular soils have also been successfully modeled using a combined DEM and PNM approach [9,10].
This study proposes an integrated DEM and PNM numerical framework to simulate SWRC for granular soils under various stress conditions. DEM models the granular soil sample subjected to stresses to replicate the corresponding soil packing. The medial axis-based network extraction algorithm proposed by Raeini et al. [12] is used to extract the network of pores and throats from the packing. Air-water flow is simulated in the pore networks representing the void structure of granular soils at various stress states to predict drying and wetting SWRCs. The organization of the paper includes four sections; the first is the introduction. Section 2 presents a thorough discussion of the suggested numerical framework. In section 3, the predicted hysteretic SWRCs for granular soils under various stress conditions are contrasted with the experimental data from the literature. Finally, section 4 presents the conclusions and suggestions for future research.

Methodology
The pore-scale numerical framework presented for predicting hysteretic SWRC of granular soils under various stress conditions can be broadly divided into three steps: firstly, the packing corresponding to the stress state of the soil is generated. It is followed by extracting the pore network, which represents the geometry and topology of the pore space. Finally, drainage and imbibition are simulated in the pore network to get hysteretic SWRCs. The following presumptions are made during the model's formulation: 1. Primarily, soil suction is due to the combined effect of capillarity and surface adsorption. Nonetheless, the capillary effect predominates in granular soils [13], and hence, adsorption is neglected in the proposed model. 2. It is presumed that the granular soils do not undergo volume change during drying and wetting cycles. 3. It is assumed that the grain size distribution (GSD) of soil does not change on stress application. This assumption is plausible for the present model, given that the stresses generally encountered by soils under field circumstances are in the order of kilopascals, while particle breakage (which can lead to a change in soil gradation) occurs at very high stresses which are in the order of megapascals.

Constructing the pore network model of granular soils subjected to stresses
For simulating SWRC of granular soils under no external stresses, a synthetic packing with the same GSD and porosity as the soil to be modeled is first generated, then the pore network is extracted from it. Granular soil particles are typically described as spheres, and the synthetic assembly is thus made up of spherical particles [9,10]. When granular soils are subjected to stresses, they undergo a change in packing density due to particle rearrangement. The porosity of the soil sample provides a quantitative measure of the variation in packing density. The current work postulates that the representative packing of a soil sample at a given stress condition is the packing generated synthetically using the GSD of the soil and a porosity equal to the porosity that the soil sample achieves after being subjected to stresses. Let us consider a granular soil sample with porosity n when it is not subjected to external stress. Changing the stress state causes the porosity to shift from n to n' while maintaining the soil gradation under the assumption of no particle breakage. A synthetic sample of spherical particles with the same GSD as that of the granular soil and porosity of n' is generated to predict the SWRC of the sample under the present stress state (Fig. 1). In this study, PFC 3D (Itasca 2014) is used for DEM simulations. The size of the synthetic assembly is governed by the GSD and is approximately 15-20 times the average particle diameter [10]. The DEM model parameters are given in Table 1. After the packing is generated, a pore network (Fig. 2) is derived from it using the network extraction algorithm of Raeini et al. [12]. The pore network is used to model constitutive relations like SWRC and is assumed to be analogous to the void structure of granular soil under the present stress state. However, DEM can simulate the packing of granular soils with only narrow GSD. For granular soils with a wide GSD, the underlying principle remains the same, i.e., generating a packing with GSD and porosity corresponding to the soil sample after stress application. The procedure proposed by Mufti and Das [14] is employed to generate a pore network of granular soils with a wide GSD and is briefly discussed herein. The following is the framework for constructing the pore network of the test soil with a wide GSD, shown in  1. Considering that DEM can accurately model granular packing with a size ratio (ratio of biggest particle size to smallest) of ~ 10, the GSD of test soil (Fig. 3) is segmented into three intervals such that the size ratio of any interval does not exceed 10. 2. Grain size (diameter) larger than 0.075 mm falls within interval 1. The grain sizes in interval 2 are between 0.075 mm and 0.0093 mm. Interval 3 contains grain sizes less than 0.0093 mm. Using the grain size information, an equivalent GSD for each interval is computed and shown in Fig. 3. 3. A synthetic packing is obtained for each interval using the equivalent GSD of that interval. For intervals 2 and 3, the porosity at which synthetic packing is modeled is n', and for interval 1, the porosity (n1) is determined using the equation: where is the weight of soil grains in an interval i (i=1, 2, 3). 4. The Raeini et al. [12] algorithm extracts the pore networks from each packing once the packing corresponding to the intervals is modeled. 5. The integration of pore networks from all intervals is the last step in the building of the multiscale pore network.
The assumption is made that the multiscale pore network accurately represents the pore space of the soil to be modeled, and simulations of two-phase flow are conducted within this pore network. The pore networks of intervals 2 and 3 (fine networks) are placed in the pores of the pore network of interval 1(coarse network) to obtain the multiscale pore network. There are three key factors that determine the formulation of the multiscale pore network: the number of pore networks of intervals 2 and 3 that need to be added, their placement within the pore network of interval 1, and the connections between networks of different intervals. The number of fine networks ( ) from intervals 2 and 3 is obtained using the following formula: where (Ω ) is the weight of the synthetic assembly corresponding to any interval i (i = 2, 3). The number of fine networks obtained using this formulation ensures that the particle size distribution of the final domain after network integration matches the GSD of the soil sample to be modeled. After determining the appropriate number of fine networks ( , ), the next step involves placing them into the pores of the coarse network. No attempt at spatial correlation has been made in the placement of the networks within the pores of the coarse network. A pore from the coarse network is randomly chosen, and then fine networks are inserted to fill the pore entirely. 6. This is followed by providing inter-network connectivity between pore networks of different intervals to obtain the representative multiscale pore network. The boundary pores of the fine network are connected with the neighboring pores of other pore networks to establish connectivity between them. Fig. 4 shows a 2D diagrammatic representation of a fine network positioned within a coarse network and demonstrates inter-network connectivity for a fine pore. In this figure, Pf is a boundary pore of a fine network, which is connected with neighboring pores Pn of other networks. The term "neighboring pores" refers to those pores whose distance from the boundary pore of a fine network (interpore distance between Pf and Pn), as measured by the length of the new throat connecting them, does not exceed the maximum length of throats that are already present in the fine network. The interpore distance , between pores Pf and Pn with radii rf and rn respectively is given by: where , is the Euclidean distance between the two pore centers. For further information, refer to Mufti and Das [14].
The pore network obtained is morphologically equivalent to the void structure of granular soils subjected to stresses and is used to simulate two-phase flow and predict SWRC along drying and wetting paths.

Predicting hysteretic SWRC
For predicting drying SWRC, the pore network is initially assumed to be fully water saturated. The water Fig.3. GSD of the granular soil to be modeled (test soil) and the equivalent GSD for each interval.
reservoir is connected to the outlet pores, while the air reservoir is connected to the inlet pores. The air pressure is gradually increased with water pressure kept constant, which increases the capillary pressure of the network. As the capillary pressure increases, pore-scale invasions occur, and the air starts to invade the pore network elements. Air can enter the pore network elements only when the capillary pressure of the network exceeds the entry pressure (Ψ ) of the element. The entry pressure of an element is given by [16]: where is the drying contact angle, is the water-air surface tension, and is the radius of an element with shape factor . Shape factor of an element is a geometric property associated with its shape and is defined as the ratio of the cross-sectional area of an element to the square of its perimeter.
is a dimensionless correction factor calculated using the equation suggested by Patzek [16]. After every capillary pressure increment, the airwater interface position is updated, and the volumetric water content ( ) of the network at the corresponding suction is computed: where is the water saturation of an element with volume and n' is the porosity of the network.
The maximum capillary pressure obtained after drainage is decreased by progressively raising the water pressure while maintaining a constant air pressure to simulate wetting SWRC. As the capillary pressure decreases, the curvature of the air-water interfaces in the corners also changes. The interface moves when the contact angle reaches the advancing/ wetting contact angle ( ). The highest capillary pressure mechanisms are favored during imbibition, unlike drainage, and hence the displacements that occur when water enters the pore space to replace air are impeded by the pores rather than the throats. Piston-like throat filling, pore body filling, and snap-off processes regulate pore scale invasions along the wetting path. The entry pressure for piston-like throat filling is: where is the radius of curvature of the air-water interface in corners and is evaluated following Patzek [16]. The computation of the entry pressures during the wetting path is more complex as water is already present in the corners after drainage and swells as the water pressure increases. The entry pressures for all the wetting pore-scale invasion mechanisms are thoroughly discussed in Mufti and Das [9]. The simulated SWRC data points are fitted using Fredlund and Xing [17] model: where a, n, and m are the fitting parameters, e is the natural logarithmic constant, is the saturated volumetric water content, is the suction, and is the suction at residual conditions.

Results
The experimental results from Miller et al. [15] are used to assess the performance of the proposed numerical framework for predicting SWRC under stress conditions. Soil composed of 75% ground silica and 25% glass beads with a GSD close to that of fine sandy silt was used in the experiments by Miller et al. [15]. The GSD of the test soil is shown in Fig. 3. Since the test soil is non-plastic, no volume change during wetting and drying cycles and negligible adsorptive forces are valid assumptions. The soil was tested at vertical net normal stresses of 10 and 200 kPa. Given that the net normal stresses were of the order of kilopascals, it is reasonable to assume that no particle breakage occurred at such low stresses. The experimental data for primary drainage and wetting is available for each net normal stress.
While simulating drainage and imbibition, the flow properties like surface tension, drying contact angle, and wetting contact angle are required. The surface tension of water is set as 0.073 N/m. Soil is a water-wet system, and hence the drying contact angle value used is 0º. The wetting contact angle is higher than the drying contact angle and is 20º to 30º in glass beads and as high as 80º in some fine-grained soils [13]. A uniform random distribution of wetting contact angles between 20º and 40º is used in this study.  Fig. 6 and Fig. 7 show the predicted SWRC along the drying and wetting paths at net normal stress of 10 kPa and 200 kPa, respectively. The predicted   SWRCs exhibit an excellent match with the experimental results. The air entry value and residual moisture content increased with the increase in net normal stress. Further, the slope of SWRC is flatter at a net normal stress of 10 kPa compared to higher net normal stress. The predicted SWRCs are consistent with the observed experimental trends. The R 2 and root mean square error (RMSE) values for the predicted drying and wetting SWRC at 10 kPa and 200 kPa are given in Table 2. High R 2 values and low RMSE is obtained, indicating high accuracy of the proposed numerical framework.

Conclusions
A combined DEM-PNM framework to model SWRC of granular soils under stress conditions is presented in this work. Granular soils undergo a change in packing when subjected to stresses. This study proposes that a grainbased model of granular soils subjected to stresses can be obtained by generating a synthetic packing with GSD and porosity identical to that of the stressed granular soil. The model is validated using the experimental results of granular soil subjected to the net normal stress of 10 kPa and 200 kPa. Drying and wetting SWRCs are simulated and compared with the experimental values. The predictions are in good agreement with the observed trend, indicating that the numerical framework can be utilized for predicting the SWRC of granular soils when subjected to stresses. Further studies are encouraged to extend the framework to expansive soils.