Numerical Modeling of the Effect of Desaturation on Liquefaction Hazard Mitigation

Abstract

Earthquake-induced liquefaction is always a concern when the soil near the surface of a site is composed of relatively loose saturated sand. One of liquefaction mitigation methods is to induce gas bubbles into the deposit to reduce the degree of saturation. A coupled pore-scale model is presented herein to investigate liquefaction resistance of desaturated granular materials. The multiphase fluid, which mimics the behavior of air and water, is modeled using the multiphase single component lattice Boltzmann method. The solid phase is modeled using the discrete element method. The coupled framework was utilized to study the behavior of a soil deposit with the different degrees of saturation of 100%, 92%, and 82% during an earthquake loading. Based on the results of the simulations performed, liquefaction occurred in the fully saturated granular deposit and was not observed anywhere at depth in the desaturated deposits. It has also been found that reducing the saturation level from 100% to 92% significantly affects behavior. In desaturated deposits, higher average coordination number, lower pore pressure buildup, and slower effective stress decay were observed compared to fully saturated deposits. However, it turned out that a further reduction in the degree of saturation from 92% to 82% does not have a significant impact on the calculated parameters.

1. Introduction

Earthquake-induced site liquefaction has been associated with very costly damage to port facilities, bridges, dams, buried pipes, and buildings of all types. Liquefaction occurs when soil has lost its strength and behaves like a liquid. In this state, the resistance of the ground drops below the level required to remain stable. During cyclic loading, saturated loose sand tends to compact overall, resulting in pore pressure buildup as void space decreases. The increase in pore pressure causes the grain contact forces to relax and eventually disappear, indicating a state of zero effective stress and complete loss of soil strength. Some of the most dramatic cases of liquefaction were observed after the 1964 earthquake in Niigata, Japan, and the 1964 earthquake in Prince William Sound, Alaska that defined the phenomenon of liquefaction. Examples of these events include the 1989 Loma Prieta (California), the 1995 Kobe (Japan), the 1999 Kocaeli (Turkey), the 2011 Christchurch (New Zealand), the 2011 Tohoku Oki (Japan), the 2012 Emilia (Italy) and the 2018 Sulawesi-Palu (Indonesia) earthquakes.

Different site mitigation techniques have been developed to decrease the risk of liquefaction in saturated sands. Dynamic compaction (e.g., [1,2,3,4,5]), different grouting methods (e.g., [6,7,8,9]), deep soil mixing (e.g., [10,11,12]), and gravel drains (e.g., [13,14,15,16]) are some of the widely used liquefaction mitigation techniques. More recently, biologically induced cementation (e.g., [17]), as well as desaturation methods, have also been presented as new techniques to mitigate liquefaction.

One of the advantages of using the induced desaturation methods to mitigate liquefaction risk is its feasibility to be used beneath existing facilities. Martin et al. [18] observed that a very small reduction in the degree of saturation decreased the risk of liquefaction effectively. Based on their studies, it was observed that one percent degree of saturation decrement of a fully saturated sand specimen, where the porosity was 40%, would cause a 28% decrease in the pore water pressure increase during a cycle. Chaney [19] and Yoshimi et al. [20] conducted studies which confirmed the significant effect of small degradation in degree of saturation on liquefaction mitigation. Based on their observation, when the degree of saturation was reduced to 90%, liquefaction resistance was found to be twice that of fully saturated conditions.

Okamura et al. [21] found that the induced desaturation technique is an effective method to mitigate the risk of liquefaction and that the desaturated state can persist for a long time, typically more than a decade. Yegian et al. [22] used electrolysis to confine gas molecules in saturated sand by ionizing hydrogen and oxygen. The other method used in their study was to trap air bubbles by draining water from the bottom and reintroducing it through the top of the sample. The studies conducted by Yegian et al. [22] confirmed that the desaturation technique is not only an efficient method to increase liquefaction resistance, but also a cost-effective technique compared to other methods. Okamura et al. [23] used plastic pipes, connected to a pressure source, to inject air into the lower part of the site through several small holes drilled. They concluded that the liquefaction resistance of the unsaturated sample with a saturation degree of about 88–94% is nearly double that of the fully saturated sample. Eseller-Bayat et al. [24] utilized the chemical compound sodium perborate monohydrate to induce oxygen bubbles in saturated sand samples by reacting with water. Based on their study, an empirical model was presented to predict the accumulation of pore pressure ratio of desaturated sands during dynamic loading.

Besides the experimental studies mentioned above, there are also numerical techniques presented by the researchers. Several constitutive models have been generated to solve the liquefaction problems in saturated soils (e.g., [25,26,27,28,29,30,31,32]). These constitutive relationships can be used in continuous models, such as the finite element method or finite difference formulation, to predict the response of a saturated or unsaturated granular deposit. Bian and Shahrour [33] used a cyclic elastoplastic constitutive model to study the effect of saturation on the liquefaction potential in partially saturated soils. Buscarnera and Di Prisco [34] used the theory of stability of materials to study unsaturated slopes. They used a coupled hydromechanical constitutive model to study the onset of instability sensed by different failure modes such as static liquefaction. Liu and Muraleetharan [35] presented a constitutive model for unsaturated sands and sandy silts. They used the model to investigate the effects of saturation level, relative density, initial state, and stress state on liquefaction during stirring. Zhang and Muraleetharan [36] used a hydromechanical elastoplastic material model to study the liquefaction potential of unsaturated sands with different degrees of saturation and confining pressures. Other techniques such as the particle finite element method was also used to characterize the behavior of non-Newtonian Bingham fluids and reproduce the essential features of liquified sand [37].

