Development of a Reactive Transport Model for Field‐Scale Simulation of Microbially Induced Carbonate Precipitation

Microbially induced carbonate precipitation (MICP) is a promising technique that could be used for soil stabilization, for permeability control in porous and fractured media, for sealing leaky hydrocarbon wells, and for immobilizing contaminants. Many further field trials are required before optimum treatment strategies can be established. These field trials will be costly and time consuming to \carry out and are currently a barrier to transitioning MICP from a lab‐scale process to a practical field‐scale deployable technology. To narrow down the range of potential treatment options into a manageable number, we present a field‐scale reactive transport model of MICP that captures the key processes of bacteria transport and attachment, urea hydrolysis, tractable CaCO3 precipitation, and modification to the porous media in terms of porosity and permeability. The model, named biogroutFoam, is implemented in OpenFOAM, and results are presented for MICP treatment in a planar fracture, three‐dimensional sand media at pore scale, and at continuum scale for an array of nine injection/abstraction wells. Results indicate that it is necessary to model bacterial attachment, that bacterial attachment should be a function of fluid velocity, and that phased injection strategies may lead to the most uniform precipitation in a porous media.


Introduction
Microbially induced carbonate precipitation (MICP) has been proposed for liquefaction reduction (DeJong et al., 2006), soil strengthening (Whiffin et al., 2007), erosion control (Hamdan & Kavazanjian, 2016;Salifu et al., 2016), permeability control (Chu et al., 2013), sealing fractured rock (Minto et al., 2016), reducing well leakage (Phillips et al., 2013), enhanced oil recovery (Wu et al., 2017), and trapping contaminants (Fujita et al., 2004). The well-studied MICP pathway of urea hydrolysis by the soil bacteria Sporosarcina pasteurii produces carbonate and pH conditions that, provided there is a source of calcium, are conducive to calcium carbonate precipitation through the following series of reactions: which, when simplified to a single overall reaction, gives To date, MICP experiments have largely focused on the pore scale and core/column scale, with some notable exceptions for loose sand (Gomez et al., 2015(Gomez et al., , 2016van Paassen et al., 2009) and fractured rock (Cuthbert et al., 2013;Minto et al., 2016;Phillips et al., 2016). Likewise, numerical modeling of MICP has also focused on the pore/core/column scale with only Cuthbert et al. (2013) Past studies have helped shape our understanding of the fundamental processes that occur during MICP, yet there remains a gap between our laboratory scale understanding of these processes and how to implement MICP in a controlled and optimized manner in the field. Given the expense of large-scale field trials and the limited number of variables that can feasibly be investigated, it would appear prudent to use numerical modeling to narrow down the large number of possible injection strategies into a handful of the most promising that can then be experimentally verified.
Current models of MICP have mostly been developed to reproduce individual experimental results and are each well-suited for their intended use. However, none (see Table 1) include all of the following features essential for field-scale modeling: 1. Model domain: a the domain should be two dimensional (2D) or 3D, b large enough to include several injection and abstraction wells, and c allow for heterogeneity in initial conditions (e.g., porosity distribution) that can be informed by geophysical survey.

Injections and flow:
a multiple injection cycles should be possible as this is often necessary to achieve the desired increase in strength or reduction in permeability, b simultaneous injections and abstractions from multiple wells should be possible to allow evaluation of promising injection strategies (van Paassen et al., 2010), c porosity and permeability should be modified by precipitating CaCO 3 so as to capture the change in reagent transport in subsequent cycles (Harbottle et al., 2016;Minto et al., 2016;Wu et al., 2017), d for MICP by bioaugmentation, bacterial transport and attachment should be explicitly modeled and attachment should be a function of velocity (El Mountassir et al., 2014;Minto et al., 2016;Tobler et al., 2014). 3. Reactions: a bacteria should become encapsulated and inactivated by precipitating CaCO 3 (Cuthbert et al., 2012), b ureolysis should be based on Michaelis-Menten kinetics (or a similar first-order function of urea concentration), c finally, while integration with a geochemistry model (such as PHREEQC) has advantages, the incorporation of complex chemical reactions imposes greater restrictions on model tractability than the transport component (Leal et al., 2017) and simplification is necessary for simulating large 2-D and 3-D domains .
The aim of this paper is to present a reactive transport model that can be used to inform injection strategies for field-scale treatment via MICP under real-world conditions. This is done through developing a flexible model that captures the key processes of MICP, validating this model on laboratory experiments, and then simulating a range of real-world injection scenarios to characterize the behavior of MICP treatment and to identify promising injection strategies for future experimental testing. This same modeling approach may be applicable to other areas of water resources research where the interactions between groundwater flow, contaminant transport, and biogeochemical processes are of interest, for example, in groundwater remediation and in soil bioremediation.

The OpenFOAM Framework
OpenFOAM was originally developed as a set of C++ libraries to solve the linked partial differential equations necessary for computational fluid dynamics (Weller et al., 1998). OpenFOAM was chosen as the base model for this work due to its robust handling of the Navier-Stokes equations of fluid flow and the ease with which reactive-transport equations could be added and solved in parallel within the open-source framework. OpenFOAM version 4.1 provided by the OpenFOAM Foundation (www.openfoam.org) was used; however, we note that use of the OpenCFD Ltd. version (www.openfoam.com) and foam-extend fork (www.sourceforge.net/projects/foam-extend) would have been equally appropriate.

Microcontinuum Flow Model: biogroutFoam
Fluid flow in porous media can be solved using the Navier-Stokes equations of continuity (equation (6)) and momentum (equation (7))) modified by the addition of a Darcy term for resistance of the porous media and a Brinkman term to account for the dissipative viscous force (Soulaine & Tchelepi, 2016), a force that only becomes significant at higher flow rates. Equations (6) and (7) are presented for an incompressible fluid neglecting gravity (hence no density-driven flow), allowing the equations to be simplified by dividing all terms by fluid density: where U s is the superficial velocity (m/s), ε is the porosity (−), t is the time (s), p is the pressure (Pa), ν is the kinematic viscosity (m 2 /s), and k is the permeability (m 2 ).
The standard OpenFOAM solver named pisoFoam was adapted to include the Darcy and Brinkman terms to account for flow through porous media. This solver utilizes the Pressure-Implicit with Splitting of Operators (PISO) algorithm to iteratively solve the Navier-Stokes equations. The modified solver has been named biogroutFoam following standard OpenFOAM naming conventions. A full list of model parameters and values used can be found in Table S1 in the supporting information, along with details of geometry and boundary conditions in Tables S2 and S3 and Figure S1 for all simulations presented.

