Practical approaches to study microbially induced calcite precipitation at the field scale

Microbially induced calcite precipitation (MICP) is a new and sustainable technology which utilizes biochemical processes to create barriers by calcium carbonate cementation; therefore, this technology has a potential to be used for sealing leakage zones in geological formations. The complexity of current MICP models and present computer power limit the size of numerical simulations. We describe a mathematical model for MICP suitable for field-scale studies. The model is simple enough to have a few equations but still retain the essence of the processes. The main mechanisms in the conceptual model are as follow: microbes attach themselves to the pore walls to form biofilm, the biofilm consumes substrate to grow, the biofilm use urea for production of calcite, and the calcite reduce the pore space which in turn decreases the rock permeability. We apply the model to study the MICP technology in two sets of reservoir properties including a well-established field-scale benchmark system for CO2 leakage. A two-phase flow model for CO2 and water is used to assess the leakage prior to and with MICP treatment. Based on the numerical results, this study confirms the potential for this technology to seal leakage paths in reservoir-caprock systems.


Introduction
Negative emissions technologies and carbon storage must be implemented to avoid dangerous climate changes (Haszeldine et al., 2018;Tong et al., 2019). Carbon capture and storage (CCS) is one of the promising scalable technologies for storing large amounts of CO 2 . Indeed, large amounts of CO 2 have already been stored in geological formations on the Norwegian continental shelf, e.g., in the Sleipner field, where more than 16 Mt CO 2 has been stored since 1996 (Furre et al., 2017). Caprocks in reservoirs provide the main trapping mechanism for CO 2 sequestration (Bentham and Kirby, 2005). The existence of faults, fractures, and abandoned wells in the primary sealing caprock of a CO 2 storage reservoir can create pathways for CO 2 to migrate back to the surface (Fang et al., 2010). Fig. 1 shows a schematic representation of CO 2 sequestration, where fractures in the caprock are a risk to leak CO 2 back to the atmosphere and to fresh water. It is therefore necessary to develop methods for mitigating CO 2 leakage to ensure its long-term storability. One of the proposed remediation measures to seal leakage zones is the use of microbes to induce precipitation of calcium carbonate (Phillips et al., 2016). Microbially induced calcite precipitation (MICP) is a new and sustainable technology which utilizes biochemical processes to create barriers by calcium carbonate cementation. The MICP technology involves the injection of diverse components into a reservoir such as microbes, substrate, and chemicals. As calcite permeability is very low, then the formation of calcite decreases the rock permeability. Thus, MICP technology has a potential to be used for sealing leakage zones in geological formations. These barriers can significantly reduce CO 2 leakage even when the leakage channels are not fully plugged (Li et al., 2019). Among other applications of MICP besides as a leakage prevention tool in CO 2 sequestration are in biomineralized concrete (Lee et al., 2018), improvement in the stiffness and strength of granular soils (Jalili et al., 2018;Whiffin et al., 2007), wastewater treatment (Torres-Aravena et al., 2018), and erosion control (Jiang and Soga, 2017).
MICP as a leakage mitigation technology is intended for use on the field scale, but performing field-scale experiments is expensive. Experiments on microsystems allow us to observe processes in more detail, which leads to improvements in core-scale experiments prior to field applications. We mention some notable works in this direction. Bai et al. (2017) performed MICP experiments in microfluidic cells to study the distribution of calcite precipitation at the pore scale and observed that calcite precipitation occurs mainly on the bottom of biofilms. Core samples from reservoirs can be used to study changes in permeability due to biofilm growth and calcite precipitation. For example, Whiffin et al. (2007) conducted a core-scale experiment to evaluate MICP as a soil strengthening process. Since the laboratory experiment was conducted under field conditions and a significant improvement of strength was observed along the column, the authors concluded that MICP can be used for largescale applications. Ebigbo et al. (2012) performed core-scale experiments under controlled conditions for studying the effect of calcite precipitation in porous media. The authors tested different injection strategies to obtain a homogeneous distribution of calcite precipitation along sand-packed columns. Their work provides a successful injection strategy for this purpose and experimental data of four columns. Mitchell et al. (2013) investigated the MICP processes in a core sample inside a high pressure flow reactor including supercritical CO 2 to simulate field conditions. Their experimental results show that MICP can be applied in the presence of supercritical CO 2 . Gomez et al. (2017) performed experiments in tanks of 1.7 m diameter and with three wells to study the reactive transport of MICP. Their results show that indigenous microorganisms could be stimulated for MICP in field-scale applications. Based on these and more experimental work reported in literature, mathematical models of this technology can be built for further studies.
Mathematical models of MICP are important as they help to predict the applicability of this technology and to optimize its benefits. Zhang and Klapper (2010) introduced a comprehensive pore-scale model for MICP which includes chemistry, mechanics, thermodynamics, fluid, and electrodiffusion transport effects. The authors performed simulations under different conditions of flow rates, concluding that the flow significantly impact the calcite distribution. Hommel et al. (2015) introduced a core-scale mathematical model for MICP which includes chemistry, mechanics, and fluid transport effects. The authors also calibrated some of the model parameters with experimental data. Minto et al. (2019) proposed a mathematical model for MICP and performed numerical studies of calcite precipitation around a production well using eight surrounding injection wells. The authors concluded that uniform calcite precipitation could be achieved by splitting the injection into phases, where different number of wells are used in each of the injection phases. Note that "phases" is used to denote both physical phases and repeatable steps in the injections strategies, and the meaning will be clear from the context.
Despite advances in modeling, simulation of the MICP process at the field-scale is challenging as current mathematical models involve the solution of large systems of highly coupled partial differential equations. In Cunningham et al. (2019) the authors suggest different approaches to handle this issue such as refinement of the grid locally, multi-scale methods, improving the time stepping, or reducing the coupling of the model equations. Tveit et al. (2018) proposed a simplified version of the MICP model presented in Hommel et al. (2015) to perform field-scale simulations. The authors study two different approaches for inducing calcite precipitation at a given distance of an injection well. Since the complexity of current MICP models and present computer power limit the size of numerical simulations, then simplified models are needed to perform field-scale studies. In this work, we build a single-phase field-scale model of MICP technology. This model includes the transport of dissolved components (microbes, substrate, and urea), biofilm activity (microbial attachment, death, detachment, and growth), and production of calcite which reduce the rock porosity and hence the effective permeability. We use the model to investigate the prevention and sealing of leakage paths located at a certain distance away from the injection well.
A simple two-phase flow model for CO 2 and water is used to assess the leakage prior to and with MICP treatment.
Our motivation to develop the mathematical model and numerical tools is as follows. We aim to have a model that captures the key processes and quantities involved in the MICP process. At the same time, we aim to have a model which is simple enough so that computational costs are less. Our main reason for the latter is that the field-scale processes require running the model on a large scale and also require multiple simulations to perform optimization studies. All these imply a heavy computational burden unless we simplify the model. Needless to say, the simplified model should still retain the essence of the processes so that it is useful. Our work is therefore a step in this direction of achieving the twin objectives. The paper is structured as follows. In Section 2 we explain in detail the MICP mathematical model, the model parameters, the computer implementation of the model, and the injection strategy. Diverse field-scale numerical experiments to prevent CO 2 leakage using MICP are presented in Section 3. A discussion on the numerical results and findings is given in Section 4. Finally, we present the conclusions in Section 5.

