A comparative study of using two numerical strategies to simulate the biochemical processes in microbially induced calcite precipitation

Dianlei Feng, Xuerui Wang, Uo Nakenhorst, Xuming Zhang, Pengzhi Pan

a Department of Hydraulic Engineering, Tongji University, Shanghai, 200092, China

b Institute of Fluid Mechanics and Environmental Physics in Civil Engineering, Leibniz University Hannover, Hannover, 30167, Germany

c Institute of Mechanics and Computational Mechanics, Leibniz University Hannover, Hannover, 30167, Germany

d State Key Laboratory of Geomechanics and Geotechnical Engineering,Institute of Rock and Soil Mechanics,Chinese Academy of Sciences,Wuhan,430071,China

Keywords:Microbially induced calcite precipitation(MICP)Advection-diffusion-reaction (ADR)equation Time discontinuous Galerkin (TDG) method OpenGeoSys-PHREEQC

ABSTRACT In this study, we carried out a comparative study of two different numerical strategies for the modeling of the biogeochemical processes in microbially induced calcite precipitation(MICP)process.A simplified MICP model was used, which is based on the mass transport theory.Two numerical strategies, namely the operator splitting (OS) and the global implicit (GI) strategies, were adopted to solve the coupled reactive mass transport problems.These two strategies were compared in the aspects of numerical accuracy,convergence property and computational efficiency by solving the presented MICP model.To look more into the details of the model,sensitivity analysis of some important modeling parameters was also carried out in this paper.

1.Introduction

Recently, microbially induced calcite precipitation (MICP) has received increasing interest in geoenvironmental engineering,as this bio-logically controlled mineralization can significantly reduce the permeability and improve the mechanical properties of soil.Potentially,MICP could be applied for the remediation of hydraulic system and reinforcement of soft underground.In some studies(Ivanov and Chu, 2008; Phillips et al., 2013; DeJong et al., 2014; Umar et al.,2016; Krajewska, 2018), the potential engineering applications of MICP are presented,especially the role of MICP in soil is highlighted.Compared to the traditional techniques to improve the engineering properties of soil, for instance chemical grouting or mechanical consolidation,the utilization of thisinnovative bio-induced method is expected to be more eco-friendly.Besides,the calcite-cemented soil exhibits much stable chemical and mechanical properties compared to the chemical-cemented soil.Although,in general the overall MICP processes have received significant attentions, some involved mechanisms, e.g.bacterial influences on reaction kinetics, and coupling effects,are currently not well understood.Thus,the control and prediction of MICP treatments remain as challenges.To get a better insight into the MICP involved processes,in the recent years,many investigations have been carried out based on either experiments or numerical modeling.A certain understanding has been gained through those investigations.To date,among most MICP numerical models,the biochemical issues involved in MICP have been recognized as one of the most important processes, which are significantly complex and often strongly coupled (Barkouki et al.,2011; van Wijngaarden et al., 2011, 2013, 2016; Ebigbo et al., 2012;Fauriel and Laloui,2012;Cuthbert et al.,2013;Hommel et al.,2015;Cunningham et al.,2019;Wang and Nackenhorst,2019,2020).Typically,such processes include the growth and decay of biomass,urea hydrolysis, hydrolysis of hydrous complex and calcite precipitation.Modeling such complex systems requires to consider, for instance,how the urea hydrolysis or calcite precipitation is affected by the temperature, the pH value, the bacteria concentration and so forth.Furthermore,one needs to consider the porosity reduction,thus the dynamic heterogeneous distribution of the permeability in the porous media caused by calcite precipitation.Despite the complexity of the physical processes,generally in the numerical modeling,MICP induced biochemical issues in soil are described by mass transport processes coupled with simplified reaction model in the porousmedia,inwhich a set of advection-diffusion-reaction(ADR)equations have to be solved.Apart from the challenges raised during the development of the MICP models,solving the mathematical models with proper numerical methods/strategies is another key issue for achieving efficient, accurate and stable simulation results.Recently,the most commonly used solution strategies for subsurface environmental reactive mass transport problems are the direct approach based on global implicit (GI) method and the sequential approach based on operator splitting (OS) method (Steefel et al., 2015).The sequential approach can be further categorized in sequential iterative approach (SIA) and sequential non-iterative approach (SNIA), of which the lattercategory ismore commonly used(Steefelet al.,2015).In general,OS strategy has the advantages of allowing the possibility of using different methods to solve chemical and transport equations,easy implementation (Carrayrou et al., 2004; Fahs et al.,2008),and providingbetteropportunities for parallelcomputation(Hundsdorfer and Verwer,1995).Especially,due to its flexible treatment of the reaction term, the stiff problem involved can be easily handled.In an early study,based on a theoretical analysis,Yeh and Tripathi(1989)concluded that in two- (2D) and three-dimensional (3D) cases, using the direct GI approach may lead to excessive memory allocation and computational time.Then the OS strategy becomes more attractive in the simulation of reactive mass transport and has been investigated and utilized in many studies(Miller and Rabideau,1993;Shao et al.,2009;Centler et al.,2010;He et al.,2015).