Chemical Species Transport: General Form
The general form of scalar transport is solved with the advection-dispersion equation (equation ( (8)). Dispersivity is represented as a tensor, thus allowing different longitudinal and transverse dispersion coefficients when appropriate (equation (9)):

10.1029/2019WR025153
Water Resources Research ∂εC ∂t þ ∇:U s C−∇:εD C ∇C ¼ 0 (8) where C is the species concentration (kg/m 3 ), t is the time (s), D C is the dispersivity for species C (m 2 /s), α C is the dynamic dispersivity coefficient (m), and d C is the effective diffusion coefficient (m 2 /s). Values for α C and d C were taken from Chrysikopoulos and Katzourakis (2015).

Model Parameterization: Porosity and Permeability
Porosity of the porous media is included in the model as a spatially varying volume scalar field. Initial porosity can be either 1-D, 2-D, or 3-D data arrays representing real measured porosities (e.g., an X-ray CT scan of a rock core or a geophysical survey), or representing an idealized porosity distribution. We use the FIJI version of ImageJ (Schindelin et al., 2015) to both create idealized porosity distributions and work with real porosity data from X-ray CT scans (e.g., Minto et al., 2018). We provide the following custom ImageJ script to translate any image data into a format that can be read by OpenFOAM as a spatially varying scalar field: https://doi.org/10.15129/7f87bcfa-e7aa-4396-ab29-7721b860345e.
Permeability is inferred from porosity using a modified Kozeny style relationship. It should be noted that Kozeny relationships can be subject to error when predicting the permeability distributions (e.g., in rock cores (Krause et al., 2011)), and may not be accurate for media cemented with MICP when the CaCO 3 formation preferentially occurs in the pore throats resulting in a large permeability reduction for little porosity change (Minto et al., 2017). However, given its simplicity, wide use, and the lack of a more suitable alternative, the modified Kozeny relationship was employed to relate permeability with porosity: where k i is the permeability in a given voxel (m 2 ) while k 0 is the base permeability (m 2 ), used as a fitting parameter, k n is a dimensionless scaling parameter often set at 3 (Costa, 2006;Dullien, 1992;Mavko & Nur, 1997), ε i is the porosity in a given voxel (−), and ε c is the critical porosity below which the media is considered impermeable (often termed the percolation threshold).

Solids and Attachment Surface Definition
In order for bacteria to attach there needs to be a solid surface; therefore, cells adjacent to a solid surface need to be identified. For pore and microcontinuum scale modeling, all cells with a porosity less than a threshold porosity are designated a solid cell. Any cell that neighbors one or more solid cells is designated an attachment cell and given the S attach value of 1. Remaining cells have an S attach value of 0.
This behavior can be overridden when performing a simulation that is entirely continuum scale by designating all cells as attachment cells so as to avoid the neighbor search and reduce the computational burden.

Bacterial Transport, Attachment, Decay, and Encapsulation
Bacteria in suspension were modeled as a chemical species with irreversible attachment to surfaces. Attachment is a function of velocity with no attachment above an upper critical velocity, attachment occurring at a maximum user defined rate below a lower critical velocity, and linear scaling of attachment between these limits. Additional bacteria retention occurs due to straining and wedging in crevices and constrictions (e.g., Herzig et al., 1970) which is here considered to be independent of velocity and bacteria concentration. 10.1029/2019WR025153

Water Resources Research
where C bact is the concentration of bacteria in suspension (OD 600 ), K attach is the attachment rate constant (1/s), U attach is the velocity dependence (−), K straining is the straining rate constant (1/s), and|U int |is the magnitude of the interstitial velocity (m/s).
Bacteria attached to a surface decay over time due to cell death and become encapsulated in precipitated CaCO 3 .
where C surfaceBact is the concentration of bacteria attached to a surface (OD 600 ), K decay is the cell death and decay rate constant (1/s), ΔCaCO 3 is the density of CaCO 3 precipitated since the last time step (kg/m 3 ), and K encapsulation is the encapsulation rate constant (kg·m 3 ·s).
Note that for simplicity, the units for C bact and C surfaceBact are presented as OD 600 which is a proxy measure for bacterial cell concentration that becomes nonlinear at higher concentrations. If modeling bacterial injections with a concentration higher than approximately 1 OD 600 , it would be more appropriate to use colony forming units. Given OpenFOAM's unit handling system, OD 600 and colony forming units are directly interchangeable provided that the appropriate specific rate of ureolysis (K u ) is used.

Urea Transport, Ureolysis, and Ammonium Production
Ureolysis is assumed to occur at the same rate for suspended bacteria and surface attached bacteria alike, following Michaelis-Menten kinetics with an ammonium inhibition term that can be turned off. Ammonium production is a function of ureolysis with two moles of ammonium produced for each mole of urea hydrolyzed. Ammonium transport follows equations (8) and (9) while loss of ammonium to the atmosphere as ammonia is not considered in this model.
where C urea is the concentration of urea (M), K u is the specific rate of ureolysis (1/s) analogous to V max in Michaelis-Menten kinetics, K m is the half-saturation rate constant (M), C NH4 is the ammonium concentration, and K NH4 is the optional ammonium inhibition (M).

Ca 2+ Transport, Carbonate Production, and CaCO 3 Precipitation
Carbonate production follows the rate of ureolysis with 1 mole CO 3 2− produced for each mole of urea hydrolyzed. Ca 2+ and CO 3 2− are incorporated into CaCO 3 in a reaction that proceeds at a predetermined rate until either Ca 2+ or CO 3 2− become limiting.
Precipitation can only occur next to a surface. In the pore-scale model, if Ca 2+ and CO 3 2− mix at a distance from a surface then they must wait until transported to the vicinity of a surface.
∂ε:C ca ∂t þ ∇:U s :C ca −∇:ε:D ca ∇C ca þ r precipitation ¼ 0 (18) where C ca is the concentration of calcium ions (M), K p is the precipitation rate constant (1/s) that delays the onset of precipitation, C CO3 is the concentration of produced carbonate (M), and MW CaCO3 is the molecular weight of CaCO 3 (kg·m 3 ·M).

Porosity Reduction and Permeability Alteration
The porosity is updated at each time step to reflect the decrease due to the volume of CaCO 3 precipitated in the previous time step, assuming that the precipitate has a single density that does not vary over space or time. Permeability is then updated to account for the new porosity using the Kozeny relationship of section 2.4.
where ρ CaCO3 is the bulk density of precipitated CaCO 3 .

Batch Reaction Calibration
Batch tests from previous experiments (MacLachlan, 2017) were used to calibrate the model under no-flow conditions and in the absence of porous media for a range of volumes and initial concentrations of each reagent. In each batch test, the ratio of bacteria to urea was varied and calcium was provided in excess so as not to limit the production of CaCO 3 . Ammonium concentration had been measured at regular intervals by ion chromatography; the final mass of CaCO 3 produced was measured by acid dissolution after 96 hr.
To model batch tests, we create a computational domain consisting of a single cell with a volume equal to that of the experimental batch, set all boundaries as impermeable walls, treat bacteria as attached to a surface, and solve the full set of equations. The output of two contrasting batches are shown in Figure 1: batch (a) starts with a relatively high bacteria concentration and low urea concentration while batch (b) starts with less bacteria and a higher urea concentration. All the urea in batch (a) was quickly consumed resulting in a plateau in CaCO 3 and NH 4 + production between one and two days. In batch (b) CaCO 3 and NH 4 + production continued up until 96 hr when the experiment was terminated, but the rate of these reactions slowed down as bacteria became encapsulated in CaCO 3 . In both batches, the rate of calcite precipitation (K p ) was sufficiently high that CO 3 2− did not accumulate.
Ammonium inhibition is considered to be minimal for whole-cell ureolysis kinetics  and so, to avoid interference with the fitting of other parameters, K NH4 was set to a very large number, effectively turning off ammonium inhibition. Likewise, bacterial decay (K decay ) was turned off as there were no experimental data to quantify the rate of decay as a distinct process from encapsulation. An assumed CaCO 3 precipitation rate constant of 1 × 10 −2 s was used. The remaining fitting parameters were ureolytic activity (K u ) and bacterial encapsulation within calcite (K encapsulation ). Unlike most of the model parameters, which have previously been reported in the literature, the rate of bacterial encapsulation was not previously known and was found to have a significant effect on the rate of ureolysis.
K u was not measured in these batch experiments (MacLachlan, 2017) but it is known from experience to typically lie in the region of 5 to 8 mM urea·min·OD for this bacterial strain (S. pasteurii, DSM-33) grown in this way (37-g/L BHI broth with 20-g/L urea, 24-hr incubation at 30°C, 100 rpm) and to vary depending on temperature and age of the bacterial stock. For this reason, K u was used as a fitting parameter for each batch (with batches prepared on the same day having the same value of K u ) while K encapsulation was simultaneously used as a global fitting parameter for all batches. Figure 2 shows the results of fitting the model to the batch test data. The calibrated values of the ureolytic activity, K u , vary from 5 to 8 mM urea·min·OD 600 . The optimum global encapsulation constant that best fit all batches was found to be 12 kg·m 3 ·s; that is, 12 kg/m 3 CaCO 3 would need to be precipitated in one second to completely encapsulate all bacteria currently on a surface in that model cell. The maximum model-calculated rate of precipitation for any of the batch tests was 9 × 10 −4 kg·m 3 ·s; hence, as typical rates of precipitation are far lower than 12 kg·m 3 ·s, only a small proportion of bacteria will therefore be encapsulated in each second. The rate of NH 4 + production decreases over the 96-hr-long experiment due to both encapsulation and reduction in urea availability ( Figure 2).

Flow Component Validation
To determine the validity of the velocity dependent bacterial attachment model (section 2.6), a series of four test models in simple planar fractures were simulated and compared with one of the experimental data sets of El The key finding of this experiment (El Mountassir et al., 2014) was that a feedback mechanism existed between CaCO 3 precipitation and local fluid velocity that resulted in the formation of preferential flow pathways. These pathways were self-organizing and stable because the constant flow rate injection resulted in increasing velocities as fracture aperture diminished due to CaCO 3 precipitation.
To simulate the experiment in Figure 3a, four different numerical models (A-D) were evaluated, each of which used the model boundary conditions and mesh element geometry shown in Figures 3b and 3c. Model A represents a smooth fracture in which bacteria can attach to the entire fracture surface, and in which there is no velocity dependent component to the attachment mechanism. The result of modeling the experiment assuming a smooth fracture (see Figure 4a) is relatively uniform CaCO 3 precipitation with a gradient in precipitate mass from the inlet to the outlet and little precipitation directly surrounding the inlet boundary, as here the CO 3 2− has been transported away before precipitation can occur (determined by flow rate and parameter K p ). This model is reflective of the majority of existing MICP models (e.g., all previous models listed in Table 1).
Model A does not truly represent the experimental conditions in which the bacteria were flocculated prior to injection. Flocculation can be expected to result in clusters of bacteria behaving as discrete particles (bacterial flocs were observed to settle under gravity). Additionally, microscale charge variations in the polycarbonate surface likely lead to preferential attachment at certain points. Model B includes some spatial variation in the initial bacteria distribution by assigning random cells in the fracture as preferential locations for initial bacterial attachment. Once 20% of the fracture aperture in these cells are occupied by CaCO 3 , the neighboring cells are then also considered as attachment surfaces. The results of model B (see Figure 4b) superficially resemble the experimental results; however, there are important differences: (1) flow paths are uniformly distributed with no major channel forming and (2) a gradient in CaCO 3 from inlet to outlet still exists. Visualizing CaCO 3 precipitation in Model B over time (see Movie S1 in the supporting information) shows that modeled precipitation occurred first in the central flow path as this was where most bacteria had attached. This precipitation then acted as an obstruction to flow in subsequent cycles which re-directed flow paths, ultimately resulting in more uniform CaCO 3 precipitation, that is, hydrodynamic feedback that seals preferential flow paths.
Models C and D include the velocity-dependent bacterial attachment proposed in equations (11) and ((12), without and with the inclusion of spatially distributed preferential initial attachment locations, respectively. In Model C, a channel with little CaCO 3 precipitate quickly forms where the velocity is highest. This channel reinforces itself as CaCO 3 precipitated on either side reduces flow through these outer regions and begins to create a wall that limits reagent access into the bulk of the fracture. Where CaCO 3 does precipitate in the channel, it appears to be due to the diffusion of CO 3 2− produced elsewhere.

Water Resources Research
For Model D ( Figure 4d) the velocity-dependent attachment limits the bacterial attachment in the central channel, and so instead of blocking the channel (as in Model B), flow in the central channel is reinforced and the surrounding fracture is filled with CaCO 3 . Hence, Model D can be characterized by formation of a single dominant channel that branches as it approaches the outlet boundary. Further, a larger mass of precipitate is found adjacent to (but outside of) the channel, while there is less precipitate in stagnant areas such as the corners.
Models C and D that include velocity-dependent bacterial attachment better represent the experimentally observed sealing behavior (El Mountassir et al., 2014) in which a single main channel forms and this channel is reinforced by hydrodynamic feedback. The inclusion of velocity-dependent bacterial attachment in the modeling of MICP appears, therefore, to be necessary for reproduction of experimental results.

Conceptual Extension to 3-D Pore-Scale Model
To conceptually verify model performance in 3-D porous media, model simulation results were compared to the spatial distribution of CaCO 3 observed with an X-ray CT scan for an MICP-treated beach sand. The beach sand was treated with five MICP cycles in a radial flow system 40 cm in diameter. Each cycle consisted of one pore volume of 0.5 OD 600 S. pasteurii, a water flush to clear the injection point, a static period of 2 hr to aid bacterial attachment, one pore volume of 1.0 M CaCl 2 /urea cementing solution, a 16-hr static reaction period, and a final water flush before the next cycle.
In Figure 4 of van Paassen et al. (2009) there was visibly less CaCO 3 within 10 cm of the injection point where we calculate based on an assumption of homogeneous spherical flow that the velocity during injections would likely have exceeded 40 cm/hr. We therefore collected a subsample of treated sand from our experiment cut approximately 15 cm from our injection point where flow velocity during treatment would have been in the region of 40 cm/hr.
The subsample was trimmed to a core 6 mm in diameter, and scanned with a Nikon Metrology X-ray CT system. Relevant X-CT scan parameters were 3,141 projections at an angular step of 0.1146°, scan energy of 120 kV, 24 μA, 708-ms exposure, and no prefiltration. The projections were reconstructed into a 32-bit 3-D volume with voxel size of 3.5 μm using Nikon XTekCT software version 4.3.4, processed in the FIJI version of ImageJ, then visualized in 3D using ParaView (Ahrens et al., 2005). The experimental X-ray CT scan results are shown in Figure 5a (see Movie S2), in which the precipitated CaCO 3 is clearly visible within the pore network. Model simulation results of the MICP experiment ( Figure 5a) predict that in much of the modeled region, velocity would exceed the attachment threshold ( Figure 5b) and attachment would be limited to stagnant areas between sand grain contact points (Figure 5c). Final precipitation results (Figure 5d) show that the CaCO 3 is mostly located at grain contact points and closely resembles the X-CT data. For comparison, Figure 5e presents a model simulation without the velocity-dependent attachment; this shows nearuniform CaCO 3 precipitation on all available surfaces and does not compare well to the experimental data in Figure 5a.
Conceptually, from comparison with the observed spatial distribution of calcium carbonate precipitate, it would appear that the model including velocity-based attachment better captures the pore-scale processes of bacteria transport, attachment, and CaCO 3 precipitation. Further sensitivity analysis and parameter optimization were not carried out in this pore-scale model for two reasons: (1) pore-scale experimental data are not available to investigate the relationship between temporally evolving precipitation patterns and pore-scale velocity distribution and (2) these pore-scale simulations are extremely time consuming to run (see Table S2 for runtimes). These long runtimes arise in the pore-scale model because mesh cells must be small enough to depict pore throat geometry and the numerical time step must be small enough to satisfy the Courant stability criteria for reliable transport of the chemical species. Mesh cells in continuum-scale models can be far larger; hence, larger time steps can be used without violating the Courant criteria.

Field-Scale Simulation
To test differing field-scale injection strategies, simulations were performed in an idealized array of wells (shown in Figure 6) with an assumed initial porosity distribution, such as might be collected from a geophysical survey of a site. Individual boundary conditions were set for each well, allowing independent control of injection and abstraction and multiple injection cycles with static periods for bacterial attachment and CaCO 3 precipitation. The four sides of the bounding box were set as pressure-controlled boundaries allowing groundwater to both exit and enter the modeled domain, depending on the local pressure gradient. Pressure, flow rate, and reagent flux were monitored individually at all boundaries to determine what entered and exited the model during treatment. Average attached bacteria concentration, mass of CaCO 3 precipitate, and porosity were output for each time step, throughout the model domain.
Wells in the test array ( Figure 6) were spaced 3 m apart and the total number of cells in the modeled domain (18 × 18 m) was approximately 15,000. Cells were polyhedral with a size ranging from 0.1 m near the wells up to a maximum of 0.5 m at the outer edges of the bounding box.
During simulation, the Courant number was set to <1 to ensure stability. As the Courant number is a function of velocity and cell size, model run duration was largely determined by the single ring of cells that define the wells, as these were both the smallest in size and had the highest fluid velocity (due to radial injection).

Effect of Injection Rate
As bacteria attachment has been shown to be significantly affected by fluid velocity, injection rates were varied to investigate its effect on the final spatial distribution of the precipitated calcium carbonate. It is expected that if velocities are too low adjacent to the wellbore, then the injection well will become blocked. Volumetric injection rates were selected to inject one pore volume in 1, 2, or 3 hr, designated as high, middle, and low injection rates, respectively. The desired area to be treated extended 3 m from the injection wells (shown as a black box in Figure 7) and this region, together with the domain average porosity of 0.23, was used to define one pore volume. The central well was used for fluid extraction by applying a constant pressure at the well, which was fixed at atmospheric pressure for all three flow rates.

Water Resources Research
For each of the three injection rates, 10 MICP injection cycles were simulated. Each injection cycle consisted of one pore volume of bacteria followed by a 30-min static period for bacterial attachment, then one pore volume of cementing solution followed by a 4-hr static reaction period. In all simulations, regardless of the injection rates, the same quantities of reagents and bacteria were injected. Figure 7 shows the posttreatment distribution of CaCO 3 at each flow rate: lower flow rates result in greater CaCO 3 concentrations close to the well, while higher flow rates better distribute the CaCO 3 within the target treatment area. Additionally, Table 2 shows that the higher flow rate treatment took less time, precipitated a greater total mass of CaCO 3 (albeit with some outside the target treatment area), and made more efficient use of the urea injected. The reason efficiency increased was that the bacteria were distributed over a larger area; hence, the ureolytic activity at any one location was decreased. This decrease in ureolytic activity results in a decreased rate of CaCO 3 precipitation, and so the rate of bacterial encapsulation decreases (see equation (14)). As a consequence, the same initial amount of bacteria was able to hydrolyze urea for longer.
An animation of the model shown in Figure 7a is included in Movie S3. This animation shows the following fields and their changes over the 10 cycle treatment duration: porosity, permeability, velocity and pressure, bacterial attachment coefficient, bacteria in suspension, bacteria attached to a surface, urea and calcium concentration, ammonium produced, and CaCO 3 precipitated.

Ammonium Recovery
Ammonium is considered a ground and surface water pollutant and recovery of produced NH 4 + may be a regulatory requirement before MICP is permitted, as was the case in the field trial of Cuthbert et al. (2013). To investigate the potential for ammonium recovery, abstraction rates from the central well were increased (controlled via changes to the pressure head in the central well) and cumulative transport of NH 4 + across all boundaries was monitored. Lowering the pressure head resulted in more fluid leaving the model domain through the central well (from which it could hypothetically be collected and treated) compared with leaving via the four-side boundaries (which could result in a release of NH 4 + to the   Table 2), but slightly reduced the efficiency of the MICP process due to the abstraction of reagents before they had an opportunity to react. Treatment extent (Figure 8) was also slightly reduced, but homogeneity was increased largely due to smaller stagnant areas between the injection wells.

Phased Injection Strategy
In all of the field-scale injection strategies presented above, it is apparent that little CaCO 3 is precipitated in the static regions between the wells, denoted by the flow glyphs in Figure 6, right. These static regions arise due to the simultaneous injection through multiple wells creating areas with no pressure gradient and hence with no flow and no possibility to deliver the reagents necessary for CaCO 3 precipitation. In industry, soil stabilization or the creation of hydraulic barriers based on the injection of cement and chemical grouts often proceeds by injecting into one well at a time, with a pump injecting reagents at a controlled flow rate and with a safe pressure cutoff above which the pump will shut off. In order to create a more uniform treatment, it may be preferable to use this approach, which results in few regions of no-flow. However, this would take far longer to treat the same volume compared with the simultaneous injection.
An alternative strategy would be to split the injection into phases, for example, the three phases shown in Figure 9 where treatment begins with simultaneous injection through the outer four wells (1,3,7,9) with abstraction from the central well (5), then treatment of the remaining four outer wells (2,4,6,8), again with abstraction from the central well, then treatment solely through the central well (5). This may have the advantage of more uniform treatment than simultaneous injection, while reducing treatment time in comparison to entirely sequential injection.
A simulation was conducted for MICP treatment injected at a constant flow rate of 8.93 × 10 −4 m 3 /s per well (Figure 9, top row). Each well received the same volume of reagents, as the flow rate is constant, and injection occurred in three phases as discussed. Due to the heterogeneous initial porosity distribution, variations in the local velocity at each well result in corresponding variations in the spatial distribution of the precipitated CaCO 3 in each of the injection phases. The pressure rose in each well as the amount of CaCO 3 increased with each treatment cycle (Figure 9, top right). Due to the low viscosity of MICP fluids well pressures remain far lower than would occur if injecting cement or chemical grouts (as used in conventional ground improvement techniques), and the pressure response is more comparable to a bioremediation treatment where nutrient solutions are injected and growth of a biofilm is stimulated. The final distribution of CaCO 3 suggests that this phased approach results in a more even treatment process than simultaneous injection, an observation which may also hold true for bioremediation through stimulation of native microorganisms. Total treatment time for the phased approach was longer; however, urea utilization efficiency was also higher (Table 2).
To make efficient use of time on site and minimize equipment costs, it has previously been proposed that a single pressure-controlled pump be used to inject into multiple wells during commercial MICP treatment (e.g., Whiffin, 2004, Figure 7.2); hence, a pressure-controlled system was also simulated. The pressurecontrolled model (Figure 9, bottom row) started with a fixed pressure of 3.91 kPa and delivered on average 58.44 m 3 of reagents (bacteria + cementing solution) to each well, making it comparable to the 58.14 m 3 delivered to each well in the flow rate-controlled model. However, due to the imposition of pressure control, the volume injected per well varied from 49.99 m 3 for Well 1 located in a low-porosity region to 110.948 m 3 for Well 3 in a high-porosity region. The flow rate per cycle also decreases markedly, particularly in Well 3, as the precipitation of CaCO 3 decreases the in situ permeability. The permeability of Well 1, however, remained unchanged as the injection rate here was severely limited from the outset by the low--permeability region around the well.

Attachment Mechanisms
Velocity-dependent bacterial attachment has not previously been included in models of the MICP process, yet we find it essential to reproduce both the hydrodynamic coupling observed by El Mountassir et al. (2014) for fractures and the pore-scale distribution of CaCO 3 precipitate observed in experimental micro X-ray CT data. When moving to the 3-D microcontinuum-scale porous media models, it was also necessary to add an additional nonvelocity-dependent removal mechanism that borrows from colloid filtration theory and accounts for straining and wedging between media grains.
At the continuum scale, the inclusion of straining and wedging likely also accounts for bacteria attachment in stagnant zones within the pore structure that are not well represented by a cell average velocity. This suggests some scale dependence of the attachment parameters which, together with the well-established porosity/permeability scale dependence (Ehrenberg, 2007;Heap & Kennedy, 2015;Whitaker & Smart, 2000), may influence the choice of model parameters.
All model parameters relating to bacterial attachment were inferred from a limited number of experiments in the literature, some of which were not intended for this purpose. Looking to the future, experiments to explicitly parameterize S. pasteurii attachment mechanisms in a range of soil types, ionic strengths, and flow conditions are required.

Encapsulation of Bacteria
Feedback between the rate of precipitation and bacteria encapsulation is one of the main restraints on the MICP process and determines both how much CaCO 3 can be precipitated for a given amount of bacteria, the constant pressure model. Wells were selected to have contrasting initial porosities: Well 1 is located in a region of low initial porosity, while Well 3 is in a highporosity region (see Figure 6). and how long it will take. For S. pasteurii in batch tests without porous media, the encapsulation constant was estimated to be 12 kg·m 3 ·s and this was taken as an indicative value for use in the field-scale simulations. Encapsulation was treated as a function of the amount of bacteria present on a surface and of the rate of CaCO 3 precipitation. This is likely a simplification of a complex process, as the function does not consider the arrangement of bacteria on the surface (e.g. uniformly distributed, aggregated, or in a biofilm), nor the morphology of the precipitating CaCO 3 (e.g., large rhombohedral crystals with low internal porosity, or many small crystals with a higher internal porosity).
Additionally, it is likely that during the 96-hr batch test there was natural cell death and possibly some growth (although no growth media were provided, hence this would be minimal). What was measured as the rate of encapsulation is actually the balance between encapsulation, death, and growth.

Injection Strategy
In radial flow systems the combined attachment processes give rise to a CaCO 3 distribution resembling that observed by van Paassen et al. (2009) with less precipitate immediately surrounding the injection point. Based on this, we propose that injection rate can be used to limit attachment near the well in both rock fractures (as was observed by Minto et al. (2016)) and in porous media. High flow rates during the injection of bacteria may allow treatment to be spread over a larger volume from a single well before blocking of the injection point occurs. This creates the opportunity for more injection cycles and hence greater soil improvement.
Simultaneous injection through multiple wells was found to greatly reduce the time required to treat a large area; however, it resulted in areas of low treatment between wells which may be undesirable. Phased injection strategies resulted in more homogeneous treatment. This strategy differs to that employed when using as conventional cement and chemical grouts for soil stabilization or the creation of hydraulic barriers as they typically create near-impermeable regions after the first injection, whereas the phased strategy relies on being able to deliver reagents through a region that has already been partially treated.
Both pressure-controlled and flow rate-controlled injection strategies appear to be valid options. Pressure-controlled injection may be simpler when treating near-surface zones where there is a risk of ground heave as the starting pressure will never be exceeded. With pressure-controlled injection there is also the advantage of delivering more reagents to areas with a high initial porosity around the well, potentially resulting in a more uniform final strength across the site. However, the volume injected is different for each well and treatment volume per cycle will also change, potentially making the handling and preparation of reagents on site more complex.

Conclusions
We present the biogroutFoam model of microbially induced carbonate precipitation suitable for field-scale modeling. Advances in this model are the inclusion of bacterial attachment mechanisms, flexible injection/abstraction boundaries necessary for modeling complex MICP treatment strategies, and a simple tool for incorporating real-world spatially varying ground properties such as those that can be derived from a geophysical survey. The model has been conceptually validated against experimental data and used to narrow down the range of possible injection strategies to determine those most promising for implementation at a field scale.
In addition to modeling MICP at field scale, this tool may be of interest for simulating reactive transport processes in lab-scale experiments where 3-D distributions of porosity and permeability can be measured by X-CT and other imaging tools. Finally, the modeling approach taken with biogroutFoam could be adapted to other biogeochemical processes, particularly those of interest in groundwater remediation and soil bioremediation.