Gaining insight into the effective stress parameter using pore-scale numerical modelling with the multiphase lattice Boltzmann method

. We investigate the factors influencing the effective stress parameter, 𝜒 , in unsaturated granular soils using pore-scale numerical modelling. The Shan-Chen multiphase lattice Boltzmann method (LBM) is used for this purpose. A procedure is outlined for measuring 𝜒 using LBM data. A full hydraulic cycle is simulated for a granular soil column and individual contributions of suction and surface tension forces to 𝜒 are measured as a function of the degree of saturation, S r . It is observed that while 𝜒 = 0 at S r = 0 and 𝜒 = 1 at S r = 1 , 𝜒 > S r for 0 < S r < 1 . The trend of the 𝜒 - S r curve is explained based on the trend in the individual components. Schematic plots are utilized to facilitate the understanding of the trends in the data. It is also observed that the 𝜒 - S r relationship is hysteretic; particularly, 𝜒 is larger during imbibition compared to drainage, for a wide range of S r , due to the more pronounced contribution of surface tension forces.


Introduction
Effective stress in unsaturated soils was first introduced by Bishop [1] as ′ = − + ( − ), where is the total stress, and are pore air and pore water pressures respectively and is the effective stress parameter, controlling the contribution of suction, − , to the effective stress, ′. For many years, it was debated whether effective stress can be used as a single stress variable for characterizing unsaturated soils similar to saturated soil [2][3][4]. Eventually it was concluded that while unsaturated soils require two independent stress state variables for a full description of their mechanical and hydraulic behaviour [5], many problems can be solved using only the effective stress [6]. As a result, in recent years, attention has been given back to effective stress for unsaturated soils. The particular difficulty with using effective stress in unsaturated soils is determining the effective stress parameter. Experimental investigations have shown many different variations of with the degree of saturation, Sr, and have provided empirical relationships [7], while theoretical solutions have suggested = Sr, if the interfacial effects are ignored [8], and χ as a function of Sr and air-water interfacial area, if the interfacial effects are accounted for [9]. Simply using = Sr is convenient but it does not properly align with the experimental data, while using the more complete formulation requires knowing the interfacial area which is not easy to measure.
In this study, we take a closer look at the effective stress parameter using pore-scale numerical modelling. Using numerical modelling, we are able to gain insight * Corresponding author: reihos@vt.edu into the factors influencing and the extent of their influence. We are also able to monitor the change of during drainage and imbibition (i.e. drying and wetting) to investigate the effect of hydraulic history.

Measuring χ using multiphase LBM
We use the multiphase lattice Boltzmann method (LBM) as the pore-scale numerical modelling technique [10]. Particularly, we use the Shan-Chen method [11,12] for single-component multiphase flow, meaning that two phases, i.e., gas and liquid, of the same component are modelled, e.g., water and water vapor. Only the fluid domain is simulated in this study and the soil grains are kept immobile, implying the assumption that the capillary forces are not large enough to cause grain movements.
In an LBM simulation, the domain is discretized with a grid called the lattice, consisting of lattice nodes; see for instance the simulation of a liquid bridge between two solid grains in Figure 1. At every timestep, fluid properties, including density, velocity and pressure, are updated for each lattice node in the fluid region, based on the imposed conditions. Using the nodal fluid densities, each fluid node is categorized as either liquid (high density) or the gas (low density). Refer to [13] for a detailed description of the numerical method. The goal is to use this information to find the effective stress parameter.
The first step is to find the force applied by the fluid to each grain. The total fluid force applied to a grain consists of two components, force due to suction, ∆ , and force due to surface tension, . These forces are shown schematically in Figure 2. The suction is applied over the surface of the grain, while surface tension is applied at the contact line between the three phases, i.e., solid, liquid and gas. We implement the following strategy to find these forces using the lattice data.  Schematic of fluid forces applied to a solid grain. is the total fluid force, ∆ is the force due to suction and is the force due to surface tension.
For calculating the suction force, the nodes that are on the surface of the grain are identified; see purple nodes in Figure 3a (note that only half of the spherical grain is shown for visualization purposes). Using these nodes, a Voronoi tessellation is performed on the surface of the grain to find the associated increment of area for each node, Figure 3b. The pressure at each node is then multiplied by the associated area to find the magnitude of the increment of force at that node. The direction of each force increment is found as the direction normal to the surface. The force increments are then integrated, while accounting for their directions, over the entire surface of the grain, to find the total force applied to the grain due to suction.
For calculating the surface tension force, the nodes that are on the contact line between the three phases are identified; see green nodes in Figure 4. For each contact line node, half the distance between its two contact line neighbours is taken as its associated length. The associated length is multiplied by the surface tension (force per unit length), which is known based on the parameters used in the numerical simulation, to find the magnitude of the increment of surface tension force at that node. The direction of the force is found by taking the cross product of the vector tangent to the contact line, , and the vector normal to the surface of the grain, , to get the vector perpendicular to both, , and then rotating by the known contact angle, θ, in a plane perpendicular to the contact line, to get the final direction, . The force increments are then integrated over the entire contact line, considering their directions, to get the total surface tension force applied to the grain.  The next step is to use the fluid forces applied to each grain to find the change of stress in the soil skeleton. Consider the unsaturated granular soil column in Figure  5, where the fluid force applied to each grain is known. For the case where there is no change in overburden pressure on the soil, the change of effective stress on plane XX′ will be only due to change of fluid forces with saturation. Therefore, we can write .
(1) On the other hand, based on the original effective stress formulation we have (2) Since is the total stress which already includes the air pressure from the top, − is basically the stress due to the weight of the soil, which for this case will be constant. Therefore, (3) Comparing equations (1) and (3) we get .
(4) If we take Sr = 1, at which there is no suction, as the baseline case and ignore buoyant forces to only focus on capillary forces, we can write at any other Sr and find and each step of the simulation using