However,the splitting error introduced in the OS method often reduces the numerical accuracy of the scheme(LeVeque,2002)and can even result in unacceptable results(Zalesak,1979;Valocchi and Malmstead,1992).Unlike OS strategy,GI strategy is expected to be more robust.With regard to the computational efficiency, some recent studies (Saaltink et al., 2001; Fahs et al., 2008) have shown that GI strategies perform even better in some cases, e.g.coupled chemical system of high nonlinearity and strong interaction between solid phase and transport processes.Besides, in Saaltink et al.(2001), it was reported that using OS method gives rise to problems associated with large reaction rates and effective porosity.By using an implicit A-stable scheme, the stiff problem induced by the reaction term can also be handled with the GI method in combination with large time steps.However, in many cases,the time discretization in such GI method can only be of 1st order, e.g.with the backward Euler scheme.Only few GI methods can fulfill the properties of being implicit, A-stable and high order in time.Among these methods, the time-discontinuous Galerkin(TDG)method has attracted many attentions in the recent years for various applications (Tezduyar et al., 2010; Sapotnick and Nackenhorst, 2012; Yang et al., 2012; Feng et al., 2017, 2019).

Applying the discontinuous Galerkin method for time discretization,namely the TDG method,allows us to easily reach the desired high time order and it has been proved to be A-stable at least for the 3rd order scheme.The method isalsoimplicit,which allows a Courant number larger than 1(Sapotnick and Nackenhorst,2012).The drawback of the TDG method is the increase of the degree of freedoms because the unknowns span a space-time domain.However, the space-time formulation involved in the TDG method shows good potentials for parallelization and implementation (Gander and Neumüller, 2016).With regard to the solution of the reactive mass transport problems in MICP,both GI and OS strategies are used in the literature.For instance,Ebigbo et al.(2012),Hommel et al.(2015),and Wang and Nackenhorst(2020)solved reactions and mass transport together using GI method; Fauriel and Laloui (2012) and van Wijngaarden et al.(2016) applied the OS strategies.However, the performance of these two solution strategies for modeling the MICP under different conditions is rarely investigated.As in practice,MICP could be applied in different scales and under flow conditions, the performances of using these two solution strategies for modeling different MICP cases are worth to be discussed.As mentioned in Hommel et al.(2016), there is a need to optimize the solution strategies to improve the computational efficiency of MICP simulation,especially in large scale modeling.This paper aims to compare two widely used numerical strategies for solving theMICPmodels.For this purpose,we developed our study based on a one-dimensional(1D)simplified two-component MICP model, where urea and calcite are considered.In this model,biomass is assumed to be constant in both space and time.The MICP model that we applied in this studyconsists of a coupled partial differential equation(PDE)-ordinary differential equation (ODE) system.The PDE in the model is a time-dependent advection-diffusion-reaction (ADR) equation with nonlinear reaction term.Both GI and OS standard SNIA approaches are adopted to solve the ADR equation.For time discretization,the DG method and Crank-Nicolson method are used in the GI approach and the OS-SNIA strategy,respectively.For spatial discretization,the standard Galerkin method is applied in both strategies.In this paper,before modeling the MICP,we firstly study a linear time-dependent ADR equation for the verification of these two numerical strategies.Thereafter, the MICP simulation is carried out.The model parameters are calibrated using the data published in Fauriel and Laloui (2012).Based on the calibrated model,the two numerical strategies are compared in terms of convergence behavior and computational efficiency.One of the main crucial points in solving the mass transport equations numerically is the control of numerical errors,namely the numerical dissipation and dispersion.Especially, the numerical dispersion usually results in nonphysical oscillations when solving the mass transport problems with large grid Péclet number.For this reason, we also compare the presented numerical strategies in dealing with the numerical oscillations.Moreover, sensitivity analysis is carried out to consider the possible MICP applications under different operation conditions.

2.Mathematical model

2.1.Chemical reactions

