Simulation of hysteresis in soil–water retention with a correlation between the invading pressure to the governing effective flow area

In order to determine the soil–water retention (SWR) behavior of a particulate medium, the invading phase pressure in the inter-particle level is correlated to the governing effective pore area in the wetting and drying paths. In a three-phase medium that consists of air, wetting fluid and solids, the invading phase on the drying path is air, whereas on the wetting path the wetting fluid advances into the cavities. On a drying path where the area of a cavity is minimum, the air-entry pressure (AEP) of a pore throat is determined by numerically solving the Young–Laplace curvature equation. This can be done using the finite difference method and Newton–Raphson (Jacobian) approximation technique. Next, a relation between the pore area and the value of AEP is developed by varying the distance between solids around the pore throat. Similarly, the water-entry pressure (WEP) is correlated to a maximum pore area of cavity. After packing the particulate domain with the given particle size distribution (PSD) and void ratio values, the primary/main drying and wetting paths of the wetting fluid are simulated and the effect of hysteresis in SWR is shown. It is considered that the total suction equals to matric suction value and the water bridges between two adjacent particles are formed in the form of pendular rings. In this study, the considered material is non-plastic and the shrinkage and swelling during the drying and wetting phases or any change in pore structure are neglected. The simulation results are compared to experimentally determined as well as estimated data from the literature and a great agreement between the results is found, which offers a reliable way around conducting tedious and expensive SWR tests.


Introduction
The known assumption that soil will be fully saturated or dry leads to unrealistic results which are not truly valid in many geotechnical applications. In many applications, such as tunnels, retaining walls, shallow and deep foundations, the stabilization of excavations and slopes, linings and covers of waste containment facilities, the soil is generally unsaturated and is composed of a three or multiphase system. Therefore, the soil suction (negative pore water pressure due to the soil, water and air phase) should be implemented in the calculations of soil mechanics. The basic relationship between the applied suction and water content, which has been termed the soil-water retention curve (SWRC) or also known as the soil-water characteristic curve (SWCC), is used to define the behavior of unsaturated soils [27,35]. The following parameters are used to describe the drying path of SWRC: (1) the volumetric or gravimetric water content at full saturation, (2) the air-entry value or pressure (AEP) at the suction point where air enters the pores in the drying process [9], and (3) the residual water content, when the lower water content threshold has no significant or extremely slow change in the soil water content.
Different experimental techniques have been developed for determining the SWRC, such as the hanging column, pressure chamber, capillary column [3,10,20,34] or SWRC monitoring methods [15]. However, these tests are timeconsuming and dependent on operator skills. In order to overcome these obstacles, numerous empirical methods have been developed. These empirical methods can be separated into two main approaches: (1) estimation methods, which use the data from conventional laboratory tests, such as the particle size distribution (PSD), and form the pedo-transfer 1 3 44 Page 2 of 15 functions (PTF) to predict the SWRC [4,12,28,37] and, (2) curve fitting methods, which interpolate or extrapolate the entire curve by fitting a predefined or assumed equation form on limited experimental results [6,11,13,26,27,36,38]. The estimation techniques don't require any experimental SWR points and solely depend on physical properties of soils, such as PSD, void ratio, and density. With the implementation of genetic programming and neural network analysis and training the algorithm with the existing experimental SWR data, the governing general relations to predict the SWR are proposed [18,19,25]. Toker [33] and Sakaki et al. [29] have suggested empirical equations to determine the AEP of the particulate medium using D 10 , D 30 , and D 50 . However, the outcome of these relations is not accurate and cannot encompass the wide range of possible particle size distributions and properties. Zhou et al. [39] and Suh et al. [32] determined the capillary pressure in irregularly shaped pore throats, whereas Mason and Morrow [21] unveiled the capillary behavior of perfectly wetting liquid in irregular triangular tubes. For non-wetting fluid, the capillary pressure for different liquid-solid contact angles and packing arrangements of spherical particles in 2D and 3D is calculated [2,16,22]. When the energy is minimized, a relation between the invasion pressure and contact angle as well with the packing factor is developed and shows that the capillary pressure for main meniscus is different than that for the water bridges [23]. In a novel apporach and while using X-Ray computed tomography (XRCT) images, Gueven et al. [14] evaluated the porous medium, particle sizes, pore throats as well as the permeability properties of sintered glass beads.
An attempt by Sattari and Toker [30] to model the drying path of SWRC of non-plastic particulate geomaterial resulted in a precise estimation of experimental data. In their method, after packing a mesoscale medium with spherical particles that match a given PSD and void ratio (e), the porous medium is identified and transformed into a pore network model. The AEP of pore throats at their narrowest cross-section between three particles are calculated by fitting a finite difference mesh and solving the Young-Laplace equation for the air-water interface by using the Newton-Raphson (Jacobian) approximation.

