A multi-phase biogeochemical model for mitigating earthquake-induced liquefaction via microbially induced desaturation and calcium carbonate precipitation

A next-generation biogeochemical model was developed to explore the impact of the native water source on microbially induced desaturation and precipitation (MIDP) via denitrification. MIDP is a non-disruptive, nature-based ground improvement technique that offers the promise of cost-effective mitigation of earthquake-induced soil liquefaction under and adjacent to existing structures. MIDP leverages native soil bacteria to reduce the potential for liquefaction triggering in the short term through biogenic gas generation (treatment completed within hours to days) and over the longer term through calcium carbonate precipitation (treatment completed in weeks to months). This next-generation biogeochemical model expands earlier modeling to consider multi-phase speciation, bacterial competition, inhibition, and precipitation. The biogeochemical model was used to explore the impact of varying treatment recipes on MIDP products and by-products in a natural seawater environment. The case study presented herein demonstrates the importance of optimizing treatment recipes to minimize unwanted by-products (e.g., H₂S production) or incomplete denitrification (e.g., nitrate and nitrite accumulation).

1 Introduction

Microbially induced desaturation and precipitation (MIDP) is a biogeotechnical technique that takes advantage of native subsurface denitrifying bacteria to mitigate earthquake-induced soil liquefaction (O'Donnell et al., 2017a, b; Pham et al., 2018). MIDP mitigates liquefaction in two ways: generation of nitrogen gas (N₂) and mineral precipitation (usually calcium carbonate, CaCO₃). The impact of MIDP on the soil system is illustrated in Fig. 1. The generated N₂ desaturates the soil, increasing its compressibility and reducing the increase in pore water pressure during cyclic loading, which is the root cause of earthquake-induced liquefaction. Carbonate precipitation increases the soil strength, thereby increasing the intensity of earthquake sharing necessary to trigger liquefaction. A primary benefit of MIDP for liquefaction mitigation is that, being non-disruptive, it can be used underneath existing structures (O'Donnell et al., 2017a; Hall, 2021). Trillions of dollars of existing infrastructure is at risk due to the potential for liquefaction, and currently that risk cannot be mitigated in a cost-effective way. MIDP is currently being evaluated at different experimental scales as a solution to this problem (O'Donnell et al., 2017a, b; Moug et al., 2022).

Figure 1 Pore-scale illustration of MIDP in water (dark blue)-saturated sand (light yellow), detailing the gas (light blue), calcium carbonate (white), and biomass (green) production.

For environmental and economic reasons, we assume that local source water will be used to prepare the MIDP treatment solution (i.e., dissolve the substrates needed to induce denitrification) in practice. Therefore, in developing this technique we must consider how naturally occurring constituents in the source water may affect denitrification. For example, a competing electron acceptor (e.g., sulfate) may consume the electron donor (e.g., acetate), leading to incomplete denitrification and the formation of unwanted products (e.g., hydrogen sulfide). In addition, the biogeochemical reactions during MIDP result in changes to alkalinity and pH that may alter process kinetics and final MIDP products.

To quantify the impact of source-water composition and to investigate how complex biogeochemical interactions occur during field application, we developed a next-generation biogeochemical model of MIDP. This mathematical model expands upon previous modeling of MIDP, which did not consider the impact of source water on MIDP or the impact of MIDP on the aqueous subsurface environment (Pham, 2017; O'Donnell et al., 2019). Our next-generation MIDP model includes all essential biogeochemical processes based on the constituents commonly observed in the natural groundwater environments, substrates added to stimulate MIDP, and mechanisms that lead to desaturation and precipitation: e.g., N₂-gas formation, acid–base speciation, and CaCO₃ precipitation. Since MIDP often is deployed in coastal areas (due to the prevalence of liquefiable soil deposits in this environment), we include conditions typical for coastal seawater in our model.

2 Model foundation

The next-generation model builds upon previous MIDP models (Pham, 2017; O'Donnell et al., 2019) but broadens the range of processes considered by the model. Our next-generation model considers microbial growth and decay, alternative microbial metabolic processes, gas production, mineral-solids' production, alkalinity and pH, microbial inhibition, and desaturation and precipitation in both fresh water and coastal environments. A comparison of the components and processes considered by the two earlier MIDP models and our next-generation model is provided in Table S1 in the Supplement.

The next-generation model was constructed in MATLAB (Little and Moler, 2017), and the code and necessary files are publicly available online at https://doi.org/10.5281/zenodo.7410676. The modeling equations (e.g., microbial growth, CaCO₃ precipitation, and biogenic gas evolution) were programmed within the original, publicly available van Turnhout toolbox, a general-form mechanistic batch model for environmental systems that considers species in the gas, liquid, and solid phases (van Turnhout et al., 2016). The van Turnhout toolbox is a program that includes a system of ordinary differential equations that model biogeochemical reactions. The van Turnhout toolbox is coupled with ORCHESTRA, an extensive database of established geochemical equilibria based on MINTEQ, to simulate chemical speciation during said modeling (Meeussen, 2003). The MIDP-specific biogeochemical model components (i.e., stoichiometry, type of inhibition and kinetics, potential chemical species) were specified in an input spreadsheet that the program accesses. The degree of saturation and percent (by weight) of mineral precipitation were calculated outside of the van Turnhout toolbox using model results, as discussed in Sect. 3.2 of this paper.

The toolbox's logic flow and calculation sequence are as follows (Meeussen, 2003; van Turnhout et al., 2016), using H₂CO₃, HCO${}_{\mathrm{3}}^{-}$, CO${}_{\mathrm{3}}^{\mathrm{2}-}$, H⁺, and OH⁻ to illustrate the process for the carbonate system.

At t=0, the program loads the input concentrations' file, which includes the concentration of all total species (e.g., H₂CO₃ representing DIC, H⁺) and the stoichiometry for metabolic and kinetic reactions: e.g.,

$$\begin{array}{rl}& \mathrm{0.222}{\mathrm{NO}}_{\mathrm{3}}^{-}+\mathrm{0.125}{\mathrm{C}}_{\mathrm{2}}{\mathrm{H}}_{\mathrm{3}}{\mathrm{O}}_{\mathrm{2}}^{-}+\mathrm{0.146}{\mathrm{H}}^{+}\\ & \to \mathrm{0.202}{\mathrm{NO}}_{\mathrm{2}}^{-}+\mathrm{0.147}{\mathrm{H}}_{\mathrm{2}}{\mathrm{CO}}_{\mathrm{3}}\\ & +\mathrm{0.103}{\mathrm{CH}}_{\mathrm{1.8}}{\mathrm{O}}_{\mathrm{0.5}}{\mathrm{N}}_{\mathrm{0.2}}+\mathrm{0.021}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}.\end{array}$$

Ordinary differential equations are used to determine compound consumption and production based on the reaction stoichiometry and kinetic equations (e.g., precipitation, biotransformation, and mass transfer) at each time step. For example, the graphic in Fig. 2 illustrates that, as C₂H₃O${}_{\mathrm{2}}^{-}$ is consumed from microbial consumption, H₂CO₃ is produced.

Figure 2 Illustration of consumption of the consumption of C₂H₃O${}_{\mathrm{2}}^{-}$ and production of H₂CO₃ during MIDP.

At each time step, the following set of linear equations is solved to determine the relative derived concentrations of [H₂CO₃], $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$, $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$, [H⁺], and [OH⁻] from H₂CO₃ produced in the previous steps. This is done in the ORCHESTRA biochemical module.

For the mass balance equations, the left side of the equation is the total dissolved inorganic carbon, H₂CO₃, from the stoichiometry described in steps 1 and 2. The right side shows the derived concentrations of species as a result of speciation, and these are indicated with square brackets.

$${\mathrm{H}}_{\mathrm{2}}{\mathrm{CO}}_{\mathrm{3}}=\left[{\mathrm{H}}_{\mathrm{2}}{\mathrm{CO}}_{\mathrm{3}}\right]+\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]+\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$$

Electroneutrality principles were considered such that all potentially produced charged species related to this balance were included.

$$\left[{\mathrm{H}}^{+}\right]=\left[{\mathrm{OH}}^{-}\right]+\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]+\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$$