In the most MICP applications, the ureolytic bacteria, typically Sporosarcina pasteurii, are injected into the soil.They are usually supplied with urea ((NH2)2CO) and calcium (Ca2+).The enzyme urease expressed by S.pasteurii catalyzes the hydrolysis of urea into ammonium() and carbonate () ions:

In the presence of calcium, the calcium carbonate (CaCO3) can precipitate:

Combining Eqs.(1)and(2),the overall MICP reaction is given by

2.2.Reactive mass transport

For the considered MICP problem,we assume that the calcium is sufficient in the system.According to the overall MICP reaction equation (Eq.(3)), urea is expected to be the most important component driving the reaction kinetics.Thus, we only consider the transport process of urea, and the calcite is assumed to be immobile.In the considered system, the porous medium is assumed to be fully saturated with pore water, and the biomass isassumed to be constant in space and time.The injected urea mass is assumed to be totally dissolved in the pore fluid,and the sorption of the urea mass onto the solid matrix is neglected.The transport of the urea mass in the pore fluid is described by the ADR equation:

where Cureais the concentration of urea(mol/L);φ is the porosity;v is the advective velocity(m/s),which is assumed to be constant in this study;and D is the diffusion-dispersion tensor(m2/s),which is composed of diffusion(Ddiff)and dispersion(Ddisp)parts,of which the latter depends on the advective velocity.In the 1D case, Ddispcan be expressed as

where αLis the dispersion coefficient(m),and I is the unit matrix.In Eq.(4),qureais a sink term(mol/(L s))caused by the reaction.Similar to many previous studies(van Wijngaarden et al.,2011,2013,2016;Fauriel and Laloui, 2012; Wang and Nackenhorst, 2019), the sink term is determined by an overall kinetic rate kurea(mol/(L s)) as follows:

Further, the Michaelis-Menten type kinetics is adopted for the overall reaction rate:

where U is the constant (mol/(L s)) for the maximum urease rate;and Kmis the half saturation constant(mol/L),at which the reaction rate reduces to half of the maximum value.

Substituting the reaction term into Eq.(4), the final ADR equation for the urea transport is derived:

As aforementioned, the produced calcite is assumed to be immobile.Thus, the increase of calcite concentration (mol/L) is calculated locally based on the overall kinetic rate:

The final system equations are composed of an ADR equation(Eq.(8))and an ODE(Eq.(9)).By applying the corresponding initial and boundary conditions, these equations are solved using finite element method (FEM).

2.3.Dimensionless form of the governing equations

Introducing the following dimensionless variables:

the dimensionless form of the governing equations corresponding to Eqs.(8) and (9) can be derived as

In Eq.(10), all these introduced dimensionless variables are indexed with“*”and the reference(or the characteristic)values are marked with the index “0”.For example, l0and t0denote the characteristic length and time,respectively;and C0is the reference concentration value.

For 1D problems, the dimensionless parameters in Eq.(10)reduce to scalars, e.g.α*= α*(v = v) and β*= β*(D = D)∙The variables α*is usually known as the Courant number and β*is the Fourier number or Neumann number.In this case, the 1D dimensionless equation corresponding to Eq.(11a) can be written in an alternative manner as

where X is the dimensionless coordinate;and S,PeIIand DaIdenote the Strouhal number, the second Péclet number, and the first Damköhler number, respectively.These three dimensionless parameters are defined as

3.Numerical methods

We adopt two types of numerical methods, namely the OS and GI methods,to solve the coupled PDE-ODE problem.Since the main numerical challenges in this model arise from solving the PDE(Eqs.(13a)-(c)), we only present the numerical schemes of these two methods for solving the model equation:

where C denotes the variable in the spatial and temporal domains Ω and T, respectively; and R(C) represents the nonlinear reaction term.

In this study,the solution of the time-dependent ADR equation via the OS strategy is based on the finite element code OpenGeoSys(OGS) (Kolditz et al.,2012) and the simulator PHREEQC (Parkhurst and Appelo,1999;Jang et al.,2018).The transport process is solved in OGS, while the reaction process is solved in PHREEQC.These solutions in the two numerical codes are coupled with OS strategy,which is conducted in two fractional steps.Such a two-step strategy is the simplest OS method,however,it is still widely used in many engineering applications.We also apply the TDG method,which is a GI method, for solving the governing equation in this study.Implicit schemes with high-order time accuracy can be developed with the TDG method in a straightforward manner.The GI strategy is implemented in our in-house Matlab code.

3.1.OS method

