Micro-feature-motivated numerical analysis of the coupled bio-chemo-hydro-mechanical behaviour in MICP

A coupled bio-chemo-hydro-mechanical model (BCHM) is developed to investigate the permeability reduction and stiffness improvement in soil by microbially induced calcite precipitation (MICP). Specifically, in our model based on the geometric method a link between the micro- and macroscopic features is generated. This allows the model to capture the macroscopic material property changes caused by variations in the microstructure during MICP. The developed model was calibrated and validated with the experimental data from different literature sources. Besides, the model was applied in a scenario simulation to predict the hydro-mechanical response of MICP-soil under continuous biochemical, hydraulic and mechanical treatments. Our modelling study indicates that for a reasonable prediction of the permeability reduction and stiffness improvement by MICP in both space and time, the coupled BCHM processes and the influences from the microstructural aspects should be considered. Due to its capability to capture the dynamic BCHM interactions in flexible settings, this model could potentially be adopted as a designing tool for real MICP applications.


Introduction
In recent years microbially induced calcite precipitation (MICP) opened new horizons in geotechnical engineering, as it provides the potential of development of novel costeffective and sustainable solutions to a wide range of engineering problems, e.g. improvement in soft underground, contamination control, erosion protection, remediation of liquid migration and so forth [12,48]. In the MICP treatment usually either local bacteria are stimulated, or urease producing bacteria, typically Sporosarcina pasteurii are injected into the target area in soil. The enzyme urease produced by S. pasteurii catalyses the hydrolysis reaction of urea (ðNH 2 Þ 2 CO), where ammonium (NH þ 4 ) and carbonate ion (CO 2À 3 ) are produced, Once calcium (Ca 2þ ) is present in the system calcium carbonate (CaCO 3 ) can be formed and precipitates, If MICP occurs in the porous media such as in soil, the precipitated calcite occupies the pore space reducing the connections of flow paths and building bonds among the solid grains. As consequences, the soil permeability can be remarkably reduced, and the mechanical stiffness of soil can be noticeably improved. Compared to the traditional method for soil improvement, such as chemical grouting, treatment through MICP is believed to be more environmentally friendly and the MICP-treated soil is expected to exhibit more stable properties. Despite these promised advantages MICP has not yet been widely applied in soil treatment. The main reasons are: (1) there is lack of methods to ensure the uniform treatment and to control the treatment quality, especially in the MICP application in the large field; (2) the cost is not competitive compared to other conventional grouting methods; (3) the ammonium by-product could lead to groundwater contamination [22]. Furthermore, the complicated processes, especially the couplings effects involved in MICP are currently not well understood. To get a better insight into MICP, in the past decade researchers have carried out numerous experimental and numerical work. A detailed review of the state of the art of MICP investigations can be found in [42,48]. Generally, two different categories of MICP tests can be found in the literature. One is the biochemical test which focuses on monitoring the concentration evolution of bacteria and chemical components under different treatment conditions (e.g. [13,29,32,36,45,56]). In some of these experiments the final distribution or the temporal evolution of the hydraulic conductivity was additionally measured (e.g. [32,45]). The other one is the pure mechanical test in which the strain-stress behaviour of the MICP-treated soils with different calcite contents was tested under different loading conditions (e.g. [1,4,16,30,47,49]). From the year 2016, a new perspective in MICP has drawn attention which focuses on investigating the MICP-caused variations of microstructure in soil and their impacts on the macroscopic properties. For example, in Dadda et al. [9] and Dadda et al. [10], using 3D X-ray micro-tomography, the MICP-induced changes in the microstructural properties with the corresponding changes in the macroscopic physical properties were investigated. In Terzis and Laloui [47] the critical contact surface area in the MICP-treated soil is analysed through X-ray micro-tomography and 3D volume reconstruction. Within this work a series of relevant micro-scale parameters were derived. Lin et al. [31] proposed a series of analytical models to investigate the evolution of permeability and shear/pressure (S/P) waves velocity in MICP soil with different idealised pore-scale calcite distributions.
Except for the experimental work numerical modelling plays an essential role in the conceptual analysis of the coupled BCHM mechanisms in MICP. Similar to the experimental work, the mainstreams of the current modelling work can be divided into two categories: biochemical modelling and mechanical modelling. In the biochemical modelling numerous studies have been attempted to understand the reactive mass transport and the hydrodynamic processes in MICP at different scales. In [52] a mathematical model was proposed to describe the advection, diffusion and reaction processes in MICP. In this model a single overall kinetic rate which is dependent on the bacteria and urea concentration is introduced to consider the entire reaction processes. Kozeny-Carman equation is adopted to calculate the permeability reduction following the decrease of porosity. Later, the expression of this overall kinetic rate is modified by van Wijngaarden et al. [57] to simulate a MICP test on a 5-m column. To quantify the MICP performance in the reduction of fracture aperture in a field test, Cuthbert et al. [8] coupled the reactive mass transport model with a density-dependent groundwater model where the MICP-induced urea hydrolysis is described by a first-order kinetic model. In [33] a more complex reaction model is proposed in which the urea hydrolysis rate is calculated by a mixed equilibrium kinetics, whereas the calcite precipitation is considered as a pure kinetic reaction. Ebigbo et al. [13] developed a Darcy-scale MICP model in which both liquid and gas phases are considered. In this model the rates of urea hydrolysis and calcite precipitation are calculated separately by two different kinetic models, and the porosity permeability relation is described by a modified Kozeny-Carman equation. This Darcy-scale model was later revised in Hommel et al. [24] with a modified expression for the urea hydrolysis rate. More recently, large-scale 3D MICP modellings were conducted in Cunningham et al. [5] and Nassar et al. [37]. Based on the numerical code Open-FOAM, Minto et al. [35] developed a reactive transport model which allows the simulation of MICP processes from the pore scale to the large continuum scale. In our previous study [53] a coupled biochemical-hydraulic model was developed to predict the temporal and spatial permeability reduction in soil during the MICP treatment. This coupled model has been validated against two different experimental data sets. As one of the most current works, Hommel et al. [25] proposed a reactive mass transport model for enzymatically induced calcium carbonate precipitation in which the influences of pH value and temperature on reactions are considered. Except for the above-mentioned macro-continuum model, Qin et al. [41] attempted to analyse the reactive transport process of MICP in micro-scale based on the pore network model.
In most continuum biochemical models mentioned above, the initial porosity and permeability are considered as the only two intrinsic properties of the host soil in the calculation of the reaction process and the porosity/permeability variations. In this way, the different MICP behaviours in soils with similar initial porosity and permeability but different microstructures are challenging to be reproduced.
Compared to the biochemical modelling work limited attention has been paid to the modelling of the mechanical behaviour in MICP. Recently, an elasto-plastic model for MICP-treated soil was developed by Gai and Sánchez [19] based on the critical state yield surface and sub-loading concepts. In this model the increase of the mechanical strength is considered in a yield surface whose envelope is dependent on the calcite contents. The softening behaviour due to the degradation of the calcite bonds is considered by reducing the yield surface area following the increase of the plastic strain. By modifying an existing model for cemented soil Nweke and Pestana [38] proposed a constitutive model for MICP soil in which the degradation of the stiffness and strength caused by mechanical failure is considered. In this model the properties of the bio-cemented material are assumed to be dependent on the confining stress and the intrinsic microstructure. Later, this model was modified by Nweke and Pestana [39] through adding a cement content-dependent parameter. More recent, Xiao et al. [59] developed a constitutive model to describe the behaviour of the MICP-treated calcareous sand under cyclic loading, where a hyperelastical model is coupled with thermodynamic plasticity. In this model, a reactant index is defined to consider the increase of dry density, cohesion and bonding effects during MICP treatment. Besides, a degradation law is introduced to describe the decrease of cohesion and degradation of bond structure during cycle loading. Except for the modelling work on the macro-scale mechanical constitutive relations, several studies have been attempted to understand the micro-scale mechanical behaviour of the bio-cemented soil using the discrete element method (DEM) (e.g. [14,17,60]).
In the most current MICP mechanical models the mechanical properties are assumed to only rely on the average calcite contents. However, in such a way the dynamic changes of the material properties and their spatial variations cannot be captured. One exception is the numerical model proposed in Gajo et al. [21], in which the temporal evolution of the mechanical properties is captured in a coupled chemo-mechanical scheme. Besides, in this model several micro-scale-related parameters are introduced and incorporated in a series of cross-scale functions to consider the effects of micro-scale structure on the macro-scale properties. However, in this model the reaction process is simplified as a mass increase or decrease rate, and the hydraulic process is not considered. Therefore in this model the spatial distribution of the bacteria, chemical components and their influences on the reaction behaviour are not addressed. As far as we know, in the current literature only one modelling work that presented by Fauriel and Laloui [15] crosses the above-mentioned two major modelling categories and addresses the fully coupled BCHM processes in MICP. However, this model has not been validated against any experiments. Besides, the micro-features are not considered.
The current state of the MICP modelling work reveals two major gaps. First, the coupled BCHM behaviour of MICP is not well understood. Second, there is a scarcity of models that link the variations in the micro-features with the macro-properties changes. Regarding the current gaps in MICP modelling the main objectives of this study are to (1) develop a comprehensive model to analyse the relevant coupled BCHM processes in MICP, (2) address the influence of microstructural properties on the macroscopic MICP behaviour. The paper is structured as follows: in Sect. 2, the relevant physical processes and modelling concept are introduced. In Sect. 3, the developed model is applied in three numerical experiments, where different coupling effects are analysed and discussed based on the computational results. Finally, conclusions are given in Sect. 4.

