Critical occlusion via biofilm induced calcite precipitation in porous media

A model for biofilm induced calcite precipitation with pressure driven flow is presented at the scale of a single pore within a porous medium. The system, an extension of previous work (Zhang and Klapper 2010 Water Sci. Technol. 61 2957–64, Zhang and Klapper 2011 Int. J. Non-Linear Mech. 46 657–66), is based on a mixture model including biomaterial, mineral, and water with dissolved components. Computational results suggest the possibility of critical occlusion in the sense that there is a distinguished trans-pore pressure head such that for pressure drops below this level, pore clogging occurs relatively quickly while for pressure drops above, clogging occurs after much longer times if at all. Beyond its relevance to engineered biofilm applications, this phenomenon is suggestive of the subtleties of embedding simple biofilm models in larger media.


Introduction
Biofilms are collections of microorganisms held together by sticky, self-secreted matrices of polymers and other substances [12,13]. Wherever reliable (even if scarce) sources of water, nutrients, and chemical or light energy are present, biofilms are also often to be found. This is in large part due to the wide metabolic functionality distributed through the microbial world and the consequent ability for microbial communities to use this potential to efficiently exploit many different environments. Biofilms seem to be especially effective structures for housing these microbial factories.
The same flexibility that allows microbial communities to function efficiently in many different environments also presents interesting opportunities for engineered systems. Indeed, in fact, engineered microbial communities including biofilms have long been used as, essentially, chemical catalysts in for example waste treatment and remediation contexts. Less common as yet is application of microbes to solving mechanical problems. One particularly intriguing possibility is use of biofilm induced mineralization as a form of structural stabilization or as a permeability reduction method. The underlying idea, described in more detail below, is to induce pH increases via microbial activity with the consequent result of mineral precipitation, chiefly carbonates, see figure 1 for an example. Microbes can penetrate into small-scale spaces of porous materials and then initiate biofilm growth, and if these biofilms can be coaxed into inducing mineral precipitation then the resulting deposits can effectively fill pores, interrupting through-flow and/or cementing the given material. Possible applications include soil stabilization [1], construction of subsurface barriers against pollutant spread or reservoir leakage (including release of sequestered carbon dioxide) [3], and repair of concrete and stone structures (including damage due to subsurface hydraulic fracturing) [5].
Biofilm-induced precipitation within a porous material is a complicated phenomenon, however, involving advective and diffusive transport, biofilm biology and mechanics, in addition to the precipitation process itself, all of these occurring within the complicated pore structure of an opaque material. Observational studies are challenging, even at the lab scale [16]. Thus a computational model capable of testing and generating hypotheses, particularly at the pore scale, can be useful. The authors have developed such a model previously [21,22] and here use this model to investigate pore occlusion. In particular, we consider the trade-offs resulting from increased flow rates that, on the one hand, increase transport in of constituents of mineral formation but, on the other hand, at the same time increase the transport out of mineral products themselves. The results of this competition are not at all clear a priori, as in addition to the advective transport, there are the additional factors of differing diffusion rates for different chemical species, electrodiffusion, and changes in shear stresses induced by precipitated mineral, among other things.
In addition to the model used here, there are other models of biofilm-induced pore clogging available [7,17]. These models do not however study the precipitation process at the single pore scale where it is possible to include the fundamental chemistry and physics in fair detail. Our aim here is to do so in order to study the actual individual pore occlusion process and, at the same time, to provide a framework for upscaling.

