A Bio-Chemo-Hydro-Mechanical Model for the Simulation of Biocementation in Soils: One-Dimensional Finite Element Simulations

Abstract

Microbially induced calcite precipitation is a soil improvement technique in which bacteria are used to produce calcium carbonate (biocement), precipitated after the hydrolysis of urea by the urease enzyme present in the microorganisms. This technique is becoming popular, and there have been several real cases of its use; however, the dosages and reaction times used to attain a required percentage of biocement mainly stem from previous experimental tests, and calculations are not performed. Thus, it is fundamental to have more robust tools and the existence of numerical models able to compute the amount precipitated, such as the one proposed in this paper, can be an important contribution. A two-phase porous medium model is created to analyse the precipitation process. The solid phase contains soil particles, bacteria and biocement, while the fluid phase contains water, urea and other dissolved species. A coupled bio-chemo-hydro-mechanical finite element formulation is defined, embodying the biochemical reaction, water seepage, the diffusion of species and soil deformation. The main novelties of this study are as follows: (i) porosity changes are computed considering the generation of solid mass due to biocement precipitation, and, therefore, soil permeability is updated during the calculation, with these highly coupled equations being integrated in time simultaneously and not sequentially; and (ii) the model is calibrated with experimental tests conceived especially for this purpose. The model is then used to compute the biocement precipitated in a sand column simulating a real experimental test. The results of the simulations present a distribution of biocement along the column closer to that observed in the experimental tests, validating the model.

1. Introduction

Biocementation, also referred as microbially induced calcite precipitation (MICP), is a ground improvement technique that involves the use of microorganisms (particularly bacteria) to generate calcium carbonate (CaCO₃) predominantly in the form of calcite. The used microorganisms are usually non-pathogenic and can be found in different types of soils.

Figure 1 presents a detailed representation of the mechanism of MICP. Basically, a feeding solution with urea and calcium chloride is added to a medium of urease-producing bacteria. The urease enzymes hydrolyse urea, producing ammonium (${\mathrm{N}\mathrm{H}}_4^{+}$) and bicarbonate (${\mathrm{C}\mathrm{O}}_3^{2-})$ ions, according to the following chemical reaction (1):

$${\mathrm{C}\mathrm{O}\left({\mathrm{N}\mathrm{H}}_2\right)}_{(\mathrm{a}\mathrm{q}.)}+{2{\mathrm{H}}_2\mathrm{O}}_{\left(\mathrm{l}\right)}\to {{2\mathrm{N}\mathrm{H}}_4^{+}}_{(\mathrm{a}\mathrm{q}.)}+{{\mathrm{C}\mathrm{O}}_3^{2-}}_{(\mathrm{a}\mathrm{q}.)}.$$

(1)

Calcium chloride is usually supplied in the feeding solution as a calcium source, which dissolves into calcium (Ca²⁺) and chloride ions. The bacteria cells serve as a nucleation site for calcium ions which react with the bicarbonate ions, resulting from the urea hydrolysis, precipitating solid calcium carbonate (CaCO₃), according to the following chemical reaction (2):

$${{\mathrm{C}\mathrm{a}}^{2+}}_{(\mathrm{a}\mathrm{q}.)}+{{\mathrm{C}\mathrm{O}}_3^{2-}}_{(\mathrm{a}\mathrm{q}.)}\to {\mathrm{C}\mathrm{a}{\mathrm{C}\mathrm{O}}_3}_{\left(\mathrm{s}\right)},$$

(2)

This calcium carbonate, usually referred to as calcite, may lead to bacteria encapsulation when it accumulates, forming an envelope and trapping bacteria.

The precipitated CaCO₃ clogs the pores and bonds soil grains, thus decreasing soil permeability and increasing its stiffness and strength [1,2,3], sealing rock fractures [4] and increasing its resistance to erosive processes [5,6,7,8]. However, in the case of large-scale field applications, one of the main limitations of this technique is to achieve a homogenous biocement distribution in the treated volume [1,9,10].

Treatment dosages, the number of injections and the reaction time to attain a required percentage of calcium carbonate are determined through experimental tests performed in the laboratory before the field work is carried out. However, it is fundamental to have more robust tools. Numerical models able to compute the amount precipitated, such as the one proposed in this paper, would be an important contribution. The majority of the works concerning the numerical modelling of biocementation in soils are reactive transport models (e.g., [11,12,13,14,15,16]), which only focus on the fluid transport along the porous matrix by modelling water flow, chemical species transport and reactions related to MICP processes. There are also a few hydro-mechanical models (e.g., [17,18]) that consider soil as a deformable medium.

This research presents a chemo-hydro-mechanical coupled finite element model to predict precipitated calcite content in soils during biocementation and subsequent changes in the hydro-mechanical properties of this porous medium. The main purpose of this model is to be able to predict the amount of precipitated biocement for a given treatment protocol, i.e., by knowing the dosages of urea and calcium and the duration of the treatment. This new bio-chemo-hydro-mechanical model has some differences to the above-mentioned ones, listed as follows:

(i)

The model considers soil grains and attached bacteria as a single species in the solid phase for the sake of simplicity. In this way, it is assumed that the attached bacteria are already homogeneously distributed in the solid matrix, and no bacterial injections are simulated;

(ii)

The species present are solid grains with attached bacteria, solid calcite, water, urea, ammonium, calcium and chloride ions;

(iii)

All the couplings are integrated in the formulation. Moreover, for the usual hydro-mechanical coupled behaviour in saturated soils, in which the deformation caused by mechanical actions affects the water flow (consolidation) or changes in water pore pressure affect the effective stresses, biochemical coupling is introduced. The chemical reactions result in the precipitation of solid calcium carbonate (referred to as calcite from now on) in the pores, reducing the soil porosity. At every time step, permeability changes are calculated as a function of the new porosity values (the pore-clogging effect), which may affect the flow velocity. The elastic stiffness of the medium is considered to be constant at this stage but could be updated to consider the presence of calcite forming bonds between the soil particles.

The most important novelty of the presented model is that the biochemical reaction rate parameters are calibrated experimentally with small-scale sand column tests.

Moreover, and contrary to other numerical techniques proposed in the literature, the highly coupled equations of the model, which considers porosity and permeability changes due to biocement precipitation, are simultaneously and not sequentially integrated in time.

2. Model Formulation

The formulation developed encompasses the hydro-geologic point of view (e.g., [11,12,14]), which only addresses diffusion couplings that involve water seepage and water advective diffusion/dispersion of ionic species, and the geo-mechanical point of view (e.g., [17,18]), which addresses hydro-mechanical couplings and accounts for soil deformation.

A two-phase saturated porous medium is considered, in which each phase is composed of several species, as presented by the schematic pore-scale representation in Figure 2. The solid phase (S) contains soil particles with attached bacteria, homogeneously distributed in the whole medium and jointly denoted by the symbol (s) as well as calcite, denoted by the symbol, (c). It is assumed that calcite is the only calcium carbonate morphotype that is considered as a resulting product of the ureolysis-induced calcium carbonate precipitation process. The fluid phase (W) contains water (w), urea (u), calcium (Ca), chloride (Cl) and ammonium (a) ions.

The main hypotheses of this model follow the strongly interacting model by Bataille and Kestin [19]:

• Mass balance is required for each species;
• Momentum balance is required for the porous medium as a whole;
• Species in the fluid phase are endowed with their own velocities; the velocity of the species in the solid phase is equal to the velocity of the phase (${\mathbf{v}}_{k\mathrm{S}}={\mathbf{v}}_{\mathrm{S}}, k\in \mathrm{S}$);
• The urea hydrolysis reaction, as shown in Equation (1), produces bicarbonate and ammonium ions;
• Calcium and bicarbonate ions in the fluid phase precipitate, forming calcite in the solid phase, according to Equation (2);
• As a simplification, the reactions in iv. and v. have the same rate. Therefore, the mass of bicarbonate formed in the hydrolysis of urea immediately precipitates when the calcite is formed, and bicarbonate ions are not a part of the model.