Acid equilibrium for H₂CO₃:

$$K_{\mathrm{acid}}={\displaystyle \frac{\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]\left[\mathrm{2}{\mathrm{H}}^{+}\right]}{\left[{\mathrm{H}}_{\mathrm{2}}{\mathrm{CO}}_{\mathrm{3}}\right]}}.$$

Acid equilibrium for $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$:

$$K_{\mathrm{acid}}={\displaystyle \frac{\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]\left[{\mathrm{H}}^{+}\right]}{\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]}}.$$

Water equilibrium:

$$K_{\mathrm{w}}=\left[{\mathrm{OH}}^{-}\right]\left[{\mathrm{H}}^{+}\right]=\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{14}}.$$

pH is calculated based on the derived H⁺ concentration based on solving simultaneously.

While the carbonate system is used here as an example, this stepwise process is used for all acid–base species and considers the total system set of reactions and species to achieve equilibrium. The total system electroneutrality considered in the model for all considered species is as follows:

$$\begin{array}{rl}& \left[{\mathrm{H}}^{+}\right]+\left[{\mathrm{NH}}_{\mathrm{4}}^{+}\right]+\left[{\mathrm{Ca}}^{\mathrm{2}+}\right]+\left[{\mathrm{CaHCO}}_{\mathrm{3}}^{+}\right]+\left[{\mathrm{CaOH}}^{+}\right]\\ & +\left[{\mathrm{CaC}}_{\mathrm{2}}{\mathrm{H}}_{\mathrm{3}}{\mathrm{O}}_{\mathrm{2}}^{+}\right]+\left[{\mathrm{Fe}}^{\mathrm{3}+}\right]+\left[{\mathrm{Fe}}^{\mathrm{2}+}\right]+\left[{\mathrm{FeOH}}^{+}\right]\\ & +\left[\mathrm{Fe}(\mathrm{OH}{)}_{\mathrm{2}}^{+}\right]+\left[{\mathrm{Na}}^{+}\right]\\ & =\left[{\mathrm{OH}}^{-}\right]+\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]+\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]+\left[{\mathrm{NO}}_{\mathrm{3}}^{-}\right]+\left[{\mathrm{NO}}_{\mathrm{2}}^{-}\right]\\ & +\left[{\mathrm{C}}_{\mathrm{2}}{\mathrm{H}}_{\mathrm{3}}{\mathrm{O}}_{\mathrm{2}}^{-}\right]+\left[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\right]+\left[{\mathrm{HSO}}_{\mathrm{4}}^{-}\right]+\left[{\mathrm{HS}}^{-}\right]+\left[{\mathrm{Cl}}^{-}\right].\end{array}$$