System definition and coupling processes
The present model is developed based on our previous modelling work [53], upon which the coupled deformation process and the relevant micro-structural effects are additionally considered. In the model soil is assumed to be homogeneous and isotropic porous medium which is fully saturated with a single fluid phase. The computed domain is on the representative elementary volume (REV) which is composed of a solid phase (s) and a liquid phase (l). In the REV transport of multi-species and multi-components is considered. Solid grains (g), calcium carbonate (CaCO 3 ) and the attached bacteria (b s ) are assumed to be immobile and considered as the entire components of the solid phase. In the liquid phase five components are considered: pore water (w), urea (urea), calcium (Ca 2þ ), ammonium (NH4 þ ) and the suspended bacteria (b l ). As shown in Fig. 1 the following processes and coupling effects are addressed: (1) multi-dimensional flow in porous media; (2) bacterial behaviour including transport, decay and attachment; (3) mass variations in liquid and solid phases caused by biochemical reactions; (4) mass transport of aqueous components governed by advection and dispersion dynamics; (5) deformation of the solid skeleton; (6) reaction kinetics depending on the temporal and spatial concentration distribution of reactants and bacteria; (7) reduction of porosity and permeability due to the filling of calcite in pore space; (8) increase of mechanical stiffness due to the calcite cementation in the soil matrix; (9) change of the effective stress due to the variations in the pore pressure field; and (10) change of pore pressure due to the deformation of solid skeleton. Except for the above-mentioned macroscopic couplings the microstructure-related change in the reaction behaviour and the macroscopic properties are considered.

Mass transport
The transport of the suspended bacteria and aqueous chemical components is governed by advection, diffusion and dispersion, where / is the porosity, C i denotes the concentration of the aqueous components i (i 2 urea, Ca 2þ , NH þ 4 and suspended bacteria b l ). q i is the sink/source term caused by reactions and q is the vector of Darcy velocity. D is the diffusiondispersion tensor which is composed of diffusion D i f and mechanical dispersion D p . D p is further given by where a L and a T are the longitudinal and transversal dispersion coefficients, respectively. v ¼ q=/ is the vector of pore water velocity.