Fig. 5.
Schematic soil column. is the force applied by the fluid to grain j due to suction and surface tension. is the contact force applied to grain i at the base.

Simulation Setup
We use the granular soil column shown in Figure 6a, consisting of 1155 spherical grains with a grain size distribution shown in Figure 6b and a porosity of 0.35. The boundary conditions are closed at the bottom, open to air at the top, and periodic horizontally. A full hydraulic cycle is simulated for this soil column; starting from full saturation the liquid is drained from the bottom until an almost dry state and subsequently the liquid is injected back until full saturation. The simulation is performed as a volume-controlled test, where the saturation is changed by a certain increment and the resultant suction is measured after equilibrium is reached. The resulting soil-water characteristic curve (SWCC) for this specimen is presented in Figure 7. At each step of the drainage/imbibition process, is calculated using equation (6). In addition, individual components of due to suction, , and surface tension, lg , are calculated by separately summing the suction and surface tension forces as follows:

X as a function of degree of saturation
The results of measured , , and lg as a function of degree of saturation for both drainage and imbibition are shown in Figure 8. Note that the degree of saturation used is the saturation of the plane at the base of the column rather than the saturation of the entire soil column. The reason is that the method used in this study for measuring only applies to the effective stress at the base of the soil column not the average effective stress throughout the specimen, therefore, it correlates with the saturation at the base not the overall saturation. The results indicate that = 0 at Sr = 0 and = 1 at Sr = 1, as expected theoretically, while > Sr for 0 < Sr < 1 during both drainage and imbibition. The results also suggest that the -Sr relationship has hysteresis, independent of the suction hysteresis shown in Figure 7. Particularly, during imbibition is higher than drainage, for a wide range of Sr. To understand the trend of the -Sr curve and find the reason for its hysteresis, individual components that constitute are investigated.  Looking back at the schematic in Figure 2, we see that the suction force on each grain is proportional to the projected area on which the suction acts ( 2 ). Using this knowledge, we can deduce that for an assembly of grains, the contribution of suction forces to , , represents the ratio of the projected area to the total cross-sectional area; see derivation in Figure 9.
for both drainage and imbibition starts at a high value at Sr = 1 and decreases as Sr decreases due to the projected areas decreasing, as shown schematically in Figures 9   and 10. However, during imbibition is lower compared to drainage due to the differences in fluid phase distributions. During imbibition, the gas phase is separated into multiple individual gas zones (clusters), while during drainage there is only one continuous gas phase from the top to the bottom. The gas clusters at the bottom of the specimen during imbibition exert an upward suction force to the grains, corresponding to the projected areas taking a negative sign in the summation and, therefore, the being smaller compared to drainage. Now if we focus on the surface tension force in Figure 2, we see that it is affected by the length of the contact line (2 ) as well as the orientation of the interfacial surface ( ). The contribution of surface tension forces to , lg , is similarly affected by the total length of the contact line and orientation of the interfacial surface, as well as the magnitude of suction itself which appears in the denominator of equation (8). Referring back to Figure 8, lg for both drainage and imbibition starts at a high value at Sr = 1 and decreases as Sr decreases, mainly due to the length of contact line decreasing. There is also a clear change of slope happening at Sr of about 0.1 during drainage and Sr of about 0.2 during imbibition. This is due to the additional contribution of the orientation of the forces. Consider the grain in Figure 11 to be a grain at the bottom of the soil column. At high Sr, the phase distribution will be similar to Figure 11a. As Sr decreases, the liquid level drops and a liquid bridge is formed underneath the grain, Figure 11b. In this transition, the length of the contact line decreases while increases, and the counteracting effects of the two results in the variation of lg happening at a low rate. Eventually, the liquid level drops enough for all the grains at the bottom to have a bridge underneath, which is referred to as the pendular regime. At this point, further reduction of Sr would require the bridges to shrink, Figure 11c. In this transition both the contact line length and decrease, therefore, the decrease of lg happens at a higher rate. This change of slope occurs at a higher Sr during imbibition compared to drainage (0.2 versus 0.1) due to the particular phase distribution during imbibition, characterized by a multitude of gas clusters and liquid bridges [13]. Another important observation from Figure 8 is that lg is much higher during imbibition compared to drainage. This is partly due to surface tension forces being larger during imbibition and the other part due to suction being smaller. The lower suction during imbibition is due to the well-known hysteresis of the SWCC, the source of which is discussed in [13]. The reason for the larger surface tension forces can be explained again using the phase distributions. Referring back to Figure 9b and 10b, we see that during drainage, the surface tension forces, shown in green arrows, mainly push outwards from the gas finger that is connecting the top to the bottom, whereas during imbibition, the gas clusters at the bottom pull the grains down and therefore, result in a larger downward total surface tension force.
As a result of summation of and lg , deviates from = Sr as Sr decreases, reaches the maximum deviation when the number of liquid bridges at the bottom are maximized, and converges back to = Sr as Sr approaches zero. This general trend is similar for both drainage and imbibition, however the magnitudes differ. While imbibition is slightly smaller compared to the drainage , imbibition lg is much larger compared to drainage lg for Sr > 0.1, resulting in the overall to take a larger value during imbibition in this Sr range. At Sr < 0.1 there is no hysteresis in -Sr because the base plane is in a similar pendular regime during both drainage and imbibition.
To ensure that our results are consistent with what is expected theoretically, compare the measured with the theoretical solution by Nikooee et al [9], = + Δ (10) where is a material coefficient and is the wetting-nonwetting (liquid-gas) interfacial area per volume. We use = lg , per the suggestion of Likos [14] and measure based on the simulation results of the phase distributions. The comparison is shown in Figure 12. We have used both for the entire soil column and for the plane at the base of the column (which corresponds to an interfacial length rather than area). We see a great match between the measured and calculated results for the latter case.