Methodology
In the implemented model, the pore water volume is divided into two main volumes: (1) bulk water that are trapped between particles and surrounding pore throats, which are the connections between the adjacent pores of the network, and (2) water bridges in the form of pendular rings which are formed between pairs adjacent particles. In the drying path, the wetting fluid that had a liquid-solid contact angle of less than 90° is driven out of pore network with the assumption of equal suction pressure on meniscus and water bridges throughout the model domain. Finally, the hysteresis in SWRC of non-plastic soil with no change in pore structure during the drying and wetting path is simulated and compared to the existing experimental data from Pham [27]. The water retention behavior of a particulate medium is simulated by applying a correlated relation between governing flow area and the invading pressure of the SWR.
In order to determine the SWRC behavior of the particulate medium, initially, the spherical particles are packed in the domain based on the PSD and are then validated according to the target void ratio. Afterward, the pore network is identified by using the Delaunay Triangulation and subsequently the pore water volume is divided into water bridges (pendular rings) and bulk water volumes. The nonlinearities in the convex surface of pendular rings under applied external suction pressure is neglected and the volume of water bridges with applied suction pressures in the main meniscus is determined. By solving the Young-Laplace curvature equation for incrementally increasing pressure differences using the Finite Difference Method (FDM) and Newton-Raphson approximation the AEP of pore throats is determined. A relation between the AEP and the governing area of cavities is developed, which later is utilized in the simulation of drainage and wetting processes, this diverges from Sattari and Toker [30] which solved the AEP of every pore throat separately.

Packing algorithm
The packing of spherical particles is done in the mesoscale by discretizing the PSD into smaller segments and then determining the number of the particles (N i ) and their diameter (D i ) in each segment.
where i = 2 to n (number of the divided segments) and M is the mass fraction, or percentage of particles corresponding to the sieve diameter (D). N i is the number of particles retained on sieve size D i (i.e. particles between sizes D i-1 and D i ). N 1 is an input of the algorithm and is the number of particles in the coarsest segment of the PSD. By using the periodic boundary condition (PBC) in the boundaries of the medium, the effect of the boundaries is eliminated [17,24]. With the assumption of PBC, the packed mesoscale medium is extended by the maximum grain diameter in X and Y directions. The medium is packed by dropping each particle from the upper boundary and attaching it to the previously packed particles. A search for the number of contact possibilities is carried out and a new particle can reach a stable condition while forming three contacts with already packed particles. One of the inputs of the algorithm is the void ratio of the soil, which has a large influence on final SWRC and AEP results. The density of the packed medium is controlled mainly by defining the contact possibility limit for each particle. When each newly arranged particle is packed with a higher contact numbers, a compact medium is formed, and vice versa. By iterations of the packing process with different contact numbers, target void ratio is reached. In order to pack a dense medium, more iterations for the particle arrangement should be performed and which increases the computation time. A loose packing for 2000 particles (uniform) will take less than 5 min on an average PC. However, for a dense packing with the same number of particles, the packing time increases by a factor of two.