Mixture model
We model a single pore within a larger lab or environmental pore system. For computational feasibility we consider a two dimensional (2D) domain with parallel walls and inflow on one end and outflow on the other. Biofilm, containing ureolysis-capable microbes, is present on the walls, and is supplied with urea and calcium ions from the flow. Through processes described below, calcite is produced and collects into a mineral phase.
We describe the entire system using a mixture model platform and refer to [21,22] for details. The model, briefly outlined here, consists of three phases-solvent, biomaterial, and calcite-with corresponding volume fractions ϕ ϕ , s b , and ϕ c respectively, as well as a number of dissolved (in the solvent) chemical concentrations. The latter satisfy advection-electrodiffusion-reaction equations stated below. Equations for the former, the material phases, together with those for a limiting substrate concentration c (dissolved in the solvent), an advective velocity v, and hydrostatic pressure p are , so equations are needed for only two of the volume fractions. Also, for time intervals of interest here, little biofilm growth is supposed so that we can set c = 0 and omit the limiting substrate equation (3).
In equations (1)-(5), the, the functional f is a chemical free energy density described below, Λ c and Λ b are mobility parameters for the calcite and biofilm phases respectively, and calcite, biofilm and solvent phases. The solid calcite phase is approximated to be an extremely viscous fluid, sufficiently viscous so that it would not deform significantly under the flow within the computation times used. Biofilm is also modeled as a very viscous fluid; in fact, for sufficiently slow flow time scales (as is the case here), this is reasonable [18]. The source (4), describing calcite precipitation rate, is defined by a commonly used quadratic rate law [19,23] where k p is the calcite precipitation rate constant and C vol is the inverse molar density of mineralized calcite [14]. S, defined below, is the calcite saturation state and S crit is a critical supersaturation parameter. The energy density 1 ln ln ln , indicates that biomaterial prefers to be hydrated, i.e. its contribution is negative for ϕ ϕ . Other terms and parameters in (6) as well as the extra parameters in the equations for ϕ b and c do not play an important role here; discussion of them can be found in [20].