By using a two-step OS method, the numerical solution of Eq.(14) can be obtained by solving the following two subproblems in sequence:

Problem I:

Problem II:

We use the method of lines for solving Problem I where linear elements are applied for the spatial discretization and the Crank-Nicoson method is used for the time discretization.Discretizing the unknown variable C in space yields

By using the FEM, the spatial unknown variable is numerically approximated by Chwhich is a linear combination of the nodal variables c={ci}and the corresponding shape functions φ={φi}.In Eq.(16), the index i varies from 1 to the total spatial degree of freedoms Nn.

Using Eq.(16) for the discretization of the weak form of Eq.(15a), we have

where ψ = {ψi} is the test function, which is chosen in the same space as the shape function φ.

The discretization eventually ends up with an ordinary matrix system:

In this paper, we apply a Neumann boundary condition(D∇CH)·n = 0(see numerical studies in later sections).Therefore,Eq.(18) reduces to

Using the Crank-Nicolson method for the time discretization of Eq.(21) yields a fully discretized system of equations:

where Δt = tn+1-tndenotes the time step between two time points tnand tn+1∙

3.2.GI method

We also apply the TDG method, which is a GI method, for solving the governing equations of the MICP model.Different from the OS method,the time-dependent ADR equation can be solved in one step within a space-time domain.As shown in Fig.1, a timedependent 2D problem defined in Ω can be numerically solved within a 3D space-time slice Qn= Ω×Tn, which is a union of all time/space elements of the slice Qn= Ueand Tn= [tn-1,tn] In other words,the temporal degree of freedoms are treated the same as the spatial degree of freedoms.

Using the discontinuous Galerkin method for the time discretization, the weak form of the governing equation (Eq.(14))within Qnreads

where M = {Mij} and K = {Kij} represent the mass and system matrices, respectively, which can be written in detail as

In the above equations, we define

Fig.1.A TDG space-time slice and element(corresponding to a 3rd time order scheme)for a 2D (in space) problem (Feng et al., 2019).

It should be noticed that the reaction term in Eq.(14) is nonlinear,hence we use the Newton-Raphson method for solving the PDE by introducing an iteration procedure:

where β=0,1,…is an iteration index.Substituting Eq.(25)into the weak form of the governing equation(Eq.(23))yields the linearized space-time weak form.Interpolating the unknown variable field Chover the space-time slice Qnyields

where Ntand Nndenote the numbers of degree of freedoms in time and space in Qn, respectively.Different from the OS method, the nodal variables c = {cjk} has a dimension not only in space(denoted by the index “k”) but also in time (denoted by the index“j”), and the corresponding spatial and temporal shape functions are denoted asand, respectively.Using Eq.(26) to discretize the linearized space-time weak form of Eq.(14)eventually yields a fully discretized space-time system of equations as

where ⊗denotes the Kronecker product, and vec(·) is the vectorization transformation of a matrix.The mass matrix M is defined in the same form as in Eq.(19a),however,the shape functions in the equation are replaced withand, respectively.The tangent system matrixand the system matrixread

Fig.2.Temporal approximation with linear shape function (of 3rd order at the discontinuous nodes) in the TDG method.

4.Comparison of the two numerical strategies

In this section, we compare the OS and GI methods in two numerical cases.The first one is a simple 1D reactive mass transport problem with constant reaction rate.Since the analytical solution of such simple problem is available, this case is also used for verification of the numerical code for both numerical strategies.The second one is the MICP modeling which is governed by the equations described in Section 2.2.

In a typical MICP test in the laboratory,the injection rate ranges from 0.14 mm/s to 0.85 mm/s (Martinez et al., 2013; van Wijngaarden et al., 2016; Dadda et al., 2017; Tan et al., 2017).The diffusion coefficient D is usually determined in the order of 10-9m2(Cussler and Cussler,2009)and the dispersion coefficient αLvaries in the range of 10-4-10-2m in the laboratory scale(Frippiat et al.,2008).The measured maximum urease activity is found in the range of 0.0024-0.009 mol/(m3s)(Wang and Nackenhorst,2019).According to Eq.(13c),depending on the model settings,the values of the above mentioned dimensionless numbers could be quite different in the modeled case considering different MICP treatment conditions.Thus,by using different strategies to solve the different MICP application cases, one would expect, on one hand, different computational performances, and on the other hand, different computational results.These two aspects are discussed in Sections 4.2.5 and 5, respectively.

4.1.1D problem with constant reaction rate

4.1.1.Problem description