During the formulation, several mass and volume quantities (defined in

Table 1. Definition of several mass and volume quantities. The suffix $kK$ is used in this research as a designation for the species $k$ in phase $K$.

Quantity (Unit)	Definition
Current volume of solid phase—$V_{\mathrm{S}}$ (m3)	$V_{\mathrm{S}}={\displaystyle \sum _{kϵ\mathrm{S}}}{V}_{k\mathrm{S}}$
Current volume of fluid phase—$V_{\mathrm{W}}$ (m3)	$V_{\mathrm{W}}={\displaystyle \sum _{kϵ\mathrm{W}}}{V}_{k\mathrm{W}}$
Current total volume—$V$ (m3)	$V=V_{\mathrm{S}}+V_{\mathrm{W}}$
Volume fraction of the species—$n_{kK}$	$n_{kK}=V_{kK}/V$
Volume content—$v_{kK}$	$v_{kK}=V_{kK}/V_{\mathbf{0}}$
Mass content—$m_{kK}$ (kg m−3)	$m_{kK}=M_{kK}/V_{\mathbf{0}}$
Molar concentration—$c_{k\mathrm{W}}$ (mol m−3)	$c_{k\mathrm{W}}=N_{k\mathrm{W}}/V_{\mathrm{W}}$
Intrinsic density of the species k—${\rho }_k$ (kg m−3)	${\rho }_k=m_k^{\left(M\right)}/v_k^{\left(M\right)}=M_{kK}/V_{kK}$
Apparent density of the species—${\rho }^{kK}$ (kg m−3)	${\rho }^{kK}=M_{kK}/V=n_{kK}{\rho }_k$
Volume fraction (porosity) of the fluid phase—$n_{\mathrm{W}}$	$n_{\mathrm{W}}={\displaystyle \sum _{k\in \mathrm{W}}}{n}_{k\mathrm{W}}$
Volume content of the fluid phase—$v_{\mathrm{W}}$	$v_{\mathrm{W}}={\displaystyle \sum _{k\in \mathrm{W}}}{v}_{k\mathrm{W}}$
Volume content of the solid phase—$v_{\mathrm{S}}$	$v_{\mathrm{S}}={\displaystyle \sum _{k\in \mathrm{S}}}{v}_{k\mathrm{S}}.$

) are used. The current volume, mass and number of moles of the species k in phase K are denoted by $V_{kK}$, $M_{kK}$ and $N_{kK}$, respectively. The molar mass, the molar volume and the intrinsic density of the species k are denoted by $m_k^{\left(M\right)}$, $v_k^{\left(M\right)}$ and ${\rho }_k$, respectively. The initial volume of the porous medium is denoted by $V_{\mathbf{0}}$, and the current solid, fluid and total volumes are denoted by $V_{\mathrm{S}}$, $V_{\mathrm{W}}$ and $V$, respectively. Several other quantities may be defined for species k in phase K, like the volume fraction $\left(n_{kK}\right)$, the volume content $\left(v_{kK}\right)$, the mass content $\left(m_{kK}\right)$, the molar concentration $\left(c_{k\mathrm{W}}\right)$ and the apparent density $\left({\rho }^{kK}\right)$. Summing the contributions of the species, it is also possible to define the volume fraction (or porosity) of the fluid phase $\left(n_{\mathrm{W}}\right)$, the volume content of the fluid phase $\left(v_{\mathrm{W}}\right)$ and the volume content of the solid phase $\left(v_{\mathrm{S}}\right)$.

2.1. Mass Balance of Species k in Phase K

The species mass balance can be written, in the current configuration, as presented by Loret and Simões in [20,21]:

$${\displaystyle \frac{{\mathrm{d}}^k{\rho }^{kK}}{\mathrm{d}t}}+{\rho }^{kK}\mathrm{d}\mathrm{i}\mathrm{v} {\mathbf{v}}_{kK}={\widehat{\rho }}^{kK}$$

(3)

where ${\mathrm{d}}^k( )/\mathrm{d}t$ is the total time derivative as seen by an observer moving with the particles of species k, ${\mathbf{v}}_{kK}$ is the velocity of species k in phase K and ${\widehat{\rho }}^{kK}$ is the rate of the mass of the species per unit of volume that is created or lost during the chemical reactions. Upon linearization, the mass balance equation in the reference configuration is as follows (see Equation (A4) in Appendix A):

$${\displaystyle \frac{\mathrm{d}{v}_{kK}}{\mathrm{d}t}}+\mathrm{d}\mathrm{i}\mathrm{v} {\mathbf{J}}_{kK}={\displaystyle \frac{{\widehat{m}}^{kK}}{{\rho }_k}}$$

(4)

where ${\mathbf{J}}_{kK}=n_{kK}\left({\mathbf{v}}_{kK}-{\mathbf{v}}_{\mathrm{S}}\right)$ is the flux of the species through the solid skeleton and ${\widehat{m}}^{kK}$ is the rate of the mass of the species per unit of volume that is created or lost during the chemical reactions, in the reference configuration.

By adding the contribution of all species $k\in \mathrm{W}$ and all species $k\in \mathrm{S}$, we obtain the following global equation of mass balance (see Equation (A10) in Appendix A):

$$\mathrm{d}\mathrm{i}\mathrm{v} {\mathbf{v}}_{\mathrm{S}}+\mathrm{d}\mathrm{i}\mathrm{v} {\mathbf{J}}_{\mathrm{W}}={\sum }_{k\in \mathrm{W}}{\displaystyle \frac{{\widehat{m}}^{k\mathrm{W}}}{{\rho }_k}}+{\sum }_{k\in \mathrm{S}}{\displaystyle \frac{{\widehat{m}}^{k\mathrm{S}}}{{\rho }_k}}$$

(5)

where ${\mathbf{J}}_{\mathrm{W}}={\sum }_{k\in \mathrm{W}}{\mathbf{J}}_{k\mathrm{W}}$ represents the flux of water.

In the case of a species $k\in \mathrm{W}$, the diffusive flux is defined as follows:

$${\mathbf{J}}_{k\mathrm{W}}^{\mathrm{d}}=n_{k\mathrm{W}}\left({\mathbf{v}}_{k\mathrm{W}}-{\mathbf{v}}_{\mathrm{w}\mathrm{W}}\right)={\mathbf{J}}_{k\mathrm{W}}-{\displaystyle \frac{n_{k\mathrm{W}}}{n_{\mathrm{w}\mathrm{W}}}}{\mathbf{J}}_{\mathrm{w}\mathrm{W}}$$

(6)

Therefore, the mass balance in Equation (4) for a species in the fluid phase can be rewritten as follows:

$$n_{\mathrm{w}}{v}_k^{\left(M\right)}{\displaystyle \frac{\mathrm{d}{c}_{k\mathrm{W}}}{\mathrm{d}t}}-\mathrm{d}\mathrm{i}\mathrm{v} {\mathsf{\alpha }}_k+v_k^{\left(M\right)}{\mathbf{J}}_{\mathrm{w}}·\mathsf{\nabla }{c}_{k\mathrm{W}}+c_{k\mathrm{W}}{v}_k^{\left(M\right)}\left({\sum }_{l\in \mathrm{W}}{\displaystyle \frac{{\widehat{m}}^{l\mathrm{W}}}{{\rho }_k}}\right)-{\displaystyle \frac{{\widehat{m}}^{k\mathrm{W}}}{{\rho }_k}}=\mathbf{0},$$

(7)

where ${\mathsf{\alpha }}_k=-{\sum }_{l\in \mathrm{W}}\left({\mathrm{I}}_{kl}-c_{k\mathrm{W}}{v}_k^{\left(M\right)}\right){\mathbf{J}}_{l\mathrm{w}}^d,$ and $\mathbf{I}$ is the identity matrix.

On the other hand, for calcite, Equation (4) simply results in the following:

$${\displaystyle \frac{\mathrm{d}{m}_{c\mathrm{S}}}{\mathrm{d}t}}={\widehat{m}}^{c\mathrm{S}}.$$

(8)

In the reaction of urea hydrolysis, one mole of urea is consumed when two moles of ammonium are formed at the rate r₁. Therefore, using the notation by van Wijngaarden et al. [11], the rate of the mass change of the species involved in this reaction is defined as follows:

$${\displaystyle \frac{{\widehat{m}}^{u\mathrm{W}}}{{\rho }_u}}={{-v}_u^{\left(M\right)}n}_{\mathrm{w}}{r}_1$$

(9)

$${\displaystyle \frac{{\widehat{m}}^{a\mathrm{W}}}{{\rho }_a}}={{2v}_a^{\left(M\right)}n}_{\mathrm{w}}{r}_1$$

(10)

In the reaction of calcite precipitation, one mole of Ca²⁺ is consumed when one mole of calcite is formed at the rate r₂. Therefore, the rate of the mass change of the species involved in this reaction is as follows:

$${\displaystyle \frac{{\widehat{m}}^{c\mathrm{S}}}{{\rho }_c}}={v_c^{\left(M\right)}n}_{\mathrm{w}}{r}_2$$

(11)

$${\displaystyle \frac{{\widehat{m}}^{Ca\mathrm{W}}}{{\rho }_{Ca}}}={{-v}_{Ca}^{\left(M\right)}n}_{\mathrm{w}}{r}_2$$

(12)