Biochemical reactions
Overall controlled kinetic reactions. Combining Eqs. 1 with 2 an overall MICP reaction could be obtained, According to Eq. 5, in the reaction model we considered the relevant reactants urea, Ca 2þ and the end products NH þ 4 , CaCO 3 . The intermediate products, e.g. dissolved hydrates and bicarbonate, are neglected. van Paassen [50] has confirmed in a series of batch experiments and numerical studies that in a system with sufficient calcium, once calcite has nucleated the precipitation rate is approximately consistent with the hydrolysis rate (r c % r h ) in the most reaction duration. Following this assumption, we introduce an overall kinetic rate to describe the rate of both urea hydrolysis and calcite precipitation. Similar to the modelling work by van Paassen [50], Lauchnor et al. [29], Fauriel and Laloui [15], van Wijngaarden et al.
[51], Minto et al. [35] and Wang and Nackenhorst [53], this overall kinetic rate is assumed to follow the Michaelis-Menten kinetics, where u sp is the maximum urease capacity of bacteria, C b t is the total bacterial concentration and equal to the sum of the suspended C b l and attached bacterial C b s concentration. k m is the half-saturation constant. Substituting the overall reaction rate into the mass transport equation (Eq. 3), the sink/source term for each chemical component can be explicitly expressed as where n i denotes the chemical stoichiometry. According to Eq. 5, n urea ¼ À1, n Ca 2þ ¼ À1 and n NH þ 4 ¼ 2. Bacterial behaviour. Except for the transport of suspended bacteria (Eq. 3), bacterial decay and attachment are considered. Usually, before MICP treatment bacteria are accumulated in the laboratory to get the optimal density thus activity. During the MICP treatment due to the possible lack of oxygen, encapsulation of the bacteria by calcite, bacteria decay is expected to be dominated over the bacterial growth [57]. Thus, in our model bacterial growth is neglected. Both bacterial decay and attachment are described by first-order kinetics. Thus, the change rate of the total bacterial concentration is given by The change rate of the suspended and attached bacterial concentration due to decay and attachment is expressed as where k d and k att are constants for decay and attachment rate, respectively. By substituting Eq. 10 for decay and attachment of the suspended bacteria, the sink term in the mass transport equation for the suspended bacteria (Eq. 3) can be derived, Cross-scale relations. The work in [23,26] suggested that as a microstructure-related feature the availability of nucleation sites is an important factor that can affect calcite precipitation. To consider the effects of the nucleation sites availability we introduce a specific reaction surface a r in the reaction model. The specific reaction surface is defined as for calcite bonds deposition and growth available area per unit volume of soil. We assume that the nucleation sites of calcite are only available around the contact areas of solid grains. This holds the fact that during injection flow velocity around the contact area is usually much lower than in other areas, which favours the retaining and attachment of the bacterial mass, and thus the calcite precipitation. These phenomena have been observed in several experimental studies under continuous flow condition (e.g. in [10,47]). Noticed that, the assumptions made above may not be valid for MICP treatment under stop-flow condition, where calcite could also precipitate on the entire particle surface and coating the solid grains. Further, the assumed reduction of deposition sites may not be valid for the MICP treatments, where biomass is continuously supplied or stimulated. In those treatments, new deposition sites could be formed due to the attachment of the newly injected bacteria on the surface of solid grains and precipitated calcite.
To incorporate the microstructural properties in the defined specific reaction surface the geometric method (GM) [44] is used, by which all solid grains are idealised as spheres with identical average diameter, and all calcite bonds are assumed to be cylindrical with identical average length and radius. According to the study in [20,21], two different types of reaction surface in large and small porosity configurations are considered (Fig. 2).
In the case of a large porosity configuration where solid grains have no direct contacts (Fig. 2a), it is assumed that the initially available contacts surface area at each contact point can be fully filled by cylindrical bodies with a radius of R bm . With the definition of N b as the number of initially available contact points per unit volume, the initially available contacts surface area is given by, With precipitation and covering of calcite on the solid surface, the specific reaction surface reduces gradually. If calcite bonds are assumed to be cylindrical with average radius of R b . The available specific reactive surface during MICP treatment can be expressed as where N ab denotes the number of the deposition sites which are occupied by calcite bonds. Noticed that in Eq. 13 the values for N b and R bm are dependent on the micro-properties of the host soil. In the model we defined a critical porosity / c , below which the pore space is no longer interconnected. Thus, the void volume around the contact area can be approximated by (/ 0 À / c ), where / 0 is the initial porosity. According to the definition of R bm , this initial void volume can also be approximated by (N b pR 2 bm L b ), with L b as the average distance between two neighbour grain. Assuming that these two approximations are identical the R bm can be calculated by In the case of a small porosity configuration where solid grains have direct contacts (Fig. 2b), the specific reaction surface becomes very small, thus the geometric assumptions can no longer be held. According to [20,21] an asymptotic relationship is adopted to describe the reaction surface, where a r0 is a constitutive constant. More detailed derivation of this asymptotic relationship can be found in [21]. According to [20,21] the expressions for the specific reaction surface for large and small porosity configurations can be combined by using porosity related terms / and 1 À / as weighting factors, Noticed that, MICP induces gradual reduction of porosity. As a consequence, the weighting term for large porosity configuration in Eq. 16 decreases continuously, which represents the variation of the MICP-treated soil from a large porosity configuration to a small one. Similar as in [21], the growth rate of N ab is assumed to depend on the current available deposition sites (N b À N ab ) and the increase rate of the calcite mass _ M CaCO 3 , where k 2 is a kinetic constant. Further, we assume that the calcite bonds grow isotropically and homogeneously in all three orthogonal directions around the deposited nucleation sites. Thus, the growth rate of the radius of the calcite bond can be calculated according to its volume increase rate _ V CaCO 3 , where k rb is a kinetic constant. Under the assumption of a linear relationship between the increase of overall kinetic rate and the increase of the specific reaction surface, Eq. 6 can be further given by, where k 1 is a reaction rate constant.
With the introduction of the specific reactive surface a r in the reaction model the following changes of the reaction behaviour due to the variations of the microscopic properties can be addressed: (1) decrease of the rate of urea hydrolysis and calcite precipitation due to the decrease of deposition sites as more and more contact points are occupied by calcite; (2) decrease of precipitation rate due to the reduction of the surface area of solid grains as they are gradually covered by calcite.
According to Gajo et al. [21], based on the similar approximation for the specific reaction surface a r another quantity a b to quantify the cementation level of soil by intact bonds can be derived, where h and g are constitutive constants. The value of a b varies in the range of [0, 1], where a b ¼ 0 represents the non-cemented soil and a b ¼ 1 represents fully cemented state. The constant a b is further adopted to describe the mechanical behaviour of the MICP-treated soil, which will be discussed in Sect. 2.2.4.