MICP model
In this section we describe the mathematical model for MICP, introducing first concepts and definitions related to this technology. Fig. 2 shows a schematic representation of the sealing mechanism using MICP. Here, we observe a fractured zone in the caprock being remediated by calcite. MICP can be defined as a bio-geochemical process which results in precipitation of calcite (the low-pressure, hexagonal form of CaCO 3 ). Calcium carbonate (CaCO 3 ) is a mineral that naturally precipitate as a result of microbial metabolic activities. Biofilm formation is a process whereby microorganisms attach themselves to a surface and produce an adhesive matrix of extracellular polymeric substances (EPS). A biofilm requires different factors for survival, e.g., nutrients and oxygen (for aerobic microbes). In this work, we denote the solution which includes all the different components a biofilm needs for growing (e.g., nutrients, oxygen) as substrate. In addition, we denote suspended biomass in water as microbes. Urea [CO(NH 2 ) 2 ] is a water-soluble compound found in the urine and other bodily fluids of mammals or produced synthetically. Urease is an enzyme catalyzing the hydrolysis of urea to ammonium (NH + 4 ) and carbonate (CO 2− 3 ). Sporosarcina pasteurii is a non-pathogenic bacterium commonly used for the MICP process as it shows a high urease enzyme activity (Bhaduri et al., 2016). Calcium (Ca 2+ ) is a mineral which contributes to processes of organisms cell. With the main concepts introduced, we proceed to describe the conceptual and mathematical model of MICP used in this paper.

Conceptual model
We consider a constant-temperature reservoir saturated with water, where calcite and biofilm only occur on the rock walls, i.e., in the space domain there are one liquid phase (water) and three solid phases (biofilm, calcite, and rock matrix). The microbes, substrate, and urea are dissolved in water prior to injection and they are transported only in the water phase by advection and dispersion. The biofilm and calcite are assumed to be impermeable and incompressible. The governing processes in the biofilm are growth, death, attachment, and detachment.
The most studied MICP process is urea hydrolysis (ureolysis) via the enzyme urease produced by special microbes, in a calcium-rich environment (Rong et al., 2012;Whiffin et al., 2007): In general, CaCO 3 precipitation is governed by four main factors (Hammes and Verstraete, 2002): calcium concentration, carbonate concentration, pH, and availability of nucleation sites. Lauchnor et al. (2015) performed experiments on S. pasteurii, showing that urea and microbial concentration have a more significant impact on the ureolysis rate than pH variations. In addition, we consider that the amount of urease is only related to the amount of biofilm, neglecting the microbes in the liquid phase as its contribution is minor (Ebigbo et al., 2012). Assuming enough calcium concentration in the water, we model the calcite formation as a function only dependent on urea and biofilm. This assumption can be justified since calcium and urea can be injected together, and would thus distribute in a similar manner, ensuring that both concentrations are high in the location where calcite precipitation is aimed.
To summarize, the system of interest consists of a 3D reservoir (porous medium), one source (injection well), one fluid phase (water), two solid phases (biofilm and calcite), and three transported components (microbes, substrate, and urea).