The efficiency of the reactions of urea hydrolysis and calcium carbonate precipitation depends on several factors, such as temperature and pH [11], bacterial activity and the mineralogy of the clay particles in the soil [22]. Encapsulation by calcium carbonate crystals can also occur, creating a diffusion barrier around the bacteria. In addition, aerobic bacteria injected into an anaerobic soil may eventually die due to the lack of oxygenation [11]. All these factors can reduce the efficiency of the reaction.

Regarding the reaction of urea hydrolysis, the Michaelis–Menten kinetics mathematical model is usually used to define the reaction rate. Although some authors (e.g., [12]) defined different velocities for the reactions of urea hydrolysis and calcite precipitation, in this work, the two reactions have the same rate (as in [11]):

$$r_1=r_2=r=v_{max}{\displaystyle \frac{c_{u\mathrm{W}}}{K_m+c_{u\mathrm{W}}}}{e}^{\left(-{\displaystyle \frac{t}{t_{max}}}\right)}$$

(13)

where $v_{m\mathrm{a}x},$ $K_m$ and $t_{max}$ are constants to be obtained. Since the bacteria distribution is assumed to be homogeneous, parameter $v_{m\mathrm{a}x}$ is constant, and, therefore, Equation (13) presents an exponential time decay for the reaction rate, where $t_{max}$ is related to the lifespan or the activity of the bacteria.

Due to the simplification adopted in Equation (13), even though bicarbonate ions are involved in the reactions, their mass balance is not a part of the model since it is considered that the mass of the bicarbonate that is formed immediately precipitates when the calcite is formed.

2.2. Momentum Balance of the Species

The global equation for the conservation of linear momentum is as follows:

$$\mathrm{d}\mathrm{i}\mathrm{v} \mathsf{\sigma }+\rho \mathbf{g}=\mathbf{0}$$

(14)

where $\mathsf{\sigma }$ is the total stress tensor, which for soils, according to Terzaghi’s principle, is defined as follows:

$$\mathsf{\sigma }={\mathsf{\sigma }}^{\prime }-p_{\mathrm{W}}\mathbf{I}$$

(15)

where ${\mathsf{\sigma }}^{\prime }$ is the effective stress tensor, $p_{\mathrm{W}}$ is the pore–water pressure (the negative sign is introduced as it is a general convention to consider the tensile components of stress as positive [23]), $\mathbf{I}$ is the identity matrix, $\mathbf{g}$ is the gravity acceleration (9.8 m/s²) and $\rho ={\sum }_{k\in \mathrm{W},\mathrm{S}}{\rho }^{kK}$.

Assuming a linear elastic behaviour, the effective stress tensor is related to the deformation tensor $\mathsf{\epsilon }$ by Hooke’s law as follows:

$${\mathsf{\sigma }}^{\prime }=\mathbf{E}:\mathsf{\epsilon }$$

(16)

where $\mathbf{E}$ is the fourth-rank elasticity tensor, and the colon denotes the double dot product. In the 1-D case, only one material constant, the Young’s modulus E, is needed.

2.3. Diffusion Equation for the Ionic Species (Fick’s Law)

The diffusive/dispersive flux of the ionic species is related to the concentration gradients according to Fick’s law as follows:

$${\mathbf{J}}_{k\mathrm{W}}^{\mathrm{d}}=-{v_k^{\left(M\right)}n}_{\mathrm{w}}\mathbf{D}·\mathsf{\nabla }{c}_{k\mathrm{W}}$$

(17)

where D is the matrix of the diffusion constants. In the 1-D case, the diffusive/dispersive flux has the following form [11]:

$${\mathbf{J}}_{k\mathrm{W}}^{\mathrm{d}}=-{v_k^{\left(M\right)}n}_{\mathrm{w}}{\mathsf{\alpha }}_L‖{\mathit{v}}_{\mathrm{w}\mathrm{W}}‖\mathsf{\nabla }{c}_{k\mathrm{W}}$$

(18)

where ${\mathsf{\alpha }}_L$ is the longitudinal diffusion coefficient, which is assumed to be constant for all the ionic species.

2.4. Flow Equation for the Fluid Phase (Darcy’s Law)

It is assumed that the water flow in the porous medium is governed by Darcy’s law as follows:

$${\mathbf{J}}_{\mathrm{w}}=-\mathit{k}·\left(\mathsf{\nabla }{p}_{\mathrm{W}}-{\displaystyle \frac{{\rho }^{\mathrm{W}}\mathit{g}}{n_{\mathrm{w}}}}\right),$$

(19)

where k is the water permeability tensor (an isotropic medium is assumed) and ${\rho }^{\mathrm{W}}={\sum }_{kϵ\mathrm{W}}{\rho }^{k\mathrm{W}}$. Particularly for the 1-D case, the following can be written:

$$J_{\mathrm{w}}=-k_{\mathit{x}}\left({\displaystyle \frac{{\mathrm{d}p}_{\mathrm{W}}}{\mathrm{d}x}}+{\displaystyle \frac{{\rho }^{\mathrm{W}}g}{n_{\mathrm{w}}}}\right),$$

(20)

where $k_x=k_h/{\rho }_{w}g$ is the longitudinal permeability obtained in ${\mathrm{m}}^3\mathrm{s}{\mathrm{k}\mathrm{g}}^{-1}$, and $k_h$ is the hydraulic conductivity, related to the porosity $n_{\mathrm{w}}$ (see Equation (A14) in Appendix B for its update) according to the Kozeny–Carman equation as follows:

$$k_h={\overline{k}}_h{\displaystyle \frac{{\left(n_{\mathrm{W}}\right)}^3}{{\left(1-n_{\mathrm{W}}\right)}^2}}.$$

(21)

3. Finite Element Formulation

3.1. The Semi-Discrete Field Equations

The 1-D finite element formulation of the problem in its weak form is obtained by multiplying the equations of the global momentum balance (14), global mass balance (5), aqueous species mass balance (7) and calcium carbonate mass balance (8) by the fields of the virtual displacement (${\mathit{u}}^{*}$), water pressure ($p^{*}$), concentration ($c^{*}$) and mass per unit of volume ($m^{*}$), respectively. Their integration over the element volume $\Omega$ using the Gauss theorem results in the following:

$${\int }_{\Omega }\mathsf{\nabla }{\mathit{u}}^{\mathit{*}}:\mathsf{\sigma }\mathrm{d}\mathrm{V}-{\int }_{\Omega }\mathsf{\rho }{\mathit{u}}^{\mathit{*}}·\mathbf{g}\mathrm{d}\mathrm{V}={\int }_{\partial \Omega }{\mathit{u}}^{\mathit{*}}·\mathsf{\sigma }·\mathbf{n}\mathrm{d}\mathrm{S}$$

(22)

$$\begin{gathered} {\int }_{\Omega }{p}^{\mathit{*}}\mathrm{d}\mathrm{i}\mathrm{v} {\mathbf{v}}_{\mathrm{s}}\mathrm{d}\mathrm{V}-{\int }_{\Omega }\mathsf{\nabla }{p}^{\mathit{*}}{·\mathbf{J}}_{\mathrm{w}}\mathrm{d}\mathrm{V}-{\int }_{\Omega }{p}^{*}\left({\sum }_{k\in \mathrm{W}}{\displaystyle \frac{{\widehat{m}}^{k\mathrm{W}}}{{\rho }_k}}+{\sum }_{k\in \mathrm{S}}{\displaystyle \frac{{\widehat{m}}^{k\mathrm{S}}}{{\rho }_k}}\right)\mathrm{d}\mathrm{V}= \\ -{\int }_{\partial \Omega }{p}^{\mathit{*}}{\mathbf{J}}_{\mathrm{w}}·\mathbf{n}\mathrm{d}\mathrm{S} \end{gathered}$$

(23)

$$\begin{gathered} {\int }_{\Omega }{c}^{\mathit{*}}{n}_{\mathrm{w}}{\displaystyle \frac{\mathrm{d}{c}_{k\mathrm{W}}}{\mathrm{d}t}}{v}_k^{\left(M\right)}dV+{\int }_{\Omega }{c}^{\mathit{*}}\left[\left({\sum }_{l\in \mathrm{W}}{\displaystyle \frac{{\widehat{m}}^{l\mathrm{W}}}{{\rho }_l}}\right)c_{k\mathrm{W}}{v}_k^{\left(M\right)}-{\displaystyle \frac{{\widehat{m}}^{k\mathrm{W}}}{{\rho }_k}}\right]\mathrm{d}\mathrm{V}+ \\ {\int }_{\Omega }{\mathsf{\nabla }c}^{\mathit{*}}·{\mathsf{\alpha }}_k\mathrm{d}\mathrm{V}+{\int }_{\Omega }{c}^{*}{\mathbf{J}}_{\mathrm{w}}·{\mathsf{\nabla }c}_{k\mathrm{W}}{v}_k^{\left(M\right)}\mathrm{d}\mathrm{V}={\int }_{\partial \Omega }{c}^{\mathit{*}}{\mathsf{\alpha }}_k·\mathbf{n}\mathrm{d}\mathrm{S}, k=\mathrm{u},\mathrm{a},\mathrm{C}\mathrm{a},\mathrm{C}\mathrm{l} \end{gathered}$$