The discrete element method (DEM) effectively accounts for the discontinuous behavior of granular soils. DEM was introduced by Cundall [38] for rock and granular problems. DEM was later applied to soil mechanical problems by Cundall and Strack [39]. DEM uses Newton’s second law and contact relations to calculate motion and contact forces between particles. The Lattice-Boltzmann method (LBM) is a powerful technique for simulating the behavior of fluid systems. The LBM can be used to simulate a variety of fluid behaviors, including complex boundary conditions, buoyancy effect, and phase separation. The key concept in LBM considers the behavior of the group of molecules or atoms rather than of any single individual. LBM simplifies the Boltzmann equation by reducing the number of fluid particles and constraining the position of the fluid particles and their directions of motion. LBM was utilized to model the features of liquefied sand, such as the decreasing viscosity with increasing shear strain rate [40]. One of the strengths of LBM is its ability to simulate multiphase fluids with single and multiple components. This means that surface tension, evaporation, cavitation and many other problems can be simulated by LBM. Component refers to a chemically independent component of a system while phase refers to the different states (gas, liquid, solid) of a component.

Multiphase multi-component LBM can be used to model the behavior of a system with different chemical components like oil and water system. The interaction force between the particles in multiphase multi-component LBM is repulsive. Various researches were conducted to study the multiphase multi-component LBM (e.g., [41,42,43,44]). Sukop et al. [45] used multiphase multi-component LBM to model the equilibrium distributions of two immiscible fluids in porous media. Galindo-Torres et al. [46] used the multiphase multi-component LBM model to investigate the dynamic effects in imbibition-drainage cycles within a porous medium. Galindo-Torres et al. [47] conducted a series of simulations using multiphase multi-component LBM to study the capillary pressure and saturation evolve in a porous medium during drainage and imbibition cycles. Richefeu et al. [48] utilized multiphase multi-component LBM to simulate capillary condensation to generate homogeneous distributions of liquid clusters in a 2-D porous medium. Multiphase single-component LBM models the behavior of a system with a single component and separated phases such as a liquid and its vapor in porous media. Unlike the multiphase multi-component models, the interaction force between the particles in multiphase single-component models is attractive. Grunau et al. [49] developed multiphase single-component LBM, where red and blue colored particles represent the behavior of different fluid phases. Their simulations of two-phase fluid flows were confirmed with the experimental studies and analytical solutions. Subsequently, Shan and Chen [50] developed a single component polyphase model in which an additional forcing term is added to the velocity field to produce phase separation. Huang et al. [51] used the one-component polyphase LBM to simulate the viscous coupling effects of two-phase flow in porous media. Eshghinejadfard et al. [52] used multiphase single-component LBM to calculate permeability by simulating laminar flows in porous media.

The coupled LBM-DEM framework was used in different geotechnical studies. Han and Cundall [53] used LBM and DEM to simulate three different porous media flow problems. Their conducted study confirmed the coupled LBM-DEM as a powerful technique to investigate the flow in porous media. Abdelhamid and El Shamy [54] generated a two-dimensional fully coupled fluid-particle model to study the fundamental mechanisms of flood-induced erosion. The solid phase was modeled by using DEM, and the fluid phase was modeled at a meso-scale level using LBM. El Shamy and Abdelhamid [55] introduced a three-dimensional coupled pore-scale model of pore-fluid interacting with discrete solid particles. This approach was utilized to model liquefaction in a fully saturated soil deposit subjected to low and large amplitude seismic excitations. Abdelhamid and El Shamy [56] presented a three-dimensional fully coupled pore-scale model to study the mechanism of fine-particle migration in granular filters using LBM and DEM. Qiu [57] used a two-dimensional model of LBM and DEM to investigate the flow in porous media. El Shamy and Abdelhamid [58] utilized coupled LBM-DEM to study the behavior of saturated granular soils during a bidirectional shaking.

In this paper, the results of soil-water systems with different degrees of saturation are presented on the basis of a coupled pore-scale idealization of the pore fluid and a discrete description of the solid particles. The microscale representation of the solid phase is obtained by DEM, while the multiphase liquid is modeled at the pore scale by single-component LBM. The fluid forces exerted on the granules are determined based on the momentum exchange between the fluid and the solid particles. The introduced procedure is used to model the response of a soil deposit with different degrees of saturation subjected to seismic excitation. The simulation results show that a small decrease in the degree of saturation of a fully saturated deposit has a significant impact on the magnitude of the excess pore pressure ratio.

2. Coupled Fluid-Particle Model

A fully coupled model is presented here to simulate granular deposits with different degrees of saturation subject to a dynamic base excitation. The fluid is modeled on the pore scale using the LBM. Here the method of generating a multiphase LBM model introduced by Shan and Chen [50] is used. The microscale response of the solid particle packet is captured using DEM. In the coupled LBM-DEM model, there is a different time step for each of the phases. In this study, the DEM time step is defined as smaller than the LBM time step. Therefore, an integration of the subcycle time is required for the calculation in the fixed phase. In the coupled model, the fluid and momentum equations are solved for each particle using an explicit time integration scheme. The rebound liquid boundary condition is applied between the surface of the solid particle and the liquid.

2.1. Multiphase Single-Component Lattice Boltzmann Method