Mathematical model
We build a mathematical model based on the assumptions laid out in the conceptual model. We adopt a continuum approach, where the processes in the system are described by conservation laws and coupling relationships. We use the subscripts {b, c, m, s, u, w} to refer to biofilm, calcite, microbes, substrate, urea, and water respectively.

Flow equations
The mass conservation and Darcy's law equations for the water phase are given by: where φ is the rock porosity, p w the reservoir pressure, v v v w the flow velocity, ρ w the fluid density, K the absolute permeability, g g g the gravity, µ w the water viscosity, and q w the source/sink term.

Leakage paths
We adopt a common approach found in Class et al. (2009), where the leakage paths in the caprock are modeled as a porous medium with higher permeability than the formation. An advantage of this approach is that the model equations do not need further modification for implementation while a drawback is that it requires a fine grid to represent explicitly the leakage paths. This may be contrasted with the widely-used approach of discrete fracture networks (DFN) where one uses a mixed dimensional setting and represents the fractures as a n − 1 dimensional objects embedded in a n dimensional porous geometry. We refer to Berre et al. (2019) for a recent review of different conceptual models for fracture and Kumar et al. (2020) and Martin et al. (2005) for a formal derivation of some of these models. We also mention that our model here can be easily adapted for different conceptual models including fractures being modeled as DFNs.

Injection well
In this work we consider only one injection well, where we use the following linear law to relate the injection rate Q w (Lie, 2019;Peaceman, 1978): Here, p I is the pressure inside the wellbore, h is the height of the well, r I the well radius, z bh a reference depth (perforation), K is the permeability in the direction of the injection, and r e an equivalent radius at which the cell pressure is equal to the average of the exact pressure.

Transport equations
To describe the transport of microbes, substrate, and urea, we consider the following advection-dispersion-reaction transport equations: Here, c ξ is the concentration of component ξ in water, J J J ξ the flux of ξ, R ξ the reaction term of ξ, and D ξ the dispersion tensor. Here, we assume that the aqueous phase density does not depend on the component concentrations.

Dispersion effects
When the components are transported throughout the reservoir, two different mechanisms affect their movement: mechanical dispersion and molecular diffusion. The former is an effect arising out of mixing due to flow and heterogeneities while the latter accounts for movement of the components from a region of higher to lower concentration.
We adopt the following model for the dispersion of components (Bear, 1972) where α L and α T are the longitudinal and transverse dispersion coefficients, v v v = v v v w /φ is the effective velocity of the aqueous phase, and D ξ the effective diffusion coefficient of component ξ.

Solid-phase equations
As previously mentioned, we consider biofilm formation and calcite precipitation fixed in space (at the pore scale, it represents the biofilm and calcite precipitate at the rock surface). Thus, the following mass balance equations describe the evolution of biofilm and precipitation of calcite where ρ χ are densities and R χ accounts for the calcite and biofilm processes and are being described later in this section.

Substrate consumption
The substrate consumption rate R s is modeled by the Monod equation where µ s is the maximum rate of substrate utilization and k s is the half-velocity coefficient for the substrate.

Urea utilization
The urea conversion rate R u is modeled by the Monod equation (Hommel et al., 2015;Lauchnor et al., 2015) where µ u is the maximum rate of urea utilization and k u is the half-velocity coefficient for urea. This model for ureolysis was introduced in Hommel et al. (2015) based on the work by Lauchnor et al. (2015), where µ u is split into maximum urease activity and mass ratio of urease to biofilm.

Calcite precipitation
The calcite precipitation is the result of a complex geochemical process. Qin et al. (2016) stated that in a calciumrich environment, when only the final calcite distribution is of interest, an approximation of the calcite precipitation rate can be given by the negative value of the urea utilization rate (i.e., R c = −R u ). Since the molar mass of urea is different from calcite, we add a yield coefficient Y uc (units of produced calcite over units of urea utilization) to account for this in the mathematical model. Then, we write the calcite production rate as We note that R c only depends on the amount of biofilm and urea, which significantly reduces the computational cost compared to more complex formulations [e.g., Ebigbo et al. (2012), Hommel et al. (2015), and Minto et al. (2019)].

Microbes
Two opposing processes determine the evolution of microbes: growth and loss. The growth term comprises of two contributions. First, the consumption of substrate by the microbes lead to its growth. This is modeled by a Monod equation c m φY sb µ s c s /(k s + c s ) where Y sb is a yield coefficient (units of produced biomass over units of consumed substrate). Second, its growth taking place via detachment or erosion of biofilm due to flow. Microbes detach from the biofilm back to the water phase due to shear forces on the interface by the water flow. The erosion is modeled by φ b ρ b k str φ||∇p w − ρ w g g g|| 0.58 where k str is the detachment rate (Rittmann, 1982). The loss term also has two contributions. First, the death of the microbes as a result of aging, which is modeled by a linear death rate −c m φk d where k d is the microbial death coefficient. Second, the microbes attach themselves to the pore wall and biofilm. This is modeled by a linear attachment rate −c m φk a where k a is the microbial attachment coefficient. In sum, the rate for the microbes R m is given by