Pore structure
The water volume ( V w ) in a pore network is divided into two separate volumes: (1) Bulk water volume ( V w b ) and (2) volume of pendular rings ( V w p ). The volume of adsorbed water has been neglected in this research, as it has an insignificant effect on both the final SWR results and the residual water content. The bulk water pores are determined by using the Delaunay Triangulation in 3D. Therefore, a bulk pore is defined between every 4 particles and with this identification process, the connectivity and volume of bulk pores is determined (Fig. 1).
As an alternative method for the packing and porous medium algorithm, the open source software Yade [31] can be used for packing and C++ library CGAL [5] can be used for the Delaunay triangulation procedure and pore network model. After the drainage of the bulk pore, the water bridges between two adjacent particles will be formed due to capillary action. The volume of the water bridges between two particles is determined using the model of Chen et al. [8]. In this work, the water bridges (i.e. pendular rings) are considered to have a torus shape, which is formed between two unidentical spherical particles (Fig. 2).
The volume of pendular rings ( V w p ) is, where v 1 and v 2 are, where R 1 and R 2 are principal radii of curvature at the pendular ring and 1 and 2 are the filling angles on each particle, depending on the applied suction. According to the Young-Laplace equation (Eq. 5), the value of filling angles can be determined through geometry.
where R 1 and R 2 are, where d is the distance between two particles shown in Fig. 2. When the bulk water of a cavity is drained the pendular rings will be formed. According to Eq. 2 through Eq. 7, the volume of pendular rings depend on surface tension, applied suction, particles' radii, inter-particle distance and particle-water-air contact angle. Pendular rings between non-contacting particles will be drained when applied suction is greater than what their geometry requires according to the Young-Laplace equation (R 1 < 0 or R 2 < 0), whereas rings between particles in contact will shrink to match increased suction. For a drying path, the packed mesoscale medium is considered to be a fully saturated pore network and the water-air interface is located at the upper boundary of the packed medium. With porous medium identification, the connectivity of all pores to each other is determined and a pore network is formed. When an applied suction is higher than the AEP of pore throat, the connected bulk pore is drained and new water-air interfaces are formed and identified. The applied suction ( ) will increase when the AEP of all pore throats are higher than the current suction.

Assessing the AEP with finite difference method
In order to determine the drying and wetting path of the SWRC, the AEP and water-entry pressure (WEP) of the invading phase in the micro pore scale should be determined. The process for assessing the AEP through a single pore throat is explained here and similar steps can be taken to determine the WEP. According to Sattari and Toker [30], the AEP between three packed particles can be calculated by solving the Young-Laplace equation, using the Newton-Raphson (Jacobian) approximation and FDM. A grid is laid between three particles through their center points (Fig. 3a). The area of pore throat (A p ) is illustrated in Fig. 3b as the checkered, blue surface trapped between the triangle passing through the centers of the particles. This is the narrowest cross-section of a pore throat between three particles. The grid density, which is defined as the minimum number of grid points located on a grid line inside the pore throat area in Fig. 3b, is 100. This parameter controls the mesh size of the grid and its effect on calculated AEP is given in Sect. 2.4. Three types of grid points, with coordinates (x G ,y G ,z G ), have been identified: (1) curvature points, representing the Young-Laplace curvature equation (Eq. 5) located in pore throat, (2) boundary points, which are located inside the particles along the contact of water-air interface with particles and should satisfy the particle-water-air contact angle, and (3) the points which are located outside of the triangle formed between three particles, which will be removed from calculation process. The only unknowns are the z coordinates of the grid points (z G ). For the boundary points (Eq. 8) in an equilibrium condition f z m,n G = 0 , which means the contact angle between the air-water surface and the particles is satisfied. where m, n = 1: number of grid points, (x, y, z) j is the center coordinates of particle i, (x G ,y G ,z G ) m,n are the coordinates of grid points representing the air-water interface, z m,n G X and z m,n G Y are first-order partial derivatives in X and Y directions determined from finite difference approach. For curvature points representing the Young-Laplace equation (Eq. 5) the following condition should be satisfied.
where z m,n G XX , z m,n G YY are second-order partial derivatives in X and Y direction and z m,n G XY is a mixed partial derivative. The effect that the particle-water-air contact angle has on final SWR results is illustrated in Fig. 4 [30]. In this drying simulation, the particles diameter range is from 0.125 to 0.150 mm and the porosity of the representative volume is 0.62. The gravimetric and volumetric water contents can be easily calculated based on specific gravity of solids, dry density and void ratio based on experimental data.
With the application of the Newton-Raphson (Jacobian) method, where t is the Newton-Raphson iteration number and J −1 is the inverse of Jacobian matrix. Solving this equation allows us to have z coordinates of grid points that correspond to applied suction (i.e. it satisfies the Young-Laplace equation at all points of the air-water interface). Note that the Dirichlet boundary condition, which grants the continuous curvature of the air-water interface on sides of the triangular mesh at the pore throats throughout the packed bed medium, is not implemented in this paper. In the case of asymmetric pores and particles as are presented here, the implementation of Dirichlet boundary condition would be too demanding and would require an immense amount of computational power.