Chemistry
The system chemistry is Ureolysis, reaction (7), is, effectively, catalyzed by microbes in the biofilm supplied with urea ( ions. Then, in the presence of soluble calcium ions, calcite (CaCO 3 ) precipitation can occur. Reactions (7) and (10) occur on the 'experimental' time scale and are included in the reaction-electrodiffusion-advection equations for their constituents in a standard manner, see below. Reactions (8) and (9), however, equilibrate rapidly compared to time scales of interest and can thus be assumed to be at quasi-equilibrium. The methodology for computing their equilibrium distribution is described in appendix A.  are in quasi-equilibrium with values determined using the methodology provided in appendix A. The other chemical species obey mass transport equations including effects of convection, electrodiffusion and chemical reactions (7) and (10). We assume that the dissolved chemical species are only present in the solvent phase. In that case, [ ] Here R i is the reaction term for species i, D i is the diffusivity of species i, ζ is the electric potential (discussed below), and a is a units-related conversion factor [10].
The reaction functions R i are defined as follows.
u r e a 1 1 where k urea is the urea hydrolysis rate coefficient and γ 1 is the urea activity coefficient (computation of activity coefficients and their importance for results is discussed in appendix A); the presence of ϕ b is a consequence of the fact that ureolysis is catalyzed by bacteria and only occurs in the presence of biofilm. Next, 1 due to reaction stoichiometry. Species 4, 5, 6, 9, 10 are at quasi-equilibrium, so that where K SO is the equilibrium calcite solubility product, S is the calcite saturation state, and k p is the calcite precipitation rate constant. Variability of diffusion constants between differently charged species leads to charge separation which in turn induces an electric field with potential designated ζ.   5, 6, 7, 8, 9, 10 E is the index set for the charged chemical species. See [22] for details.

Numerical methods
Details of the numerical methods including parameter values can be found in [20,22] as well as in appendix B. Briefly summarizing, the equations are non-dimensionalized and then solved by a finite difference method on a uniform rectangular mesh in space chosen so as to satisfy the CFL condition in scaled variables.

Advection-reaction-diffusion equations for quantities
i are solved using a Crank-Nicolson scheme except that the advection term is treated by a fully implicit upwind scheme, again as previously described [8,20]. The coupled continuity equation (1) and momentum transport equation (2) are solved using a velocity-corrected projection scheme as previously described [8,20]. Unlike before, though, constant pressure head conditions (rather than constant flux) are used here following [9], see appendix B for implementation details. Constant flux conditions would not allow blocking of the flow (which would result in zero net flux) and would instead force pressure changes to reach unrealistically large values in order to force the required flux through the system.
Computations are performed on a 2D rectangular domain of size 1 mm 2 with flow going from left to right, although this computational domain is considered to be embedded within a longer domain of length = L 9 mm for purposes of maintaining an imposed pressure gradient, see appendix B. The flow is then driven by an imposed constant pressure head ΔP over the longer domain. Initial velocity is set to be Since the computational times considered here are relatively short, we set c = 0 and thus allow no growth of biomaterial.

Results
Perhaps the most fundamental issue one might consider with respect to biofilm-induced mineralization in porous media is that of expected loss of permeability. The most effective method for permeability reduction is to in fact entirely shut pores off through mineral occlusion.
We present here then a model of pore occlusion and observe, at least in the model, an unexpected and significant phenomenon, namely the existence of a pressure drop below which clogging occurs rapidly and above which clogging occurs much more slowly if at all. We consider a model pore with flow through the pore driven by a constant, imposed pressure head. Over time, biofilm-induced calcite precipitation has the potential to narrow the region of the pore available to through-flow even to the point of effectively occluding it altogether. Precipitation results from sufficient co-occurence of precipitant components; previous study has indicated that the rates of accumulation can be limited by transport of raw materials, calcium ions and urea particularly, into the pore region and then unto activity sites [22]. Transport rates are largely controlled by advection rates which are in turn largely controlled in this system by the magnitude of the imposed pressure head, at least at early times before occlusion occurs. Hence we study the effect of varying the magnitude of the imposed pressure head. Based on the importance of material transport, one might expect larger pressure heads to result in increased transport of raw materials into activity centers and hence to faster clogging, but, on the other hand, both precipitation itself, see equation (10), and mineral accumulation, see equation (4), are relatively slow processes operating on time scales that may well be slower than the advective time scale so that there is also the possibility that increased advection can remove material from the pore before it is able to precipitate and/or accumulate.
To investigate, we performed a number of computations according to the specifications outlined in section 2.  figure 3, do we see mineralization penetrate significantly across streamlines at later times and into the bulk fluid region. This is a slower process; at t = 8 h, precipitate is seen mainly in the biofilm and downstream regions, and only later, at t = 16 h, do we see noticeable precipitate migrating into the main channel and significantly diminished flow (two orders of magnitude smaller than initially). The occlusion step has characteristics of an instability: diffusion into the main channel and resulting clogging slows the flow further which in turn speeds the occlusion process with, in this case, nearly complete blockage by t = 24 h. Note that precipitate is even able to propagate in the upstream direction. . Initial biomaterial volume fraction (a) and velocity (b): note that initial velocity changes rapidly at very early times to accommodate the presence of biomaterial because of high biomaterial viscosity. Biomaterial volume fraction (c) and velocity (d) at t = 48 h: biomaterial is sheared by flow, note that maximum velocity has decreased somewhat but through flow still occurs. Calcite volume fraction (e): note that most precipitated calcite is on, near, or downstream of the biomaterial interface. Ionic strength (f): note ionic strength is largest in the vicinity of the biofilm-bulk fluid interface which, also, is where most ureolytic activity occurs. . Calcite volume fraction (a) and velocity (b) at time t = 8 h: note that some reduction in maximum velocity can be observed, though most mineral accumulation is seen within the biofilm and downstream. Calcite volume fraction (c) and velocity (d) at t = 16 h: note that the pore is mostly occluded and maximum velocity has decreased significantly. Calcite volume fraction (e) and velocity (f): occlusion is essential complete.
shows contours of the corresponding velocities u. Observe that a change from occluding to nonoccluding behavior appears to occur for an imposed pressure head producing an initial maximum velocity somewhere between = × − . Comparing to figures 4(c) and 5(c), we observe that there is somewhat more calcite precipitation at t = 240 h as compared to t = 24 h, but the maximum values of u are still of the same order of magnitude. So while it may be possible that the pore will eventually fully occlude given enough time, behavior is qualitatively different than that seen for slightly smaller imposed pressure drops. (It could be argued that biofilm growth, which was neglected here, may be significant over a 240 hr period-however, the purpose of this longer run is only to verify that solution behavior is indeed qualitatively different than for slightly smaller initial pressure head.) To further illustrate, we plot maximum values of u at t = 24 h for a number of different imposed pressure drops (again translated into maximum initial velocity u 0 ) in figure 7. The critical occlusion transition is particularly noticeable in figure 7(b) with the vertical axis in log scale.
To test the robustness, we searched for the same critical transition after varying key model parameter choices, see figure 8 for two examples.  through the tunnel increases and the ureolysis (catalyzed by the presence of biofilm) proceeds at a slower rate. We observe that critical occlusion occurs at a smaller u 0 , one which allows longer residence time for the reaction to occur. Figure 8  due to the faster reaction caused by the larger k urea . We remark that k urea is a particularly interesting parameter to vary for the following reason. According to [15], k urea can vary in a wide range (from − 0.069 day 1 to − 0.99 day 1 ) due to change of experimental conditions such as temperature, while k p (the calcite precipitation rate) and S crit (the critical supersaturation parameter) vary in a much narrower range. Overall, our computations indicate that the critical  occlusion phenomenon is robust, though the actual critical value of u 0 (or ΔP) depends on parameterization values.

Toy model
We consider the following toy model which is highly simplified but illustrative of the underlying physics. Suppose that precipitation reactions are occurring in a specific region V of space, with lengths normalized so that the volume of V is one. Flux into and out of this region, driven by a pressure drop across V of magnitude ΔP, occurs at rate D. A limiting substrate, with concentration in V denoted by S(t), is transported into V (via the fluid flux) from an outside supply held at fixed concentration S 0 and is transported out of V with concentration S(t). While in V, the substrate is converted to precipitate at rate rS with rate constant r. Precipitate is also lost by being washed out of V by the outbound fluid flux. Volume fraction of precipitate in V is denoted by X(t). Altogether, then, the model iṡ  (14), S follows from the first equation, so we need only consider X. Sketching the left-and right-hand sides of (14) as functions of X, solutions are indicated by intersections of the two curves. We see two possible diagrams, see figure 9, left and right panels, corresponding to relatively large and small ΔP, respectively. As Δ → P 0, ( ) D X r also approaches 0, so the blue curve approaches the horizontal line Y S 0 . Thus the representative diagram changes from figure 9, left panel, to figure 9, right panel, at a critical value of ΔP. Note both panels include the occlusion solution = X Y S 0 . The difference is that for the large ΔP case, a second, non-occluding solution occurs (note though that if D(X) decreases very rapidly as a function of X, then it is possible that no non-occluding equilibrium exists). It is easily checked that the non-occluding solution is the stable equilibrium in the case that it occurs. Conversely, when there is no non-occluding equilibrium, the occlusion solution is stable.
The underlying message of the toy model is that the non-dimensional quantity D r, measuring the ratio of reaction time to residence time (roughly speaking, a Damkohler number), determines the ability of the system to clear precipitate faster than it accumulates. We in particular note that there is a critical value of ΔP for which D r is small enough that precipitate cannot be cleared, in which case occlusion results. There is not a precise way to define D or r in the full system, but nevertheless the predictions of the toy model are, even beyond the ΔP dependence, qualitatively consistent with observations, see for example full model computations shown in figure 8 where certain key parameters are varied from previous computations. As described before, in figure 8 left panel, occlusion results are shown for initial biofilm colony radius size reduced from 0.2 mm previously to 0.1 mm, increasing the system flux and decreasing its overall reaction rate (note that biofilm is the effective enzyme for ureolysis) as compared to previous computations, i.e. increasing D and decreasing r, see figure 7 right panel. Note that the new system requires smaller u 0 and hence smaller ΔP to occlude. In figure 8 right panel, occlusion results are shown for urea hydrolysis rate coefficient k urea increased from − 0.2 day 1 previously to − 0.4 day 1 , increasing the reaction rate as compared to previous computations, i.e., increasing r, again compare to figure 7 right panel. Note that the new system now requires larger u 0 and hence larger ΔP to occlude.

Discussion
Clearly, biofilm-induced mineralization is a complex process involving chemistry (distribution of dissolved nitrogen and carbon species), physics (advection and electrodiffusion), biology (microbially modulated ureolysis) and thermodynamics (mineral precipitation and accumulation), all of which operate on their own time scales, yet still interact in essential ways to produce  (14)). The y-intercept of the blue curve increases with decreasing ΔP. Occluding solution (where D = 0 and both right-and left-hand sides of (14) equal Y S 0 ) is available in both cases but unstable for sufficiently large ΔP. a result. Using an improved version of a previously presented mixture model for biofilm induced calcite precipitation including now effects of ionic strength and activity on chemical reactivity, numerical simulations illustrate an important mineralization phenomenon, namely pore occlusion. In particular, results suggest that as the constant pressure head is increased, the system undergoes a qualitative behavior change from relatively rapid clogging to slow, or nonclogging behavior. That is, we seem to observe a critical occlusion phenomenon. We vary, essentially, initial flow rate, a relatively controllable system component in many engineering applications. Flow rate affects time scales for transport in and out of mineral precursors and components. When advection is the dominant transport mechanism relative to diffusion, i.e. when imposed pressure head ΔP is relatively large, then calcite precipitation mainly occurs near the biofilm-solvent interface, the location where ureolysis happens. In that case precipitation processes can carry through before necessary components are transported away by the flow. In contrast, when ΔP is relatively small, diffusive transport increases in importance leading to an instability where precipitation occurs out in the main channel, further slowing the flow and so further reducing the importance of advective transport.
This apparent instability is potentially important in the context of engineered biofilminduced mineralization barriers within porous media. In such materials, access is difficult (indeed, this is a prime motivation for using biofilms) so any exploitable control on the location and timing of permeability reduction can be useful. We are able to suggest that input flow rate may provide such an ability. Beyond such possible utility, critical clogging provides an illustration of another interesting and subtle way in which biofilms can act as reactive materials and how they can influence and be influenced by their physical environment. This is also a warning: large-scale porous media models often treat biological effects as linear or nearly so. We propose here at least one instance, and an important instance, where this may not be the case. where N is the total number of chemical species present. If I is small, i.e. the solution is sufficiently dilute, then γ ≈ 1 i so that activity and concentration are comparable. Figure  as a function of I on the left, and plots activity coefficients for monovalent and divalent ions as a function of I on the right. For concentrations observed in this study, activity coefficients for neutral and monovalent ions are close to one, but activity coefficients for bivalent ions can be noticeably smaller than 1.
Fast reactions given by (A.1) are represented by the following set of mass action equations  Combining with water dissociation balance and charge neutrality balance constraints  (5) Ifˆ< otherwise, go to step (2).We remark that special care has to be given in step (3)  numerically; details can be found in [22].