Biofilm processes
As in the case of microbes above, the biofilm development is determined by the net of its growth and loss. Consumption of substrate by the biofilm lead to its growth. This is modeled by the Monod equation The microbes in the biofilm die as a result of aging and being encapsulated by the produced calcite (De Muynck et al., 2010). The former is modeled by a linear death rate (Ebigbo et al., 2012). As described previously, the microbial attachment leading to its growth is modeled by c m φk a while the erosion leading to its loss is expressed by −φ b ρ b k str φ||∇p w − ρ w g g g|| 0.58 . In sum, the rate for the evolution of the biofilm is given by

Porosity reduction
The void space in the porous medium change in time as a function of the biofilm and calcite volume fractions φ b and φ c respectively. Using the definitions of φ b and φ c , we have the following equality

Permeability modification
Porosity-permeability relationships are used frequently in mathematical modeling to account for permeability reduction as a result of biofilm and calcite growth. Diverse porosity-permeability relationships have been proposed for the last decades. These relationships can also include the permeability of biofilm and be derived as a result of upscaling pore-scale models van Noorden et al., 2010). In this paper, we follow Thullner et al. (2002) and use a porosity-permeability relationship where significant reduction in CO 2 leakage can be achieved even when the leakage paths are not fully plugged, Here, K 0 is the initial rock permeability, φ crit is the critical porosity when the permeability becomes a minimum value K min , and η is a fitting factor.

Remarks on the MICP model
The development of the present mathematical model is inspired by previous works on the MICP technology Ebigbo et al., 2012;Hommel et al., 2015;Lauchnor et al., 2015;Qin et al., 2016). One of the most complete models for the MICP technology is presented in Hommel et al. (2015). This MICP model includes detailed chemistry reactions, mechanics, and fluid transport effects. Given the complexity of the model and the current computing power, solving simultaneously all equations would limit the size of the problem. Hence we build a simpler mathematical model so that the computational costs are less but the essence of the processes is retained. We summarize the main assumptions that we have adopted to build the simplified MICP model: only one fluid phase (water) and two solid phases (calcite and biofilm) are presented, there are only three transported components (microbes, substrate, and urea) dissolved in the fluid phase, the amount of urease is only related to the amount of biofilm, and the calcite formation only depends on urea and biofilm. The mathematical model is given by Eqs. (1-12). This model consists of six mass balance equations and six cross coupling constitutive relationships.

Implementation
The EOR module in the MATLAB ® reservoir simulation tool (MRST), a free open-source software for reservoir modeling and simulation, is modified to implement the MICP mathematical model (Bao et al., 2017;Lie, 2019).
Specifically, the polymer example (black-oil model + one transport equation) is modified (single-phase flow + three transport equations + two mass balance equations) to solve the MICP mathematical model. A comprehensive discussion of the solution of the polymer model can be found in Bao et al. (2017). The MICP mathematical model is solved on domains with cell-centered grids. Two-point flux approximation (TPFA) and backward Euler (BE) are used for the space and time discretization respectively. The resulting system of equations is linearized using the Newton-Raphson method. In contrast to the polymer model, we implement dispersion of the transported components, permeability changes due to calcite and biofilm formation, and biofilm detachment due to shear forces. The spatial discretization is performed using internal functions in MRST and the external mesh generator DistMesh (Persson and Strang, 2004). The MICP processes can be simulated over time and the simulator stops when full-plugging of at least one cell is reached (i.e., φ = 0). The links to download the corresponding code can be found above the references at the end of the manuscript.

Model parameters
Mathematical models require the numerical values of coefficients in the equations to be solved. These model  Table 1 summarizes the model parameters for the numerical simulations. We comment on the maximum rate of urea utilization µ u , yield coefficient Y uc , and minimum permeability K min . Lauchnor et al. (2015) estimated values for the kinetics of ureolysis by S. pasteurii. We consider a value of µ u = 1.61 × 10 −2 s −1 [here we use the value of mass ratio of urease to biofilm of 3.81 × 10 −4 and 0.06 kg/mol for urea multiplied by 706.7 mol/(kg s) ]. The molar mass ratio of calcite (0.1 kg/mol) to urea gives a value of 1.67 for the yield coefficient Y uc . The value of K min is set to 10 −20 m 2 which is of the order of magnitude of permeability in a caprock to retain fluids for CCS (Schlumberger, 2020).
The equivalent radius r e for the injection well depends on the grid. For a domain with rectangular grid blocks, the equivalent radius is given by r e = 0.14 ∆x 2 + ∆y 2 (Peaceman, 1978). We set the well radius to r I =0.15 m (Ebigbo et al., 2007). Regarding input concentrations, the maximum amount of urea dissolved in water is limited by its solubility, e.g., 1079 kg/m 3 at 20°C. In the MICP experiment reported in Whiffin et al. (2007) the concentration of urea corresponds to 66 kg/m 3 . The concentration of microbes is typically given in colony forming units (CFU) or in optical density of a sample at 600 nm (OD 600 ). Two values of concentrations for S. pasteurii used in experiments and reported in literature are 4×10 7 CFU/ml and 1.583 OD 600 . The former is equivalent to 0.01 kg/m 3 using a cell weight of 2.5×10 −16 kg/CFU (Norland et al., 1987) while the latter is approximately equal to 17×10 8 CFU/ml (Jin et al., 2018), which, using the cell weight, is converted to 0.425 kg/m 3 . Different substrate concentrations are reported in literature, e.g., 5 kg/m 3 of cellobiose (Linville et al., 2013). Here we consider the following concentrations: c m =0.01 kg/m 3 , c s =1 kg/m 3 , and c u =66 kg/m 3 . While the urea and microbial concentrations are taken from laboratory experiments, we choose a substrate concentration such that the substrate is totally consumed by the microbes.
Different studies can be conducted on mathematical models with a few parameters. For example, sensitivity analysis on the mathematical model allows us to identify critical model parameters. We refer to Landa-Marbán et al.  Given that the position of the well is fixed in the domain, the control variables for the injection strategy are the flux rate (water) along the height of the well, i.e., Q w (z, t) and concentrations of components (microbes, substrate, and urea). This injection strategy involves several phases where the three components are injected in the following order: microbes, substrate, and urea. First, microbes are injected for a total time t I 1 . This injection is followed by water injection (t I 2 ) to move the microbes away from the injection well. Subsequently, there is a no-flow period to facilitate attachment of microbes to the pore walls (t I 3 ). Substrate is injected (t I 4 ), followed by water displacement (t I 5 ), and subsequently there is a no-flow period (t I 6 ) to stimulate biofilm formation away from the injection well and around the sealing target. Urea is injected (t I 7 ), displaced by water (t I 8 ), and subsequently a no-flow period (t I 9 ) to precipitate calcite at the sealing target. We refer to these nine stages as phase I. Several phases can be applied to decrease the permeability in the targeted zone, see Fig. 3.