The compounds listed above are defined in the input spreadsheet as possible species. The model then uses ORCHESTRA to determine the concentration of each species based on the biogeochemical reactions and system equilibrium.

3 Model principles

Denitrification is a multi-step process of nitrogen-species reduction. During each reduction step, energy and biomass are produced when paired with oxidation of an electron donor that produces a thermodynamically favorable reduction-oxidation (redox) reaction. The four steps of nitrogen reduction in denitrification conform to the following reduction half reactions, each consuming two or one electron equivalent (e⁻) (Rittmann and McCarty, 2020).

Nitrate reduction to nitrite:

$$\mathrm{0.5}{\mathrm{NO}}_{\mathrm{3}}^{-}+\mathrm{2}{\mathrm{H}}^{+}+\mathrm{2}{e}^{-}\to \mathrm{0.5}{\mathrm{NO}}_{\mathrm{2}}^{-}+\mathrm{0.5}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}.$$

Nitrite reduction to nitric oxide:

$${\mathrm{NO}}_{\mathrm{2}}^{-}+\mathrm{2}{\mathrm{H}}^{+}+e^{-}\to \mathrm{NO}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}.$$

Nitric oxide reduction to nitrous oxide:

$$\mathrm{NO}+{\mathrm{H}}^{+}+e^{-}\to \mathrm{0.5}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}+\mathrm{0.5}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}.$$

Nitrous oxide reduction to dinitrogen:

$$\mathrm{0.5}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}+{\mathrm{H}}^{+}+e^{-}\to \mathrm{0.5}{\mathrm{N}}_{\mathrm{2}}+\mathrm{0.5}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}.$$

In the model developed herein, the four steps were simplified to two steps, nitrate to nitrite and nitrite to dinitrogen gas. The reductions of nitric oxide and nitrous oxide are assumed to occur completely, as they are thermodynamically favorable (found in the Supplement). In contrast, the accumulation of NO${}_{\mathrm{2}}^{-}$ must be considered explicitly because it is a denitrification inhibition “bottleneck” and presents a risk to human health (Pham et al., 2018; Almeida et al., 1995). The variables for the equations used in the model described in the next several sections are detailed in Table 1.

Table 1 Notation, units, and description for variables used in the biogeochemical model.

3.1 Microbial metabolism, growth, and decay

The processes within the model follow Monod kinetics, represented as multiplicative dual-substrate limitation (O'Donnell et al., 2019; Bae and Rittmann, 1996).

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}{C}_{\mathrm{d}}}{\mathrm{d}t}}=-\widehat{q}{X}_{\mathrm{a}}{\displaystyle \frac{C_{\mathrm{d}}}{K_{\mathrm{d}}+C_{\mathrm{d}}}}\cdot {\displaystyle \frac{C_{\mathrm{a}}}{K_{\mathrm{a}}+C_{\mathrm{a}}}}{I}_{\mathrm{i}}\end{array}$$

Equation (1) considers the electron-donor substrate (C_d) and three electron-acceptor substrates (C_a): NO${}_{\mathrm{3}}^{-}$ and NO${}_{\mathrm{2}}^{-}$ for denitrifying bacteria and SO${}_{\mathrm{4}}^{\mathrm{2}-}$ for sulfate-reducing bacteria. For preliminary analysis, we assumed an initial denitrifier biomass concentration of 0.5 mmol L⁻¹ and sulfate-reducing biomass concentration of 0.25 mmol L⁻¹. The values of the constants ($\widehat{q}$, K_d, and K_a) are in the Supplement, along with the derivations of said parameters based on Rittmann and McCarty (2020). The Supplement also includes preliminary model results, assuming a reduced initial denitrifier biomass concentration of 0.05 mmol L⁻¹ and sulfate-reducing biomass concentration of 0.025 mmol L⁻¹. The comparison of the findings detailed in Sect. 4, the results section of this paper, and those found in the Supplement illustrate the influence of the initial biomass concentration on the model results. The main impact of lowering the initial biomass concentration was a delayed start-up time for the microbiological processes, though the trends remained the same for all modeled scenarios. These derivations were used to determine reaction stoichiometry, true yield (Y), and μ_max (maximum specific growth rate) for all electron-donor and electron-acceptor pairs and the nitrogen source. The values of the kinetic and stoichiometric parameters are detailed in Table 2. The inhibition factor I_i is described in a later section of this paper.

Values of half-maximum-rate concentrations (K_d and K_a) in the literature show variability for each electron-donor and electron-acceptor pair due to the wide range of environments of the microorganisms (e.g., sediment, estuarine water, wastewater) and the high degree of diversity of microorganisms able to carry out these reactions (Abdul-Talib et al., 2002; Papaspyrou et al., 2014; Vavilin and Rytov, 2015). Table 3 details the constants we used as representative values for each K_d and K_a (for Eq. 1) based on relevant electron-donor and electron-acceptor pairs and sources of those values. While these values are not specific to a coastal seawater environment, they have been experimentally validated.

Table 2 Reaction stoichiometry, yield (Y), and maximum specific growth rates (μ_max) expected during MIDP, considering acetate to be the electron donor and natural electron acceptors. Units for all parameters are in Table 1.