Flow in porous media
The governing equation for flow in porous media is derived according to the macroscopic mass conservation, where X i ¼ q i m i denotes the mass flux of aqueous component i with m i is the molar mass of component i (i 2 urea; Ca 2þ ; NH þ 4 ). Using the chain rule and by introducing the fluid compressibility modulus dq l dp l , Eq. 21 can be further expressed as where p l is the pore pressure, q l is pore fluid density, v s is the vector of the solid velocity and Q l ¼ ðX i þ q w Þ=q l is treated as the sink/source term. The term of the specific discharge between the liquid and solid phase /ðv f À v s Þ is derived based on the Darcy's law, where K is the permeability tensor, l l is the fluid viscosity and g is the gravity vector. In Eq. 22 the term for porosity variation is assumed to be only dependent on the solid deformation, as the calcite induced porosity variation is updated explicitly with the deformation process. According to the assumptions of the linear poro-elasticity the deformation dependent porosity variation can be derived based on the solid mass conservation [34], where K s is the bulk modulus of solid grain and a is the Biot coefficient which is assumed to be 1 in the calculation. Substituting Eqs. 24 in 22 yields where _ e v ¼ r Á v s denotes the change rate of the volumetric strain, which can be calculated in the deformation process (sec. 2.2.4). It is emphasised that the terms containing _ e v and K s represent the most relevant effects from the mechanical field to the hydraulic field in both spatial and temporal perspectives. The term ð / K l þ aÀ/ K s Þ represents the storage capacity of the porous medium, which is assumed to be constant in the model with a specified storage coef- As the magnitude of the porosity change is about ten order lower than the value for compressibility modulus, the influence of porosity variation on the storage capacity of soil can be neglected.
Porosity and permeability variations. It has been observed in many experimental studies (e.g. [8,15,51,52,57]) that the filling of the calcite in pore space is the main reason for the pore space reduction. Thus, in the model the reduced volume of pore space is assumed to be equal to the produced calcite volume, where V bulk ; V pore ; V bulk are the bulk volume of the porous medium, the volume of pore space and calcite. q CaCO 3 and m CaCO 3 are the density and molar mass of calcite, respectively. According to [53], in a small time interval the decrease rate of the total porosity can be obtained by In our previous study [53] the performance of different porosity permeability relationships in MICP modelling has been investigated. Among those, the modified Kozeny-Carman (KC) equation based on the effective porosity concept is proven to be the most suitable one for MICP. Thus, in the present model the same modified KC equation is adopted, Compared to the original KC equation the total initial and current porosity is replaced by the effective initial and current porosity ð/ eff Þ i , ð/ eff Þ 0 . The effective porosity is defined as the fraction of pore volume which is interconnected. According to [28] the effective porosity is given by where a is a constant. / c is the critical porosity below which the referred area is no longer fluid permeable. The constants a and / c are usually dependent on the microstructure of the porous media. It is worth remarking that the introduction of the modified KC equation and the effective porosity allows the model to capture the following observations in MICP. The pore throats are clogged by calcite leading to form isolated non-permeable areas even though in those areas not all pore spaces are fully filled with calcite ( [7,40,54,55,58]). A detailed explanation of the physical meaning of the effective porosity and the modified KC equation can be found in [53].

Deformation of the solid skeleton
The effective stress concept based on the Biot theory is applied to describe the influence of pore pressure on solid deformation. In the stress field, the diagonal of the total stress tensor is composed of effective stress and pore pressure. With the convention that stresses are positive when they are tensile whereas pore pressure is positive when it is compressive, the effective stress is given by where r eff is the effective stress tensor. Using Voigt notation for the stress tensor m ¼ ½1; 1; 1; 0; 0; 0 T is the matrix form of Kronecker's delta. In the calculation the effective stress is considered as the actual contributor to the deformation of the solid skeleton. In the model only the elastic deformation is considered. The relationship between the stress tensor increment to the elastic strain rate tensor is given by where D e is the elasticity tensor.
Increase of the mechanical stiffness. A calcite bondsdependent elastic modulus is introduced to describe the MICP effects on the mechanical stiffness (Eq. 33). Through the observation via scanning electron microscopy (SEM), Terzis and Laloui [47] suggested that the intact calcite bonds predominately contribute to the improvement in the soil mechanical properties. In contrast, incomplete bonds are supposed to play a minor role. Therefore, in the model the elastic modulus is correlated with quantity for the intact calcite bonds a b (Eq. 20) instead of with the total calcite content. The expression for the elastic modulus is obtained from Gajo et al. [20], in which a linear interpolation between the uncemented and fully cemented soil is proposed, where E g and E b are, respectively, the elastic modulus of the host soil and the fully cemented soil, and is a constitutive constant.