Numerical studies
In this section, we consider several examples that are divided into two parts. In the first part we study MICP in systems where we target calcite precipitation at selected parts of the aquifer [e.g., Minto et al. (2019), Nassar et al. (2018), and Tveit et al. (2018)]. This mimics a situation where MICP technology is applied to prevent formation of leakage paths in the caprock, that is in regions with closed fractures/faults that could be opened when CO 2 is injected. In the second part we study MICP in systems where leakage paths are modeled explicitly [e.g., Cunningham et al. (2019)]. Here, we focus on the benchmark problem introduced in Ebigbo et al. (2007) and Class et al. (2009), where two aquifers are separated by a caprock with a leakage path.
For the numerical examples we consider two set of reservoir properties. The first set of properties is taken from Tveit et al. (2018), where the authors studied the MICP technology for sealing at a given distance of an injection well. One of the motivations to include the same properties in this work is to compare qualitatively the simulation results between the two different model implementations. The second set of properties is taken from Ebigbo et al. (2007). Let K A denote the permeability in the aquifer, K L the permeability of the leakage path, L the length, W the width, and H the height of the aquifer, h the height of the caprock, l the distance of the leakage zone from the well, a the aperture of the leakage path, and ω the aperture of the potential leakage zone. Table 2 summarizes the properties of both systems. In the examples, we will perform simulations on 1D, 2D, and 3D domains. On each of the domains we will study different aspects of the injection strategy (Section 2.5) and the dynamics of the MICP process. The learnings from one domain will be useful by itself, but will also inform the studies on the other domains, ultimately leading up to running the 3D benchmark problem in the second part. Since this benchmark problem involves the solution on a large domain, we will neglect the dispersion effects to decrease the computational time (only for the 2Dld and 3Dld domains since we will compare their numerical results). We remark that for the numerical simulations the 1D and 2D domains are given as 3D grids (e.g., the 1D horizontal domain is represented by a grid of dimensions L×1 m×1 m). Fig. 4 shows the four different domains we consider for the numerical experiments. In all experiments the potential leakage region is located at a distance l from the injection well. The simplest spatial domain for numerical studies is a 1D horizontal domain as shown in Fig. 4a. This domain consists of an injection well, a potential leakage region, and an open boundary. Two 2D horizontal extensions of this domain are given in Fig. 4b and Fig. 4c. The former represents a potential leakage region with a given aperture ω, while the latter represents a potential leakage region of aperture ω and width W . Fig. 4d shows a 2D vertical domain with a height H where the potential leakage region is on the top caprock.

Injection well
Open boundary Storage reservoir Potential leakage x z x y x x y

1D horizontal domain (1Dhd)
We first investigate the dynamic evolution of the model components (i.e., microbes, substrate, urea, biofilm, and calcite) during the injection of phase I on the 1D horizontal domain in Fig. 4a. The different values for the times in the injection of phase I are the following: t I 1 = 20 h, t I 2 = 40 h, t I 3 = 140 h, t I 4 = 160 h, t I 5 = 180 h, t I 6 = 230 h, t I 7 = 250 h, t I 8 = 270 h, and t I 9 = 300 h. These injection times are identical to the ones studied in Tveit et al. (2018). Given the values of injection times, the corresponding injection rate which leads to permeability reduction on the target zone is Q I w =3.8 × 10 −5 m 3 /s. The numerical results are shown in Fig. 5. We observe that after 300 hours all urea is used to produce calcite around the potential leaky zone. The pore space around the target zone is reduced significantly after injection of phase I.