Fluid can be considered as a medium made of tiny particles like atoms and molecules, where these particles stream and collide with each other. The collision between particles generates interparticle forces which can be represented by the ordinary differential equation of Newton’s second law. At each time step, location and velocity of each particle must be identified. This approach is called molecular dynamics (MD). The main drawback of MD is the large number of equations that must be solved in each step. However, it is not required to identify the behavior of every single molecule or atom. The main idea in LBM is considering the behavior of the collection of molecules or atoms as a whole instead of every single one. The scale of the collection of molecules is known as meso-scale, and it stands between micro-scale and macro-scale. In LBM, instead of Navier–Stokes equations, the Boltzmann equation is solved. The most important strengths of LBM are the ability to simulate multiphase fluid and the flow in complex boundaries. The terminology used in LBM is DnQm where n represents the dimension and m represents the number of linkages. The objective of this research is to simulate three-dimensional problems. Two popular models that can be used in LBM simulations of three-dimensional problems are the D3Q15 and D3Q19. The LBM D3Q15 (Figure 1) is used herein to reduce the computational cost of the simulations. The fictitious particle can move in the 15 directions in D3Q15 model. The central distribution function $f_0$ has no speed. There are six distribution functions which have lattice unit velocity in only one direction $f_1$, $f_2$, $f_3$, $f_4$, $f_5$ and $f_6$. The eight remaining distribution functions are at the corners of the lattice, $f_7$, $f_8$, $f_9$, $f_{10}$, $f_{11}$, $f_{12}$, $f_{13}$, $f_{14}$ and $f_{15}$ which can move in three directions and have lattice unit velocity in each direction.

Table 1. Velocity vectors for D3Q15.

Distribution	Velocity in	Velocity in	Velocity in
Function	X Direction	Y Direction	Z Direction
f${}_0$	0	0	0
f${}_1$	1	0	0
f${}_2$	0	1	0
f${}_3$	−1	0	0
f${}_4$	0	−1	0
f${}_5$	0	0	1
f${}_6$	0	0	−1
f${}_7$	1	1	1
f${}_8$	1	1	−1
f${}_9$	1	−1	−1
f${}_{10}$	1	−1	1
f${}_{11}$	−1	1	−1
f${}_{12}$	−1	1	1
f${}_{13}$	−1	−1	1
f${}_{14}$	−1	−1	−1

shows velocity vectors for any of these distribution functions.

Other than velocity vectors, the weighting factor for each direction is necessary for LBM calculation. For the D3Q15 model, the weighting factor is given as follows:

$$\begin{array}{c}\hfill w_i=\left\{\begin{array}{cc}2/9\hfill & i=0\hfill \\ 1/9\hfill & i=1–6\hfill \\ 1/72\hfill & i=7–14\hfill \end{array}\right.\end{array}$$

(1)

The Boltzmann distribution function is shown by $f(\mathbf{x},\mathbf{v},t)$, where f is the number of the collection of the molecules positioned between $\mathbf{x}$ and $\mathbf{x}+\delta \mathbf{x}$ with velocities between $\mathbf{v}$ and $\mathbf{v}+\delta \mathbf{v}$ at the time t. If the force $\mathbf{N}$ acts on a molecule with unit mass, the velocity of the molecule changes from $\mathbf{v}$ to $\mathbf{v}+\mathbf{N}\delta t$ and its position changes from $\mathbf{x}$ to $\mathbf{x}+\mathbf{v}\delta t$ after undergoing that force. After a collision occurs between particles, f will change in the $\delta \mathbf{x}\delta \mathbf{v}$ interval. The rate of change for the distribution function is known as the collision operator shown by $\mathrm{\Omega }$ (Mohammad [59]). The relationship between distribution functions before and after the collision is given by:

$$f(\mathbf{x}+\mathbf{v}\delta t,\mathbf{v}+\mathbf{N}\delta t,t+\delta t)\delta \mathbf{x}\delta \mathbf{v}-f(\mathbf{x},\mathbf{v},t)\delta \mathbf{x}\delta \mathbf{v}=\mathrm{\Omega }\delta \mathbf{x}\delta \mathbf{v}\delta t$$

(2)

Solving the collision term is very complicated. The popular Bhatnagar, Gross, and Krook (BGK) method (Bhatnagar et al. [60]) for collision operator is used herein. The collision operator is defined as:

$$\mathrm{\Omega }=\frac{(f^{eq}-f)}{\tau }$$

(3)

where $\tau$ is the relaxation factor and $f^{eq}$ is the equilibrium distribution function. The general form of the equilibrium distribution function can be written as:

$$f_i^{eq}=w_i(A+B{\mathbf{e}}_{\mathbf{i}}.\mathbf{u}+C{({\mathbf{e}}_{\mathbf{i}}.\mathbf{u})}^2+D(\mathbf{u}.\mathbf{u}))$$

(4)

where $\mathbf{u}$ is the macroscopic flow velocity vector, A, B, C, and D are constants and need to be determined based on the conservation principle (mass, momentum, and energy) and ${\mathbf{e}}_{\mathbf{i}}$ is velocity vectors for any of the distribution functions which is introduced in Table 1. The density factor, $\mathsf{\rho }$, is identical to the summation of distribution functions in all directions given by:

$$\mathsf{\rho }=\mathrm{\Sigma }{f}_i$$

(5)

The Bhatnagar et al. [60] equilibrium distribution function for fluid can be written as:

$$f_i^{eq}=\rho w_i(1+\frac{3}{c^2}{\mathbf{e}}_{\mathbf{i}}.\mathbf{u}+\frac{9}{2c^4}{({\mathbf{e}}_{\mathbf{i}}.\mathbf{u})}^2-\frac{3}{2c^2}\mathbf{u}.\mathbf{u})$$

(6)

in which c is the lattice speed, which is equivalent to the ratio between lattice spacing and lattice time step ($c=\mathsf{\delta }$x/$\mathsf{\delta }$t). In this research, c is considered equal to 1. Macroscopic properties of the fluid, velocity $\mathbf{u}$ and pressure p, are calculated as follows:

$$\begin{array}{c}\hfill \mathbf{u}=\frac{1}{\mathsf{\rho }}{\mathrm{\Sigma }}_{i=0}^n{f}_i{\mathbf{e}}_{\mathbf{i}}\end{array}$$

(7)

$$\begin{array}{c}\hfill p=C_s^2\rho \end{array}$$

(8)

where $C_s$ is the speed of sound on the lattice and equals to $c/\sqrt{3}$.

Adding an attractive force to the single-phase LBM model makes phase separation into a liquid and its vapor possible. In this study, the multiphase single-component Shan-Chen model [50] is used, which can simulate the behavior of water and air bubbles. No tracking of liquid-vapor interphase is needed in LBM. In the Shan-Chen model, the interaction strength between fluid particles is described by an interparticle potential function. The separation of the liquid and vapor phases occurs naturally in response to the attraction/repulsion force between different phases [61]. The Equation of State (EOS) allows the separation between phases. The Van der Waals EOS was used herein to describe the behavior observed in real gases as:

$$P=\frac{nRT}{V-nb}-a{\left(\frac{n}{V}\right)}^2$$

(9)

where P is the gas pressure, V is the gas volume, n is the number of moles, R is gas constant, T is temperature, a is a constant related to the average attraction force between particles, and b is the volume of one mole of particles. The second term on the right-hand side of Equation (9) is an attractive force between molecules. Since a is positive, this term reduces the pressure in comparison to that of a perfect gas. The $-nb$ term in the denominator accounts for the non-negligible volume of molecules. To simulate multiphase fluid and to consider the attraction force, long-range interactions between fluid particles are necessary. The interactions with nearest neighbor particle densities will be sufficient to simulate the basic phenomena of multiphase fluid interactions. For single component multiphase fluids, an attractive (cohesive) force, ${\mathbf{F}}^{\prime }$, between nearest-neighbor fluid particles is needed:

$${\mathbf{F}}^{\prime }(\mathbf{x},\mathbf{t})=-G\psi (\mathbf{x},t){\mathrm{\Sigma }}_{i=1}^{14}{w}_i\psi (\mathbf{x}+{\mathbf{e}}_{\mathbf{i}}\delta t,t){\mathbf{e}}_{\mathbf{i}}$$

(10)

where G is interaction strength factor and $\psi$ is the interaction potential function. A negative G value smaller than −4 indicates an attractive force that separates the fluid into two phases with different densities, including the liquid and vapor [62]. Yuan and Schaefer [63] mentioned that in the coupling of the single component model with the EOS, G becomes unimportant as its effects are canceled out in calculating the force of interaction. In this study, G was set equal to −5.2. According to Shan and Chen [50], the interaction potential function must be monotonically increasing and bounded:

$$\psi \left(\rho \right)=1-exp\left(\rho \right)$$

(11)

The formulations mentioned above incorporate the molecular attraction aspect of the Van der Waals gas model. For mimicking the attraction force between fluid molecules, ${\mathbf{F}}^{\prime }$ should be considered in the BGK operator. The LBM scheme presented for multiphase problems should consider the attractive force such that:

$$f_i(\mathbf{x}+{\mathbf{e}}_{\mathbf{i}}\delta t,t+\delta t)=\frac{\delta t}{\tau }{f}_i^{eq}(\mathbf{x},t)+f_i(\mathbf{x},t)(1-\frac{\delta t}{\tau })+{\mathbf{F}}^{\prime }\delta t(1-\frac{1}{2\tau })$$

(12)

The developed technique by Shan and Chen [50] is a popular method for modeling single-phase multi-component fluid behavior. Numerous studies that were conducted by utilizing Shan and Chen [50] method are in close agreement with analytical and experimental solutions. Some of the recent studies that were conducted by using Shan and Chen [50] method for simulating the behavior of multiphase fluid in porous media are Galindo-Torres et al. [46], Richefeu et al. [48] and Galindo-Torres et al. [47].

The interaction of the fluid and moving solid particles is considered by coupling DEM and LBM. The most simple and common technique for coupling the fluid and solid particles is the bounce-back technique (e.g., [47]). In this technique, particle surface has to be placed halfway between fluid and solid nodes. Even though the bounce-back method is of the first-order in error, lots of researchers are using this technique because of its simplicity. Galindo-Torres et al. [46,47] used the bounce-back boundary conditions for modeling multiphase LBM simulations in a porous medium. The forces and torques applied by the fluid on the solid grains are calculated based on the momentum the fluid exchanges with the solid particles. This method was originally used by Ladd [64] to calculate the forces on solid particles. Based on their study the fluid force is applied through each link connecting fluid node to solid node as:

$$\mathbf{F}({\mathbf{x}}_b,{\mathbf{e}}_{\mathbf{i}})={\mathbf{e}}_{\mathbf{i}}[f_i({\mathbf{x}}_b,t)+f_i({\mathbf{x}}_b+{\mathbf{e}}_{\mathbf{i}}\delta t,t)]$$

(13)

where ${\mathbf{x}}_b$ is the position of a solid node. The total fluid force (${\mathbf{F}}_{Fluid}$), applied on the particle is sum of the contribution of each link connecting a fluid node to a solid node:

$${\mathbf{F}}_{Fluid}={\mathrm{\Sigma }}_{all}{\mathrm{\Sigma }}_{i\ne 0}\mathbf{F}({\mathbf{x}}_b,{\mathbf{e}}_{\mathbf{i}})$$

(14)

The total torque (${\mathbf{T}}_{Fluid}$), applied on the particle is:

$${\mathbf{T}}_{Fluid}={\mathrm{\Sigma }}_{all}{\mathrm{\Sigma }}_{i\ne 0}\mathbf{r}\times \mathbf{F}({\mathbf{x}}_b,{\mathbf{e}}_{\mathbf{i}})$$

(15)

where $\mathbf{r}$ is a vector connects the particle center to its boundary. The magnitude of the resultant forces is affected by the number of lattice units covered by the solid particle.

2.2. Discrete Element Method

The micro-scale particle displacements and rotations impact the macro-scale response of a granular packing. DEM is an efficient method to model the behavior of granular materials on a micro-scale. In DEM, granular material is modeled as a collection of individual grains (Cundall and Strack [39]). Using DEM makes it possible to trace the movement of every single particle. Contact forces and displacements can be calculated by tracing the movement of individual solid particles. DEM calculations alternate between particles where Newton’s second law is applied and contacts where force-displacement law is used. In this study, other than contact and body forces, fluid forces and torques are acting on every single particle. Newton’s second law determines the motion of every particle under the mentioned loads:

$$\begin{array}{c}\hfill m\ddot{\mathbf{x}}=m\mathbf{g}+\mathrm{\Sigma }{\mathbf{f}}_{\mathbf{c}}+{\mathbf{F}}_{Fluid}\end{array}$$

(16)

$$\begin{array}{c}\hfill I\dot{\mathsf{\omega }}=\mathrm{\Sigma }{\mathbf{r}}_{\mathbf{c}}\times {\mathbf{f}}_{\mathbf{c}}+{\mathbf{T}}_{Fluid}\end{array}$$

(17)

where m is particle mass, $\mathbf{x}$ is particle position vector, I is particle moment of inertia, $\mathsf{\omega }$ is rotational velocity vector, $\mathbf{g}$ is gravitational acceleration vector, ${\mathbf{f}}_{\mathbf{c}}$ is the inter-particle force at each contact, and ${\mathbf{r}}_{\mathbf{c}}$ is a vector links particle center to the contact location. Momentum exchange between the solid particles and fluid is considered to calculate the applied force and torque. In DEM, solid particles are considered rigid; however, they can overlap. The relationship between relative displacement of two contacting particles and their contact force is generated by force-displacement law. The contact forces are acting in normal and shear directions:

$${\mathbf{f}}_{\mathbf{i}}={\mathbf{f}}_{\mathbf{i}}^{\mathbf{n}}+{\mathbf{f}}_{\mathbf{i}}^{\mathbf{s}}$$

(18)

in which ${\mathbf{f}}_{\mathbf{i}}^{\mathbf{s}}$ and ${\mathbf{f}}_{\mathbf{i}}^{\mathbf{n}}$ are normal and shear components of the force vector. The normal contact force is defined as:

$${\mathbf{f}}_{\mathbf{i}}^{\mathbf{n}}=K_n{\mathbf{x}}^{\mathbf{n}}{\mathbf{n}}_{\mathbf{i}}$$

(19)

where $K_n$ is the normal stiffness, and ${\mathbf{x}}^{\mathbf{n}}{\mathbf{n}}_{\mathbf{i}}$ represents the overlap in the normal direction. The shear contact force is computed incrementally according to:

$$\mathsf{\delta }{\mathbf{f}}_{\mathbf{i}}^{\mathbf{s}}=-K_s\mathsf{\delta }{\mathbf{x}}_{\mathbf{i}}^{\mathbf{s}}$$

(20)

where $-K_s$ is the shear stiffness at the contact, and $\mathsf{\delta }{\mathbf{x}}_{\mathbf{i}}^{\mathbf{s}}$ is the shear component of the contact displacement increment vector.

2.3. Computational Scheme

The solution of the coupled response of solid and fluid phases was achieved using an explicit algorithm. The DEM analysis was conducted using the PFC3D software developed by ITASCA [65]. The multiphase fluid code was generated by the authors in C++ and linked to the PFC environment. For solving the coupled model, one time step must be set for each of the solid and the fluid phases. The DEM time step ($\delta t$DEM) corresponds with solid phases and the LBM time step ($\delta t$LBM) corresponds with fluid phases. In order to have a constant time step in DEM, $\delta t$DEM must be chosen to be smaller than the DEM critical time step. The DEM time step is defined to be smaller than the fluid time step $\delta t$LBM. Thus, the fluid time step $\delta t$LBM is equal to $N\delta t$DEM, where N is an integer number.

By having this set of time steps, a subcycling time integration for the DEM part is necessary. For every fluid time step, the fluid code receives the location of each particle from the previous DEM time step. After completing the necessary calculations for LBM streaming and LBM collision steps, the hydrodynamic forces and torques are calculated. The DEM code obtains the forces and torques and applies them to the right-hand side of the equations of motion (Equations (16) and (17)). The obtained forces and torques remain constant throughout the subcycling. Finally, particle velocities and locations are obtained using an explicit central finite difference algorithm [65].

3. Model Verification

While there are several publications presenting the effects of desaturation on liquefaction mitigation through 1 g and geotechnical centrifuge experimental results (e.g., [18,19,20,66,67,68,69]), current computational power does not make it practical to have a one-to-one comparison with those experiments. Therefore, two series of verification simulations were performed before running the main model to ensure the LBM code and the coupled LBM-DEM code were without errors. The first series of simulations were modeled to find the relationship between surface tension ($\sigma$) and the inverse of the bubble radius ($\frac{1}{r}$). The second series of simulations were developed to capture the general form of the soil-water characteristic curve in a porous medium.

A spherical light fluid contained in a denser one, which forms a bubble, is affected by inward pressure and outward pressure. The inward pressure works to enlarge the bubble while the outward pressure squeezes it. The surface tension acts as a counterbalancing force to hold the bubble at a constant size. Based on Laplace law, the relationship between the inside and outside pressure difference ($\Delta P$) and ($\frac{1}{r}$) is a straight line with zero intercept, where the slope of the line represents surface tension. By simulating series of bubbles with different radius sizes, surface tension can be estimated. Figure 2a,b show bubbles with various sizes in the periodic domain. Three-dimensional models were generated in a 100$\delta x$× 100$\delta x$× 100$\delta x$ domain with periodic boundaries. In these models, denser fluid, which represents the surrounding phase, is 20 times heavier than the embraced fluid. Figure 3 presents the plot of $\Delta P$ vs. $\frac{1}{r}$ in dimensionless lattice units (lu), which is commonly done in LBM simulations (e.g., [46,47,48,61]). LBM lattice units could be transformed to physical units through dimensional analysis [46]. The relationship follows the linear relationship of $\Delta P$ vs. $\frac{1}{r}$ in Laplace law (Equation (21)) with a slope of 0.095 lu, which corresponds to the surface tension in analyzed fluids:

$$\Delta P=\frac{\mathsf{\sigma }}{r}$$

(21)

The relationship between the Laplace pressure difference ($\Delta P$) and volumetric degree of saturation (S) is also referred to as the water retention curve or the soil-water characteristic curve. Laplace pressure difference (suction), which is defined in Equation (22), is the difference between mean vapor pressure ($P_v$) and mean liquid pressure ($P_l$). In this paper, mean pressure of the lighter fluid is represented by $P_v$, and the mean pressure of the denser fluid is shown by $P_l$:

$$\Delta P=P_v-P_l$$

(22)

The soil-water characteristic curve follows a general form where $\Delta P$ declines in three stages. The first and last steps, where S decreases with lower slope, can be distinguished from the intermediate stage, where S declines rapidly as $\Delta P$ increases. The general shape of the soil-water characteristic curve was retrieved by the preformed LBM-DEM coupled simulations (Figure 4). A deposit with a smaller height than the one described in Section 4 was used in these simulations. The height of the deposit used herein is one fourth of the height of the deposit described in the next section. The density ratio between the two phases of fluids in the conducted simulations is 20. Here S was calculated as the ratio of the number of LBM nodes occupied by the denser phase to the number of LBM nodes occupied by both fluid phases in the voids.

4. Simulations

Three granular deposits with different degrees of saturation were generated to study the behavior of unsaturated soils subjected to applied shaking. A one-dimensional basic dynamic excitation in the x-direction was applied to the deposition shown in Figure 5. To model an infinite system, periodic boundary conditions for both liquid and solid phases are established at the lateral boundaries (El Shamy and Abdelhamid [55]). The periodic boundary conditions allow to simulate the huge problem in a limited computational domain. The lower limit is simulated as the liquid tightness limit for both phases. In the solid phase, the base wall mimics the behavior of a solid wall that can move in a horizontal direction. In the liquid phase, the rebound boundary condition is used to simulate the behavior of the solid wall. The constant pressure boundary condition is used to mimic the fluid surface. A semi-infinite deposit was created by solid particles falling under the effect of gravity in a parallelepipedal region with periodic boundaries in the lateral directions. The granular deposit was generated using the underwater bottom density to account for the buoyancy effect on the granules.

Initial sensitivity analysis simulations were performed to assess the spacing between the lattice grid. It was found that the diameter of the bubbles must be greater than 60 LBM nodes to keep the bubbles in a stable size. Considering the computational costs, the size of the generated bubbles are bigger than the solid particles (Figure 6). Essentially one bubble was generated at selected depth locations to yield the targeted degree of saturation (Figure 6). For creating the bubbles, some additional solid particles were generated between the deposit particles by PFC. The centers of the newly added solid particles were used to find the most suitable space between the deposit particles. For the bubbles, lighter fluid was introduced to the LBM model as spheres using the center of those particles. The diameter of these bubbles is covering 60 LBM nodes. After generating the bubbles, the newly added particles were deleted. Note that the location of the bubbles did not remain constant during the simulations. Additionally, the adhesive force between solid particles due to the presence of the fluids was not considered in conducted simulations since these forces tend to be negligible for the simulated range of degrees of saturation. Mei et al. [70] recommended using at least 10 lattice units across the particle diameter. In this study, the average particle diameter is covering 12.5 lattice units.

To decrease the number of particles in DEM, the number of nodes in LBM and the duration of the applied shaking, a high gravitation acceleration was applied to the model. The applied high-g level mimics the small scale geotechnical models in centrifuge testing conditions.

Table 2. Scaling laws used for the conducted simulations.

Parameter	Model Prototype Ratio
Gravity	N
Length	1/N
Acceleration	N
Velocity	1
Density	1
Mass	1/N${}^3$
Time	1/N
Force	1/N${}^2$
Stress	1

shows the scaling laws utilized to convert different parameters from prototype to model units. In the employed high g-level technique, the scaling factor is only applied to macro-scale parameters, and micro-scale parameters including the average particle diameter are considered identical to the original size in the prototype. This method was found to be a useful technique in DEM simulations to decrease the computational domain to a manageable size in several applications (e.g., [71,72,73]). A summary of computational parameter details and proprieties for solid particles and fluids are listed in

Table 3. Simulation details.

Solid particles
Diameter	4.8 mm–7.2 mm
Normal stiffness	5 $\times {10}^6$ N/m
Shear stiffness	5 $\times {10}^6$ N/m
Normal damping ratio	0.1
Shear damping ratio	0.1
Friction Coefficient	0.5
Density	1650 kg/m${}^3$
Number of particles	3100
Fluid
Density	1000 kg/m${}^3$
Viscosity	5.0 Pa.s
LBM nodes number	103 × 480 × 103
Computation parameters
g-level	30
LBM Time step	8 $\times {10}^{-6}$ s
DEM Time step	8 $\times {10}^{-7}$ s

.

The initial permeability of the deposit was estimated by the Kozeny-Carman Equation [74]:

$$k=\frac{\rho g{\varphi }^3}{5\mu (1-{\varphi }^2)}{\left(\frac{\Sigma \frac{4}{3}\pi r^3}{\Sigma 4\pi r^2}\right)}^2$$

(23)

where $\rho$ is the fluid density, r is particle radius, g is the gravitational acceleration, and $\varphi$ is the average deposit porosity which is 0.41. The deposit was saturated with a highly viscous fluid to compensate for the effects of the employed 30-g field and large particle sizes (an average diameter of about 6 mm) and comply with the scaling of permeability. For the prototype fluid viscosity of 0.166 Pa.s and the employed particle sizes, the initial permeability of the deposit was estimated to be 2.21 mm/s (same order of coarse sand permeability when saturated with water). The use of much higher viscosity than that of water (166 times higher) compensates for the employed large particle size that would not normally liquefy if it was saturated with water or lightly viscous fluid.

4.1. Computational Details

In the simulations carried out, a deposit of solid particles 190 mm high with lateral dimensions of 50 mm × 50 mm was generated. The pellets had an average diameter of 6 mm with a minimum diameter of 4.8 mm and a maximum diameter of 7.2 mm (Figure 5) and the average saturated density of the deposit was approximately 1967.4 kg/m${}^3$. The high gravitational acceleration applied to the granular deposit was 30 g. Based on the centrifuge scaling laws (Table 2), the macroscopic length and time scales in the model were 30 times smaller than those in the prototype. As a result, the amplitude of acceleration and the frequency of movement in the model are 30 times higher than in the prototype. The model looks like a prototype where the depth of the granular deposit is 5.7 m.

The one-dimensional dynamic base excitation (in the x-direction) was applied as a velocity-time history to the solid base wall on which the particle deposit rested. The time history of the corresponding acceleration of the dynamic input excitation in the prototype shell is shown in Figure 7. The sinusoidal acceleration input signal is increased gradually until it reaches the maximum amplitude at 4.5 s. It then remained constant for 7.5 s and gradually decreased to zero after 13 s. To estimate the low-stress shear modulus, mean shear wave velocity, and natural frequency of the fully saturated deposit, a sinusoidal fundamental acceleration with an amplitude of 0.01 g and a frequency of 3 was applied to the deposit. Hz under a gravity of 30 g. The acceleration and shear-strain histories were obtained. Using these stories, the low-stress mean shear modulus was found to be about 70 MPa, which corresponds to a mean shear wave velocity of about 187 m/s. Therefore, the natural frequency of the deposit (defined as the shear wave velocity divided by four times the height of the deposit) was approximately 8.2 Hz.

The results of models with 100%, 92%, and 82% of the degree of saturation were analyzed. Several parameters were monitored at a time interval of 0.000016 s in the model (0.0048 s in prototype) during the simulations. These parameters include the vertical effective stress, excess pore-pressure ratio, coordination number, shear stress–strain time histories, and averaged particle acceleration. Results are presented in prototype units unless otherwise specified.

4.2. Simulation Results

Excess pore pressure ratio (ratio of fluid pore pressure increase to initial effective vertical stress) and effective stress-time histories were used to identify liquefaction from a large-scale perspective. The accumulation of an excessive pore pressure ratio was the key factor in identifying the initiation of liquefaction. In general, liquefaction has been defined as the moment when the pore pressure ratio reaches the value one. At this time, the contact forces between the particles have approached zero and therefore the effective stress has disappeared. The temporal evolutions of the pore pressure ratio at selected depth locations for models with different degrees of saturation are shown in Figure 8. The results show that the pore pressure ratio approached values between 0.8 and 1 in the upper 2.7 m of the fully saturated deposit. In the deepest parts of the fully saturated deposit, the pore pressure ratio remained around 0.8.

Based on the simulations performed, liquefaction occurred at the fully saturated deposit. In the desaturated deposits, this ratio remained below 0.5 for all levels, despite the increase in the pore pressure ratio at shallow places (Figure 8). No liquefaction was observed at any depth in the 82% and 92% saturated deposit. It can also be concluded that the decrease in the degree of saturation from 100% to 92% has a significant impact on the decrease in the pore pressure ratio. However, further reducing the saturation level from 92% to 82% has less significant effects. This answer shows that a small percentage reduction in the degree of saturation has a significant effect on reducing the risk of liquefaction. The same behavior was observed in several experimental studies (e.g., [18,19,20]). It should be noted, however, that most of these experimental studies showed a considerable difference in behavior as the degree of saturation decreased say from ∼95% to ∼85%, which was not observed in current simulations. This could be attributed to the testing environment that imposes undrained conditions in the experiments while this condition is not enforced in the presented LBM-DEM model. On the other hand, the field study reported by [23] showed that no significant difference was observed regarding the increase in liquefaction resistance for the degrees of saturation of about 88∼94%. Similarly, the free-field results of the centrifuge experiments performed by [68] showed no significant difference between the response at the degrees of saturation of 93.5% and 79.5%. The results of the in situ testing and centrifuge experiments, which are in agreement with the results from performed simulations, highlight the importance of system level modeling as opposed to element level testing.

Liquefaction is also marked when the vertical effective stress of the soil decreases to zero. The vertical effective stress of the fully saturated deposit declined and approached zero at all depth locations. The reduction in the vertical effective stress was significantly lower in the desaturated deposits compared to the fully saturated one (Figure 9). There is not a significant distinction between the vertical effective stress time histories of the desaturated deposits as can be seen in Figure 9.

During liquefaction, the shear moduli of the granular soil approached zero as a result of the vanishing effective stress. The shear stress–strain loops for the deposits at selected depth locations are shown in Figure 10. The shear stress–strain loops show a stiffness degradation of shear moduli at different depths of the deposits with different degrees of saturation during shaking. The shear modulus of the fully saturated deposit approaches smaller values in comparison with the desaturated deposits. Shear modulus reduction curves for soil deposits with different degrees of saturation are shown in Figure 11. The reduction in shear modulus in the desaturated deposits is mainly due to the experienced large shear strains. However, the degradation of the shear modulus in the saturated deposit is due to the combined effect of pore pressure buildup and subsequent decrease in shear stress as well as the large shear strains experienced by the deposit.

Figure 12 shows the settlement time histories in the conducted simulations. There is not a significant difference between the final settlements and settlement time histories of the deposits with various degrees of saturation. Even though the reduction in porosity causes the pore-pressure buildup and, consequently lead to liquefaction, there is not a significant difference between the settlement time histories of the liquefied and desaturated granular deposits. It can be concluded that the main reason for the settlement in the simulated deposits was the applied shaking and not liquefaction of the deposit.

It was of importance to examine the averaged acceleration in the shaking direction of the solid particles at each depth location (Figure 13). The behavior of acceleration time history presents the resisting force to the motion of the solid particles during shaking. As the soil liquefies, particle contacts are lost and smaller acceleration amplitudes are transmitted throughout the liquefied portion of the deposit (Figure 13a). As with the other evaluated parameters, there is not an apparent difference between the acceleration time histories of the desaturated deposits (Figure 13b,c).

The degree of saturation variation during shaking at different depth locations of desaturated deposits is shown in Figure 14. In general, the degree of saturation was not constant with depth, and its local value changed with time as the bubbles moved up and down and were squeezed during shaking. The averaged degree of saturation variation of the deposits during the shaking is shown in Figure 15. There is a small increase in the degree of saturation for both desaturated deposits. This increase of the degree of saturation represents the volume reduction of the air bubbles during the cyclic loading. The air bubbles volume reduction is happening because they absorb the excess pore pressure generated during shaking. From the micro-scale point of view, liquefaction occurs at the locations where the coordination number of the solid particles (the average number of contacts per particle) reduces until they start floating. In order to have a stable pack of frictional spheres, the coordination number is required to remain 4 or above (Edwards [75]). Figure 16 shows the coordination number time histories of the deposits with different degrees of saturation. The top 0.6 m of the fully saturated soil deposit, coordination number dropped significantly to a value below 4. The 2.7 m and 4.8 m depth locations of the fully saturated deposit the coordination number was sustained above 3, but below the value of 4. This response shows that instability occurred in every depth locations of the fully saturated deposit. Figure 16 also shows that the coordination numbers in the 82% and 92% saturated deposits remained virtually constant. This behavior confirms that liquefaction was not observed at any depth location in the desaturated deposits. Coordination number time histories for both desaturated deposits were almost identical during the applied shaking. It can be concluded that the decrement of the degree of saturation from 92% to 82% does not have a significant effect on the coordination number time history. The coordination number time history is affected distinctly when the degree of saturation drops from 100% to 92%.

5. Conclusions

This work investigates the potential of a three-dimensional pore model developed to analyze liquefaction in fully desaturated granular deposits under seismic excitation. The multiphase single-component lattice Boltzmann method was used to model the liquid phase at the pore scale, while the solid phase was idealized at the microscale using the discrete element method. The presented coupled LBM-DEM model considers the discrete behavior of granular particles, including the actual fluid flow in the pore space. The forces exerted on the solid particles by the fluid phase were calculated by considering the exchange of momentum between the solid particles and the fluid. The model was used to simulate the response of soil deposits with different degrees of saturation during seismic excitation. Three granular deposits subjected to dynamic excitation with saturation degrees of 100%, 92% and 82% were modeled.

Based on the results of the simulations carried out, liquefaction took place in the fully saturated granular deposit. For the same set of conditions applied for the fully saturated deposit, no liquefaction was observed at any depth location of the unsaturated deposits. It has also been noted that reducing the degree of saturation from 100% to 92% primarily affects behavior. At the depth near-surface location, the excess pore pressure ratio dropped from near 1 to less than 0.4. Furthermore, at the deepest points of the 92% saturated deposit, a significant drop in the pore pressure ratio in excess of the values observed in the fully saturated deposit was observed.

A significant behavior change was observed in the effective stress time histories after inducing the bubbles to the fully saturated deposit. Acceleration amplitude increased by desaturating the fully saturated deposit during shaking in all depth locations. Shear modulus approached higher values in the 92% saturated deposit in comparison to the fully saturated one for a specific shear strain. Further reduction in the degree of saturation from 92% to 82% did not have a significant effect on the computed parameters. Similar behavior was observed for the time history of the coordination number. Coordination number time history was affected significantly by an 8% decrease of the degree of saturation in the fully saturated deposit. However, it remained close to the value of 4 in both 82% and 92% saturated deposits. No significant difference was observed between the settlement time history and the final settlement of the granular deposits with different degrees of saturation. The volume reduction due to the dynamic excitation was the leading cause of the settlement and not liquefaction.

Because of the high cost of the computational scheme, simulating the practical size of the solid particles and air bubbles was not feasible herein. However, with the current increasing trend of high performance computing platforms that utilizes parallel computing, there is high potential for the presented technique to model systems of actual particle and bubble sizes.