Solution strategies
Under the considerations above, the final system to be solved is composed of the field equations for mass transport of aqueous components (C) (Eq. 3), for reaction kinetics (B) (Eq. 19), for flow in porous media (H) (Eq. 25), for the deformation of solid skeleton (M) (Eq. 32), and the corresponding constitutive relationships describing the coupling effects among each process. In total, eight primary unknowns are contained in the system, namely concentration of urea, calcium, ammonium, suspended bacteria, attached bacteria, calcite, pore pressure and displacement. The final algebraic equations describing the entire system are solved using the finite element method (FEM), which has been implemented in the opensource FEM code OpenGeoSys (OGS) [3,27]. The standard Galerkin method is adopted for the spatial discretisation, whereas the temporal discretisation is based on the implicit backward Euler method. As shown in Fig. 1, it is expected that reactions are strongly coupled with the mass transport process. Thus, they are solved in a monolithic way (BC). To be specific, the reaction term (Eq. 19) is incorporated as the sink/source term in the mass transport equations (Eq. 3) using the global implicit method [43]. On the contrary, the flow (H), deformation (M) and the BC subsystem are solved iteratively in a staggered way. In each iteration the couplings are considered by incorporating the dependency of the parameters and primary variables from different processes. The detailed solution strategy for the entire system is illustrated Fig. 3. Within one time step, the hydraulic process is solved first with the updated density, permeability, total and effective porosity from the previous time step. The fluid velocity calculated by solving the H process is then transferred into the BC process. With the specific reactive surface calculated from the previous time step, the BC process is solved and the spatial distributions of each aqueous biochemical components are obtained. As the solid components, the concentration of attached bacteria and calcite (mass with respect to the liquid volume) is calculated at each local node, where the indices n and i refer to the nth node and the ith time step, respectively. After obtaining the calcite content the total/effective porosity and the permeability are accordingly updated (Eqs. 27,29,28). Further, based on the calculated calcite content and total porosity, the microstructure-related parameters, i.e. N ab , R b , a b (Eqs. 17,18,20) are determined. With a b calculated from the BC process and pore pressure derived from the H process, the M process is solved. Next, the change rate of the volumetric strain calculated in the M process is substituted in the H process. With these updated parameters the next iteration is started with the same sequence. Same as in our previous study [53] a criterion for stopping the iterations is defined, where X i;kþ1 and X i;k are the vectors of the primary variables of each process that are calculated in the ith time step by the (k þ 1)th and kth iterations, respectively.

Numerical case study
In this section, three simulation cases are considered to test the model capability in the prediction of the coupled processes in MICP. First, the experiment reported in Tan et al. [45] is simulated, by which the model parameters are calibrated and the microstructural effects on the biochemical behaviour are analysed. In the second case the experiment performed in Dadda et al. [9] is considered, where the MICP-induced macroscopic as well as the related microscopic property changes are estimated. The last case is a scenario modelling on a MICP treatment that consists of continuous BCHM treatment phases. In this simulation case the dynamic hydro-mechanical response of soil during MICP treatment is analysed. The first two cases are coupled BCH modelling, whereas the last one is a fully coupled BCHM modelling.

Model set-up
In Tan et al. [45] four cylindrical sand columns with a height of 7.5 cm and a diameter of 3.5 cm were treated by MICP. These four samples have different initial properties ( Table 1). The experiment was conducted with two phases, a biological injection phase followed by a chemical injection phase. During the biological injection the number of the in-sample retained bacterial cells was determined according to the number of the total injected bacterial cells and the number of bacterial cells measured in the effluent. During the chemical injection phase solution containing 0.5 mol/l urea and calcium chloride was injected at a constant rate of 0.97 mL/min. The chemical injection was continued for 72 h during which the hydraulic conductivity  in the samples was measured every 12 h. At the end of the treatment the finally produced calcite mass in each sample was measured. This experiment has already been considered in our previous work [53], where a similar expression for the kinetic rate of urea hydrolysis as Eq. 19 was adopted. However, in [53] the u sp is assumed to be constant, in which the microstructural influences are not considered. Besides, in the previous modelling an exponential term (expðÀ t t d Þ) with a time constant t d was multiplied by the kinetic rate to capture the observed decrease of reaction rate over time. In the present work the same axis-symmetric model as in [53] with a height of 7.5 cm and a radius of 1.75 cm is set up. Considering the small model size, the initial bacterial concentration is assumed to be equivalent to the measured number of bacterial cells that were retained in samples at the beginning of chemical injection (i.e. 2.15e9 in sample 1, 2.06e9 in sample 2, 1.97e9 in sample 3 and 1.61e9 in sample 4). The initial concentration of urea, calcium, calcite and ammonium is set to be zero over the whole model domain. On the injection boundary, Neumann's type boundary condition corresponding to the constant injection rate of 0:97 ml=min is imposed for the liquid flow. Dirichlet's type boundary condition with the constant urea/calcium concentration of 0.5 mol=L is applied for the mass transport. On the noninjection boundary, the atmospheric pressure (Dirichlet type) is assigned for the liquid flow. For the mass transport this boundary is considered as an open advective boundary by which the passing mass flux is automatically calculated according to the advective velocity. The model parameters are listed in Table 2. The values of the most macroscopic parameters are obtained either from the literature or through model calibration in our previous study [53]. The values of the microscopic quantities and the cross-scale functions are determined mainly based on the geometric idealisation or the related experimental measurements. By assuming that all solid grains are identical spherical grains with the measured average diameter, the total number of grains per unit volume can be determined by, where D 50 is the measured average diameter of the solid grains. The parameters remaining for calibration in the present work are k 1 , k 2 and k rb . Their values were determined by fitting the measured calcite mass and evolution of the hydraulic conductivity in sample 2. The simulations of the tests on the other samples (sample 1, 3, 4) are performed with the calibrated values from the simulation of the test on sample 2. Case3 5e-3, (-) [18] Constitutive parameters (Eqs. 6, 17) Decay rate k d (s À1 ) 3.18e-7 [46] Attachment rate k att (s À1 ) 1e-4 Assumed Mechanical parameters (case3) Elastic modulus: Non-treated soil E g (MPa) 80 [47] Fully cemented soil E b (MPa) 300 Fitted Poisson's ratio m (-) 0.25 [11] Constitutive parameter (-)