2D horizontal circular domain (2Dhcd)
We now consider the 2D horizontal circular domain in Fig. 4b studied in Tveit et al. (2018). The authors used a sequential approach to solve the mathematical model on a fine triangular grid, which was implemented using FiPY (Guyer et al., 2009). The main purpose of this example is to compare qualitatively with the results in Tveit et al. (2018). Most of the model parameters considered in Tveit et al. (2018) have the same values as in Table 1 or are of the same order of magnitude. The radius of the domain, target location of MICP, initial porosity, and permeability are the same as in the 1D experiment, which also are the mean values in the log-normal distributions in Tveit et al. (2018). We simulate the injection of one phase of MICP using the same injection times as in the previous example (1Dhd). We set the injection rate equal to the rate from the 1Dhd example multiplied by L, i.e., Q I w =2.8 × 10 −3 m 3 /s. Fig. 6 shows the grid, initial permeability, and permeability reduction for our numerical simulations. Here we remark that a scale factor is used to relate the flow rates between the domains. The corresponding scale factor from the 1Dhd to the 2Dhcd is given by the radius L which results in reduction of permeability at the same distance from the injection well in both domains as shown in Figs. 5f and 6b.
Comparing qualitatively the permeability reduction reported in Tveit et al. (2018) to the one seen in Fig. 6b, we observe that both simulations predict the reduction of permeability at the targeted distance from the injection well of 15 m. The different approaches to model the biofilm detachment [e.g., detachment from growing biofilm in Tveit et al. (2018) and detachment due to erosion in this work] and the numerical schemes [sequential in Tveit et al. (2018) and fully implicit in this work] results in the discrepancies [i.e., difference in the value and diffuse rings outside the target area in Tveit et al. (2018)] between the predicted permeability reductions. In addition, the computational cost of the present grid is lower compared to the uniform fine triangular grid studied in Tveit et al. (2018). Thus, in the subsequent experiments, we discretize the spatial domain in a similar manner as shown in Fig. 6a, where the grid around the injection well and the region where the calcite precipitation occurs is fine (order of tens of centimeter), and gradually make it coarser (order of meters) towards the domain boundaries.

2D horizontal rectangular domain (2Dhrd)
We focus on the 2D horizontal domain in Fig. 4c. We set the simulation domain size to 2L = 150 m and W = 20 m. For this example we investigate the reduction of permeability in a potential leakage zone along the width of the aquifer. We simulate the injection of one phase of MICP using the same injection times as in the previous example. We set the injection rate equal to the rate from the 1Dhd example multiplied by W , i.e., Q I w =1.5 × 10 −3 m 3 /s. Fig. 7a shows the permeability reduction after phase I of the injection. We observe that the closer to the lateral boundaries we target the calcite precipitation, the further into the aquifer the components need to be injected, due to the radial flow. Consequently, not all parts of the potential leakage region are covered by one phase of MICP injection. We apply a second phase of injection with the same concentrations and time intervals as phase I but at a higher injection rate Q II w =2Q I w . This higher injection rate for the second phase is selected after testing different injection rates. Fig. 7b shows the permeability reduction after a second phase of injection. Thus, several injection phases at different rates are needed to reduce the permeability inside the potential leakage region. Figure 7: Permeability reduction after injection of (a) phase I and (b) phase II. The potential leakage region is inside the black rectangle.

2D vertical rectangular domain (2Dvrd)
In the next example we study the 2D vertical domain shown in Fig. 4d. We set the simulation domain size following the benchmark study in Ebigbo et al. (2007), that is, L = 500 m and H = 30 m. We simulate the injection of one phase of MICP using the same injection times as in the previous example. We investigate two different injection approaches along the well for the components and water to efficiently get calcite precipitation at the potential leakage region in the caprock. For the first simulation the injection of components and water is only in the first 3 m of the well. For the last simulation we change the water injection to be along the whole height of the well.
Given that the distance to the leakage zone for this reservoir is larger than the one in the previous examples, then a larger injection rate is needed. An injection rate which leads to permeability reduction on the target zone is Q I w = 6 × 10 −3 m 3 /s. Fig. 8 shows the permeability reduction for both injection approaches. We observe that the first injection approach leads to calcite precipitation also along the vertical direction. This is not desired as it could lead to encapsulation of the injection well. With the second injection approach, we accomplish calcite precipitation around the potential leakage region located near the caprock. Then we consider the latter injection approach in the next examples where vertical wells are also simulated.

MICP to seal leakage paths
Diverse reservoir representations where leakage paths are explicitly modeled can be found in literature. In this work, we focus on the two domains shown in Fig. 9. A simple representation of a 2D domain with one leakage path between two aquifers is shown in Fig. 9a. A well-established 3D benchmark for CO 2 leakage is given by the domain in Fig. 9b (Class et al., 2009;Ebigbo et al., 2007).