(24)

$${\int }_{\Omega }{m}^{\mathit{*}}\left({\displaystyle \frac{\mathrm{d}{m}_{c\mathrm{S}}}{\mathrm{d}t}}-{\widehat{m}}^{c\mathrm{S}}\right)\mathrm{d}\mathrm{V}=\mathbf{0}$$

(25)

where $\partial \Omega$ is the element surface, and $\mathbf{n}$ is its unit normal vector.

The seven primary unknown fields (the solid displacements, water pressure, concentrations of the ionic fluid species—namely urea, ammonium, calcium and chloride—and the calcite mass per unit of volume) are approximated inside each bar finite element (e) by using a linear combination of shape functions (${\mathbf{N}}_{\mathrm{U}}$, ${\mathbf{N}}_{{\mathrm{P}}_{\mathrm{W}}}$, ${\mathbf{N}}_{\mathrm{C}}$ and ${\mathbf{N}}_{{\mathrm{M}}_c}$) and nodal values (${\mathbf{U}}^{\left(e\right)}$, ${{\mathbf{P}}_{\mathrm{w}}}^{\left(e\right)}$, ${{\mathbf{C}}_{k\mathrm{W}}}^{\left(e\right)}$ and ${{\mathbf{M}}_c}^{\left(e\right)}$):

$${\mathbf{u}}^{\left(e\right)}={\mathbf{N}}_{\mathrm{U}}·{\mathbf{U}}^{\left(e\right)}$$

(26)

$${p_{\mathrm{W}}}^{\left(e\right)}={\mathbf{N}}_{{\mathrm{P}}_{\mathrm{W}}}·{{\mathbf{P}}_{\mathrm{w}}}^{\left(e\right)}$$

(27)

$${c_{k\mathrm{W}}}^{\left(e\right)}={\mathbf{N}}_{\mathrm{C}}·{{\mathbf{C}}_{k\mathrm{W}}}^{\left(e\right)}, k=\mathrm{u},\mathrm{a},\mathrm{C}\mathrm{a},\mathrm{C}\mathrm{l}$$

(28)

$${m_{\mathrm{c}\mathrm{S}}}^{\left(e\right)}={\mathbf{N}}_{{\mathrm{M}}_c}·{{\mathbf{M}}_c}^{\left(e\right)}$$

(29)

The same approximations are made for the virtual fields. Quadratic interpolation is employed for solid displacements (three nodes per element), and linear interpolation is employed for the other primary unknowns (two nodes per element).

By substituting the virtual fields’ approximation into (22)–(25), we obtain the following element first-order semi-discrete equation:

$${\mathbf{F}}_{\mathrm{i}\mathrm{n}\mathrm{t}}^{\left(\mathrm{e}\right)}={\mathbf{F}}_{\mathrm{e}\mathrm{x}\mathrm{t}}^{\left(\mathrm{e}\right)},$$

(30)

where the element external forces vector is as follows:

$${\mathbf{F}}_{\mathrm{e}\mathrm{x}\mathrm{t}}^{\left(\mathrm{e}\right)}=\left\{\begin{array}{c}\begin{array}{c}{\int }_{\partial \Omega }{\mathbf{N}}_{\mathrm{U}}^{\mathrm{T}}\mathsf{\sigma }·\mathbf{n}\mathrm{d}\mathrm{S}\\ -{\int }_{\partial \Omega }{\mathbf{N}}_{{\mathrm{P}}_{\mathrm{W}}}^{\mathrm{T}}{\mathbf{J}}_{\mathrm{w}}·\mathbf{n}\mathrm{d}\mathrm{S}\end{array}\\ \begin{array}{c}\begin{array}{c}{\int }_{\partial \Omega }{\mathbf{N}}_{\mathrm{C}}^{\mathrm{T}}{\mathsf{\alpha }}_{\mathrm{u}}·\mathbf{n}\mathrm{d}\mathrm{S}\\ \mathbf{0}\end{array}\\ {\int }_{\partial \Omega }{\mathbf{N}}_{\mathrm{C}}^{\mathrm{T}}{\mathsf{\alpha }}_{\mathrm{C}\mathrm{a}}·\mathbf{n}\mathrm{d}\mathrm{S}\\ {\int }_{\partial \Omega }{\mathbf{N}}_{\mathrm{C}}^{\mathrm{T}}{\mathsf{\alpha }}_{\mathrm{C}\mathrm{l}}·\mathbf{n}\mathrm{d}\mathrm{S}\end{array}\\ \mathbf{0}\end{array}\right\}$$

(31)

and, since ${\sum }_{k\in \mathrm{W}}{\displaystyle \frac{{\widehat{m}}^{k\mathrm{W}}}{{\rho }_k}}=\left({-v}_{\mathrm{u}}^{\left(M\right)}+2v_{\mathrm{a}}^{\left(M\right)}-v_{\mathrm{C}\mathrm{a}}^{\left(M\right)}\right){\displaystyle \frac{1}{m_{\mathrm{c}}^{\left(M\right)}}}{\displaystyle \frac{\mathrm{d}{m}_{c\mathrm{S}}}{\mathrm{d}t}}$ and ${\sum }_{k\in \mathrm{W}}{\displaystyle \frac{{\widehat{m}}^{k\mathrm{W}}}{{\rho }_k}}+{\sum }_{k\in \mathrm{S}}{\displaystyle \frac{{\widehat{m}}^{k\mathrm{S}}}{{\rho }_k}}=\left({-v}_{\mathrm{u}}^{\left(M\right)}+2v_{\mathrm{a}}^{\left(M\right)}-v_{\mathrm{C}\mathrm{a}}^{\left(M\right)}+v_{\mathrm{c}}^{\left(M\right)}\right){\displaystyle \frac{1}{m_{\mathrm{c}}^{\left(M\right)}}}{\displaystyle \frac{\mathrm{d}{m}_{c\mathrm{S}}}{\mathrm{d}t}}$, the element internal forces vector is as follows:

$${\mathbf{F}}_{\mathrm{i}\mathrm{n}\mathrm{t}}^{\left(\mathrm{e}\right)}=\left\{\begin{array}{c}\begin{array}{c}\begin{array}{c}{\int }_{\Omega }{\mathbf{B}}_{\mathrm{U}}^{\mathrm{T}}\mathsf{\sigma }\mathrm{d}\mathrm{V}-{\int }_{\Omega }{\mathbf{N}}_{\mathrm{U}}^{\mathrm{T}}\rho \mathbf{g}\mathrm{d}\mathrm{V}\\ \_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\end{array}\\ \begin{array}{c}{\int }_{\Omega }{\mathbf{N}}_{{\mathrm{P}}_{\mathrm{W}}}^{\mathrm{T}}\mathrm{d}\mathrm{i}\mathrm{v} {\mathit{v}}_{\mathrm{s}}\mathrm{d}\mathrm{V}-{\int }_{\Omega }{\mathsf{\nabla }\mathbf{N}}_{{\mathrm{P}}_{\mathrm{W}}}^{\mathrm{T}}{\mathbf{J}}_{\mathrm{w}}\mathrm{d}\mathrm{V}\\ \begin{array}{c}-\left({-v}_{\mathrm{u}}^{\left(M\right)}+2v_{\mathrm{a}}^{\left(M\right)}-v_{\mathrm{C}\mathrm{a}}^{\left(M\right)}+v_{\mathrm{c}}^{\left(M\right)}\right){\displaystyle \frac{1}{m_{\mathrm{c}}^{\left(M\right)}}}{\int }_{\Omega }{\mathbf{N}}_{{\mathrm{P}}_{\mathrm{W}}}^{\mathrm{T}}{\displaystyle \frac{\mathrm{d}{m}_{c\mathrm{S}}}{\mathrm{d}t}}\mathrm{d}\mathrm{V}\\ \_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\end{array}\\ {\int }_{\Omega }{\mathbf{N}}_{\mathrm{C}}^{\mathrm{T}}{n}_{\mathrm{w}}{\displaystyle \frac{\mathrm{d}{c}_{\mathrm{u}\mathrm{W}}}{\mathrm{d}t}}{v}_{\mathrm{u}}^{\left(M\right)}\mathrm{d}\mathrm{V}+{\int }_{\Omega }{\mathsf{\nabla }\mathbf{N}}_{\mathrm{C}}^{\mathrm{T}}{\mathsf{\alpha }}_{\mathrm{u}}\mathrm{d}\mathrm{V}+{\int }_{\Omega }{\mathbf{N}}_{\mathrm{C}}^{\mathrm{T}}{\mathbf{J}}_{\mathrm{w}}·{\mathsf{\nabla }c}_{\mathrm{u}\mathrm{W}}{v}_{\mathrm{u}}^{\left(M\right)}\mathrm{d}\mathrm{V}\end{array}\end{array}\\ \begin{array}{c}\begin{array}{c}\begin{array}{c}+{\int }_{\Omega }{\mathbf{N}}_{\mathrm{C}}^{\mathrm{T}}{v}_{\mathrm{u}}^{\left(M\right)}\left[c_{\mathrm{u}\mathrm{W}}\left({-v}_{\mathrm{u}}^{\left(M\right)}+2v_{\mathrm{a}}^{\left(M\right)}-v_{\mathrm{C}\mathrm{a}}^{\left(M\right)}\right)+1\right]{\displaystyle \frac{1}{m_{\mathrm{c}}^{\left(M\right)}}}{\displaystyle \frac{\mathrm{d}{m}_{c\mathrm{S}}}{\mathrm{d}t}}\mathrm{d}\mathrm{V}\\ \_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\end{array}\\ \begin{array}{c}{\int }_{\Omega }{\mathbf{N}}_{\mathrm{C}}^{\mathrm{T}}{n}_{\mathrm{w}}{\displaystyle \frac{\mathrm{d}{c}_{\mathrm{a}\mathrm{W}}}{\mathrm{d}t}}{v}_{\mathrm{a}}^{\left(M\right)}\mathrm{d}\mathrm{V}+{\int }_{\Omega }{\mathsf{\nabla }\mathbf{N}}_{\mathrm{C}}^{\mathrm{T}}{\mathsf{\alpha }}_{\mathrm{a}}\mathrm{d}\mathrm{V}+{\int }_{\Omega }{\mathbf{N}}_{\mathrm{C}}^{\mathrm{T}}{\mathbf{J}}_{\mathrm{w}}·{\mathsf{\nabla }c}_{\mathrm{a}\mathrm{W}}{v}_{\mathrm{a}}^{\left(M\right)}\mathrm{d}\mathrm{V}\\ \begin{array}{c}+{\int }_{\Omega }{\mathbf{N}}_{\mathrm{C}}^{\mathrm{T}}{v}_{\mathrm{a}}^{\left(M\right)}\left[c_{\mathrm{a}\mathrm{W}}\left({-v}_{\mathrm{u}}^{\left(M\right)}+2v_{\mathrm{a}}^{\left(M\right)}-v_{\mathrm{C}\mathrm{a}}^{\left(M\right)}\right)-2\right]{\displaystyle \frac{1}{m_{\mathrm{c}}^{\left(M\right)}}}{\displaystyle \frac{\mathrm{d}{m}_{c\mathrm{S}}}{\mathrm{d}t}}\mathrm{d}\mathrm{V}\\ \_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\end{array}\\ {\int }_{\Omega }{\mathbf{N}}_{\mathrm{C}}^{\mathrm{T}}{n}_{\mathrm{w}}{\displaystyle \frac{\mathrm{d}{c}_{\mathrm{C}\mathrm{a}\mathrm{W}}}{\mathrm{d}t}}{v}_{\mathrm{C}\mathrm{a}}^{\left(M\right)}\mathrm{d}\mathrm{V}+{\int }_{\Omega }{\mathsf{\nabla }\mathbf{N}}_{\mathrm{C}}^{\mathrm{T}}{\mathsf{\alpha }}_{\mathrm{C}\mathrm{a}}\mathrm{d}\mathrm{V}+{\int }_{\Omega }{\mathbf{N}}_{\mathrm{C}}^{\mathrm{T}}{\mathbf{J}}_{\mathrm{w}}·{\mathsf{\nabla }c}_{\mathrm{C}\mathrm{a}\mathrm{W}}{v}_{\mathrm{C}\mathrm{a}}^{\left(M\right)}\mathrm{d}\mathrm{V}\end{array}\end{array}\\ \begin{array}{c}+{\int }_{\Omega }{\mathbf{N}}_{\mathrm{C}}^{\mathrm{T}}{v}_{\mathrm{C}\mathrm{a}}^{\left(M\right)}\left[c_{\mathrm{C}\mathrm{a}\mathrm{W}}\left({-v}_{\mathrm{u}}^{\left(M\right)}+2v_{\mathrm{a}}^{\left(M\right)}-v_{\mathrm{C}\mathrm{a}}^{\left(M\right)}\right)+1\right]{\displaystyle \frac{1}{m_{\mathrm{c}}^{\left(M\right)}}}{\displaystyle \frac{\mathrm{d}{m}_{c\mathrm{S}}}{\mathrm{d}t}}\mathrm{d}\mathrm{V}\\ \_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\end{array}\\ \begin{array}{c}{\int }_{\Omega }{\mathbf{N}}_{\mathrm{C}}^{\mathrm{T}}{n}_{\mathrm{w}}{\displaystyle \frac{\mathrm{d}{c}_{\mathrm{C}\mathrm{l}\mathrm{W}}}{\mathrm{d}t}}{v}_{\mathrm{C}\mathrm{l}}^{\left(M\right)}\mathrm{d}\mathrm{V}+{\int }_{\Omega }{\mathsf{\nabla }\mathbf{N}}_{\mathrm{C}}^{\mathrm{T}}{\mathsf{\alpha }}_{\mathrm{C}\mathrm{l}}\mathrm{d}\mathrm{V}+{\int }_{\Omega }{\mathbf{N}}_{\mathrm{C}}^{\mathrm{T}}{\mathbf{J}}_{\mathrm{w}}·{\mathsf{\nabla }c}_{\mathrm{C}\mathrm{l}\mathrm{W}}{v}_{\mathrm{C}\mathrm{l}}^{\left(M\right)}\mathrm{d}\mathrm{V}\\ \begin{array}{c}+{\int }_{\Omega }{\mathbf{N}}_{\mathrm{C}}^{\mathrm{T}}{v}_{\mathrm{C}\mathrm{l}}^{\left(M\right)}{c}_{\mathrm{C}\mathrm{l}\mathrm{W}}\left({-v}_{\mathrm{u}}^{\left(M\right)}+2v_{\mathrm{a}}^{\left(M\right)}-v_{\mathrm{C}\mathrm{a}}^{\left(M\right)}\right){\displaystyle \frac{1}{m_{\mathrm{c}}^{\left(M\right)}}}{\displaystyle \frac{\mathrm{d}{m}_{c\mathrm{S}}}{\mathrm{d}t}}\mathrm{d}\mathrm{V}\\ \_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\end{array}\end{array}\end{array}\\ {\int }_{\Omega }{\mathbf{N}}_{M_{\mathrm{C}}}^{\mathrm{T}}\left({\displaystyle \frac{\mathrm{d}{m}_{c\mathrm{S}}}{\mathrm{d}t}}-m_{\mathrm{c}}^{\left(M\right)}{n}_{\mathrm{w}}r\right)\mathrm{d}\mathrm{V}\end{array}\right\}.$$

(32)

From the element internal force vector, the element tangent stiffness and tangent diffusion matrices may be retrieved as ${\mathbf{K}}^{\left(\mathrm{e}\right)}=\mathrm{d}{\mathbf{F}}_{\mathrm{i}\mathrm{n}\mathrm{t}}^{\left(\mathrm{e}\right)}/\mathrm{d}{\mathbf{X}}^{\left(\mathrm{e}\right)}$ and ${\mathbf{C}}^{\left(\mathrm{e}\right)}=\mathrm{d}{\mathbf{F}}_{\mathrm{i}\mathrm{n}\mathrm{t}}^{\left(\mathrm{e}\right)}/\mathrm{d}{\mathbf{V}}^{\left(\mathrm{e}\right)}$, where ${\mathbf{X}}^{\left(\mathrm{e}\right)}$ is the unknown vector containing the nodal primary variables (displacements, water pressure, concentrations of ionic species and calcite mass content) as follows:

$${\mathbf{X}}^{\left(\mathrm{e}\right)}=\left[\begin{array}{ccc}{\mathbf{U}}^{\left(e\right)}& {{\mathbf{P}}_{\mathrm{w}}}^{\left(e\right)}& \begin{array}{ccc}{{\mathit{C}}_{\mathrm{u}\mathrm{W}}}^{\left(e\right)}& {{\mathit{C}}_{\mathrm{a}\mathrm{W}}}^{\left(e\right)}& \begin{array}{ccc}{{\mathit{C}}_{\mathrm{C}\mathrm{a}\mathrm{W}}}^{\left(e\right)}& {{\mathit{C}}_{\mathrm{C}\mathrm{l}\mathrm{W}}}^{\left(e\right)}& {{\mathbf{M}}_c}^{\left(e\right)}\end{array}\end{array}\end{array}\right]$$

(33)

and ${\mathbf{V}}^{\left(\mathrm{e}\right)}={\dot{\mathbf{X}}}^{\left(\mathrm{e}\right)}$.

The Gaussian quadrature rule is applied to obtain the components of the internal forces vector (32), as well as the components of the stiffness and diffusion matrices. Due to the presence of highly nonlinear terms in this vector and in these matrices, the number of integration points per element is two, for all vectors and matrices, despite quadratic interpolation being used for the displacements and linear interpolation being used for the remaining primary variables. A list with the formulas of additional dependent variables is presented in Appendix B, while the block structure of the element diffusion and stiffness matrices can be found in Appendix C.

In this 1-D formulation, essential boundary conditions may be imposed at the top and the bottom of the soil column to specify the values of the primary variables: the displacements, water pressure or ionic concentrations. Natural boundary conditions may be imposed at the top and at the bottom of the soil column to specify forces per unit of area, the water flux or the diffusive fluxes of the urea, calcium and chloride ionic species (see (31)). On the other hand, the extension of this formulation to 2-D or 3-D biocementation problems is straightforward.

3.2. Time Integration

The formulated semi-discrete equations are integrated in time. A midpoint scheme is employed for the integration, signifying that, for each increment, the equations are evaluated at time ${\mathrm{t}}_{n+\mathsf{\alpha }}={\mathrm{t}}_n+\mathsf{\alpha }\mathsf{\Delta }\mathrm{t}$, with $\mathsf{\alpha }=0.5$ and $\mathsf{\Delta }\mathrm{t}={\mathrm{t}}_{n+1}-{\mathrm{t}}_n$.

Since the internal forces are a non-linear function of the nodal primary unknowns, within each step, several iterations are performed until the system’s convergence is reached, that is, the residual global force $\mathbf{R}$ is zero:

$${\mathbf{R}}_{n+\mathsf{\alpha }}={\mathbf{F}}_{n+\mathsf{\alpha }}^{\mathrm{e}\mathrm{x}\mathrm{t}}-{\mathbf{F}}_{n+\mathsf{\alpha }}^{\mathrm{i}\mathrm{n}\mathrm{t}}\simeq \mathbf{0}$$

(34)

where ${\mathbf{F}}^{\mathrm{e}\mathrm{x}\mathrm{t}}$ is the global vector of external forces, dependent on the nodal “loads” $\mathbf{S},$ while the global vector of internal forces depends on $\mathbf{X}$ and $\mathbf{V}$. The quantities $\mathbf{A}=\mathbf{S},\mathbf{X},\mathbf{V}$, evaluated at time ${\mathrm{t}}_{n+\mathsf{\alpha }}$ are defined as ${\mathbf{A}}_{n+\mathsf{\alpha }}=\left(1-\mathsf{\alpha }\right){\mathbf{A}}_n+\mathsf{\alpha }{\mathbf{A}}_{n+1}$.

The Newton–Raphson method is applied, with the previous linearization of the equations around an approximation of the solution through the development of the internal forces vector in Taylor series neglecting the higher-order terms, requiring the residual vector at the next iteration to be zero:

$${\mathbf{R}}_{n+\mathsf{\alpha }}^{i+1}={\mathbf{F}}_{n+\mathsf{\alpha }}^{\mathrm{e}\mathrm{x}\mathrm{t}}-{\mathbf{F}}_{n+\mathsf{\alpha }}^{\mathrm{i}\mathrm{n}\mathrm{t}}\simeq {\mathbf{R}}_{n+\mathsf{\alpha }}^i-{{\mathsf{\alpha }\mathbf{C}}^{\mathrm{*}}}_{n+\mathsf{\alpha }}^i\Delta \mathbf{V}=\mathbf{0},$$

(35)

where ${\mathbf{C}}^{\mathrm{*}}=\mathbf{C}+\mathsf{\alpha }\Delta \mathrm{t}\mathbf{K}$ is the effective tangent diffusion matrix. Equation (35) is solved in order to obtain the increment of velocity $\Delta \mathbf{V}$ that is used to update $\mathbf{V}$ and $\mathbf{X}$, according to the following:

$$\left\{\begin{array}{c}for i=0 { \mathbf{V}}_{n+1}^0={\mathbf{V}}_n, {\mathbf{X}}_{n+1}^0={\mathbf{X}}_n+\left(1-\alpha \right)\Delta t{\mathbf{V}}_n\\ for i\ge 1 {\mathbf{V}}_{n+1}^i={\mathbf{V}}_{n+1}^{i-1}+\Delta \mathbf{V}, {\mathbf{X}}_{n+1}^i={\mathbf{X}}_{n+1}^0+\alpha \Delta t{\mathbf{V}}_{n+1}^i.\end{array}\right.$$

(36)

While the choice of the steps’ number used in the incremental process may be variable and free, the number of iterations within each step depends on a predefined criterion of convergence. As the iteration process progresses, returning in each iteration $\Delta \mathbf{V}$ that is used to update the velocity vector, the convergence of the process is achieved when the residual, and consequently $\Delta \mathbf{V}$, vanish. Thus, the proposed convergence criterion is as follows:

$${\displaystyle \frac{{‖\overline{\Delta \mathbf{V}}‖}^{\mathit{i}}}{{‖\overline{\Delta \mathbf{V}}‖}^{\mathbf{0}}}}\le TOL,$$

(37)

where ${‖\overline{\Delta \mathbf{V}}‖}^{\mathit{i}}$ is the norm of a nondimensional $\Delta \mathbf{V}$ vector at step $n$ and iteration $i$, and $TOL$ is the tolerance, which is a very small number.

The nondimensionalisation of $\Delta \mathbf{V}$ is required since both $\mathbf{X}$ and $\mathbf{V}$ encompass nodal values with different physical dimensions and orders of magnitude (e.g., displacements, water pressure, concentrations of ionic species and the mass content of calcite). Consequently, the positions related to the displacement variables in the computed $\Delta \mathbf{V}$ vectors, are divided by the corresponding maximum value occurring at the initial (0) iteration. The same is done for the positions in $\Delta \mathbf{V}$ related to the water pressure, to the concentrations of the ionic species and to the calcite mass content:

$$\begin{gathered} {\left({\overline{\Delta \mathbf{V}}}^i\right)}_{\mathbf{U}}={\displaystyle \frac{{\left({\Delta \mathbf{V}}^{\mathit{i}}\right)}_{\mathbf{U}}}{{\left({\Delta \mathbf{V}}^{\mathbf{0}}\right)}_{\mathrm{m}\mathrm{a}\mathrm{x},\mathbf{U}}}}, {\left({\overline{\Delta \mathbf{V}}}^i\right)}_{{\mathbf{P}}_{\mathrm{w}}}={\displaystyle \frac{{\left({\Delta \mathbf{V}}^{\mathit{i}}\right)}_{{\mathbf{P}}_{\mathrm{w}}}}{{\left({\Delta \mathbf{V}}^{\mathbf{0}}\right)}_{\mathrm{m}\mathrm{a}\mathrm{x},{\mathbf{P}}_{\mathrm{w}}}}}, {\left({\overline{\Delta \mathbf{V}}}^{\mathrm{i}}\right)}_{{\mathbf{C}}_{k\mathrm{W}}}={\displaystyle \frac{{\left({\Delta \mathbf{V}}^{\mathit{i}}\right)}_{{\mathbf{C}}_{k\mathrm{W}}}}{{\left({\Delta \mathbf{V}}^{\mathbf{0}}\right)}_{\mathrm{m}\mathrm{a}\mathrm{x},{\mathbf{C}}_{k\mathrm{W}}}}}, \mathrm{and} {\left({\overline{\Delta \mathbf{V}}}^{\mathrm{i}}\right)}_{{\mathbf{M}}_{\mathrm{c}}} \\ ={\displaystyle \frac{{\left({\Delta \mathbf{V}}^{\mathbf{i}}\right)}_{{\mathbf{M}}_{\mathrm{c}}}}{{\left({\Delta \mathbf{V}}^{\mathbf{0}}\right)}_{\mathrm{m}\mathrm{a}\mathrm{x},{\mathbf{M}}_{\mathrm{c}}}}}. \end{gathered}$$

(38)

When this is performed, the same tolerance $TOL$ may be used for all variables.

4. Experimental Procedure and Results

4.1. Soil

A uniform grading-size commercial river sand named APAS 30 was used in the experimental tests. This soil is composed of fine-sized quartz grains with a specific gravity ${\mathrm{G}}_{\mathrm{S}}\approx 2.6.$ The samples were prepared, aiming for a bulk density around $1.50~1.52\mathrm{g}/{\mathrm{c}\mathrm{m}}^3,$ which corresponds to an initial void ratio of $e_{\mathbf{0}}\approx 0.71$ and an initial porosity ($n_{\mathrm{W}\mathbf{0}}$) of $41.5\mathrm{\%}$.

4.2. Bacteria and Feeding Solution

Sporosarcina pasteurii are the bacteria selected for this study. This type of bacteria is commonly found in nature and is harmless to humans. Bacteria were supplied by the American Type Culture Collection (ATCC) entity and were lyophilized upon receipt. The growing process of Sporosarcina pasteurii was carried out until reaching a concentration of $~{10}^8$ bacteria/mL. More detailed explanation on bacteria growing stages and process is presented in [3].

To produce calcium carbonate, a feeding solution had to be prepared with the nutrients that are necessary to supply to the bacteria. This solution was prepared with 0.5 M of urea and 0.5 M of calcium chloride (CaCl₂).

4.3. Experimental Treatment Procedure

In each experimental test, a solution with the concentration of approximately $~{10}^8$ bacteria/mL was initially added to the sand, in order to saturate the soil (Figure 3a). The sand column tests were then carried out in small-scale specimens inside syringes with 2 cm diameter and 20 mL capacity (Figure 3b). A draining layer, composed of a geotextile and filter paper, was placed at the bottom of the syringe, followed by a layer-by-layer compaction of the sand previously saturated with bacteria. The sand column, with an assumed homogeneous distribution of bacteria, had a final height of 6 cm (Figure 3c).

The sand column tests started by irrigating 15.7 mL (the equivalent to two void volumes) of feeding solution on the top of the syringe. The upper part of the syringe was filled with feeding solution, while a downward fluid flow was established by opening the bottom of the syringe (Figure 3d). The application of feeding solution took 60 s. The exit at the bottom of the syringe was closed when the fluid on the top of the soil column achieved a constant height of 2 cm (Figure 3e), stopping fluid flow through the soil. In this manner, it was possible to investigate the precipitation of biocement with time and depth caused by diffusion of calcium and urea, assuming the bacteria would be immobile because they were attached to the soil particles. Several cases were investigated varying the time the treatment fluid was immobile at the top of the sand column, with some duplicates for estimating the experimental error, as summarized in

Table 2. Summary of the sand column tests performed.

Duration of Treatment	Quantity of Tests Performed
1 h	4
2 h	3
3 h	4
5 h	3
12 h	3
24 h	3

.

For each case, at the end of the treatment, 7.85 mL (one void volume) of distilled water was applied at the top of the column to wash away the feeding solution. The soil column was then removed from the syringe and divided in three approximately equal-sized fractions (top, middle and bottom of the specimen), which were subjected to dissolution tests in hydrochloric acid to measure the calcium carbonate content and its distribution along the height of the column.

The CaCO₃ content was determined through acid leaching tests performed in each fraction. Prior to the leaching tests, each fraction was dried in oven for 24 h at 105 °C to determine the mass of solids ($M_1$). Then, the fraction was washed with hydrochloric acid (0.5 M), which dissolved the calcium carbonate. The occurrence of the reaction was verified by gas liberation and pH measurements. After that, the sample was washed with distilled water, filtered and dried again in the oven at 105 °C for 24 h to measure the final mass ($M_2$). The percentage of calcium carbonate content ($\mathrm{\%}\mathrm{c}\mathrm{c}$) was then calculated as follows:

$$\mathrm{\%}\mathrm{c}\mathrm{c}=\left(M_1-M_2\right)/M_1$$

(39)

while the calcium carbonate mass content (referred to as calcite content, m_c) of each fraction was calculated as follows:

$${\mathrm{m}}_{c\mathrm{S}}=\left(M_1^i-M_2^i\right)/\left(V{M}_1^i/\sum M_1^i\right),$$

(40)

where $i$ refers to the bottom, middle and top fractions, respectively, and $V$ is the total volume of the specimen.

4.4. Experimental Results

The percentages of calcium carbonate content, %cc, measured in all tests are presented in

Table 3. Calcium carbonate content obtained from leaching tests.

Duration of Treatment (h)	Average %cc Measured in Each Column	%cc (Average)
1	0.92	0.84
	0.91
	0.78
	0.76
2	0.65	0.74
	0.94
	0.65
3	0.62	0.68
	0.74
	0.72
	0.64
5	0.79	0.69
	0.49
	0.78
12	0.84	0.75
	0.67
	0.75
24	0.95	0.83
	0.69
	0.84

(including average values for each soil column) considering the period of time that the treatment fluid was kept immobilized at the top of the columns (duration of treatment).

Despite some small variations, the precipitation of calcium carbonate remained practically constant after 1 h of treatment, probably because bacterial encapsulation had occurred. Since treatment was conducted under no-flow conditions after the first 60 s, the transport of ionic species after that time was restricted to purely diffusive flow, which occurred at very low velocities. This could also have reduced feeding fluid circulation and caused the accumulation of chemical compounds in the pores, therefore creating difficulties for the survival of bacteria [3]. On the other hand, the decrease in calcite content that sometimes was observed after 1 h of treatment could be attributed to the dissolution of irregularly shaped and less stable CaCO₃ minerals [24].

The first observation is in line with the findings of other authors. For example, Cuthbert et al. [25] evaluated the effects of high initial concentrations of feeding solution (0.5 M, the same concentration applied in this study) on calcite formation. They found precipitation of large calcium carbonate crystals immediately after the treatment due to higher ureolysis rate at that stage. According to the authors, these large calcium carbonate crystals that were attached to the bacteria blocked further access of the feeding solution. Xiao et al. [24] also studied the distribution of calcite in microchannels with different concentrations of the feeding solution. For larger concentrations, they observed a fast precipitation of calcium carbonate crystals without further development along the remaining treatment time. They explained this stabilization by the early encapsulation of bacteria. This type of behaviour can be mathematically described by Equation (23) for the ureolysis reaction, and therefore, experimental data could be used to calibrate the parameters in this equation. Because the amounts measured were stable after 1 h of treatment, the data found in these tests were used as reference.

Figure 4 presents the average values of the mass of calcite per unit of volume precipitated after 1 h of treatment, at the top, middle and bottom fractions. In general, more precipitated calcite was formed at the top, while less calcite mass per unit volume was formed at the bottom. This same distribution was observed in other experimental sand column studies ([26,27]) and could be explained by pore clogging of the upper fraction, occurring after the injection of the feeding solution.

5. Numerical Simulations

5.1. Model Parameters

The model was used to simulate the experiments; basically, a 1-D behaviour was observed in the soil columns. The adopted constants are presented in

Table 4. Values of the constants used in the simulation.

Constant	Value
Soil constants
Specific gravity ${(G}_{\mathrm{S}}={\rho }_{\mathrm{S}}/{\rho }_{\mathrm{W}})$	2.60
Young’s modulus (E)	50 × 103 kPa
Hydraulic conductivity parameter (${\overline{k}}_h)$	3.5 × 10−5 m/s
Initial porosity ($n_{\mathrm{W}\mathbf{0}}$)	0.415
Chemical constants
Calcite molar mass ($m_{\mathrm{c}}^{\left(M\right)}$)	100.1 × 10−3 kg/mol
Calcite molar volume ${(v}_{\mathrm{c}}^{\left(M\right)})$	36.9 × 10−6 m3/mol
Urea molar mass ${(m}_{\mathrm{u}}^{\left(M\right)}$)	60.1 × 10−3 kg/mol
Urea molar volume ${(v}_{\mathrm{u}}^{\left(M\right)})$	45.4 × 10−6 m3/mol
Ammonium molar mass ${(m}_{\mathrm{a}}^{\left(M\right)})$	18.0 × 10−3 kg/mol
Ammonium molar volume $\left(v_{\mathrm{a}}^{\left(M\right)}\right)$	20.0 × 10−6 m3/mol
Calcium molar mass ${(m}_{\mathrm{C}\mathrm{a}}^{\left(M\right)})$	40.1 × 10−3 kg/mol
Calcium molar volume ($v_{\mathrm{C}\mathrm{a}}^{\left(M\right)})$	1.75 × 10−6 m3/mol
Chloride molar mass $\left(m_{\mathrm{C}\mathrm{l}}^{\left(M\right)}\right)$	35.5 × 10−3 kg/mol
Chloride molar volume ($v_{\mathrm{C}\mathrm{l}}^{\left(M\right)}$)	15.42 × 10−6 m3/mol
Longitudinal diffusion constant (${\mathsf{\alpha }}_L$)	0.02 m
Reaction rate parameters
$v_{max}$	2 mol/(m3s)
$K_m$	125 mol/m3
$t_{max}$	260 s

. The reaction rate parameters ($v_{max},K_m,t_{max})$ and the diffusion constant (${\mathsf{\alpha }}_L$) were calibrated to better fit the calcite mass values measured experimentally after 1 h of treatment. This calibration process will be explained later in Section 5.3.

Bar (1-D) finite elements were employed with quadratic interpolation for the solid displacements (with three nodes per element) and linear interpolation for the other primary unknowns (with two nodes per element). Several simulations were performed varying the sizes of the elements, the time steps and tolerance values, until the convergence of the results was found. In

Table 5. Convergence study: percentage of total calcium carbonate obtained in numerical simulations after 1 h of treatment, as a function of mesh refinement and integration time step.

Number of Elements	Time Step
	0.05 s	0.01 s	0.005 s	0.001 s
30	0.75208%	0.75225%	0.75229%	0.75225%
60	0.75205%	0.75221%	0.75222%	0.75222%
90	0.75205%	0.75221%	0.75222%	0.75222%
120	0.75204%	0.75221%	0.75222%	0.75222%

, the percentage of total calcium carbonate obtained in the numerical simulations after 1 h of treatment is shown as a function of the mesh refinement and integration time step. The chosen model had a total of 60 elements, comprising 121 nodes and 487 global degrees of freedom. A constant time step of 0.01 s and a 0.001 tolerance were used in the time integration.

5.2. Initial and Boundary Conditions

For the calculation, it was assumed that the soil column was saturated with a homogeneous bacteria distribution for the studied depth, without the presence of urea, calcium, chloride or calcite. Two time stages were considered, as described in Figure 5. The first one (for t ≤ 60 s) corresponded to the initial addition of the feeding solution by establishing a constant fluid flow through the column. In the second one (for t > 60 s), the fluid flow was stopped and the treatment progressed by diffusion over different periods of time. The boundary conditions corresponding to both stages are also presented in Figure 5:

• Concerning the mechanical model, the bottom boundary was restrained, and therefore, all the displacements were null during the entire calculation;
• Concerning the hydraulic model, the flow is free during the first 60 s at the bottom boundary, and the pore pressures remained constant at the top boundary after that;
• Concerning the biochemical model, at the top boundary, during the first 60 s when the feeding solution was applied, there was a sudden increase in the concentration of urea and calcium from 0 to 0.5 M and in the concentration of chloride from 0 to 1 M; these concentrations remained constant until the end of the analysis. These values were set to mimic the application of the feeding solution in the experimental tests, in which the fluid flow applied at the top of the columns was equal to 0.2 mm/s during the first 60 s of the analysis and in which a layer of feeding solution remained, with a constant height and constant concentration, at the top of the syringe until the end of the analysis.

5.3. Calibration of the Reaction Rate Parameters and of the Diffusion Constant

The calibration of the reaction rate parameters in (13) and of the diffusion constant in (18) was performed after several numerical simulations, carried out with values varying among 1~10 mol/(m³s) for the maximum reaction velocity ($v_{max}$), 100~300 mol/m³ for the half-saturation constant ($K_m$), 200~3600 s for $t_{max}$ and 0.01~0.2 m for the diffusion constant (${\mathsf{\alpha }}_L$). The values presented in Table 4 were the ones that fitted the calcite mass values measured experimentally for treatments lasting 1 h (3600 s) better.

In Figure 6a, it may be seen that the numerical profile of the mass of calcite per unit of volume obtained with the calibrated parameters presented in Table 4 fits the minimal and maximal experimental values for the top and bottom fractions, with both experimental and numerical values presenting the same distribution trend. Moreover, it may be seen in Figure 6b that, for the calibrated parameters, the reactions stopped around t = 1200 s.

5.4. Results and Discussion

Figure 7 presents the profiles along the column height of the chemical species involved in the biochemical reactions ($c_{\mathrm{u}\mathrm{W}},c_{\mathrm{C}\mathrm{a}\mathrm{W}},c_{\mathrm{C}\mathrm{l}\mathrm{W}},c_{\mathrm{a}\mathrm{W}}$) computed at different time instants. It may be observed in Figure 7a–c that, as expected, the concentrations of urea, calcium and chloride ions remained constant at the top boundary because this was the point of application of the feeding solution. Along the sand column, it can be confirmed that the concentrations of urea (Figure 7a) and calcium (Figure 7b) decreased at the same rate with time, which is expected because urea and calcium are consumed in the reactions. The concentrations along the column of chloride ions increased only during the first 60 s, while the feeding solution was added, and then remained constant because these ions do not participate in the MICP process (Figure 7c).

The progression of the biochemical reaction was identified by the reduction in the concentrations of urea and calcium, and consequently by the increase in the ammonium concentrations (Figure 7d). Accordingly, there was an increase in the mass of calcite per unit volume with time (Figure 8a), which was always more noticeable at the top of the sand column. This calcite distribution trend is confirmed by other authors for similar boundary conditions ([11,27]). The ammonium concentrations and calcite mass remained constant after 1200 s because then the reaction rate became practically zero (see Figure 6b). The precipitation of solid calcite into the voids reduced the soil porosity, which is confirmed in Figure 8b.

The evolution with the time of fluid flow is presented in Figure 9a. During the first 60 s of the simulation, the fluid flow was practically constant along the specimen and equal to the flow applied at the top of the column; when the bottom exit was closed, the flow was reduced to zero. The evolution with time of the pore–water pressure is presented in Figure 9b. The water pressure profile started as hydrostatic before applying the boundary condition, and then it increased at the top of the column, creating a constant hydraulic gradient in the first 60 s. Then it became hydrostatic again, once the flow at the bottom was interrupted and the feeding solution accumulated above the soil column, becoming equal to the weight of the liquid column above this horizontal section. Due to the hydro-mechanical coupled behaviour, this slight pore pressure increment implies a reduction in effective stresses, and for this reason, very small swelling deformations (about 0.0067%) were computed (Figure 9c). The observed very small deformations justify the linear elastic effective stress–strain relationship that was assumed in this work.

Finally, Figure 10 presents a comparison between numerical and experimental values of the total calcium carbonate content for all tested columns and all treatment durations. The experimental average results are represented together with the minimal and maximal experimental values (see Table 3). The numerical simulation curve presents a constant value that intercepts all the experimental data range, therefore validating the model. The numerical values shown in Figure 10 are independent of the treatment duration because, according to Figure 6b, the reactions stopped around t = 1200 s. This can be explained by bacteria encapsulation, as already discussed ([24,25]). The reduction in calcite content that is sometimes observed in the experimental values in Figure 10 could possibly be attributed to the dissolution of irregular-shaped and less stable CaCO₃ polymorphs ([24]).

6. Conclusions

A finite element numerical model to perform a coupled biochemical-hydro-mechanical analysis was developed for the simulation of the biocementation process in a sand column. This was intended to be a tool to compute the amount of biocement precipitated for a given set of treatment dosages and hydraulic and mechanical boundary conditions. With the hydraulic model, it was possible to simulate the transport by conduction and diffusion of the ions, with the hydro-mechanical coupled model being related to changes in the soil’s effective stresses. The precipitation of solid mass (calcite, i.e., the biocement) determined porosity changes, affecting permeability (bio-chemo-hydraulic coupling). The model is prepared to be modified to include this pore-clogging effect into stiffness and strength (bio-chemo-mechanical coupling).

The parameters of the model were calibrated with the results of experimental small-scale sand column tests, in which the biocementation treatment had a duration of 1 h. The distribution of calcite obtained in the numerical simulations, with a noticeable concentration at the column’s top, was in good agreement with the experimental results for all the other tested treatment periods, therefore validating the model.

The reaction rate parameters determine the biochemical reaction and are of paramount importance in this analysis. Nevertheless, the imposed boundary conditions dictated a single application of the feeding solution and established a no-flow condition after the first 60 s of the treatment. This possibly limited the circulation of ionic specimens and lowered the bacterial activity, which explained the constancy of the calcite content value obtained for several treatment durations.

Finally, the approach adopted in this study can be used for other applications. Other bio-chemical reactions can be simulated if the chemical species and stoichiometric calculations are adapted for the specific case at hand, as well as the bacterial activity parameters. For example, it may be possible to simulate methane formation by methanogenic bacteria or corrosion by iron-oxidizing bacteria.