Numerical study of a desaturation process by air injection in sandy soil

This work aims to numerically study the desaturation process in a coarse-grained sandy deposit by means of air injection. It is well-known that the soil cyclic strength to liquefaction of a saturated, sandy deposit is positively affected by the presence of gas in the void space in either a dissolved or a free form. A numerical study is performed to investigate some of the factors affecting the desaturation advance and controlling the soil desaturation process by air injection. Among these factors are: a) soil-water characteristic curve and intrinsic permeability (hydraulic parameters), b) injection time, and c) standard approach to two injection wells placement. Only the standard mechanisms of biphasic flow are investigated (i.e., incompressible and isothermal biphasic flow in an isotropic homogeneous porous media with capillary effects).


Introduction
The liquefaction phenomenon has been associated with the generation of massive damage to diverse categories of buildings during preceding earthquakes, in particular structures supported in shallow foundations. This massive damage is perceived from the soil liquefaction around vast areas, like Niigata's in 1964, where almost all the two to four level buildings founded on shallow foundations were affected [1]. Also, the 1990 Philippines earthquake, where almost all affected structures were two to four levels supported on shallow foundations [2] as well as the Darfield earthquake in 2010, where soil liquefaction produced damage in two-level structures founded on spread footings that suffered lateral spreading [3]. Furthermore, the Christchurch earthquake in 2011 that occasioned ground bulk upward and foundations inclination in small to medium height edifications [4], among others. These case studies exposed the strong association between support behavior and potential damage of existing structures founded on shallow foundations associated with liquefaction.
These conditions prompted the development of countermeasures that have been employed in situ to mitigate deformations associated to liquefaction. They are accompanied by procedures that consider foundation behaviour, structural response, and laboratory-scale models in well-controlled environments. Information on available liquefaction soil improvement techniques is provided by Seed et al. [5]; however, these techniques are only applied to projects involving essential structures [6]. It is recognized that most soil improvement techniques use large amounts of energy due to the materials production and the construction process. There is an apparent necessity to develop cost-effective liquefaction mitigation techniques available for existing structures [7].
The effectiveness of inducing partial saturation as mitigation of liquefiable soils has been discussed in the literature [8]. Methods have been proposed and tested to desaturate the soil, decreasing bulk modulus and increasing the air-water mixture compressibility inside pores, because the modulus of the wetting phase decreases dramatically with a minimal volume of air [9]. During cyclic loads, the gas absorbs the excess pore pressure generated by reducing its volume [10]. Methods that include air injection [8], water electrolysis [7], chemical reactions [11], and microbiological processes [12], among others, have been proposed. These methods offer a less invasive, cost-effective, and practical alternative when mitigation can affect soil permissible load capacity during the process [13].
This work aims to present an overview of liquefiable soil improvement methods through the generation of gas bubbles due to the necessity to promote the application of new methodologies and to identify and collect technical references. This work also emphasizes practical aspects to include considerations on geotechnical soil properties. Additionally, a mathematical inspection of desaturation by air injection is carried out. It is done since a slight degree of saturation changes have substantial effects on the soil cyclic shear strength. A parametric analysis is performed, evaluating parameters that actively control the soil desaturation process, behaviour of partially saturated soils, injection time, and desaturation radius using to parallel injectors. In order to develop the predictions, the lite (free) version of the PDE Solutions Inc [14] software was implemented.