2D leaky domain (2Dld)
We focus on the domain shown in Fig. 9a. In Ebigbo et al. (2007) the leakage is given as a result of a well which is modeled as a porous medium with higher permeability than the formation. To asses the leakage rate of CO 2 before and after application of MICP, we solve a simple two-phase flow model for CO 2 and water (see Appendix A). We simulate the injection of one phase of MICP using the same injection times and injection rate as in the previous example (2Dvrd). After simulation the observed permeability reduction in the leakage path was minor.
One alternative is to apply additional phases of MICP at different flow rates keeping the same time intervals to aim for the sealing of the leakage path. Since performing simulations on this 2D domain is computationally cheap, we proceed to design an injection strategy such that a minimum number of phases are needed for the sealing of the leakage path. It is beyond the scope of this paper to perform optimization studies. Here we use an ad-hoc approach where we keep the same order of injection of the components, values of concentrations, height of injection along the  Fig. 10 shows the numerical results of this injection strategy on the 2D leaky domain. For a better visualization of the different MICP processes, in Fig. 11a we plot the average value normalized by its maximum value achieve in phase I or II for water velocity, microbial, substrate, and urea concentrations, biofilm and calcite volume fractions, and permeability reduction in the leakage path. We observe a remarkably increase of calcite after injection of urea in phase II which in turn decreases significantly the volume fraction of biofilm. Fig. 11b shows the leakage rate without and after MICP injection of phase I and II after 100 days of CO 2 injection. In the numerical results, we calculate the leakage as the CO 2 flux at the middle of the leaky well (z = 80 m) (Ebigbo et al., 2007). We observe that the leakage rate is practically stopped after two phases.

3D leaky domain (3Dld)
We consider the 3D benchmark reservoir as described in Class et al. (2009) andEbigbo et al. (2007) shown in Fig.   9b. Since the properties of the previous domain (2Dld) are also equal to the ones in the benchmark, we expect to obtain similar results after applying the same injection strategy. Thus, we simulate the injection of two phases of MICP using identical time intervals as in the previous example. We set the injection rate equal to the rate from the 2Dld example multiplied by L, i.e., Q I w =4 m 3 /s. Fig. 12 shows the numerical results after applying phase I and II of MICP. Fig. 13 shows the different MICP processes at the leaky well and the leakage rate before and after MICP treatments. We observe that the dynamics of the processes are similar to the ones plotted for the 2D leaky domain. We also observe that the curve without MICP injection is in good agreement with the ones presented in the benchmark study for CO 2 leakage in Class et al. (2009). Then, as observed in the 2D leaky domain, the leakage stops after applying two phases of MICP treatment.