Effect of contact angle
The contact angle effect on air-entry pressure (AEP) values is small and the main difference observed is in the volumes of water bridges calculated during the drainage process. So, in this research, the contact angle effect for determining a relationship between the AEP and the area of an inter-particle pore throat is neglected due to the numerical difficulties which are discussed in Sattari and Toker [30]. However, the calculated volume of water bridges as described in Sect. 2.2 are still dependent on the liquid-solid contact angle values.

The correlation between AEP and area of pore throat
The AEP of a pore throat not only depends on its area size (as defined in Fig. 3b) but also on the diameter of the particles forming the pore throat. The first step in finding the correlation between AEP and the area of the pore throat is to determine the particle range where the assumption of spherical particle shapes is valid. This range is considered to include coarse to silt-sized particles, avoiding clay particles due to their plate-like shape and surface chemistry. Therefore, the diameter of the particles used in this simulation is limited to being between 1 to 0.002 mm. Three particles are packed together and by varying the distance between them, different pore throat areas are formed. The AEP of the pore throat between the particles is determined as described in Sect. 2.3. Figure 5 depicts the effect of the mesh size on the AEP values which are determined for a pore throat formed between three uniform particles with radii of 0.5 mm.
According to the results, a grid density of 50 is considered to be adequate for eliminating the mesh size effect.
The z coordinate of each point that satisfies the equilibrium at each suction increment is determined according to Eq. 10. Figure 6 illustrates the calculated geometry of water-air interface at four suction values for a pore throat which is formed between particles with the radius (R) of 0.010 mm. Figures 6a-c depict the air-water interface at equilibrium for the indicated suction values up to 15 kPa. Figure 6d shows the solutions convergence failure at 16.1 kPa suction, which is then inferred to be the AEP. The procedure is repeated for different particle sizes ranging from 1 mm to 0.002 mm. When a pore throat is composed of three particles with different diameters, it is assumed that the particle with the smallest diameter controls the AEP value. Figure 7 shows the correlation between the AEP and the area of pore throats formed between different particle radii. The smaller the pore throat area is, the larger the AEP. With an increasing area of a pore throat, the effect of particle radius on final AEP value increases. Both AEP and area of pore throats are normalized with the smallest particle's radius, as presented in Fig. 8 as a double logarithmic plot. The double linear version of the same plot, which is not provided here for the lack of visual clarity, hints at an exponential relationship between these two normalized parameters.
where R is the minimum radius among the three particles and A p is the area of the pore throat shown in Fig. 3b. This general equation can be implemented in any pore network model with bedded particles to calculate the AEP of pore throats. This exponential relation between the AEP and area of pore throats (Eq. 11) is implemented in the simulation of the drying and wetting path of SWR in this work. Fig. 4 The effect of contact angle (θ) varying from 0° to 30° on final drying SWR results, where the diameter of particles range from 0.125 to 0.150 mm and porosity is 0.62

