Investigating the effect of porosity on the soil water retention curve using the multiphase Lattice Boltzmann Method
aa r X i v : . [ phy s i c s . g e o - ph ] F e b Investigating the effect of porosity on the soil water retention curve using themultiphase Lattice Boltzmann Method
Reihaneh
Hosseini , ∗ , Krishna
Kumar , ∗∗ , and Jean-Yves
Delenne The University of Texas at Austin, 301 E Dean Keeton St, Austin, TX 78712 IATE, Univ. Montpellier, CIRAD, INRAE, Montpellier SupAgro, Montpellier, France
Abstract.
The soil water retention curve (SWRC) is the most commonly used relationship in the study ofunsaturated soil. In this paper, the e ff ect of porosity on the SWRC is investigated by numerically modelingunsaturated soil using the Shan-Chen multiphase Lattice Boltzmann Method. The shape of simulated SWRCsare compared against that predicted by the van Genuchten model, demonstrating a good fit except at low degreesof saturation. The simulated SWRCs show an increase in the air-entry value as porosity decreases. The Soil Water Retention Curve (SWRC) correlates thedegree of saturation or water content of unsaturated soilto its matric suction or capillary pressure. Idealized func-tional forms of SWRC [1, 2] are often used to estimatemechanical properties of unsaturated soil such as e ff ec-tive stress and shear strength [3, 4] as well as hydraulicproperties such as relative permeability [5, 6]. As the soilundergoes volumetric deformation during shearing, its wa-ter retention behavior changes with the change of porosity.Many researchers have proposed models that consider thee ff ect of porosity (or void ratio) on SWRC [7–11]. In thisstudy, we investigate the e ff ect of porosity on the suctioncharacteristics of unsaturated soil from a micromechanicsperspective [12], by simulating the multiphase system us-ing the Shan-Chen Lattice Boltzmann Method (LBM). An in-house 3D multiphase LBM code is developed formodelling the multiphase fluid domain of unsaturatedgranular material. The D3Q19 scheme [13] for velocitydiscretization and the BGK collision operator are used.The interaction between the three phases, i.e. liquid, gasand solids, is considered by introducing Shan-Chen (SC)-type interaction forces into the model [14, 15]. The SCinteraction forces per unit volume are calculated using F SC ( x ) = − ψ ( x ) G X i w i ψ ( x + c i ∆ t ) c i ∆ t (1)where ψ is an e ff ective density, G is a parameter that con-trols the strength of the interaction (negative for attrac-tion), { c i } is the discrete velocity set with { w i } as the corre-sponding weights, and ∆ t is the time step which is set to 1 ∗ e-mail: [email protected] ∗∗ e-mail: [email protected] dimensionless lattice unit (lu). ψ is defined as ψ ( ρ ) = s c s ∆ t G ( p ( ρ ) − ρ c s ) (2)where c s is lattice sound speed equal to √ ∆ x ∆ t for theD3Q19 model [13, 16] with both ∆ x and ∆ t set to 1 lu, ρ isdensity at point x , and p is pressure. Equation (2) is a rear-rangement of the Equation of State (EOS) of the SC-typemultiphase system, which allows incorporating di ff erentEOS by simply redefining p [17]. The Carnahan-Starling(C-S) EOS is used to increase the numerical stability forsystems with large liquid-gas density ratios: p = ρ RT + b ρ/ + ( b ρ/ − ( b ρ/ (1 − b ρ/ − a ρ (3)with a = . R T c / p c and b = . RT c / p c [17]. T c is the temperature below which phase separation occursand p c is the pressure at which the first and second deriva-tives of the EOS at T c are zero. The C-S EOS curves fordi ff erent temperatures are shown in Figure 1.The developed code is validated by comparing the co-existence densities, the densities at which the liquid andgas phases coexist in a multiphase system, from the nu-merical simulations against the theoretical solution, at thedi ff erent T / T c values shown in Figure 1. The density ateach node is initialized to a value between the liquid andgas densities, with a small random perturbation to allowphase separation, and the system is allowed to equilibrate.Once a liquid droplet is formed, the average liquid and gasdensities are measured, and the corresponding pressuresare calculated using Equation (3). For this study a = b = R = igure 1. Comparison of the numerical (multiphase LBM sim-ulation) and theoretical coexistence curves for the C-S EOS atdi ff erent T / T c . In this study, we use four di ff erent grain configurationswith initial porosities ( η = V v oid / V total ) of 0.50, 0.40, 0.30,0.15 (corresponding to void ratios of 1.00, 0.67, 0.43,0.18) to investigate the e ff ect of porosity on the SWRC.These configurations are created by randomly positioninggrains inside a fixed-size domain until the required poros-ity is achieved. As the objective of this study is to developSWRCs at constant porosities, the grain positions are fixedand grain-grain interactions are not considered. The gran-ular assembly comprises of grain diameters between 10 luand 30 lu, with a uniform distribution. The lattice spac-ing, ∆ x , is set to 1 lu, which corresponds to a minimum of10 grid points per grain diameter. The LBM domain is setto 200 × ×
200 lu , with periodic boundaries in everydirection to eliminate boundary e ff ects.In all simulations, the relaxation parameter, τ , for theBGK collision operator is set to 1 lu. The parameters forthe C-S EOS are chosen as a = b = R = T c = . p c = . ρ c = . T / T c of 0.7 is selected for these sim-ulations. G is set to -1 (by using Equation (2) for ψ , theparameter G cancels out when this equation is substitutedin Equation (1), therefore, the magnitude of G is irrelevantand only its sign is of importance). The density of solids, ρ s , in a multiphase simulation only controls the contactangle between the solids and the liquid and does not cor-respond to the physical density of the solids; ρ s closer tothe liquid density, ρ l , results in a more hydrophilic surface,whereas ρ s closer to the gas density, ρ g , creates a more hy-drophobic surface [18]. In these simulations, ρ s is set to0.98 ρ l , corresponding to a contact angle of 14 ◦ .The dimensionless lattice units, lu, used for all the pa-rameters in this study can be converted to physical unitsby choosing appropriate conversion factors. It can beshown that the fluid kinematic viscosity in an LBM simu-lation is ν = c s ( τ − .
5) in lu [13]; the conversion factors for length, C l , and time, C t , should be chosen such that c s ( τ − . C l / C t corresponds to the correct physical kine-matic viscosity. In addition to this criteria, accuracy andstability should be taken into account when choosing theseconversion factors. For this qualitative study, all values aresimply reported in lu.In all four configurations of the granular assembly, thedensity of the fluid domain is initialized to a value be-tween the liquid and gas coexistence densities (0.358 luand 0.009 lu, respectively), with a random perturbation of ± .
1, to allow phase separation. After an adequate numberof steps, when the fluid has phase separated and reachedequilibrium, liquid is injected into the system by increas-ing the density of all fluid nodes by 0 .
005 lu every 1000steps. It is ensured that 1000 steps is enough for the sys-tem to reach equilibrium before the next injection. TheSWRC is generated by calculating the average suction, δ p , and degrees of saturation, S r , at the end of each in-jection stage. δ p is defined as the di ff erence between theaverage gas pressure and the average liquid pressure in theentire fluid domain. S r is measured by dividing the totalnumber of liquid nodes by the total number of fluid nodes(voids). In addition, the number of liquid clusters at eachinjection stage is computed using the Depth-First Searchalgorithm [19]. Figure 2 shows a 3D view of the grain configuration with η = S r = ff erent states of liquid clustering: pendular (3a), fu-nicular (3b), capillary (3c), and droplet states [20]. In thependular state, pairs of grains are connected by binary liq-uid bridges (e.g. dashed-circle zones). In the funicularstate, a combination of binary liquid bridges (e.g. dashed-circle zone) and liquid clusters connecting multiple grains(e.g. dashed-triangle zone) is present. In the capillarystate, the liquid almost fills the entire pore space betweengrains, however, its surface forms concave menisci (e.g.dashed-square zones) exerting a net cohesive force on thegrains. In the droplet state, the liquid immerses the entiredomain and exerts no capillary forces on the grains.The pore-scale evolution of clusters with increase insaturation is investigate by examining the evolution ofnumber of liquid clusters. Figure 4 shows the number ofliquid clusters normalized by the number of disconnectedpores in the soil structure ( N cluster / N pores ) as a function of S r , for di ff erent porosities. The circles indicate points atwhich snapshots of the simulation are provided in Fig-ure 3. As S r increases from zero, N cluster / N pores initiallyincreases rapidly up to a peak value, due to the formationof a large number of small liquid droplets on the surface ofthe grains. With further increase in saturation, these smalldroplets aggregate together to form liquid films aroundthe grains and subsequently coalesce to binary bridges be-tween grains, thereby decreasing the number of individualclusters at a steep rate. The number of distinct clusters igure 2.
3D view of the multiphase domain with η = S r = gradually decreases as the system moves from the pendu-lar state to the funicular state, and finally, to the capillarystate. Eventually, all the pores are filled with liquid, andthe cluster count becomes equal to the number of discon-nected pore spaces.As the porosity decreases, the peak of the curvesin Figure 4 shift to the right, indicating that a higher S r is required for the initial clusters to connect and form liq-uid films on the grain surfaces. This is expected as lowerporosity corresponds to a larger number of tiny pore spacesand a higher surface area of the granular assembly. Inaddition, in the case of the lowest porosity, the numberof distinct clusters decreases more gradually throughoutthe transition from the pendular state to the capillary state,compared to the other porosity cases.The SWRCs simulated for di ff erent porosities areshown in Figure 5. The simulated SWRCs are com-pared against fitted curves based on the simplified vanGenuchten model [2, 8] S r = ( 11 + ( αδ p ) n ) m (4)where α , n and m are fitting parameters. The fitting pa-rameters for the curves in this study are reported in Ta-ble 1. Parameter α controls the shift of the curve to left andright and is inversely proportional to the air-expulsion / air-entry value (AEV). It can be deduced from the shift ofthe simulated curves in Figure 5 or the value of α in Ta-ble 1, that the AEV for the simulated SWRCs increaseswith decreasing porosity. The inverse relation between theAEV and porosity is consistent with all proposed modelsfor the SWRC which take void ratio into account [7–11].Parameter n controls the slope of the curve and is relatedto the pore size distribution, with higher n correspondingto a more uniform distribution [21]. Comparing the n pa-rameter for the di ff erent curves shows that the curves withthe highest and lowest porosities have steeper slopes at themiddle section of the SWRC, compared to the curves withintermediate porosities. This could be explained by as-suming that the lowest and highest porosity models haveuniform small and large pores, respectively, while the in- aaaabc Figure 3.
Snapshots of 2D slices of the multiphase domain with η = S r = S r = S r = termediate porosity models have a combination of smalland large pores; however, the verification of this assump-tion is left for future studies. Parameter m controls thesymmetry of the curve when δ p is plotted in log space. Asseen in Figure 5, the maximum suction of the simulatedcurves is constrained, creating an unsymmetrical shape forthe SWRC at lower porosities. Therefore, m is higher forlower porosities.The SWRC for the simulation with the highest poros-ity ( η = igure 4. Evolution of normalized number of clusters as a func-tion of degree of saturation, for di ff erent porosities. The circlesindicate points at which snapshots of the simulation are providedin Figure 3 Figure 5.
SWRCs from the LBM simulations for di ff erentporosities. The dashed lines represent the fitted curves based onthe van Genuchten model. Table 1. van Genuchten fitting parameters. η α n m S r range of 0.1 to about 0.8. However, at higher satura-tion ( S r > S r = η = S r occurs as expected, however,at S r below 0.2, the suction remains unchanged. Interest-ingly, all simulated curves converge to the same value ofmaximum suction, indicating that the maximum suction inthese multiphase LBM simulations may be constrained bythe modeling parameters. Further research is required tounderstand the cause of this behavior. The e ff ect of porosity on the Soil Water Retention Curve(SWRC) was studied by simulating unsaturated soil us-ing the 3D multiphase Shan-Chen-type Lattice BoltzmannMethod. Four di ff erent grain configurations with porosi-ties ranging from 0.15 to 0.50 were used. The four di ff er-ent states of liquid clustering, namely pendular, funicular,capillary, and droplet states, were observed in the multi-phase simulation for the di ff erent porosity cases. The sim-ulations were able to correctly capture the increase of theair-entry value with decreasing porosity and the shift ofthe SWRC towards higher suction. However, the resultsindicated that the maximum suction is constrained at lowdegrees of saturation and further research is needed to ex-plain this behavior. References [1] R.H. Brooks, A.T. Corey, Hydrology Papers 3 (1964)[2] M.T. van Genuchten, Soil Sci. Soc. Am. J. , 892(1980)[3] N. Khalili, M.H. Khabbaz, Géotechnique , 681(1998)[4] N. Lu, W.J. Likos, J. Geotech. Geoenviron. Eng. ,131 (2006)[5] Y. Mualem, Water Resour. Res. , 513 (1976)[6] D.G. Fredlund, A. Xing, S. Huang, Can. Geotech. J. , 533 (1994)[7] M. Aubertin, J.F. Richard, R.P. Chapuis, Can.Geotech. J. , 55 (1998)[8] D. Gallipoli, S.J. Wheeler, M. Karstunen, Géotech-nique , 105 (2003)[9] A. Tarantino, Géotechnique , 751 (2009)[10] D. Masin, Int. J. Numer. Anal. Meth. Geomech. ,73 (2009)[11] A.N. Zhou, D. Sheng, J.P. Carter, Géotechnique ,669 (2012)[12] V. Richefeu, F. Radjai, J.Y. Delenne, Comput.Geotech. , 353 (2016)[13] Y.H. Qian, D. D’Humières, P. Lallemand, EPL ,479 (1992)[14] X. Shan, H. Chen, Phys. Rev. E , 1815 (1993)[15] X. Shan, H. Chen, Phys. Rev. E , 2941 (1994)[16] X. He, L.S. Luo, Phys. Rev. E , 6811 (1997)[17] P. Yuan, L. Schaefer, Phys. Fluids18