A horizontal column is continuously flushed with constant injection rate from the left side.Such a problem can be modeled as a 1D transient ADR equation:

where r is the constant reaction rate.By adopting the initial condition as illustrated in Fig.3,the analytical solution of Eq.(30)reads

Fig.3.Initial distribution of the reactant concentration, and the concentration distribution after 21,000 s of transport according to the numerical solution and that computed using both the OS and GI methods.

where erf is the error function, and the diffusion-dispersion coefficient D is defined as D = Ddiff+αLv∙

This 1D verification problem has been solved using both the OS and the 3rd order TDG method GI-TDG(3).The computational domain is uniformly discretized into 200 elements.Totally 21,000 s of transport is simulated with uniform time step of Δt=100 s.The modeling parameters are listed in Table 1.In this study, the grid Péclet number and Courant number are 1.11 and 0.316,respectively.

Table 1Parameters and setup used for the 1D constant reaction case.

4.1.2.Convergence analysis

The numerical and the corresponding analytical solutions of the problem at the end time point of the simulation are plotted in Fig.3.The numerical results computed by both numerical strategies represent the analytical solution well.To quantify the performance of both numerical methods, we use the following error as a measure of the numerical accuracy:

where Ndofdenotes the number of degrees of freedom;and Cjandare the analytical and numerical solutions at the jth node,respectively.The numerical error corresponding to the OS method is 2.3162 ×10-4, whereas the GI-TDG(3) method results are in an error of 1.6041 ×10-6.

We further present a convergence analysis for the solution of the verification linear problem with these two numerical methods.The Courant number is fixed to a constant value as Cr=0.316 in the convergence analysis, and the grid size Δx (dimensionless) varies from 6.25 ×10-4to 10-1accordingly.

On one hand, in the OS strategy, we use the Crank-Nicoson method (2nd order) for the discretization of Problem I (Eq.(15a))and an ODE solver implemented in PHREEQC for solving Problem II(Eq.(15b)).By default,a third order Runge-Kutta method is used in PHREEQC.On the other hand,both the 1st and 3rd time order TDG methods are investigated in the GI strategy.

The results of the convergence analysis are shown in Fig.4.It is demonstrated that both numerical strategies converge when solving this verification problem.However, the OS strategy can only reach a 1st order convergence rate.On the other hand, the GITDG(3) method shows a 2nd order convergence rate.

4.2.MICP modeling

4.2.1.Model setup

In this section, we further compare the numerical solutions of the MICP model using the OS and the GI-TDG(3) strategies.Unlike the model presented in Section 4.1,in the MICP model,the reaction rate is a nonlinear function of the urea concentration(Eq.(7)).The numerical investigations in this section are carried out in a horizontal column of 50 cm long (Fig.5).Urea with a constant concentration of 0.5 mol/L is injected from the left side of the domain.Since the length of the column is much larger than its diameter,we model the whole process as a 1D problem.The computational domain has been uniformly discretized into 200 elements in space,and totally 5400 s of MICP treatment is simulated with a time step of Δt=100 s.Most of the modeling parameters,as listed in Table 2,are taken from Fauriel and Laloui (2012) except for the flow velocity, which is obtained through calibration.

Table 2Parameters of the MICP reference model.

Fig.4.Convergence analysis of the 1D linear problem.

Fig.5.1D MICP injection test.

4.2.2.Simulation results

The calculated profiles of urea and calcite along the axis after 5400 s of injection by using both OS and GI methods are illustrated in Fig.6.Similar to the computational results in the linear reaction case (Section 4.1), no significant difference in the numerical solutions of the MICP problem is found by using these two different numerical strategies.Due to the continuous injection of urea at the left point,the maximum urea concentration is observed there,and the urea concentration decreases in the domain.Since the reaction rate is dependent on the urea, correspondingly, the maximum calcite is produced at the injection source and decreases over the domain.As shown in Fig.6,through merely calibration of the flow velocity, the urea profile of the presented MICP model can well reproduce the corresponding profile properties presented in Fauriel and Laloui (2012).As a remark, much more complicated mechanisms, e.g.bacterial transport, decay, absorption of chemical components, variation of porosity, permeability, and mechanical stiffness increase are considered in Fauriel and Laloui (2012).Although in their study, constant injection rate and thus constant Darcy velocity (φv in Eq.(8)) are assigned, because of the consideration of porosity changes by MICP,and variations of the pore fluid velocity(v)are expected in space and time.In contrast,in our study,porosity is assumed to be constant, thus a constant pore fluid velocity is adopted as advection velocity in Eq.(14).Therefore, the calibrated constant advection velocity(7.9×10-6m/s)in our study is different from the injection velocity in Fauriel and Laloui(2012).