The hysteresis in SWR
The hysteresis in SWRC during the drying (desorption) and wetting (sorption) path of the granular medium with wetting fluid is due to 4 reasons as mentioned by Tuller and Or [35]: (1) the ink-bottle effect, which means that the AEP of a pore throat controls the drying path whereas the Water Entry Pressure (WEP) of a bulk pore or cavity controls the filling of pores in the wetting path (Fig. 9a), (2) the trapped air forms cavities during the wetting path (Fig. 9b), where all the pore throats connected to the pore have air-water interfaces, (3) the swelling and shrinkage of the medium during the drying and wetting path, which has a higher influence on plastic minerals, such as clay, and (4) the difference between the air-water-particle contact angle in advancing and receding water menisci.
The simulation of a drying path of SWRC (Fig. 10a) is carried out with the assumption of a fully saturated granular medium, where the air-water surfaces are assigned to the upper boundary of the representative volume. When  . 7 The correlation between the AEP and area of pore throat for different particle radii varying from 0.5 to 0.01 mm the simulated/calculated AEP of a pore throat which has an interface is lower than the applied suction, the connected bulk pore is drained and spontaneously the volumes of water bridges that remain between adjacent particles due to capillary pressure are calculated. Focusing on the first and second explanations for hysteresis that are listed above, the wetting path and hysteresis of SWRC is analyzed by neglecting the swelling and shrinkage of the medium as well as eliminating the effect of contact angle. Therefore, for a dry condition where the only water volumes are the pendular rings that are developed during the initial drying path, the wetting path of SWRC is simulated (Fig. 10b). Initially, the bottom boundary is assigned as water-air interface and the applied suction value is controlled by the WEP of equivalent pore areas. The equivalent pore area is determined by the area of a circle which has the same radius as a sphere that has the same volume as the bulk pore. The equivalent pore area is substituted in the place of (A p ) in order to calculate WEP by Eq. 11.
During the simulation an algorithm controls for the possibility of there being trapped air cavities, which happens when a dry pore is surrounded by saturated pores. A cavity can only be fully saturated or fully drained and the water rise within a cavity is neglected. Figure 11. is a conceptual illustration of the hysteresis during initial drying, main wetting and main drying curves under applied suction pressures. Fig. 8 Normalized AEP versus normalized pore throat area, in log-log scale, for 5 sets of particle radii ranging from 0.5 to 0.01 mm Fig. 9 The hysteresis in wetting and drying paths of SWR due to a ink-bottle effect, and b trapped air cavity during the wetting path

Simulation of the drying path of SWR
In order to validate the model, experimental data from literature is used. The simulation results are compared not only with the experimental data, but also with two estimation methods from the literature [4,12], where physical properties of soil, such as PSD, void ratio, and density, are the inputs of the methods and no fitting algorithm is employed.

Material
By utilizing the literature and the three PSD given in Fig. 12, an attempt at validating the model is made. Figure 13 depicts the packing arrangements of these soils with different void ratios. Table 1 presents all the relevant soil parameters and inputs from the literature which are required for the simulation of SWR. The first set of data is taken from Fredlund et al. [12], where drying SWR for a sand ranging from 1 to 0.1 mm is considered. The target void ratio of the simulation is 0.61. The second set of data is taken from Toker [33], where an experimental drying SWR of New Jersey fine sand (NJFS) (D = 0.6 to 0.02 mm) with sub-rounded particle shapes is given. The void ratio, in this case, is 0.66. The last set of data displays quartz sand ranging diameter from 0.6 to 0.06 mm which is taken from Ahmadi-Adli [1]. The shape of the particles is sub-angular and the void ratio of the soil is 0.94, which is a loose packing. Fig. 10 The flowchart of the simulation steps for a drying path and b wetting path of SWR  Fig. 11 The hysteresis in SWR during initial drying, main wetting and main drying path