Methods to induce degree of saturation 2.1 Air-injection
This method consists in inserting a pipe to liquefiable soils and exhaust air from the tip. Partial desaturation of soil by air injection has gained attention in recent years as an innovative technique to decrease soil liquefaction susceptibility. Some experimental models are developed in shaking tables [15], centrifugal [16] and in-field [17], demonstrating that the desaturation advance front is affected by the soil type, stratification, injection pressure, and depth [18]. Zeybek and Madabhushi [19]studied the desaturated soil deformation trough a sample column with and without shallow foundations support, finding that bearing capacity failure was avoided. On the other hand, the induced liquefaction depth, deviatoric strains, and surface ground settlements were reduced with desaturation. Zeybek and Madadbushi [19] stated that during gas generation, under monotonic loading, volumetric contractions occurred due to a decrease in effective stresses and an increase in soil compressibility. Also, expansive volumetric strains linked to coagulation and ascendant vertical flow of the non-wetting phase were observed, triggering a partial bearing capacity failure condition. Then, after desaturation, volumetric deformations arose due to seismic acceleration applied. However, these deformations as well as the bearing capacity failure were limited during reconsolidation due to the dissipation of excess pore pressure. In a similar study, Zeybek and Madabhushi [20] stated that the reconsolidation and foundation settlement were minimized for high foundations pressure, increasing the post shaking bearing capacity, but intensifying accelerations transferred to the structure.
To study the longevity of desaturation, laboratory tests performed in a zone desaturated 26 years ago demonstrated that none of the samples had more than 90% of the degree of saturation [21]. In a similar study, made by Okamura and Teraoka [22], it is inferred that generating degrees of saturation below 90% would produce desaturation conditions for more than ten years. Zeybek and Madabhushib [16] studied the durability of air bubbles under hydrostatic conditions, at low and high pore pressure, vertical upward and downward flow, variable pore pressure, and lateral shaking.

Electrolysis of water
Electrolysis introduces gas inside the saturated skeleton pores without using a non-wetting phase pressure generating decomposition of water into oxygen and hydrogen by a continuous electric current [23]. Yegian, Eseller, and Alshawabkeh [24] used an electric current to produce gas in a fully saturated sample. After electrolysis, desaturated samples were evaluated cyclically, and no liquefaction was generated. Yegian, Eseller, and Alshawabkeh [8] defined the amount of gas produced due to electrolysis and indicate that for a reduction in 3% in degree of saturation, liquefaction was avoided. Eseller et al. [25] analyzed the permanence of gas bubbles during long-term and groundwater flow in a sample to which the water electrolysis process was carried out. Ended the process, the presence of gas was perceived after a week without showing possible diffusion. Also, Yegian, Eseller, and Alshawabkeh [7] observed the diffusion of trapped gas after 442 days of water electrolysis, and desaturation was maintained.

Chemical method
This technique uses sodium perborate monohydrate that makes a reaction when it comes into contact with water and produces a large fount of oxygen [11]. The sodium percarbonate is environmentally friendly and easy to apply without health and security concerns. Moreover, sodium percarbonate injection is usually applied in soil bioremediation to remove organic contaminants [26].
Eseller-Bayat [27] desaturate soil with efferdent powder mixtures, avoiding liquefaction during dynamic tests. In a similar work, Eseller-Bayat et al. [11] stated that samples with a degree of saturation of less than 90% do not reach liquefaction during cyclic tests. Complementary results were found in Yegian et al. [7], and it was concluded that for defined degrees of saturation and a shear strain amplitude, the higher the relative density, the lesser the maximum pore pressure. Furthermore, Gokyer [28] developed a code-named SUTRA-Bubble, characterizing three-dimensional chemical kinetics, oxygen gas generation, and degrees of saturation decrease, to calculate the injection pressure and the time interval of a chemical solution concentration to generate the required degrees of saturation plume and concentration of sodium percarbonate. Yegian et al. [7] indicate that desaturation can be prolonged in time below hydrostatic conditions, vertical flow gradients, and dynamic load, showing desaturation invariability.

Microbial method
The bio desaturation is implemented by forming tiny insitu nitrogen gas bubbles, using microbial denitrification [29]. The nutrients are dissolved in the water and can travel similar to the water in the sand [12]. The partial induction of the degrees of saturation through biogas' generation reduces the excess pore pressure and increases the ground´s shear strength, which benefits the foundation's response directly [30]. It is recognized that the liquefaction strength of the soils is substantially increased when those are slightly desaturated by the inclusion of gaseous nitrogen produced by the denitrifying bacteria in their voids [31]. However, He [32] perceived that for the generation of in-situ biogas in a sand column assembly by denitrifying bacteria, the nitrogen bubbles were not consistent during water flows with a hydraulic gradient equal to 0.1.
As has been shown, the generation of air bubbles in the soil considerably improves the soil's strength to liquefaction. Next, a numerical process will be presented in which the desaturation characteristics achieved through an air injection process in sandy soil are highlighted.