Discussion
The first part of the numerical examples includes MICP simulations to prevent the formation of leakage paths in an aquifer. We first study the spatial distribution of the diverse MICP variables (namely microbes, substrate, urea, biofilm, and calcite) on a simple 1D horizontal domain. Simulations on this domain are suitable for testing large numbers of injection strategies as it requires the lowest running time. In the 2D horizontal circular domain the potential leaky zone is given by a small area at a given distance from the injection well. However, the fluid injection through the well is in all horizontal directions leading to calcite precipitation around the injection well. To only target the potential leaky zone then the direction of injection could be controlled to decrease the cost of injected components and to not reduce the aquifer storage capability in unnecessary regions. Whereas this is an evident observation, we could not find experimental nor numerical studies involving directional wells in the literature. Thus, including directional wells in the simulators could be useful for MICP studies. An important observation of the simulations on the 2D horizontal rectangular domain is that several phases of injections might be needed to precipitate calcite along the width of the aquifer as a result of the radial flow from the well. The last example in the first part is a 2D vertical rectangular domain. Here we study the injection of water and components along the well. When we inject only at the top of the well, the calcite "encapsulates" the injection well along the vertical direction. We have observed that calcite precipitation only near the caprock is achieved injecting the components only at the top and water along the well. In the numerical studies we designate an arbitrary fixed part of the well for the injection of the components. However, the choice of this control variable is likely to have a significant impact in the simulations (e.g., due to transversal dispersion of the components). Hence, diverse injection heights should be studied for different systems to cut the injection times/cost.
The second part of the numerical examples focus on MICP simulations on domains where the leakage paths are modeled explicitly. We first study the sealing of a leakage path in a caprock between two aquifers on a simple 2D domain. We proceed to find values of time intervals and injection rates to achieve sealing in the leakage path.
The designed injection strategy leads to sealing of the leakage path after two phases of MICP injection. The last numerical experiment is performed on a 3D benchmark leaky well. Since reservoir properties of both domains are equal, we use the same injection strategy as in the 2D domain. After simulations, we observe that the leakage path is blocked successfully. Despite the satisfactory and straightforward application of the injection strategy from the 2D domain to the 3D leaky well, this is ultimately restricted by the simplicity of the 3D domain. For instance, a different injection strategy should be designed for a 3D domain with a fracture across the width of the caprock.
As observed in the simulations on the 2D horizontal rectangular domain in Fig. 7, we expect that several phases of MICP injection are required to significantly reduce the permeability in the leakage path.
The MICP model presented in this paper is built from a simplified description of the underlying processes based on previous publications as described in Section 2. The mathematical model is simple enough to have a few equations (six mass balance equations and six cross coupling constitutive relationships) and input parameters (twenty-one parameters) but still retain the essence of the processes. Extending the model to two-phase flow (water+MICP and CO 2 ) lets us consider additional processes, e.g., dissolution of calcite due to the presence of CO 2 ; however, solving simultaneously all equations would limit the size of the problem as discussed in Cunningham et al. (2019).
As proposed by the authors, this could be solved by applying multi-physics methods, e.g., using analytical solutions for the flow.
Differences between the leakage rate curves on Figs. 11 and 13 arise from different issues. Prior to MICP Figure 12: (a) Initial permeability, (b-c) permeability reductions after phase I and II, and (d-f) amount of CO 2 in the three different scenarios. treatment, the percentage of leakage rate on the 2Dld is nearly 50% while for the 3Dld is lower than 0.25%. The reason for this difference is that on the 2D system all CO 2 injected either flows under or through the leakage path while for the 3D system only a small portion of CO 2 flows under and through the leakage path (since the flow is radial in the 3D system). Nevertheless, we would expect to have a similar curve shape for both systems. For the 3D system we observe a sharp rise of CO 2 leakage when the CO 2 reaches the leakage zone (approximately after 10 days of injection) and then drops (approximately after 50 days of injection). This is attributed to boundary effects (Ebigbo et al., 2007), as the CO 2 leakage is expected to continue increasing slowly after the initial sharp rise as shown in Fig. 11 (Nordbotten et al., 2005).
In addition, we observe a different decrease of CO 2 leakage after one phase of MICP treatment in both systems. From Fig. 11a we observe that the permeability reduction on the leakage path for the 2Dld system is around 25% while for the 3Dld system is around 80% (Fig. 13a). In this work we describe the 2Dld as a general simple system to study MICP in a leakage path. This domain is commonly related to an approximation of a system with a fault across the caprock [e.g., see Tavassoli et al. (2018)]. Thus, numerical simulations in both systems (2Dld and 3Dld domains) lead to different flux rates and velocities in the leakage path which in turn result in different permeability reduction after application of phase I (this explains the differences between velocities, concentrations, volume fractions, and permeability reduction values in Figs. 11a and 13a). For these examples we observe that the flow velocity through the leakage path is slower in the 2Dld domain than in the 3Dld. From the mathematical model (Eq. 8) we observe that the calcite precipitation is very sensitive to urea availability and biofilm formation. As a consequence of the difference between velocities, the whole MICP process is slow down in the 2Dld domain in comparison to the 3Dld domain. In addition, from Fig. 10 we observe calcite precipitation between the injection and leaky well while from Fig. 12 we observe calcite precipitation outside this region, i.e., between the injection well and outer boundary, since the flow field in the 3D domain also transport the different components around the leakage zone. Then the difference between both flow fields have an impact on the MICP processes in the leakage path as observe in Figs. 11a and 13a. Notwithstanding these differences, the leakage path is successfully remediated in both systems (2Dld and 3Dld domains) after application of the second MICP treatment. A comprehensive investigation of boundary and grid effects, in addition to studies in more complex 3D domains, is not feasible using the current implementation of the mathematical model. Here we have used MRST for the testing of the model as implementation on this framework is not difficult; however, there are computation time limitations using this toolbox. Our current plan is to implement this mathematical model using the open porous media (OPM) initiative which allows to perform more computationally challenging simulations (Rasmussen et al., 2020). Hence we will use the OPM simulator to study these effects and MICP in more complex 3D domains.
In this work we only study one injection strategy for the MICP application consisting of several periodical phases. As described above, each phase involves the injection of three components, injection of only water, and periods of no flow, given a total of 18 control variables: three concentrations, six injection rates, and nine period times. To enable other researchers to benefit from our work and test different injection strategies, we have made the code available through a repository. Optimization could be applied to this problem, e.g., considering the costs of injected components, required energy for injection, and total time of injection. Moreover, changing the sequence of components for the subsequent phases could lead to a better injection strategy, e.g., injection of substrate after the first phase for microbial resuscitation. Then, studying different sequence of components plus different heights of the well adds more complexity to the optimization. Though beyond the scope of the current work, further study in this direction is needed for optimization of the MICP technology.

Conclusions
In this work, we present a mathematical model for the MICP technology including the transport of injected components (microbes, substrate, and urea), biofilm formation, and calcite precipitation. We study an injection strategy involving several phases where each phase includes microbes, substrate, and urea with periods of only injection of water and no flow. We conduct diverse numerical experiments for various one-and two-dimensional domains to study the MICP process. Finally, we show the results of applying MICP to a well-established 3D benchmark problem for CO 2 leakage.
Based on this work conclusions are as follows: • Learnings from 1D and 2D studies help us to develop practical injection approaches for 3D simulations.
• Only using a part of the well for injection of components and water leads to calcite precipitation along the whole vertical direction; however, using the top part of the well for injection of components and the rest of the well for water injection leads to calcite precipitation on the top of the aquifer.
• Several phases of injection might be needed for decreasing the permeability in (potential) leakage regions as a result of the radial flow by the injection well.
• Numerical simulations show a stop in CO 2 leakage after MICP treatment.
Appendix A Two-phase flow mathematical model for CO 2 and water We describe the simplified two-phase flow model used in Class et al. (2009) for a benchmark study in problems related to CO 2 storage. CO 2 and water are assumed immiscible and incompressible. We denote by s w the water saturation and s CO2 the CO 2 saturation (s w +s CO2 = 1). We write Darcy's law and the mass conservation equations for each α phase (α = w, CO 2 ) v v v α = − k r,α µ α K(∇p α − ρ α g g g), φ ∂s α ∂t where k r,α is relative permeability. The relative permeabilities are set as a linear function of the saturations (k r,α = s α ), and capillary pressure is neglected (p w = p CO2 ). The model parameters are summarized in Table A-1.