Results comparison
The computed results in the current model are compared to both the experimental measurements and the results from our previous modelling work [53] (Fig. 4). Due to the similar initial porosity in all four samples, as shown in our previous modelling results, without addressing the different initial microstructures, only slight differences of the final calcite mass are computed among the four samples. Thus, our previous results have certain discrepancy from the experimental observations (Fig. 4a). In the current approach with the consideration of the microstructural effects on the reaction behaviour a much better agreement between the measured and computed results is achieved. Correspondingly, the observed evolution of the hydraulic conductivity, especially in sample 4, is better reproduced in the current model (Fig. 4b). Moreover, from the measurements it is seen that the hydraulic conductivity decreased intensively only at the early injection stage (the first 24 hours). After that the decrease rate of hydraulic conductivity reduces gradually which indicates a decrease in the overall reaction rate over time. This decrease of reaction rate could be caused by, on the one hand, the bacterial decay, on the other hand the reduction of the specific reactive surface due to the filling of calcite in the pore space and covering of surface of solid grains by calcite. To capture this observed decrease of reaction rate over the time, in our previous work [53] we defined a time constant in the reaction term in which the influence of the reduction in specific reactive surface on reaction rate is also incorporated. However, the introduction of a time constant is phenomenological rather than physical. In the current model, with the adoption of the micro-related parameters and the related cross-scale functions the changes of the reaction surface due to the variation of microstructure and their influence on the overall reaction behaviour can be physically represented.

Evolution of the microscopic and macroscopic properties (case2)
To get a deeper insight into the relationships between the evolution of micro-and macroscopic properties in MICP we performed numerical modelling of the experiments reported in [9]. In [9] eight columns of Fontainebleau sand (D 50 ¼ 0:21 mm) with the same diameter of 68 mm but different heights (four samples have a height of 540 mm and four samples have a height of 300 mm) were treated by MICP. Two treatments were conducted by each sample following the same procedure. In each treatment first one pore volume of the bacterial solution with 1OD 600 of Sporosarcina Pasteurii was injected from the bottom of the sample at a constant flow rate of 0:2 mm=s. After 1 h two injections of one pore volume of the chemical solution containing 1.4 mol/l of urea and calcium were carried out at a constant injection rate of 0.14 mm/s. A retention time of 10 h was set between these two injections. More treatment details are found in [9]. As the length of the sample is much larger than its diameter a one-dimensional (1D) model with a length of 560 mm is considered. Initially, the concentration of bacteria and all chemical components are set to be zero over the model domain. During injection phases, similar boundary conditions to those prescribed in the simulation case 1 for the non-injection boundary (top) and the injection boundary (bottom) (Sect. 3.1.1) are assigned with the specified injection rate and concentration in the experiment. During the retention phase, for flow both top and bottom boundaries are assigned with constant atmospheric pressure. For mass transport, open advective boundaries are prescribed at both boundaries. The model parameters are listed in Table 2. The values of the most parameters are consistent with those used in numerical case 1. Only a L , k 1 and / c are assigned with different values, as the first two Fig. 4 Comparison of the computed results with the experimental measurements and our previous results [53]: a final calcite mass versus initial porosity of each sample, b temporal evolution of the hydraulic conductivity parameters are scale-dependent and the last one is related to the microstructure.