Conclusion
We measure the effective stress parameter, , as a function of degree of saturation, Sr, at the base of a granular soil column during a full hydraulic cycle, using numerical simulations with the multiphase lattice Boltzmann method (LBM). We find that is equal to one at Sr of one, gradually deviates from = Sr as Sr decreases, with the maximum deviation happening where the base plane enters the pendular regime, and converges to zero as Sr approaches zero. We also see that the -Sr relationship is hysteretic, with being larger during imbibition compared to drainage for Sr > 0.1. We find the reason for this hysteresis in the individual contributions of suction and surface tension forces to . While the contribution of suction forces to is smaller during imbibition compared to drainage, the contribution of surface tension forces is much larger, resulting in total to be larger during imbibition. The larger surface tension forces during imbibition are due to the differences in phase distributions between drainage and imbibition, particularly the presence of individual gas clusters during imbibition which do not exist during drainage. At Sr < 0.1 no hysteresis is observed since during both drainage and imbibition the base plane is in a similar pendular regime.
This study shows the potential of using pore-scale numerical modelling for finding a better understating of the effective stress in unsaturated soils. Future work will focus on improving the simulation capabilities to overcome the current limitations. Particularly, in this study, the grains are immobile, measurements are made only at the base of the soil column, and the soil is granular. Coupling the current LBM code with a discrete element method (DEM) code will allow grain movements to be modelled as well as measurements of effective stress to be made at any location and any direction in the specimen. It will also allow improvements to be made to the shape of the grains and the interparticle forces to be able to simulate finegrained soils.