Validation according to the experimental data
The medium is packed according to the PSDs given in Fig. 12. The packing algorithm is validated by controlling the simulated void ratio and its comparison to the experimental void ratio (e) given in Table 1. If the void ratio is not within the tolerance of 0.02 of the target void ratio, the packing process is repeated until the tolerance criterion is satisfied. The number of generated particles and resulting void ratios are given in Table 2. Note that due to the application of PBC, the actual number of the particles in the medium is more than 7762. With porous medium identification, the pore network is constructed. Afterwards, the drainage process is applied (Fig. 10a). A suction ( ) is applied and according to Eq. 11 the AEP Fig. 12 The particle size distribution (PSD) of three soils from the literature Fig. 13 The packed domain for soils from a Fredlund et al. [12] (N Particles = 2833), b Ahmadi-Adli, [1] (N Particles = 8846), and c Toker [33] (N Particles = 7762), where void ratios are 0.583, 0.942 and 0.667, respectively of all pore throats that contain a water-air interface is checked for a breach. When > AEP of any pore throat, the air is invaded into the connected pore and the volume of water remaining in the pore in the form of pendular rings is calculated according to Eqs. 2 and 5. The drainage process is repeated until the air-entry pressures of all pore throats with an air-water interface are greater than the applied suction (AEP > ). Then, the suction value is increased and the whole process is repeated. The minimum applied suction value is 0.1 kPa, and suction is increased incrementally. The calculation time for a drainage process is given in Table 2. Considering that in traditional SWR tests, a single suction value in an experimental test requires around 1 day until the equilibrium is reached, the SWR simulation is significantly more efficient and time-saving. The calculations are done with the MAT-LAB platform, on a personal computer with Intel Core i7, 3.60 GHz. The drying SWRC of sandy soil given by Fredlund et al. [12] is given in Fig. 14a. The comparison of the simulation  results to the experimental data shows the accuracy of the model. The model is able to capture the AEP, slope of drying SWR and residual water content with high accuracy. The result of the NJFS soil given by Toker [33] is shown in Fig. 14b. The simulation outcome is in a great agreement with the experimental data. The Fredlund et al. [12] estimation method is unable to simulate the slope and residual water content value. The final set of drying SWR results is taken from Ahmadi-Adli [1] on quartz sand (Fig. 14c). The packed medium is very loose and according to the results, the model is able to capture the slope of SWR and residual water content with high accuracy. The existing error in the determination of AEP in comparison to other estimation methods is negligible. The air-entry is defined at the intersection of the tangent lines at the initiation point of the invading air pressure. The residual water content and suction correspond to a water volume and suction at which all bulk water is drained and the only remaining water volume is due to the developed pendular rings in the inter-particle medium. It is defined at the intersection of the tangent lines in the bulk water volume of zero. The quantitative comparison of experimental, simulation and estimation results are presented in Tables 3 and 4.

Simulation of the hysteresis in SWRC
The simulation of hysteresis in SWRC is carried out by using the experimental data of two different soil types. The experimental data is taken from Pham [27], where the hysteresis on sandy and silty samples is investigated. First, the initial drying path of SWRC is simulated by using the steps shown in Fig. 10a. Afterward, the main wetting path of SWRC is simulated according to the flowchart in Fig. 10b, and finally the main drying path is plotted.

Material
The PSD results of Beaver Creek Sand and the Processed Silt samples [7] are shown in Fig. 15. The processed silt is obtained from a natural silt after the reduction of its sand and clay fractions. However, a small fraction of clay minerals (around 7% of total mass) still remains in the sample which causes swelling and shrinkage during the drying and wetting phase. The particle and dry density of the sandy sample are 2.67 and 1.652 (g/cm 3 ) and for the silty sample are 2.67 and 1.726 (g/cm 3 ), respectively.