4.2.3.Convergence analysis

Fig.6.Computed axial profiles of calcite and urea after 5400 s of MICP injection using the OS and GI strategies, and the computed urea profile in Fauriel and Laloui (2012).The spatial element size is Δx = 2.5 ×10-3 m.

In order to compare the numerical accuracy of the OS and GI strategies for modeling the MICP problem, we carry out a convergence analysis in this section.The Courant number is kept to be constant as Cr = 0.316 and the gird size Δx varies from 0.000625 m to 0.1 m.We use the integrated urea mass Mureaover the whole domain as a measure for the relative error:

where the index j is employed to distinguish different numerical solutions calculated by the adoption of different grid sizes and time steps(to ensure a constant Courant number),and Murearefers to the calculated urea mass by using a very fine mesh with Δx = 2.5 × 10-4m with the GI-TDG(3) method.Fig.7 shows the convergence behaviors of the different numerical strategies.Unlike solving the linear problem, the OS strategy fails to present converged solutions.This could be due to the accumulation of the splitting error when solving nonlinear problems.Although it is known that the splitting error will only reduce the numerical accuracy to the 1st order, in this study, we obtain an even worse solution.In contrast, the GI strategy results in converged numerical solutions, and a 2nd order convergence rate is observed with the GI-TDG(3) method.

4.2.4.Computational efficiency

We further compare the computational time of solving the MICP problem using different numerical strategies.The numerical setups are the same as that presented in Section 4.2.1.As shown in Fig.8,the computational time increases with increasing number of degree of freedom.In this study, the GI strategies are implemented using the MATLAB code which is supposed to be less computational efficiency than the implementation of OS strategy-based C language.Moreover, the problem is solved within a time-space domain by using the TDG method, which involves more degrees of freedom than using the method of lines.The presented GI strategy’s computation cost shall be higher than the OS.However,as shown in Fig.8, the GI strategies generally result in shorter computational time than the OS strategy for the simulations of the MICP problem in this study.On one hand,this might be due to the fast convergence of the GI-TDG methods for solving the nonlinear ADR problems.On the other hand, the data transfer between OGS and PHREEQC causes additional computational time in each computation step.

The GI-TDG method benefits from its A-stability property and can tolerate large time steps even with a Courant-Friedrichs-Lewy(CFL) number larger than one.Such a property can lead to more efficient computations than using the OS strategy.In thiscomparison study, we limit ourselves to the time steps for which both numerical strategies are stable.The biochemical processes involved in the MICP problems are complex.With more understanding of the MICP process mechanisms, we can foresee that a coupled multiple time scale model will be needed to describe the physical problems in the future.The time step for solving such problems is limited to the most crucial process (or equation).For this reason,the GI strategies presented in this study are potentially a better choice for modeling the MICP problems.

Fig.7.Results of the convergence analysis for MICP model using different numerical strategies.

Fig.8.Comparison of computational cost using OS and GI methods.

4.2.5.Numerical oscillation

Even though the numerical solutions of the MICP model are almost identical by using the presented OS and GI strategies, one can still observe differences when the mass transport process becomes more advection dominated than the previously introduced scenario in Section 4.2.1.In this section, we investigate a MICP problem with a large grid Péclet number as Peg= 12.99.To set up such a numerical example,we change some model parameters.We increase the injection flow velocity to 8.5 × 10-5m/s which corresponds to the injection rate in the MICP test 3A presented in Martinez et al.(2013).Moreover,according to the measured range of values of laboratory-scale dispersivity reported in Frippiat et al.(2008), we apply a smaller dispersivity here as 1.83 ×10-4m.

The numerical solutions of the urea concentration at the early time t = 116 s are shown in Fig.9.Either the OS or the GI strategy provides reasonable results.However, one can observe a small difference at the upstream of the transport front.The OS strategy produces a more significant overshoot than the GI strategy.The maximum urea concentrations simulated by using the OS and GI methods are 0.508 mol/L and 0.501 mol/L, respectively.The generation of such numerical oscillation is due to the sharp initial gradient of the urea concentration in combination with a large grid Péclet number.The results presented in Fig.9 demonstrate that the GI strategy is superior in controlling the numerical oscillation.However, the numerical oscillation observed from both strategies are small, and they all disappear at a later time.For most MICP problems with slow injections, both of these strategies are expected to produce smooth enough results.However, if a large injection rate is considered,one should be careful on the choice of a proper numerical method to obtain stable solutions.