Table 3 Half-maximum-rate concentrations, K_d and K_a, used for each electron-donor and electron-acceptor pair.

Microbial growth within the model is represented via reaction kinetics and stoichiometry expressed in Eq. (2):

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\mathrm{d}{X}_{\mathrm{a}}}{\mathrm{d}t}}=X_{\mathrm{a}}Y\widehat{q}-b.\end{array}$$

Biomass yields (Y) are listed in Table 2. For sulfate-reducing bacteria, b was set to 0.03 d⁻¹, whereas it was set to 0.05 d⁻¹ for denitrifiers (Rittmann and McCarty, 2020). As a result of decay, NH${}_{\mathrm{4}}^{+}$ is released and can serve as a nitrogen source for denitrification. Since NH${}_{\mathrm{4}}^{+}$ is thermodynamically favorable over NO${}_{\mathrm{3}}^{-}$ as a nitrogen source, it is used first before NO${}_{\mathrm{3}}^{-}$ during denitrification using a user-defined switch. We used the inhibition function described in Sect. 3.2 as the switch to interrupt biomass from using NO${}_{\mathrm{3}}^{-}$ as the nitrogen source in the presence of NH${}_{\mathrm{4}}^{+}$. Decay involves endogenous respiration, and we assumed that 80 % of decayed biomass is available as an acetate for metabolism, while 20 % becomes inert biomass (Rittmann and McCarty, 2020). The stoichiometry for decay is

$$\begin{array}{rl}& \mathrm{0.238}{\mathrm{CH}}_{\mathrm{1.8}}{\mathrm{O}}_{\mathrm{0.5}}{\mathrm{N}}_{\mathrm{0.2}}+\mathrm{0.012}{\mathrm{H}}_{\mathrm{2}}{\mathrm{CO}}_{\mathrm{3}}+\mathrm{0.095}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\\ & \to \mathrm{0.125}{\mathrm{C}}_{\mathrm{2}}{\mathrm{H}}_{\mathrm{3}}{\mathrm{O}}_{\mathrm{2}}^{-}+\mathrm{0.048}{\mathrm{NH}}_{\mathrm{4}}^{+}+\mathrm{0.077}{\mathrm{H}}^{+}.\end{array}$$

3.2 Inhibition

The van Turnhout toolbox has the capability to model different inhibition mechanisms, but we only used non-competitive inhibition during denitrification because the enzymes that perform nitrate, nitrite, and sulfate reduction are different and not self-inhibitory (Glass and Silverstein, 1998). Denitrification inhibition, which slows nitrate and nitrite reduction rates (Glass et al., 1997), was included for the reduction of nitrate to nitrite and nitrite to N₂ gas. I_i is a general term for inhibition of either step, with i indicating which reaction. The form of I_i, shown in Eq. (3) is for non-competitive inhibition, and the inhibition coefficients for each inhibitor are found in Table 4:

$$\begin{array}{}\text{(3)}& I_{\mathrm{i}}={\displaystyle \frac{K_{\mathrm{i}}}{K_{\mathrm{i}}+C_{\mathrm{i}}}}.\end{array}$$

Table 4 Non-competitive inhibition coefficients (K_i).

^a Unacclimated environments were DI and drinking water. ^b Acclimated environments were groundwater and seawater.

Although several inhibitors could affect MIDP, HNO₂ is the most important inhibitor of the MIDP process (Lilja and Johnson, 2016). Significant inhibition (95 % rate reduction) of overall denitrification has been reported at 0.04 mg HNO₂ L⁻¹, which primarily impacted the intermediate NO${}_{\mathrm{2}}^{-}$ reduction step (Glass et al., 1997; Abeling and Seyfried, 1992), and a 60 % decrease in NO${}_{\mathrm{3}}^{-}$ reduction at 0.08 mg HNO₂ L⁻¹ was also reported (Ma et al., 2010). Within the model, HNO₂ inhibits NO${}_{\mathrm{3}}^{-}$ and NO${}_{\mathrm{2}}^{-}$ reduction using different inhibition coefficients (Table 4). The inhibition by HNO₂ is driven by pH speciation because NO${}_{\mathrm{2}}^{-}$ is dominant at a pH of 3.4 and higher and HNO₂ is negligible for pH ≥ 7.6. However, only a small concentration of HNO₂ can have a significant impact on denitrification, which underscores the importance of pH and the accumulation of the intermediate NO${}_{\mathrm{2}}^{-}$.

Inhibition between nitrate and nitrite reductions has been identified, with the presence of nitrate having a larger effect on nitrite reduction than nitrite on nitrate reduction (Lilja and Johnson, 2016; Glass et al., 1997; Almeida et al., 1995; Soto et al., 2007). Nitrite accumulation increases in the presence of nitrate until nitrate is depleted, such that nitrite reduction becomes the dominant process (Glass and Silverstein, 1998). When only nitrite remains, the rate of nitrite reduction increases. However, others have described that, as long as the electrons are adequately provided by the electron donor, competitive inhibition between nitrate and nitrite reductions is not significant (Soto et al., 2007; Ma et al., 2010; van den Berg et al., 2017). Therefore, the model does not include competitive inhibition, although it naturally includes competition for the electron donor between nitrate and nitrite reductions through thermodynamic favorability of nitrate reduction over nitrite reduction.