Validation according to the experimental data
The packed domains for silty and sandy samples are shown in Fig. 16. For the Beaver Creek Sand domain, a total of 15,360 particles are generated which resulted in 96,047 pores and a void ratio of 0.511. Similarly, for the processed silt sample, a total number of 34,720 particles and 217,399 pores are generated and the final simulated void ratio is 0.601.
The experimental results from Pham [27] on Beaver Creek Sand and the Processed Silt samples are given in Fig. 17. In these plots, the hysteresis in SWR during the initial drying, main wetting and main drying paths is observed. In comparison to the initial drying path, the effect of trapped air-bubbles during the main wetting path causes the  . 15 The particle size distribution of Beaver Creek Sand and the Processed Silt samples Fig. 16 The packed medium for a Beaver Creek Sand and b processed silt difference between the water content value in zero suction. However, during the main drying path, where the medium is unsaturated, the difference from main wetting path is primarily due to the ink-bottle effect. In the simulation, this is transferred to the pore throat area or equivalent pore area for the draining or rising water volume, respectively. The simulation results are shown as well in Fig. 17. According to the results of the Beaver Creek Sand: • In initial drying path, the AEP of the domain is well captured. However, the suction pressure corresponding to residual water content is underestimated. The main reason for this is due to the unknown size distribution of 4% mass of particles smaller than 0.06 mm (Fig. 15).
The fine solids will change the pore network structure significantly and due to the smaller pore throat area, the AEPs of these pore throats will be higher (Eq. 11). • According to the wetting path, the WEP in the medium is overestimated. This may also be due to the unknown 4% mass of fine particles as well as the assumption that the pores are fully saturated or completely drained. Therefore, the water rise within a cavity is not calculated unless it is fully saturated. The number of trapped air bubbles are also slightly underestimated as the entrapment of air volumes that occupy multiple adjacent pores is not simulated. The saturation value in zero suction is 0.97. Fig. 17 The hysteresis in a Beaver Creek Sand, and b Processed Silt samples • Similar to the initial drying path, the residual water content in the main drying path is underestimated. Table 5 depicts the quantitative comparison between the experimental and simulation outcomes. Due to a lack of experimental points it is not possible to determine the exact experimental AEP value. • The simulation consistently predicted desaturation rates that are too high (too steep of a slope in the bulk pore drainage regime) for all of initial, main drying and wetting paths.
Similarly, the hysteresis in the processed silt sample is simulated. It is observed that: • The simulation captured the AEP of the domain, slope of drying path and the residual water content well. The liquid-solid contact angle is assumed to be 45° and as it can be interoperated from the result, neglecting the unknown behavior of clay particles resulted in a slightly lower water content in the high suction range compared to the experimental data. • The simulated WEP is in agreement with the experimental data. In contrast to the previous result, the number of trapped air bubbles is slightly overestimated, which can again be due to the pore structure change during the drying and wetting path, especially when considering the swelling and shrinkage in clay minerals. The simulated saturation value in zero suction is 0.903. • According to Table 5, the residual pressure is underestimated, because the presence of clay minerals is neglected in the simulation.
The proposed simulation approach can estimate the AEP, WEP and residual water content, residual pressure as well as the slope of the wetting and drying path of SWR with sufficient accuracy. The MATLAB code for simulating hysteretic SWR is available in *. p format for free download: http:// www.geote chnic s.ifg.uni-kiel.de/en/publi catio ns/downl oads.

Conclusion
In this study, the simulation of the hysteresis in soil-water retention (SWR) behavior of non-plastic granular material is presented. The three-phase domain undergoes drying and wetting cycles, where the wetting fluid with a liquid-solid contact angle of less than 90° advances or recedes in the medium. The solids are considered to be spherical, non-plastic, and swelling and shrinkage is neglected. The domain packing algorithm is based on particle size distribution (PSD) and is calibrated according to the target void ratio. The pore network is extracted by Delaunay Triangulation and the connectivity of pores are determined. The fluid in the multi-phase permeable medium is divided into two main volumes: (1) bulk water trapped in cavities surrounded by 4 irregularly packed particles, and (2) the water bridges between two adjacent particles which are in form of pendular rings. In the drying path of SWR, the air-entry pressure (AEP), air as the invading phase, of pore throats is correlated to the minimum cavity area. On the other hand, the water-entry pressure (WEP) in the wetting path is correlated to the maximum area of a cavity. The AEP and WEP are calculated based on solving the Young-Laplace curvature equation using the finite difference method and Newton-Raphson (Jacobian) approximation. It is shown that the invading pressure of cavity not only depends on its area but also on the size of the particles forming it. After the drainage of a cavity in the drying path, and due to capillary pressure, the water bridges between two particles are formed. Its volume is determined based on applied suction magnitude, the distance between particles and the liquid-solid contact angle. The applied suction on the meniscus is assumed to be equal to that in the pendular rings. The normalized relation between invading pressures and governing area of pores is utilized to determine the drying and wetting paths of SWR and depict their hysteretic nature due to the ink-bottle effect and trapped air bubbles in the wetting path of SWR. The simulation results are compared with experimental and estimation methods, such as Arya and Paris [4] and Fredlund et al. [12], and the accuracy of the proposed method is granted. The algorithm is able to estimate the AEP, WEP, slope of SWR and residual water content of a three-phase packed domain with great accuracy, while eliminating the need for conducting tedious experiments as well as reducing the computational time needed for SWR simulations. The correlation between invading pressure and governing flow area can be implemented in any discrete element method (DEM) to determine the soil-water retention (SWR) behavior undergoing pore structure change.