5.Sensitivity analysis

Fig.9.Comparison of the urea concentration profile at early time (t = 116 s) under a large grid Péclet number condition.The profile in blue corresponds to using the 3rd time order TDG method.

In the sensitivity analysis,two reaction-relevant parameters are considered.One is the first Damköhler number DaI(Eq.(13c)),which represents the ratio of the reaction rate with respect to the advective transport rate.Based on the DaI, the system can be characterized as either reaction-dominated(high DaI)or advectiondominated(low DaI).The other one is the half-saturation constant Km(Eq.(9)), whose value is dependent on the species of the microorganism and the environmental conditions,e.g.temperature and pH value.For the sensitivity analysis,several simulations with different values of DaIand Kmare carried out.The characteristic length l0and concentration C0for the calculation of DaI(Eq.(13c))in each simulation case are set to be 0.5 m and 1 mol/L,respectively.Except for the values of DaIand Km, all the other model settings remain the same as the MICP model presented in Section 4.2.1.The integrated urea mass over the model domain in each case after 20,000 s continuous urea injection from the left boundary,and urea consumption by the MICP reaction are finally compared.

5.1.First Damköhler number

Four simulations with different DaIvalues of 1.27,2.53,25.3,and 50.6,respectively,are carried out.According to the definition of DaI(Eq.(13c)),in the real MICP application,a variation in DaIrepresents that different amounts of biomass is injected into the considered system during the specific MICP treatment duration.As a consequence,different urea masses should be consumed by reaction.The calculated total urea masses are illustrated in Fig.10 as a function of DaI.Consistent results are obtained by both GI and OS methods.With the increase of DaI,less urea remains in the system.The high DaIvalue indicates a dominant effective reaction rate over the transport rate.Thus,in the system with high DaIvalues within the same reaction time, more urea can be consumed, therefore, less urea mass remains than that in the system with low DaI.However,the reduction in urea mass can conversely lead to the decrease of the reaction rate (Eq.(9)).Therefore, as shown in Fig.10, the decrease rate in the urea mass reduces gradually with the increase of DaI.Moreover, it can be seen from the calculated results that at DaI=50,despite the continuous supplementation of urea,the total urea mass tends to 0 which indicates that the system is absolutely dominated by reaction and the injected urea can be fully consumed at the injection source before it can get further transported.In practice, such a reaction-dominated system could be found in the MICP application with low injection rate where the reaction could be much faster than the advective transport.In such systems,usually calcite precipitation takes place preferentially in the area near the injection source that leads to a non-uniform calcitedistribution.If large amounts of calcite accumulates around the injection source,the inlet would be sealed,thus the transport of the urea/calcium solution can be hindered.In contrast,in an advectiondominated system, such as that treated by the MICP with high injection rate, a more homogeneous spatial distribution of calcite is expected.However, in this case, the transport is much faster than the reactions and the reactants could get out of the system before they are sufficiently consumed.A sufficient MICP reaction as well as a desired optimal calcite distribution can be obtained only with an optimized relation between the reaction rate and transport rate.Due to the flexible settings, numerical simulation could play an important role in the optimization of such MICP treatment conditions.

Fig.10.Computed total urea mass after 20,000 s using different first Darmöhler numbers in both GI and OS methods.

5.2.Half-saturation constant

Similar to the sensitivity analysis of DaI, four simulation cases with different Kmvalues varying from 0 mol/L to 0.5 mol/L are carried out.These values belong to a reasonable range of the halfsaturation constant by the specific microorganisms.The computed total urea masses after 20,000 s of injection by using different Kmvalues are illustrated in Fig.11.The urea masses computed using the GI method are almost identical with those computed using the OS method.By the large Kmvalue, more urea mass retains in the system, which indicates that the reaction rate reduces with the increases of the Kmvalues.Physically,Kmcould be considered as a measure to quantify the relationship between the behavior of microorganisms and their growth-limiting substrate.To be specific, the low Kmrepresents that a microorganism, on one hand, has a great capacity of rapid growth with low growthlimiting substrate, and on the other hand, its maximum specific growth rate or reaction activity can be reached by low growthlimiting substrate.In other words, low Kmrepresents generally a high bacterial activity.Regarding the modeling aspect, lower Kmshould lead to a higher reaction rate and vice versa.Since the bacterial behavior is not specifically considered in our model, by changing the value of Km, variations of bacterial activity-related reaction rate due to e.g.the change of ambient environmental factors could be captured in a simple manner.In the computed results, as the Kmvalue tends to zero, according to Eq.(7), the specified maximum reaction rate is reached.It can be seen from Fig.11 that with the specified maximum reaction rate of 2 ×10-5mol/(L s), still a minimum urea mass of 64 mol/m2after 20,000 s of MICP reaction remains in the domain.