The model applies different non-competitive inhibition constants for salinity (as NaCl) for nitrate and nitrite reduction because nitrite reduction is more sensitive to salinity than nitrate. Because the magnitude of inhibition depends on experimental conditions and adaptation of the microorganisms, the value of K_i may differ for local conditions (Krishna Rao and Gnanam, 1990).

Hydrogen sulfide (H₂S) also can be inhibitory to denitrification (Pan et al., 2019). Nitrate, nitrite, and N₂O reductions have been inhibited by H₂S, though the extent and sensitivity of reduction in the presence of H₂S was experiment-dependent (Senga et al., 2006; Pan et al., 2013; Tugtas and Pavlostathis, 2007; Liang et al., 2020; Cardoso et al., 2006). Within the model, one aqueous-phase H₂S-inhibition constant was used for both NO${}_{\mathrm{3}}^{-}$ and NO${}_{\mathrm{2}}^{-}$ reduction steps.

A pH < 6 can significantly slow denitrification (Glass and Silverstein, 1998) by inhibiting enzyme activity (Šimek and Cooper, 2002) and microbial growth (Estuardo et al., 2008). When the pH goes higher than 8, enzyme activity also can be impeded, leading to reduced denitrification rates or incomplete denitrification. Incidents of a high pH often are temporary, as CaCO₃ precipitation in MIDP buffers the pH (Salek et al., 2015). The benefit of including a pH-inhibition function when predicting denitrification has been demonstrated, but the values of their governing parameters are environment-specific and require fitting (Estuardo et al., 2008). Within the model, we only considered the indirect net effect of pH through HNO₂ inhibition, which does not require environment-specific parameters because the concentration of HNO₂ is automatically calculated within the model structure.

3.3 Biogenic gas production

O'Donnell et al. (2019) considered the production of N₂ and CO₂ during denitrification but did not consider the varying subsurface stresses that would influence phase transfer. The relative concentrations of the produced biogenic gas can affect the distribution of gas at depth, since the gases have different solubilities, as well as different stoichiometries for electron-donor consumption.

Our next-generation MIDP model includes mass-transfer kinetics for transfers of N₂, CO₂, and H₂S from the aqueous phase to the gas phase (or from the gas phase). We considered gas-phase transfer kinetics because assuming instantaneous gas-phase transfer clearly would be an oversimplification, based on the review on mass transfer of biologically driven gas production completed by Kraakman et al. (2011). N₂, CO₂, and H₂S concentrations were modeled in the aqueous and gas phases. The rate of transfer of a gaseous compound from the aqueous phase to (or from) the gas phase, v_i[g], depends on the gas's degree of super-saturation and a mass-transfer-rate coefficient (Salek et al., 2015):

$$\begin{array}{}\text{(4)}& v_{i\left[\mathrm{g}\right]}=k_{L}a\left(C_{i\left[\mathrm{g}\right]}-{\displaystyle \frac{C_{i\left[\mathrm{aq}\right]}RT}{K_{\mathrm{H}}}}\right).\end{array}$$

We assigned k_La values for N₂, CO₂, and H₂S of 1 d⁻¹ (Shin et al., 2003), though the values can vary widely based on porous-medium conditions and temperature. We did not include pore-scale kinetics. The aqueous concentrations of CO₂ and H₂S depend on the pH, as described below.

The biogenic gas volume needed to achieve a target level of desaturation (S_g) by N₂ ([N₂]_g) and CO₂ ([CO₂]_g) was determined by

$$\begin{array}{}\text{(5)}& [{\mathrm{N}}_{\mathrm{2}}{]}_{\mathrm{g}}+[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{g}}={\displaystyle \frac{p{S}_{\mathrm{g}}}{RT}}\end{array}$$

in which p was assumed to be equal to the sum of the hydraulic pressure at the treatment depth (7.6 m in an upcoming example) and the atmospheric pressure. Gas-phase H₂S was not included in the desaturation calculations because its solubility is much higher than N₂ and CO₂.

Equation (6) describes the amount of input NO${}_{\mathrm{3}}^{-}$ required for desaturation by N₂ and CO₂ ([NO${}_{{\mathrm{3}}_{\mathrm{d}}}^{-}$], mol${}_{{\mathrm{NO}}_{\mathrm{3}}}$ L${}_{\mathrm{pore}}^{-\mathrm{1}}$) at the deepest target treatment depth, which is the lowest depth of the treated zone. The depth increases the pressure ($p_{{\mathrm{N}}_{\mathrm{2}}}$ and $p_{{\mathrm{CO}}_{\mathrm{2}}}$) and the needed amount of gas production to exceed the solubility threshold ($K_{\mathrm{H},{\mathrm{N}}_{\mathrm{2}}}$ and $K_{\mathrm{H},{\mathrm{CO}}_{\mathrm{2}}}$) and enter the gas phase, according to Henry's law. The equation considers the amount of gas needed to overcome the solubility threshold to achieve the target level of desaturation (Hall et al., 2018; Pham, 2017):

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\left[{\mathrm{NO}}_{{\mathrm{3}}_{\mathrm{d}}}^{-}\right]& ={\displaystyle \frac{\left(\frac{[{\mathrm{N}}_{\mathrm{2}}{]}_{\mathrm{g}}}{l}+\frac{p_{{\mathrm{N}}_{\mathrm{2}}}}{K_{\mathrm{H},{\mathrm{N}}_{\mathrm{2}}}}\right)Y_{{\mathrm{NO}}_{\mathrm{3}}^{-}}}{Y_{{\mathrm{N}}_{\mathrm{2}}}}}\\ & +{\displaystyle \frac{\left(\frac{[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{g}}}{l}+\frac{p_{{\mathrm{CO}}_{\mathrm{2}}}}{K_{\mathrm{H},{\mathrm{CO}}_{\mathrm{2}}}}\right)Y_{{\mathrm{NO}}_{\mathrm{3}}^{-}}}{Y_{{\mathrm{CO}}_{\mathrm{2}}}}}.\end{array}\end{array}$$