Results comparison
The computed and measured profiles of calcite mass fraction along the column height are shown in Fig. 5, (left). It is reported in [9] that the mass fraction of calcite reaches 5% to 6% through the first treatment, which is well predicted in the model. Further, the computed final calcite distribution can also be well reproduced. Since both bacterial and chemical solution was injected from the bottom, the produced calcite mass is found maximum at the bottom and decreases over the column height. The higher amount of calcite precipitated at the injection source causes strong reduction in local pore space, and a consequential increase of local flow velocity. The local acceleration of pore flow leads to the decrease of residence of local reaction thus drive more reactants and biomass downstream towards the further area [56]. As a consequence, more calcite precipitates in the areas further from the injection source. Thus, as shown in both observed and computed results, the distribution of the final calcite mass is relatively homogeneous. Remark that the experimental data shown in this figure correspond to the averaged calcite mass fraction that were measured in columns 1-4 [9]. The computed final porosity ranges from 0.28 to 0.315 depending on the computed local calcite volume, which fits well with the estimated final porosity of 0.3 in the experiment [9]. The comparison of the computed and measured reduction in permeability following the increase of calcite volume fraction is shown in Fig. 5, right. It can be seen that the experimental observed permeability evolution [9] can be well captured in the model. An approximately linear relationship between permeability and calcite volume fraction is estimated in the model as the computed finial porosity (0.28-0.315) is far beyond the adopted critical porosity (0.12) (Table. 2). When more calcite is produced and leading to a reduction of porosity approaching or below the critical porosity, a much stronger and nonlinear decrease in permeability with the increase of calcite volume fraction is expected.
Except for the macroscopic properties, Dadda et al. [10] investigated the variation of contact properties in MICPtreated sand based on the microscopic observations. In their study, three different types of specific contact surface per unit volume of porous media (m 2 =m 3 ) are defined, namely the specific contact surface before (SAB), after MICP (SAA) and the contact surface of calcite (SAC). Correspondingly, these three types of specific contact surface are estimated in our model according to the related computational results. To be specific, SAB is assumed to be the specific area, where the solid grains initially have direct contacts with each other (see. Fig. 2b), Under the assumption of cylindrical calcite bonds, the contact surface of calcite (SAC) can be approximated based on the computed number of deposition sites and radius of calcite bonds, Similarly, the total specific contact surface after MICP can be estimated, The evolution of the ratio of different contact surfaces estimated by [10] are well reproduced in the model (Fig. 6). With the increase of calcite volume fraction a rapid nonlinear increase of the specific contact surface area after MICP over initial contact surface area (SAA/SAB) is observed. As the calcite volume reaches 10 % the contact surface area is more than three times the initial contact area. The increase of the ratio SAC/SAA is also nonlinear with the increase of calcite volume. When the volume fraction of calcite is higher than 6 % the total contact surface is dominated by the contact surface of calcite bonds (SAC=SAA [ 70 %).
To estimate the changes of the elastic modulus in the MICP-treated samples, Dadda et al. [11] adopted two different analytical methods which are based on the contact cement theory (CCT) and numerical homogenisation, respectively. As shown in Fig. 7, by the same calcite volume fraction, the normalised effective elastic modulus (effective elastic modulus (E eff ) of the treated sand normalised by the elastic modulus of solid grain (E s )) by using the CCT method is generally higher than that estimated by the numerical homogenisation. In our model, this evolution of the elastic modulus is estimated according to Eq. 33 with the computed cementation level of calcite bonds a b (Eq. 20). In [11], it is assumed that the elastic modulus of the solid grain is equal to that of the cement phase (calcite), namely E s ¼ E c ¼ 80 GPa. Equivalently, in our model the elastic modulus of the fully cemented sand E b (Eq. 33) is set to be 80 GPa. The predicted evolution of elastic modulus in our model falls between the estimated evolution by the CCT method and by numerical homogenisation (Fig. 7) within the observed range of calcite volume fraction. Besides, the model predicted evolution trends more to the elastic modulus estimated by the CCT method than by the numerical homogenisation, as in both CCT method and our model the microstructural effects, i.e. the contact properties on the elasticity evolution, have been taken into account.
3.3 Scenario analysis of the hydro-mechanical response (case3)

Description of the scenario simulation
In this section scenario simulation of a MICP laboratory experiment is considered. As shown in Fig. 8, in this experiment a saturated cylindrical sand column with a height of 30 cm and a radius of 5 cm is tested under five different phases. The considered biochemical and hydraulic test (BCH) phases have a similar procedure as by the continuous injection tests (e.g. test 3A) in Martinez et al. [32]. Specifically, in our simulation an additional mechanical test phase is considered at the end of the experiment. In the first phase (B), the biological solution with 1.5 times the pore volume containing one optical density (OD) bacteria is injected from the top of the column at a constant injection rate of 10 ml/min. After that an 8 hours no-flow retention phase (R1) is designed for the bacterial attachment. Following that a 60-hour chemical phase (C) is conducted with continuous injection of a chemical solution containing 50 mM Ca 2þ and urea from the bottom of the sample at a constant injection rate of 2 ml/min. After a second 8 hours retention phase (R2) the final mechanical phase (M) is carried out. During the BCH and R1/R2 phases the sample is loaded at constant compression stress of 1e5 Pa at both axial and lateral boundaries. During the mechanical phase the compression loading at the top is increased linearly from 1e5 Pa to 5e5 Pa within 1 hour at an undrained condition. An axissymmetric model (Fig. 8, left) was set up for this scenario case. Same type boundary conditions for flow and mass transport as in case 1 and case 2 are applied for the injection and non-injection boundaries with the specified injection rate and concentration considered in the experiment. The parameters for the mechanical calculation are listed in Table. 2. The values for E b and e are determined by fitting the measured evolution of elastic modulus of MICP-treated medium sand by Terzis and Laloui [47] in the extended simulations of tests on sample 2 and sample 3 in case 1 (Fig. 9).

Simulation results
The evolution of the computed calcite profiles during the chemical phase is illustrated in Fig. 10, left. During the chemical phase a gradual increase of calcite contents is calculated. The maximum calcite content is found in the area near the bacterial injection source (top) rather than near the chemical source (bottom). The distribution of calcite is mainly determined by the interactions among the bacterial activities, advective flow, diffusion and reaction rate. During the bacteria injection, as bacteria is continuously injected from the top, the maximum bacterial mass is found near the top of sample due to the bacterial attachment. The bacterial mass decreases gradually from the injection source (top) to the bottom. However, during the following chemical injection, the injection source is set to be opposite to the bacterial injection source, thus the urea/calcium concentration is expected to be minimum at the early injection stage around the bacterial injection source. Therefore, at the beginning of the chemical injection almost no calcite precipitates at the top of the sample (5h in Fig. 10). With the continuous injection of urea/calcium, more and more reactants get transported towards the bacteria injection source (top) which increases the reaction rate. As a consequence, with the MICP treatment more calcite is produced in the area near the bacterial injection source. A similar trend of the calcite distribution was observed in the experiment reported in [32]. The computed profiles of elastic modulus along the column height have a similar shape to the calcite distribution. To be specific, where more calcite is produced the elastic modulus is higher. Due to the inhomogeneous calcite distribution, the initial homogeneous mechanical property becomes inhomogeneous after treatment. Since the calcite distribution is mainly governed by the biochemical and hydraulic interplay, an accurate prediction of the MICP induced temporal as well as spatial variations of the mechanical properties is possible only when the fully coupled BCHM processes are considered. For a better demonstration of the MICP effects on the hydro-mechanical behaviour, an additional simulation with constant elastic modulus was carried out, where the MICP effects on the mechanical properties are not considered. The computed axial strain and stress during the mechanical test with and without consideration of the MICP effects are shown in Fig. 11. The increase of the axial loading at the top during the mechanical test causes a gradual increase of effective axial stress and strain (Fig. 11). Due to the specified mechanical boundary conditions, namely free displacement at top and lateral boundaries, no distinct difference was calculated in the stress profiles between the cases with MICP and without MICP (Fig. 11, left). The slight difference near the bottom is mainly caused by the pore pressure deviations between these two cases (Fig. 12c). A more evident discrepancy is found by the strain distributions (Fig. 11, right). As the elastic modulus is significantly improved by MICP, the computed strain in the case with MICP effects is much lower than that without consideration of the MICP effects. Since the system is primarily controlled by stress the trends of the strain profiles calculated in the MICP case correspond well to the trends of the elastic modulus (Fig. 10, right).
In Fig. 12a we compared the evolution of pore pressure at the bottom and the axial displacement at the top of the sample during the whole treatment in the cases with MICP and without MICP treatment. In both cases the injection of chemical solution from the bottom induces a rise of pore pressure at the injection source. Initially, the mechanical loading at the top causes a small compression in the sample. At the beginning of chemical injection phase, the generation of pore pressure in both cases causes a slight extension at the top of the sample. In the MICP case, the strong reduction in the permeability, especially in the later chemical injection stage (after 60 hours) causes a significant increase of pore pressure which leads to a remarkable increase of the displacement. As the pore pressure response during MICP can cause additional deformation, a fully Fig. 11 Axial profiles of the effective stress (left), strain (right) computed at different times during the mechanical tests in the simulation cases with and without consideration of MICP effects Fig. 12 Computed axial displacement at the top and pore pressure at the bottom of the sample: a during the whole treatment; evolution of b axial displacement and c pore pressure during the mechanical testcomparison between the cases with (solid lines) and without (dashed lines) consideration of the MICP effects coupled BCHM model is necessary especially for the simulation of MICP treatment under high injection pressure, by which large deformation, thus hydraulic fracturing in the soil could be caused by the high pore pressure. With the generation of hydraulic fracturing, the flow path and thus reactive mass transport behaviour of MICP can be changed. Contrary to the MICP-treated case, without consideration of MICP because of no permeability reduction, pore pressure and displacement remain constant during the chemical injection. In the following mechanical test a much stronger response in pore pressure and displacement is induced by the mechanical loading. The magnified evolution of pore pressure and axial displacement during the mechanical test can be seen in Fig. 12b, c. Due to the overall enhanced mechanical stiffness compared to the results without MICP effects, the axial displacement at the top of the column is remarkably reduced in the MICP case (Fig. 12b). Due to the small displacement under the same mechanical loading, a weaker consolidation process is expected in the MICP-treated sand than in the untreated sand. As a consequence, as can be seen in Fig. 12c, during the undrained mechanical loading the pore pressure increases in the MICP-treated sand that is much lower than that in the sand without MICP treatment. In the practice, it is common sense that a sudden excess of pore pressure could lead to mechanical failure. Thus, from the computational results it can be addressed that apart from the most recognised MICP advantages (e.g. reduction of permeability, improvement in mechanical stiffness and strength) one more particular benefit from MICP can be inferred regarding the perspective of hydraulic mechanical coupling. To be specific, MICP can be utilised to hinder the strong pore pressure response in the soil during consolidation thus to avoid the possible related mechanical failure.

Conclusions
The microstructural effects on the macroscopic behaviour and the fully BCHM couplings are two relevant issues of MICP in soil that are currently less investigated. In the present study, we applied a simplified GM method to describe the influences of the microstructural aspects on the reaction kinetics and mechanical behaviour of the MICPtreated soil in a macro-continuum model. To address the multi-field problems in MICP, we defined four subsystems, respectively, for the biological (B), chemical (C), hydraulic (H) and mechanical (M) fields. The couplings among the BCHM fields are considered in the corresponding constitutive equations and the additional defined coupling terms in each field equation. The final equations describing the system behaviour and the coupling effects are implemented in the open-source FEM code OpenGeoSys (OGS). The biochemical reactions and the mass transport are solved in a monolithic way, whereas they are solved with the hydraulic process and the mechanical process iteratively in a staggered scheme.
The developed model is applied in three numerical examples. The first one aims at analysing the microstructural effects on the temporal evolution of the reaction behaviour, by which the experiment conducted in [45] is simulated. With the introduction of the microscopic related parameters and cross-scale functions, the observed different MICP behaviours in samples with similar initial porosity but different microstructures can be captured. In the second case the tests reported in [9] are simulated, by which the relationship between the macroscopic and microscopic property changes is investigated. The comparison of the computational results with the experimental observations indicates that our model can well reproduce not only the macroscopic behaviour changes by MICP but also the evolution of the specific reaction surface. In the last case a scenario analysis of a fully BCHM MICP test is performed. A reasonable prediction of the relevant chemical and hydro-mechanical response in MICP during different treatment phases is obtained. The numerical analysis indicates that depending on the injection strategies the soil with initial homogeneous properties could become inhomogeneous. This inhomogeneity in space and time can be captured in the model only when the reaction kinetics, mass and fluid transport, and deformation processes are solved in a coupled manner.
In conclusion, the present model provides a step forward for better understanding of the BCHM interactions and the micro-related macroscopic property changes in MICP. Since all relevant chemical components are addressed, the present reaction model can be easily extended to consider more complex treatment conditions (e.g. insufficient calcium, unfavourable pH condition). Moreover, the microrelated parameters, i.e. the adopted simple GM method in the model, could be further refined to capture more complicated microstructure-related mechanisms, e.g. definition of randomly distributed solid grains and calcite to consider the dependency of crystal nucleation and growth on the different grain sizes, surface roughness, mineralogy texture, etc. Besides, the microstructural effects on the plastic or damage behaviour of the MICP-treated sand could be addressed with the development of proper cross-scale functions. Thus, based on the current coupling concept further development of a coupled elastoplastic BCHM model would be possible.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.