A.2. Effect of ionic strength and activity coefficient
Computation of activity coefficients complicates and slows the overall numerical procedure. To justify its inclusion, then, we show here the relative importance of its inclusion. Figure 2(f), from a panel presenting simulation results with constant pressure head ΔP value such that the initial parabolic velocity profile has a maximum value of − − 10 m s 3 1 , shows contours of the ionic show contours of activity coefficients γ i for monovalent and divalent ions at t = 48 h for the same simulation, with both sets being less than one. The activity coefficient for divalent ions in particular is significantly less than one, especially in the important active region near the biofilm-bulk fluid interface. Since calcium and carbonate ions, which are central to calcite precipitation, are both divalent, then we can anticipate that the activity versus concentration difference has significant impact on precipitation rates.
To compare, then, figure B1 shows results for the same computation as shown in figure 2 except using concentrations rather than activities, i.e., γ = 1 i for all species in the computations illustrated in figure B1. In comparison to figure 2, the pore is noticeably more occluded, see figure B1(c), and maximum velocity is approximately halved, see figure B1(d). This enhanced occlusion is not unexpected since, although the activity coefficient for neutral ions is very slightly larger than 1, those activity coefficients for charged ions, particularly divalent ones, are significantly less than 1. Hence, by assuming γ = 1 i , we effectively increase the available reaction and precipitation components so that precipitation can occur at a faster rate.