$K_{\mathrm{H},{\mathrm{N}}_{\mathrm{2}}}$ was set to 1600 and $K_{\mathrm{H},{\mathrm{CO}}_{\mathrm{2}}}$ to 29 (L atm mol⁻¹). We determine the amount of nitrate needed to achieve target desaturation levels at the lowest depth because greater depth requires a higher concentration of gas to achieve target desaturation levels, as the pressures are at their respective maxima.

Equation (7) was used to determine the biogenic gas volume (V_g, L_gas L${}_{\mathrm{tot}}^{-\mathrm{1}}$),

$$\begin{array}{}\text{(7)}& V_{\mathrm{g}}={\displaystyle \frac{S_{\mathrm{g}}RT\mathit{\phi }}{p}}.\end{array}$$

3.4 Solids' precipitation and dissolution

Calcium carbonate precipitation occurs when dissolved inorganic carbon (DIC), produced from microbial substrate conversion of the electron donor, exceeds the solubility of CaCO₃ for the concentration of Ca²⁺ present. The stoichiometry for CaCO₃ precipitation is

$${\mathrm{Ca}}^{\mathrm{2}+}+{\mathrm{H}}_{\mathrm{2}}{\mathrm{CO}}_{\mathrm{3}}\to {\mathrm{CaCO}}_{\mathrm{3}}+\mathrm{2}{\mathrm{H}}^{+}.$$

At the beginning of each run, the concentration of species available for precipitation is calculated as their dissolved form as discussed in Sect. 2, and then solids' precipitation is determined. The van Turnhout toolbox considers precipitation based on equilibrium calculations from the ORCHESTRA module (Meeussen, 2003). This assumption is valid when the rates of precipitation and dissolution of minerals are much faster than the phase transfer between the aqueous and solid phases (Salek et al., 2015). Previous MIDP modeling did not consider precipitation kinetics but assumed instantaneous equilibrium (Pham, 2017; O'Donnell et al., 2019). Instantaneous equilibrium may be an over-simplification for environmental conditions (Singurindy et al., 2004) in which the mechanisms of crystal nucleation, crystal growth, and mass transfer of reactants to the contact point of crystal growth are important (Rittmann et al., 2002). Therefore, we included precipitation and dissolution kinetics in the next-generation model.

The model considers first-order precipitation and dissolution kinetics with respect to the Ca²⁺ and CO${}_{\mathrm{3}}^{\mathrm{2}-}$ concentrations (Rittmann et al., 2002; Chou et al., 1989):

$$\begin{array}{}\text{(8)}& R_{\mathrm{p}}=\mathit{ka}\left(\mathrm{1}-{\displaystyle \frac{K_{\mathrm{sp}}}{\left[{\mathrm{Ca}}^{\mathrm{2}+}\right]\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]}}\right)\left[{\mathrm{Ca}}^{\mathrm{2}+}\right].\end{array}$$

K_sp was set to $\mathrm{1.83}\times {\mathrm{10}}^{-\mathrm{8}}$ mol² L⁻² at 25 ^∘C for CaCO₃. ka is a combined kinetic coefficient because it is difficult to separate mass transfer kinetics, crystal growth rate, and solid surface area (Rittmann et al., 2002; Spanos and Koutsoukos, 1998; Rittmann et al., 2003). ka can have a large range depending on the environment and the ease of establishing precipitation nucleation points. We assumed ka was 100 d⁻¹, though this value should be used as a fitting parameter subject to experimental validation. Precipitation was implemented using the van Turnhout toolbox's method for biochemical reactions, ka was specified as a reaction rate, and K_sp was among the governing input parameters.

Equation (9) was used to determine the amount of substrate needed to achieve a target precipitation level, which is determined by the ratio between mass of precipitated CaCO₃ and mass of the soil solids ([CaCO₃], kg CaCO₃ kg soil⁻¹).

$$\begin{array}{}\text{(9)}& [{\mathrm{NO}}_{\mathrm{3}}^{-}{]}_{\mathrm{c}}={\displaystyle \frac{\left[{\mathrm{CaCO}}_{\mathrm{3}}\right]{\mathit{\rho }}_{\mathrm{soil}}{Y}_{{\mathrm{NO}}_{\mathrm{3}}^{-}}}{{\mathrm{eu}}_{{\mathrm{CaCO}}_{\mathrm{3}}}{Y}_{{\mathrm{CaCO}}_{\mathrm{3}}}}}\end{array}$$

The stoichiometric coefficients considered the total amount of input NO${}_{\mathrm{3}}^{-}$ and produced H₂CO₃ for the total assumed two-step denitrification process. The DIC available for precipitation to provide $Y_{{\mathrm{CaCO}}_{\mathrm{3}}}$ is estimated based on pH-driven speciation at each time step.