Fig.11.Computed total urea mass after 20,000 s using different values of halfsaturation constant in both GI and OS methods.

Finally, we compare the impacts of DaIand Kmon the MICP reaction behavior.The results are shown in Fig.12.For a better comparison, the values of DaIand Kmin each simulation case are normalized with the reference values (DaIr= 1.27, and Kmr= 0.01 mol/L).The reference values are identical with the values adopted in the MICP model presented in Section 4.2.1.The increase of Kmleads to a slight increase in the integrated urea mass over the domain.However,the increase of DaIresults in a stronger reduction in the urea mass.The computed urea mass with respect to each parameter indicates that the DaIhas a much strong influence on the reaction behavior than the Km.Thus, based on the current model concept to obtain a reasonable prediction of the MICP reaction behavior, the value of DaIshould be carefully determined.

6.Conclusions and discussion

In this work, we performed numerical studies of the ADR processes in MICP, where the reaction processes are considered as nonlinear and described by the Michaelis-Menten type kinetic model.To solve the coupled problems of the reaction and transport processes,we applied two different solution strategies,namely the OS and GI methods.The simulations based on the OS and GI methods were conducted using the open-source finite element code OGS-IPHREEQC and our in-house Matlab code, respectively.The performance of these two strategies in solving such nonlinear reactive mass transport problems was compared and discussed.

Both numerical codes have been verified against the analytical solution of a 1D linear reactive mass transport problem.Some modeling parameters are calibrated.It was found that our models can well reproduce one of the urea profiles computed in Fauriel and Laloui(2012),where a more complicated reactive transport model was considered.Based on the simulations of this MICP injection test,we compared the computational performance of the OS and GI numerical strategies.Although both numerical strategies provided comparable simulation results, the OS strategy failed to converge during the convergence analysis.In contrast, converged numerical results were obtained using the GI strategy.Especially, the convergence rate increased rapidly as long as a higher order time scheme was applied, i.e.TDG-GI method with 3rd order.In the comparison of the computational efficiency,the GI strategy showed a better performance as well.It is speculated that in the simulation using the OS strategy, the data transfer between the OGS and PHREEQC are more time consuming than the solution process of the algebraic system.In the testing of the numerical oscillation,theOS strategy exhibits a slight larger overshoot behavior than the GI strategy in the urea profile at the early injection stage.However,in the later phase,the oscillation vanished in both cases.This indicates that even in the case of relative large Pegnumber, e.g.MICP treatment with high injection rate, accurate results could be delivered by both solution strategies.Further, we performed the sensitivity analysis of the first Damköhler number (DaI) and the halfsaturation constant (Km), which are considered as the most important parameters in the characterization of systematic advection and reaction behavior.The results of the sensitivity analysis indicate that compared to Km,DaIhas a stronger influence on the urea mass changes in the model domain.Thus, in the application by means of designing the system DaI, e.g.through controlling the injected bacterial concentration and the injection rate,one could get a desired MICP performance.

In this study, comparable accurate simulation results and computational performance are obtained in both OS and GI strategies, which indicate that both numerical strategies are appropriate in the simulation of the biogeochemical processes in MICP,where simplified reactive mass transport processes are considered.For the case of MICP applications under more complicated flow and reaction conditions,our current model could be modified to include a more complex reaction model.Besides the reactive mass transport of further chemical components,the bacterial behavior and its influence in the reaction could be considered.Since the presented GI method provides more accurate solutions and shows better convergence behavior in this study, we expect that the GI method performs better for more complex models than the OS method.However, the GI method, especially the TDG method we used in this study, requires more computer storage than the OS strategy.Therefore, it is still an open question on which method will solve complex 3D MICP models more efficiently.We think there are more we need to look into in this aspect in the future.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

The authors Xuerui Wang and Udo Nackenhorst would like to acknowledge the financial support from the German Research Foundation (DFG) (Grant No.NA 330/20-1).Dianlei Feng wants to thank the DFG under grant No.FE 1962/1-1(426819984) for financial support.The authors also want to thank the Open Research Fund of State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics,Chinese Academy of Sciences (Grant No.Z019002).