Appendix B. Constant head boundary conditions
We follow the approach from [9] used to compute pressure driven flow, outlined as follows. Consider a 2D domain with pressure heads at each end with an interior section of the domain partially obstructed by the biofilm, see figure B2. The computational domain consists of the interior, biofilm section only. Let L be the total length of the domain, a be the distance between the top and bottom of the domain, and ΔP be the constant pressure drop across the whole domain prescribed by the pressure heads. The domain is divided into three regions: upstream from the obstructed region, the obstructed region, and downstream from the obstructed region. Let L L , , and ΔP 1 , ΔP 2 , and ΔP 3 be the pressure drops in each region so that Δ Δ Δ Δ + + = P P P P 1 2 3 . Here we set the domain to be 9 mm in length and 1 mm in height, where the computational domain is the 1 mm by 1 mm square in the middle region. Thus L = 9 mm, = = L L 4 mm  where ( ) u x y , is the x-component velocity. Note that Q is independent of x due to the incompressibility constraint. Since the flow is pressure driven, Q need not be constant in time, though, even for constant ΔP as for example calcite precipitation can cause increased obstruction with resulting decrease in flux.  Following [9], a relationship between the flux Q and the pressure drop ΔP 2 is derived and subsequently used for setting a boundary condition for the flow equation as follows. Assume that the flow is parabolic in the regions upstream and downstream from the obstructed region the flow, in which case = ( )