While we only considered calcium carbonate precipitation, the model has the flexibility to model precipitation of other minerals. The user would need to add in separate equations to model precipitation kinetics based on the reactants, desired products, and the ka and K_sp values appropriate for the desired precipitation reaction and product.

3.5 Determining pH

Because pH governs the concentration of important aqueous species based on acid–base speciation, the pH influences many of the geochemical reactions involved in MIDP. The pH was determined using the geochemical equilibrium software ORCHESTRA, which is part of the van Turnhout toolbox. ORCHESTRA uses a mass balance on all species within the system and the products of rate-dependent processes as a function of time (i.e., kinetic, biogeochemical, and phase transfer processes). At each time step, the program performs a mass balance on all complexed species and their fate (e.g., transformed through microbial processes, precipitation, gas-phase transfer) (van Turnhout et al., 2016; Meeussen, 2003).

4 Case study MIDP behavior seawater conditions: model results and discussion

To demonstrate the capabilities of the model, we illustrate MIDP behavior when targeting desaturation for liquefaction mitigation in a coastal geochemical environment. In this case study, we demonstrate the impact of precipitation on biochemical reactions and the resulting multi-phase products and by-products resulting from MIDP and other environmental biogeochemical processes (e.g., sulfate reduction). However, we only consider desaturation as a target treatment mechanism and do not model an MIDP treatment recipe optimized for precipitation as a liquefaction-mitigation mechanism.

Table 5 details the chemical characteristics used to simulate coastal groundwater conditions, which were assumed to have the same characteristic of seawater due to intrusion (hereafter, referred to as “seawater”). The treatment substrate was added to the baseline level of these components. Table 5 does not reflect the varying concentrations of calcium acetate and calcium nitrate used in the treatment recipes, which are described later. We based the target treatment zone's soil properties on a case study of microbial desaturation via denitrification in Portland, Oregon, presented by Moug et al. (2022). The deepest target treatment depth was 7.6 m. We assumed typical values of total unit weight, dry density, and porosity for uniform clean sand for the soil (Christopher et al., 2006): total unit weight of 19.5 kN m⁻³ (dry unit weight of 15.6 kN m⁻³; bulk density of 1950 kg m⁻³) and porosity of 0.39.

Table 5 Chemical characteristics assumed for a coastal seawater environment.

^a Average of measured values (Alin et al., 2017). ^b Reference composition of “standard seawater” from and calculated for pH = 7.61 for acid–base species (Millero et al., 2008; European Commission, Directorate General for Research, 2011). ^c Bruland et al. (2001).

Table 6 MIDP treatment recipes for each modeled condition.

Figure 3 Gas volumes normalized to the total soil volume (a, c, e) and degree of saturation by depth (b, d, f) for the simulated coastal seawater conditions. The desaturation target was 10 % or a saturation ratio of 90 %.