Mathematical model description
In this work, an incompressible and isothermal biphasic flow in homogeneous porous media with the capillary effects model is included. Pore fluid´s mass conservation is solved through a generalization of Darcy´s law, biphasic flow, as presented by Pinder and Gray [33].

Mass-momentum conservation equations
The mathematical model is established for both the nonwetting phase air, nw, and the wetting phase water, w, as demarcated by the phase velocity, , Equation (1): Where kint is the intrinsic permeability of the fluid i (m 2 ); kr,i is the relative permeability (a function of saturation for a given fluid); μi is the fluid's dynamic viscosity (Pa·s); pi is the fluid pressure (Pa); ρi is the fluid density (kg/m 3 ); g is the acceleration of gravity; and D is the coordinate of vertical elevation (m). Considering a homogeneous and isotropic skeleton, the macro-scale balance conservation for each phase ݅ is presented in the Equation (2) and (3): Where θ is the total porosity; Si is the saturation and t is time. Phase saturations responds to the constitutive relation presented in Equation (4): Concerning the capillarity, it is usually defined as a macro-scale capillary pressure, pc, function of the wetting phase saturation, commonly defined as Equation (5): The capillary pressure correlation and the saturation relationship eliminates two unknowns, and the mass conservation is stated by Equation (6) and (7):

Model formulation
The mass conservation is written in the form of a pressure-saturation system by summing Equations (6) and (7) generating Equation (8). The system then reads as Equation (8) and Equation (7).
Moreover, supposing that the capillary pressure is merely governed by the wetting phase degree of saturation, the capillary term સ can be redefined as Equation (9):

Soil-water characteristic curve and capillary pressure model
The capillary pressure dependence on the saturation is assumed by the Equation (10), based on Van Genuchten [34]: Where the expression m = 1 -1/n is supposed. pc (sw,eff) is the suction and differentiating the Equation (10); the capillary term in the pressure Equation can be assumed as Equation (11):

Relative permeability relationships
The concept of effective degree of saturation is included, that is a normalization of the wetting phase degree of saturation, as presented in the Equation (12): sw,r , and snw,r are the minimum limiting values of the soil residual degree of saturation for the wetting phase and non-wetting phase degree of saturation, respectively. In the model, kr,w , and kr,nw are the non-dimensional permeability functions for the wetting and non-wetting phases. Here, and as presented in Pinder and Gray [33], The model parameters are based on the Mualem porosity distribution function and the van Genuchten capillary pressure relationship. The approach makes use of the pore-size distribution function in deriving ‫‬ − ܵ ௪ models. The relative permeabilities, as presented in Van Genuchten [34] are described by Equations (13) and (14): Where a and b are the shape parameter of wetting phase permeability and the shape parameter of non- wetting phase permeability, respectively, and m and n are the material parameters in the Van Genuchten [34] expressions.