The reported desaturation levels required to increase the cyclic shear resistance for liquefaction mitigation range between 2 % and 10 % (He and Chu, 2014; O'Donnell et al., 2017a). We chose 10 %, which is at the high end of the mitigation range but well below the desaturation level at which gas starts to migrate upward or spread laterally, reported to occur at 20 % for poorly graded (i.e., uniform) fine sands (Pham, 2017).

Figure 4 Water-quality results for 28 d of MIDP in coastal seawater conditions targeting a desaturation level of 10 % in three different treatment recipes: empirically matched, 25 % excess acetate, and 25 % reduced acetate.

Figure 5 Water-quality results for the first 3 d of MIDP in coastal seawater conditions targeting a desaturation level of 10 % in three different treatment recipes: empirically matched, 25 % excess acetate, and 25 % reduced acetate.

Figure 6 DIC concentration and CaCO₃ precipitated during the first 3 d of MIDP in coastal seawater conditions targeting a desaturation level of 10 %.

Following Eq. (5), 7.10 mmol L⁻¹ of total N₂ gas is required to meet a minimum target desaturation level of 10 % throughout the entire treatment zone (assuming only desaturation via N₂ gas). To meet the target desaturation, we estimated the treatment recipe to be 22.4 mmol L⁻¹ of nitrate (1.84 g calcium nitrate L⁻¹) and 32.1 mmol L⁻¹ of acetate (2.54 g calcium acetate L⁻¹) using Eq. (6). However, based on background levels of nitrate and nitrite and the use of released ammonium as a nitrogen source, these levels were adjusted to establish the treatment recipe detailed in Table 6. The adjusted values were set to result in complete denitrification (i.e., no residual nitrate or nitrite that relies on bacterial decay as the electron donor) and to not exceed 0.1 mmol L⁻¹ of acetate after complete denitrification. We compared the impact of varying the input levels of acetate (as calcium acetate) on the MIDP treatment to the matched treatment recipe. We tested the impact of addition of an extra 25 % of acetate over our original estimations, referred to as the “excess acetate” treatment recipe, and 25 % less acetate from our original estimation, referred to as the “reduced acetate” treatment recipe. For the excess- and reduced-acetate comparisons, we did not adjust the input levels of nitrate (as calcium nitrate) from our original estimations of 22.4 mmol L⁻¹ (1.84 g calcium nitrate L⁻¹). We assumed that each treatment recipe was added in one application (i.e., not in a continuous flow-through manner).

The results of the matched treatment recipe on the subsurface gas volume and saturation profile are shown in the top two panels of Fig. 3. For the coastal seawater conditions, the target desaturation level of 10 % at 7.6 m (or a degree of saturation of 90 %) was achieved by N₂ generation in approximately 2.0 d. The amount of CO₂ produced did not reach its saturation threshold, and CO₂ did not contribute to desaturation at any of the modeled depths. The difference in volume of gas at the different levels is due to the increase in hydrostatic pressure with depth.

The middle panels of Fig. 3 indicate that adding excess acetate increased the degree of saturation at 7.6 m, which is shown by less than 90 % saturation at 7.6 m. In contrast, the bottom two panels show that adding less acetate slowed N₂ generation so that 90 % saturation was not reached at 7.6 m in 60 d.

The impacts of the different MIDP treatment recipes on water quality are shown in Fig. 4, and the initial 5 d of treatment for each modeled scenario is highlighted in Fig. 5. With the matched-treatment recipe, almost all nitrate and nitrite were consumed by 1.6 d, and less than 10⁻⁴ mol L⁻¹ of acetate remained (Fig. 5). Nitrite accumulation was transient and modest (3.3 mM at its peak); thus, complete denitrification was achieved with this treatment recipe for coastal seawater conditions (Fig. 5). After 1.6 d, sulfate reduction began and continued to occur, driven by microbial endogenous respiration, at a small rate that resulted in the production of more total CO₂ and total H₂S than the matched treatment over time (Fig. 4).

With the excess-acetate recipe, all the nitrate was completely reduced, although small, transient accumulations of nitrite and nitrous acid occurred (Fig. 5). As expected, not all the acetate was consumed with the excess-acetate recipe, and the remaining acetate led to sulfate reduction and the highest amount of produced H₂S of the three modeled scenarios (Fig. 4). Additional N₂ was produced because of the higher amounts of nitrate in the treatment recipe, leading to a level of desaturation at 7.6 m that exceeded the target 10 % value.

For the reduced-acetate test, approximately 20 % of the input nitrate remained after all the acetate was consumed, and this residual nitrate was slowly utilized beyond 1.5 d through biomass endogenous decay (Fig. 4). The peak amount of accumulated nitrite was not as high as the other conditions because of the overall limited nitrate reduction, but some nitrite accumulation remained throughout the modeled 28 d due to the lack of acetate. The dip and quick increase in nitrous acid around 1.5 to 1.7 d (Fig. 5) were due to the shift in electron donor from input acetate to bacterial decay. Not enough electron donor was available to reduce all of the input nitrate and the accumulated nitrite after 28 d (Fig. 4), even though bacterial decay caused added denitrification. N₂ produced in the reduced-acetate condition did not meet the 10 % desaturation threshold at 7.6 m at the end of the 28 d modeled period.

Nitrite reduction produced most of the base, which is demonstrated by the spike in pH after 1 and 1.4 d in all scenarios (Fig. 5), when the rate of nitrite reduction was at its maximum in each treatment. In the matched case, the pH returned to being circumneutral after 1.8 d due to the precipitation of CaCO₃, which consumes base. This trend is reinforced by rapid Ca²⁺ consumption in the early treatment time period (Fig. 5), along with production of CaCO₃ in the first ∼ 2 d, shown in the right panel of Fig. 6. DIC production in the excess-acetate treatment lagged the matched recipe and was slightly quicker in the reduced-acetate case, which also is seen with the pH trends in Figs. 4 and 5. However, the overall consumption of DIC and subsequent precipitation were slower in the excess- and reduced-acetate treatments, which correspond to the longer time for the pH to reach approximately neutral levels. The additional CaCO₃ precipitated with the excess acetate resulted from the excess of input calcium, since acetate was added as calcium acetate.

Microbial decay coupled to sulfate reduction also produced more DIC (left panel of Fig. 6), promoting additional CaCO₃ precipitation after the completion of denitrification (after day 2 in Fig. 6). In the excess-acetate and matched treatments, an increase in DIC after 15 d was due to the increased rate of sulfate reduction because of the additional electron donor (Fig. 4).

In summary, the simulations show that the matched-acetate recipe optimized MIDP treatment targeted at desaturation for coastal seawater conditions by maximizing the desired outcome (i.e., N₂ production for desaturation) while minimizing undesired by-products (e.g., nitrite and nitrous-oxide accumulation from incomplete denitrification, residual acetate, and H₂S from sulfate reduction).

5 Conclusion

The next-generation biogeochemical model expanded our previous biogeochemical models for MIDP by considering microbial stoichiometry and kinetics for two steps of denitrification and for sulfate reduction. The next-generation model also includes gas–liquid mass-transfer kinetics for N₂ and CO₂, CaCO₃ precipitation kinetics, microbial competition, and inhibition by HNO₂, salinity, and sulfide. Model simulations demonstrated that adding nitrate and acetate using a properly matched recipe led to rapid desaturation without causing unwanted outcomes: incomplete desaturation and accumulations of nitrite and nitrous oxide with too little acetate or residual acetate and accelerated H₂S generation with excess acetate. The model can be used to optimize treatment recipes for maximizing desaturation or precipitation in most subsurface groundwater environments for liquefaction mitigation. However, field data describing the environmental biogeochemical characteristics (e.g., pH, background chemical concentrations) for the most optimized results are necessary to understand the potential biogeochemical reactions and processes that may impact MIDP and, subsequently, liquefaction mitigation.