Comparison with experimental results
In this study, to calibrate and validate the present model, the desaturation test employing air-injection and evaluating magnitudes of partial desaturation generated during and halted air injection presented by Yasuhara et al. [35] is used. Yasuhara et al. [35] developed experimental measurements of saturation developed within a 77 cm height, 25 cm width, and 10 cm depth acrylic container, as presented in Fig.1. A saturated 27 cm-height soil column with a Dr=60 % was fabricated by pouring in water layers of Toyoura sand. When the specimen was constructed, an overload confining soil skeleton, to avoid sand boiling and piping, was placed during the air injection process. A pressure gauge was fixed on top of the soil ground to analyse the free water level changes continually. As a result, during air injection with an air pressure of 12 kPa, an average degree of saturation in soil column as a function of porosity and free water level variation was quantified.  [35].
A numerical analysis of the simultaneous flow of water and air during air injection was developed. Numerically, the matric suction and fluid phase's impedance were involved. The former, implementing the relation between the suction (negative water pressure) and saturation for the Toyoura sand employed for the air injection experiments, means the laboratory soil/water retention curve, as shown in Fig. 2. The latter, including the fluid phase's relative permeability, stated by the laboratory experimentation described in Fig. 3. Together with the fitting parameters presented in Table 1. The group of the mathematical equations for the multiphase flow and the parameters for predictions were implemented. Finally, the comparison results between the experimental measurements and the numerical results are shown in Fig.  4. Fig. 4 indicates that airflow and the evolution of the desaturation by air injection were defined when the injection was conducted. In the Yasuhara experimental work, it is perceived that the desaturation monotonically occurred in time, and eventually, a steady-state with a degree of saturation of 65 % was obtained. During the steady-state, the air injection is halted in the 1800s. Consequently, the saturation progressively rises on time, reaching approximately 85 %, however without recovering the initial saturated state due to remaining gas occluded within the soil skeleton.   experimental results presented. Calibrated results are situated within the experimental work with desaturation in time and then a steady-state with 65% saturation. Differences in degree of saturation are probably related to the anisotropy in the permeability of soils. It is typically recognized that the horizontal permeability is greater than the vertical one, and this effect is not studied in this work. Finally, as it is observed, residual degrees of saturation numerically obtained are nearly equivalent to those defined in the experimental work.

Fig. 4.
Comparisons between measurements from Yasuhara et al. [35], and predictions by the proposed model in this work.

Numerical results of Air Injection
This analysis is developed to understand the effect of the hydraulic parameters and injection time on the effective radius. The hydraulic parameters, typical for homogenous deposits of sands and silty sand soils, were defined from Lu and Likos [36], and they are presented in Table 2. The first group of simulations is developed on three different soil types and analysed for defined air injection. Each soil in Table 2 has related parameters dependent on the pore size, geometry, and distribution. Consequently, a collection of intrinsic permeability values was designated to characterize deposits of clean to fine sands. Intrinsic permeability is constant for any soil during predictions irrespective of the sort of parameters of the fluid being conducted, responding to soil skeleton behavior during the injection process at low pressures, without pore structure alterations [17]. Table 3 shows the air and water physical parameters used in multiphase flow predictions. Data from Chen et al. [37]. The domain used for the analysis is formally taken to be the rectangular region in Fig. 5, being made of a homogeneous soil deposit of 30 m width and 15 m height. For the wetting phase, an impermeable wall boundary condition is given to the bottom and the air injection hole. For the non-wetting phase, airflow can be developed crossing the lateral and top soil column boundaries, and at the injection, point placed at 6.0 m depth. The predictions start from a full water-saturated state, and air advances as the air pressure exceeds the hydrostatic and capillary pressures. The initial and boundary conditions for this model are summarized as follows:

Fig. 5. Finite element domain and boundary conditions.
For soil desaturation applications, Ogata and Okamura [38] suggested that the maximum injection pressure must be less than the hydrostatic pressure plus 50% of the effective vertical stress, accordingly to avoid pore structure to be damaged or disturbed during air injection execution. In this work, the maximum injection pressure suggested, according to Ogata and Okamura [38], is approximately 87 kPa. It is related to an injector at 6.0 m, and a soil unit weight of 19.0 kN⁄m 3 , greater magnitudes may generate some counter-productive internal structure reorganization [16]. Consequently, the injection pressure implemented is Pnw,inj =85 kPa. It is essential to clarify that this injection pressure value was stated to understand the desaturation front advance as a different response among soils on time. The effective radius desaturation magnitude of snw=0.1 was designated since Chaney [39], and Yoshimi et al. [40] stated that the sand liquefaction strength doubles when a 90% partial saturation induction is reached.
The results for the first group of simulations are shown in Fig. 6. It shows an effective radius as a function of soil and time. It is observed that at lower values of intrinsic permeabilities, kint= 5.55x10 -12 m 2 (Soil 1) and kint= 1.00x10 -11 m 2 (Soil 2), effective radius grows almost linearly with increments in time. For the most permeable soil, kint= 5.55x10 -11 m 2 (Soil 3), the trend is different, achieving a similar effective radius regardless of time. Larger hydraulic conductivities allowed a larger effective radius of desaturation and in a shorter time. However, for Soil 3, it can be perceived that a steady state is reached between 2000 and 3000s. Given the importance of the air injection holes distribution in the target treatment zone, the injection hole pattern and spacing should be analysed based on the ground desaturation radius and the complete coverage of the desaturated area within the injected depth. During application, injection holes may be arranged in multiple rows. The second group of simulations was developed for soil 2 to present in situ ground desaturation technique predictions with the inclusion of a second drilled hole developing air injection. The second injection hole presents the same depth, geometry, and pressure characteristics as the first injector presented in Fig. 5. A schematic plan view of the system with the second injector is shown in Fig. 7.
The second injector's contribution to the desaturation system was analysed, and the effective desaturation of the injected depth is presented for one row of two injected holes. Predictions were made with injection hole horizontal spacing of 2.0, 2.5, 3.0, and 3.5m. The desaturation process results are presented for the four proposed spacings between injection holes and injection times of 900, 1800, 2700, and 3600s. Fig. 8 is presented, from which it can be seen that for a distance between injectors of 2.0m and an injection pressure of 85kPa, the desaturation process is effective. Observing that for injection times of 900, 1800, 2700, and 3600s a degree of saturation of 95, 90, 87, and 84% is reached, respectively.
According to Fig. 9, it is possible to perceive that for separation of 2.5m, and a recommended injection pressure of 85kPa, results were obtained for an injection period of 3600s according to what is desired, precisely achieving partial induction of the degree of saturation less than 90 % through the entire horizontal distance between injectors at 6m depth. However, for lower times of 2700s and 1800s, maximum degrees of saturation of 93 and 95% respectively were perceived, showing that somehow such injection times for the chosen pressure and spacing of injectors are not convenient. Fig. 6. The effective radius for soils with hydraulic parameters described in Table 2.   On the other hand, less promising results are presented in Fig. 10, and Fig. 11, understanding that even for injection times of 3600s at spacings between injectors of 3.0 and 3.5m, degree of saturations of around 96 and 98% were found, respectively. From this, it can be inferred that for the given separation and static pressure of 85kPa, the partial induction of saturation is not adequate. It is noted in Fig. 10, and Fig.11 a loss of effectiveness in the desaturation process precisely at medium distances between injectors, with associated degrees of saturation of 100% or close. Higher injection pressures would generate severe alterations on the soil structure, and the spacing between injectors should probably be reduced.   Fig. 11. Degree of saturation for 3.5m standard approach to two injection well placement.

Conclusions
The desaturation techniques belong to a young interdisciplinary field; however, it has been shown that are effective. Due to the implementation of readily inexpensive equipment, desaturation can offer significant savings compared to the commonly applied liquefaction mitigation methodologies. Some of the parameters that revolve around methods are controlled and would allow a flexible implementation. Few desaturation systems have been designed, operated, and monitored in the field, under different conditions. In situ remediation applications information would help as a database to identify the advantages and limitations to guide future systems. The proposed multiphase flow model presents concordance and coherence with the experimental measurements of the evolution of saturations in time and space. The two conducted analyses facilitate the process of keeping track of the local desaturation processes. The former presents an effective radius as a function of soil and time. It was observed that an effective radius grew almost linearly with increments in time at lower values of intrinsic permeabilities. For the most permeable soil, the trend was different, achieving a similar effective radius regardless of the time. The latter presents the second identical injector's contribution to the system, characterizing that for shorter distances, between injectors, desaturation was generated more rapidly than at greater distances.