An immersed boundary-lattice Boltzmann model for biofilm growth in porous media

In this paper, we present a two-dimensional pore-scale numerical model to investigate the main mechanisms governing bioﬁlm growth in porous media. The ﬂuid ﬂow and solute transport equations are coupled with a bioﬁlm evolution model. Fluid ﬂow is simulated with an immersed boundary–lattice Boltz-mann model while solute transport is described with a volume-of-ﬂuid-type approach. A cellular automaton algorithm combined with immersed boundary methods was developed to describe the spreading and distribution of biomass. Bacterial attachment and detachment mechanisms are also taken into account. The capability of this model to describe correctly the couplings involved between ﬂuid circulation, nutrient transport and bacterial growth is tested under different hydrostatic and hydrodynamic conditions (i) on a ﬂat medium and (ii) for a complex porous medium. For the second case, different regimes of bioﬁlm growth are identiﬁed and are found to be related to the dimensionless parameters of the model, Damköhler and Péclet numbers and dimensionless shear stress. Finally, the impact of bioﬁlm growth on the macroscopic properties of the porous medium is investigated and we discuss the unicity of the relationships between hydraulic conductivity and bioﬁlm volume fraction.


Introduction
Biofilm is a community of one or more species of bacteria or other microorganisms (fungi, algae, yeasts) associated irreversibly with a liquid or solid surface (water, biological tissues, solid substrates located in marine environments or freshwater) ( Potera, 1999 ) and enclosed in a matrix of polysaccharide ( Kalmokoff et al.,20 01;Prakash et al.,20 03;Smith,20 05 ).Biofilm is composed primarily of microbial cells and extracellular polymeric substance (EPS).The EPS plays various roles in the formation and structure of the biofilm and also, protects the cells by preventing the access of antimicrobial and xenobiotic compounds and confers protection against environmental stresses such as UV radiation, pH shift, osmotic shock and desiccation ( Chmielewski and Frank, 2003; Decho, 20 0 0; Mah and O'Toole, 2001;Wingender et al., 1999 ).
The study of biofilm growth in porous media is currently attracting interest for environmental applications such as bioremediation of contaminated sites ( Al-Bader et al., 2013;Singh et al., 2006 ) or development of bio-barriers for groundwater protection ( Huang et al., 2011;Seo et al., 2009 ).One of the major challenges common to these decontamination methods is due to excessive accumulation of biofilm in some pores (particularly in the vicinity of the pollutant source or within the bio-barrier).This accumulation of biofilm, indeed, leads to a decrease of the pore space which causes a reduction in porosity and permeability of the system while increasing the hydrodynamic dispersion ( Kone et al., 2014 ).Eventually, the clogging mechanism may change the transport of dissolved organic compounds and therefore strongly influence the success of the application of cleaning methods for polluted groundwaters.Accurate knowledge of the formation and evolution of the biofilm over time is crucial to control its development.
In recent decades, a large experimental effort has been devoted to study the biofilm structure and its change during growth.In parallel, there has been a growing interest in the numerical modeling of the involved processes to reproduce the observed experimental results and confirm or refute existing hypotheses.Historically, many models have been developed since the 90s.If the first numerical models of bacterial growth were limited to considering the availability of the nutrient, they become more sophisticated to include the influence of shear forces ( Alpkvist and Klapper, 2007;Morgenroth and Wilderer, 20 0 0 ), bacterial extinction ( Horn and Hempel, 1997 ) and adhesion mechanisms ( Ebigbo et al., 2010;Kapellos et al., 2007 ).Recently, they succeeded to describe more complex bio-physico-chemical processes such as quorum sensing ( Frederick et al., 2011 ), physiological state of the cells http://dx.doi.org/10.1016/j.advwatres.2017.06.009 0309-1708/© 2017 Elsevier Ltd.All rights reserved.( Kapellos et al., 2007 ) and the effects of competition between bacterial populations ( Ebigbo et al., 2013;Lackner et al., 2008;Noguera et al., 1999;Pfeiffer et al., 2001;Picioreanu et al., 2004a ).The interest for biofilm-related issues in porous media has motivated the application of these models to such complex structures.A variety of models have thus been derived since the pioneering works of Suchomel et al. (1998) and Dupin et al. (2001) based on pore networks.If sophisticated 2D and 3D pore-scale descriptions have been achieved ( Ebigbo et al., 2013;2010;Kapellos et al., 2007;Knutson et al., 2005;Peszynska et al., 2016;von der Schulenburg et al., 2009 ) and have given additional insight on the bio-chemical and physical mechanisms that drive the biofilm dynamics, important questions on how the morphology and distribution of biomass may affect the hydrodynamic properties of porous medium remain to be addressed.
The aim of this work is to present a mathematical model capable of predicting the temporal evolution of biofilm growth in porous media.The model was validated by qualitative comparison with benchmark cases from the literature.First, we considered the case of a biofilm developing on a solid flat surface in hydrostatic conditions for the limiting growth regimes, i.e., transport or kinetics limited regime.Then, we have imposed a lateral flow to verify the influence of hydrodynamic conditions on the biofilm pattern.Finally, we have applied our model to a complex porous medium in order to observe the role of bacterial growth on the hydrodynamic properties of the medium.In particular, the relationship between porosity and permeability has been discussed as a function of the biofilm growth regime.

Literature survey on biofilm growth models
Overall, there are two possible approaches to modeling biofilm growth in porous media ( Wang and Zhang, 2010 ), namely the discrete (or individual-based) approach (including multi-agent and cellular automata models) and continuous approach.In the multiagent model, the fundamental entity is the agent representing the individual bacteria ( Kreft et al., 1998 ) or bacterial colony.Each agent grows by consuming the substrate and divides when the volume reaches a critical threshold.The newly-formed agent remains in contact with the parent agent in a direction randomly chosen.Then, biomass propagation occurs by pushing agents soon as they become too close ( Alpkvist et al., 20 06;20 06;Batstone et al., 2006;Kreft et al., 2001;Picioreanu et al., 2004b ).The first multi-agent model for the growth of bacterial colonies, called Bac-Sim, was developed by Kreft et al. (1998) and extended to a multispecies simulation by Kreft et al. (2001) .Multi-agent models have been successfully applied to describe and optimize various applications on biofilm reactors ( Batstone et al., 2006;Matsumoto et al., 2010;2007;Picioreanu et al., 2005;Xavier et al., 2007;2005a ).Control strategies based on biofilm disruption were also studied ( Xavier et al., 2005b ).Using the multi-agent model is relevant if the interest is focused on interactions between individuals ( Hellweger and Bucci, 2009;Laspidou et al., 2010 ).Indeed, this model offers the possibility to work at the bacteria scale and facilitates taking into account the cellular mechanisms.Empirical knowledge of the system and of its biology may be used directly to assign basic mechanisms of multi-agent model (see a recent review on advantages and limitations of individual-based models in Hellweger et al., 2016 ).However, it requires very long computation time particularly in a 3 D configuration ( von der Schulenburg et al., 2009 ).Conversely, cellular automata model divides the space into grid cells, each cell being considered as individual entities that can take several possible states (e.g., solid, fluid or biofilm).The state of a grid cell is determined by the previous state of its neighbors.In this approach, biofilm is considered as a unit cell containing EPS, water and bacteria.The first applica-tion of cellular automata model to biofilm growth was presented by Wimpenny and Colasanti (1997) .This model was based both on the work of Colasanti (1992) and the diffusion-limited aggregation (DLA) model used by Schindler and Rataj (1992) ; Schindler and Rovensky (1994) and Fujikawa (1994) .Ever since, several models have been developed using this approach ( Hermanowicz, 1998;20 01;Khassehkhan et al., 20 09;Laspidou et al., 2012;Laspidou and Rittmann, 2004;Laspidou et al., 2005;Liao et al., 2012;Picioreanu et al., 1998b ).They were able to produce complex growth patterns.Two cellular automata algorithms are particularly used in the literature.Picioreanu et al. (1998b ); 1998c ) divide into two equal parts the biomass of the cell when the threshold value is reached to mimic cell division.This biomass is then redistributed among the four adjacent cells (for a 2 D configuration) until reaching the closest unoccupied cell.These so-called pushing algorithms can lead to a spatial discontinuity of the biomass distribution, due to the division performed.To make them continuous, another approach has been further developed by Noguera et al. (1999) and Knutson et al. (2005) .In this approach, only the excess of biomass is redistributed.Cellular automata models yielded good agreement with experimental results.They have the advantage of easily taking into account the different cellular mechanisms and facilitate implementing boundary conditions and complex system geometry.However, they are less suitable for multi-bacterial biofilms while preferential orientation of the Cartesian grid may cause unrealistic growth.Unlike both previous formulations, continuous models are purely deterministic.This approach relies on the idea that bacterial growth causes the development of a pressure field inside the biofilm ( Eberl et al., 2001 ).This pressure field constitutes the driving force of a very viscous flow and causes slow displacement of the biofilm ( Wanner and Gujer, 1986 ).A similar idea has been used for the first time by Wanner and Gujer (1986) to express the velocity of the fluid/biofilm interface due to biofilm growth in a onedimensional model.The average velocity of biomass is thus determined by mechanical laws, and more precisely by the Euler equations of fluid dynamics.Continuous models are able to predict the changes in biofilm thickness, the transient dynamics of bacterial community as well as the spatial distribution of microbial species and chemical compounds in the biofilm.Various numerical methods can be used to describe the displacement of interfaces, such as the level set method ( Alpkvist and I.Klapper, 2007 ) and Volume-Of-Fluid (VOF) method ( Ebigbo et al., 2013 ).Continuous models have the advantage of being purely deterministic, rigorous, based on conservation laws and well-known PDEs.However, this type of model gives only a macroscopic insight into the simulated domain.They are therefore particularly suited to large scale systems but they are not valid when the size of the representative elementary volume is close to the bacteria or micro-colony scale.
Based on this in-depth review of biofilm models, we opted for the most flexible approach i.e., the cellular-automaton method.We must note that a major shortcoming of this approach is related to the elementary size of the unit cell used by the cellular automaton.This cell size constrains the minimal change in biofilm thickness during biomass growth and may dramatically impact model predictions at the first stage of colonisation or when the pores are nearly plugged.Hereafter, this drawback will be offset by the use of efficient non-conforming numerical methods.We refer the reader to Benioug et al. (2015) for a comparative analysis and accuracy assessment of such methods applied to flow and reactive transport in porous media.In the present study, we will use the Lattice Boltzmann Method coupled with the Immersed Boundary Method (IB-LBM) for calculating groundwater flow and a VOF-like method for solving solute transport ( Benioug et al., 2015 ).The cellular automaton model will take advantage of these numerical efforts to work on entities smaller than the mesh size.

Model description
Simplifying assumptions are often required to simulate the main mechanisms describing the complexity of biological and physico-chemical phenomena involved in the development of bacterial populations in porous media (e.g., nutrient availability, mechanical strength of the biofilm and influence of shear forces on the detachment of sessile bacteria, competition for the substrate and symbiotic effect between different bacterial populations, quorum sensing, evolution of the physiological state of cells) ( Brackman et al., 2009;Davies et al., 1998;Vega et al., 2014;Zhang et al., 2009 ).We do not pretend here to describe all the features of this complexity and consequently, we will focus on a specific class of problem where we will more particularly investigate the bioclogging issues in porous media.In our case, the following assumptions will be made: • The solid matrix is supposed non deformable and the porous medium is fully saturated.• A single bacterial species is present in the medium and only under sessile form (mono-species biofilm).• Biofilm is considered as a continuous and homogeneous medium ( Wood and Whitaker, 1998 ), although this one is inherently characterized by a structural heterogeneity (bacterial cells, EPS, nutrients) and an intrinsic biodiversity (bacteria, protozoa).• Metabolic reactions that are involved in the bacterial growth process only depend on a single solute, i.e., the organic molecule serving as a carbon and energy source and acting as the electron donor in the redox process.Other substances (electron acceptor, mineral salts, amino acids) which are necessary for the survival of bacteria are assumed in large excess in the medium.In those circumstances, the biodegradation kinetics is described by Monod kinetics ( MeeGee et al., 1970 ).• The nutrient is transported by convection -diffusion in the fluid phase and only by diffusion in the biofilm phase.In other words, macro-porosities within the biofilm where convection could occur are assumed sufficiently large to be explicitly discretized and are therefore attached to the fluid phase.• The concentration of the chemical species dissolved in the fluid phase is sufficiently low so that solute transport in the water does not affect the initial physical properties of the fluid, e.g., its dynamic viscosity and its density.In this manner, the flow and transport equations can be uncoupled.
Finally, we will use the time-scale separation assumption for the different processes involved (flow, nutrient transport and biomass growth) by adopting a quasi-stationary approach.Indeed, the characteristic relaxation time of the velocity field is very small compared with the characteristic time of transport, which itself is generally very small compared to the characteristic time of biofilm growth ( Picioreanu et al., 1998a ).This assumption is valid in most cases ( Skowlund, 1990 ). Consequently, we solve the flow and species transport equations in steady-state conditions.The temporal evolution is performed through the bacterial growth equation (and the associated cellular automaton model).
Based on the assumptions above, the problem at the microscopic scale can finally be stated as follows.

Liquid flow
In the fluid phase ( f ), we consider the steady-state flow of an incompressible Newtonian fluid, at low Reynolds number so that the mass balance and momentum conservation equations are described using the following equations: with non-slip boundary conditions for the fluid-solid and fluidbiofilm interfaces, respectively A fs and A fb , where u , p, μ f , ρ f are respectively the velocity vector, pressure, fluid dynamic viscosity and fluid density.A ij represents the interface between the i and j phases ( i, j = f (fluid), b (biofilm), s (solid)), and n ij is the external normal vector to the interfaces between phases i and j .In the following, the phase i will be denoted by i .We must keep in mind that the interface A ib between both regions is assumed to evolve versus time but sufficiently slowly so that a non-slip condition can be adopted.Note also that the viscoelastic behaviour of biomass that may cause deformation and slip flow at the biofilm surface ( Alpkvist and Klapper, 2007;Shaw et al., 1996 ) is neglected here.These equations are solved by using the IB-LB model presented in Benioug et al. (2015) .A direct-forcing Immersed Boundary (IB) model is coupled with the Lattice Boltzmann Method (LBM) to keep the fluid-solid interface sharp and maintain a direct control on numerical accuracy.More details about these numerical methods can be found in Benioug et al. (2015) and they will be only briefly reminded below.We use here a Lattice-Boltzmann model with a single relaxation time (SRT) and a BGK (Bhatnagar-Gross-Krook) collision operator.A three-dimensional lattice with 19 velocity vectors (D3Q19 model) is considered.The discretized Boltzmann equation is classically expressed as the collision propagation equation: where τ represents the relaxation time and f (eq )   k the equilibrium distribution function associated with the velocity vector e k .As detailed in Benioug et al. (2015) , the coupling with IBM is performed by adding a boundary force density into the momentum equation.Consequently, Eq. ( 2) is modified as follows: where ν is the kinematic viscosity and g represents the forcing term used to satisfy the no-slip condition on the immersed boundaries.The component g k of the forcing term along the vector e k is given by Kang and Hassan (2011) : with: and u the velocity in the bulk fluid and Ū the interpolated velocity at the boundary of the immersed cell.As illustrated in Fig. 1 , the linear interpolation of the component Ū j along the direction normal to the immersed boundary and located at a distance x from the interface, leads to: where φ( x , t ) describes the biofilm volume fraction for the grid cell x at time t .Note that LBM equations are solved in a transient form and steady state solution is assumed to have been achieved when the pressure convergence criterion is verified between two consecutive time steps.

Solute transport
The transport of chemical species both in the fluid and biofilm phases is governed by the following conservation equations: in the biofilm phase ( 11) where c f (resp.c b ) is the concentration of the solute within f (resp.b ), u is the velocity vector field as predicted from Eqs. ( 2) -( 4) and D f and D b are the molecular diffusion coefficients in f and b .K bf is the partition coefficient between the two phases f and b (its value is usually close to 1 for a fluid/biofilm system).R b is the reaction term which is described by Monod kinetics: where ρ bio , μ max , K are, respectively, the biomass concentration, the maximum reaction rate and the half-saturation constant.
Dirichlet and Neumann conditions are classically imposed on external boundaries of the domain.As previously for the flow calculation, non-boundary conforming method is required to simulate concentration fields at the pore-scale due to growth (or reduction) of b .We consider here a volume-of-fluid type model based on the assumption that the system behavior is close to the equilibrium conditions.This implies that the characteristic time for transport has the same order of magnitude in both phases and the characteristic time of the reaction kinetics is large ( Golfier et al., 2009 ).In other words, (i) the concentration gradients in the vicinity of the interface are sufficiently small and (ii) the concentrations do not vary significantly at the scale of the immersed cells.Therefore, a single concentration C both in f and b can be defined.This approximation leads to a one-equation transport model valid everywhere expressed as: where φ( x , t ) describes the volume fraction of b over .D φ , R φ and u φ represent respectively the molecular diffusion coefficient, the reaction term and the velocity vector within .They are expressed as a function of φ so that, for the two limit cases, i.e., φ = 0 (fluid cell) and φ = 1 (biofilm cell), Eqs. ( 10) and ( 11) are recovered.Thus, we have: (20)

Biomass growth and decay
The time evolution of the biomass concentration ρ bio is governed by the following mass balance equation: F λ is the stoichiometric coefficient of the biological reaction and k d is the extinction coefficient of the bacteria.This simplified law considers only the bacterial growth produced by the substrate consumption and extinction of bacteria.
The biomass concentration ρ bio calculated at each time step from Eq. ( 21) is then used to determine the value of the volume fraction φ by cell (the biomass density is supposed to be constant in the whole domain).Note that the net change in biomass concentration is calculated over a reference time step of growth, denoted t growth and corresponding to an increase in φ of 0.01 based on the maximum reaction rate μ max .As soon as the ρ bio value exceeds locally, within a biofilm cell, an arbitrary value denoted by ρ biomax (corresponding to a value of φ = 1), the biomass density is distributed spatially by using the cellular automata model (described in more details in the next section).

Biomass spreading
The spatial distribution of biofilm is modeled via a cellular automaton.The model is activated when the biomass concentration exceeds locally the imposed threshold value.As for any biofilm model, the cellular-automata rules assume that bacteria grow toward the areas where the substrate concentration is the most important.Basically, for each grid cell ( x, y ), as soon as φ( x, y ) > 1 (i.e., ρ bio ( x, y ) > ρ biomax ), the excess of biomass ρ bio = ρ bio ( x, y ) -ρ biomax is moved towards an adjacent cell.This algorithm is repeated for each cell in turn until all the excess biomass is redistributed.
During biomass spreading, depending on the nature of neighboring cells, three different cases can be identified: 1.If there is one fluid cell or more among the neighboring cells, ρ bio is placed in one of them, randomly chosen with equal probability, and this one becomes a partially-filled cell with a volume fraction φ = ρ bio ρ biomax .
2. If all the neighboring cells are biofilms but, at least one of them contains a biofilm volume fraction less than 1, ρ bio is randomly distributed in these cells with equal probability.If the volume fraction of the latter cells becomes higher than 1, the cellular automaton intervenes again to move the surplus of biomass towards adjacent cells.
3. If all the neighboring cells are biofilms and have a biofilm volume fraction equal or greater than 1, ρ bio is placed in one of them, randomly chosen with equal probability and the process of biomass spreading continues until a fluid or partially-filled cell is encountered.
Unlike classical models in which the fluid cell becomes a biofilm cell when it receives biomass, the main improvement of this algorithm is based on the management of the cells partially occupied by biofilm.This modification is, of course, motivated by the use of IB-LB methods for the flow calculation and the VOF approach employed for solute transport.However, we must note that this algorithm does not allow to specify the shape of the fluidbiofilm interface within the cell, necessary for using the Immersed Boundary Method in computing the velocity and pressure fields.To facilitate the resolution of the flow, we assumed in the following that these fictitious interfaces are parallel to the mesh facets and perpendicular to the direction of growth.

Biomass detachment
Only the detachment by erosion, i.e., under the effect of shear stresses acting on the superficial biofilm layers (i.e., cells with fluid/biofilm interface) is taken into account in the model.We assume that the partial or total detachment of the biomass present in a cell occurs when the shear stress exceeds locally a fixed threshold value, denoted τ max ( Choi and Morgenroth, 2003;Stoodley et al., 2001;Sudarsan et al., 2005 ).Although this approach is simplified with regard to the complexity of the processes involved (e.g., changes in mechanical properties of biofilm depending on the quantity of EPS produced) ( Rupp et al., 2005;Stoodley et al., 2002;Telgmann et al., 2004 ) it presents the advantage of being easy to calibrate in comparison with laboratory data sets e.g., Paul et al. (2012) and Walter et al. (2013) .It also allows us to consider both surface erosion (for low values of φ det ) and volume erosion (the maximum size of particles that can be detached is limited by the cell size) which seems to play a major role in the evolution of the biofilm thickness ( Derlon et al., 2013 ).
To determine the quantity of eroded biomass, we calculate in a first step the shear stress exerted for each biofilm cell on the facets in contact with the fluid.The maximum stress exerted on each cell is taken as reference value for the calculation of the detachment.This value is obviously updated at each time step according to the changes of the medium geometry and of the velocity field.For the biofilm cell ( i, j ) represented in Fig. 2 , we have in 2 D for instance: where the cells east, west, south and north represent the neighboring fluid cells.The parietal velocity gradients are discretized by finite differences and calculated between the fluid/biofilm interface -where the velocity is zero -and the centre of the adjacent fluid cell distant of length l stress .Thus, for the shear stress on the upper face, i.e., the north face, we obtain: In the particular instance illustrated in Fig. 2 , l stress is given by: Note that if the cell ( i, j ) would be fully occupied by biofilm, i.e., with the fluid/biofilm interface located at ( i , j + 1 2 ) we would recover l stress = y 2 .To model the process of detachment, we are based on the experimental results of Paul et al. (2012) which indicate that the thickness of the biofilm exponentially decreases with increasing shear stress.From these observations, we have established an empirical formulation relating the detached biofilm volume fraction φ det to the shear stress τ : with C det an empirical constant.The remaining volume fraction of biomass after detachment can be derived easily as follows (we assume that the biomass keeps a constant density): such as φ t is the biofilm volume fraction at the previous time step, i.e., before detachment.The last step consists in describing the fate of the detached biofilm particles, i.e., the process of biomass attachment.

Biomass attachment
As the transport of planktonic bacteria is neglected in this model, we suppose that the detached biofilm particles are only entrained by the flow.Thus, if the shear stresses are favourable, i.e., τ ≤ τ attach the biomass can potentially reattach to the grain surface or to the surface of another biofilm encountered along its pathway.It should be noted that despite its impact on the bio-plugging of groundwater media ( Bunt et al., 1993;Cowan et al., 1991;Scheuerman et al., 1998 ), this mechanism is rarely accounted for in biofilm growth models ( Kapellos et al., 2007 ).
The detached biofilm particles are assimilated to bacterial flocs advected along streamlines.They are assumed to be sufficiently small not to be submitted to inertial effects (Stokes number less than 1) and have the same density as the fluid (a biofilm consists of more than 90 % water) so that the sedimentation effects are negligible.Moreover, we take into account the assumption that attachment conditions depend only on environmental hydrodynamic stresses.Consequently, the algorithm is designed as follows: for each fluid cell containing a non-zero biomass concentration (i.e., biomass previously detached), bacteria concentration is transported along the streamline (from cell centre to cell centre) during the calculation time step t growth .For each cell crossed along the path way, the model checks the probability of biomass attachment ( Fig. 3 a).If one or more of the neighboring cells are biofilm or solid cells and if the value of the local shear rate is low enough ( τ ≤ τ attach ), then the virtual biofilm particle is attached, i.e., the fluid cell is converted into partially-filled cell ( Fig. 3 b).Otherwise, the biofilm particle follows its path way until it leaves the domain.
Finally, the general algorithm of our model including fluid flow, chemical species transport, biofilm growth, biomass attachment and detachment is illustrated in Fig. 4 .

Results and discussion
To illustrate and evaluate the predictive capability of our numerical model, two numerical benchmarks taken from the literature ( Eberl et al., 2001;Picioreanu et al., 1998a;20 0 0;Wanner et al., 2006 ) have been initially performed.The two tests were carried out in a 2D vertical configuration but for hydrostatic and hydrodynamic conditions.For the first one, the biofilm grows on a solid flat surface and the nutrient is conveyed downward only by diffusion from the upper boundary.From this benchmark, we investigate the ability of the model to properly reproduce the biofilm growth limiting regimes, i.e., under mass transfer-limited and reaction-rate-limited conditions.In the second test case, we impose a lateral flow to examine the impact of hydrodynamic con- ditions and verify that the model captures the dynamics of the developing biofilm pattern.
Finally, we apply our model to a 2 D porous medium.The aim of this study is to assess the influence of biofilm growth on the changes of the geometrical (bioclogging) and physical properties (porosity, permeability) of the porous medium.A peculiar attention is paid to the porosity-permeability relationship and its dependency to biofilm growth regimes.
To facilitate the analysis and the interpretation of results, numerical data will be reported in a dimensionless form as described below: where τ ∞ represents the maximum shear stress value encountered in the whole domain.

Case 1. Biofilm growth under hydrostatic conditions
The computational domain with a uniform mesh size (Nx − 1) × (Nz − 1) and a resolution where Nx and Nz are the number of nodes along the longitudinal ( x ) and vertical ( z ) directions.We consider n 0 biofilm colonies randomly distributed on the solid flat surface, located at z = 0 ( Fig. 5 a).Each colony occupies initially one cell in the x -direction and two cells in the z -direction.
The initial biomass concentration is the same for all the colonies of biofilm and is ρ bio (t = 0) = 1 .A constant concentration C 0 in nutrient is imposed on the upper limit of the domain, i.e., at z = 1 .Periodic boundary conditions are imposed in the xdirection in order to avoid side effects.Finally, a wall condition is imposed on the solid border at z = 0 .Dimensionless values of simulation data are given below ( Table 1 ).Note that μ max and K values depend on the fixed value of Da .

Time evolution of the biofilm structure
In this simulation, an intermediate value of the Damköhler number corresponding to Da = 10 is considered.During the first time steps of simulation, the stratification of nutrient concentration field is observed under the effect of the diffusion.The isoconcentrations are practically horizontal along the x -direction.We observe a homogeneous development of biofilm which occupies practically the entire solid surface ( Fig. 6 a).The biofilm growth is progressively causing a change in the availability of nutrient ( Fig. 6 b).This disturbance of the concentration field intensifies progressively and destabilizes the biomass development ( Fig. 6 b and c).Indeed, the substrate loading rate being the limiting factor, the bacterial colonies near the source continue to grow producing biomass at the detriment of those which are most remote.This heterogeneous distribution, that favors the large bacterial colonies ( Fig. 6 d), leads to the formation of finger-like structures ( Fig. 6 c and d) which have also been reported in the literature ( Costerton et al., 1994;Eberl et al., 2001;Picioreanu et al., 1998c ).

Influence of Da number on biofilm pattern
In this section, we assess the Damköhler number influence on the biofilm growth.The numerical simulations are performed for a large range of Da numbers varying from 5 • 10 −2 to 10 3 to investigate biofilm growth regimes, both under mass transfer-limited and/or reaction rate-limited conditions.For low Da , the influence of mass transfer rate within the biofilm is predominant.Indeed, the low nutrient consumption by bacteria results in a large penetration of substrate into the biofilm ( Fig. 7 a -right).This observation is corroborated by Fig. 7 a-left, where the average concentration C along the x -direction remains above zero over the entire thickness of the biofilm.Such hydrostatic conditions lead to a compact and smooth biofilm ( Fig. 7 a), corresponding to a growth limited by the kinetics.In contrast, when Da increases, the reaction kinetics becomes more important relative to the diffusion rate.The biofilm growth is more and more limited by the solute transport.Consequently, the outside surface of the biofilm becomes more irregular and the biofilm evolves towards a more porous structure with many channels and voids between the colonies that grow in the form of fingers ( Fig. 7 b and d).The concentration no longer penetrates inside the biofilm and the unstable growth process is amplified (the fingers become thinner and narrower).The average bacterial concentration ρ is also more spread out over the zdirection and decreases when approaching the solid surface where the most mature bacteria layers are affected by the cellular extinction due to nutrient depletion.

Impact of Péclet number
This second test case is identical with the one above except for the boundary conditions and the values of the dimensionless numbers Pe and Da .All the other numerical data given in Table 1 remain unchanged.The condition imposed on the upper limit, i.e., at z = 1 , is now a symmetry condition and the substrate concentration C = 1 is injected at x = 0 with a parabolic velocity profile ( Fig. 5 b).These conditions are expressed as follows: A random distribution of n 0 = 25 colonies was applied at the initial time.To limit the influence of the boundary conditions, the colonies were placed between x = 0 .26 and x = 0 .74 .At t = 0 , each colony are occupying one cell of width and two cells of height.
The influence of convection on substrate availability and on the changes in bacterial population is investigated by modeling the biofilm growth under different hydrodynamic regimes from the variation of the Péclet number.Consequently, the Damköhler number is kept constant for all the simulations, i.e., Da = 1 and the values of the Péclet number are Pe = 0 . 1 , 1 and 100.The results are illustrated in Fig. 8 .It is important to note that in these simulations the mechanisms of attachment/detachment of the bacteria have been desactivated.
• For the largest Pe number ( Fig. 8 c), the biofilm growth seems not to be influenced by the changes in the velocity field.The biofilm colonies show relatively uniform growth.The high flow velocities favor a constant supply of nutrient everywhere, the isoconcentrations are thus globally uniform around the biofilm patch.• The decrease in the Péclet number ( Fig. 8 a and b), causes a differentiated biofilm growth by favoring the colonies near the nutrient source.Indeed, when Péclet number decreases and Da remains constant, the convection characteristic time becomes very large relative to the reaction characteristic time.The substrate availability then decreases progressively away from the source (here, x = 0 ).This phenomenon is amplified by the

Impact of bacterial detachment
In a second step, our attention is focused on the mechanisms of bacteria attachment/detachment and their impact on biofilm growth.However, as this mechanism implicitly depends on the velocity field, we first studied the mesh-dependency of our numerical results.Two different meshes are considered: Nx × Nz = 100 × 33 and 80 × 26.The Péclet number is Pe = 50 with a critical shear stress τ max = 10 3 and the Damköhler number is Da = 1 .
The medium is initially inseminated by a volume of bacteria corresponding to 3.8% of total volume.The bacteria were randomly distributed on the solid surface but between the same limits than above, i.e., between x = 0 .26 and x = 0 .74 .The results are illustrated in Fig. 9 .First, we observe that the biofilm patterns are very similar, independently of the mesh size ( Figs. 9 a and b).The analysis of the increase of biomass concentration versus time ( Fig. 9 d) supports this finding and confirms the capability of this IB-LB based model to predict biomass detachment.Second, comparison with numerical simulations in the absence of bacterial detachment ( Figs. 9 c and d) highlights the role of this regulatory mechanism on biomass development.Impact of detachment becomes clearly visible after t = 22 ( Fig. 9 d).After t = 32 , the biofilm height reaches a plateau and stabilizes around 0.4 ( Figs. 9 a and b) while the biomass continues to grow infinitely in the absence of bacterial detachment, at least if the substrate availability is sufficient ( Fig. 9 c).The progressive biofilm growth causes a reduction of the flow cross-section and hence, an increase of the local shear stresses.When the shear stress reaches the critical value τ max , erosion of superficial biofilm layers occurs.The onset of this steady-state plateau corresponds to a balance between growth and detachment.
As a conclusion, the numerical simulations carried out in this section indicate that the model captures properly the dynamics of biofilm growth and results are in agreement with the literature.In the next section, we apply our model to porous media colonised by biofilms.

Application to a 2D porous medium
We consider now a representative volume of a complex porous medium ( Fig. 10 ).Such a geometry gives us an overview of how the pore-scale architecture (solid distribution and pore connectivity) may impact biofilm development.The domain has a size of lx × lz = 12 × 12 mm 2 with an initial porosity 0 of 0.6 and the pore throat characteristic length, l f , is estimated to be 0.3 mm ( Wood et al., 2007 ).The medium is initially inseminated by a volume of bacteria corresponding to 5.5% of total volume.The bacteria are randomly distributed on the pore walls.A velocity field is imposed at the inlet of the domain (i.e., at x = 0 ) while the pressure is kept constant at the outlet (i.e., at x = 1 ).Periodicity conditions are introduced at the boundaries z = 0 and z = 1 .The influence of fluid flow, substrate consumption and the impact of detachment and attachment mechanisms on the structural architecture of the biofilm are then investigated through different values attributed to the dimensionless Péclet and Damköhler numbers.
Table 2 contains the dimensionless parameters used for the numerical simulations and values are in the range of those found in the literature ( Bottero et al., 2013;Knutson et al., 2005;Picioreanu et al., 1998a ).Note that the value of and F λ have been chosen to achieve a full uncoupling between the biofilm growth and nutrient transport and consumption.The extinction coefficient k d was fixed at k d = 3 .4 × 10 −8 s and the half-saturation constant K varies between 0.0032 mg/l and 0.32 mg/l.The maximum growth rate μ max was imposed via the Damköhler number.
Note that the dimensionless shear stress value, τ , implicitly related to the Péclet number (through the pore velocity magnitude) vary also with the numerical computations performed ( τ max is kept constant).From all the simulations carried out, three pairs of

Low Pe, low Da, low τ
For low values of Pe ( Pe = 1 ), there is no influence of shear stress on the biofilm growth.Moreover, as the Damköhler number ( Da = 0 .01 ) is small compared to the Péclet number, transport mechanisms (diffusion in the biofilm phase and convection and diffusion in the fluid phase) prevail relative to the reaction kinetics.Consequently, the substrate is uniformly available throughout the porous medium ( Fig. 11 a) and the biomass grows in a relative uniform manner around the solid grains ( Fig. 11 c).
However, even for these low values of Pe − Da, biomass accumulation is slightly more pronounced near the nutrient source, i.e., at x = 0 ( Fig. 11 b) due to a low but non-zero concentration gradient along the flow direction.This small discrepancy in growth is amplified with time and leads to a local decrease of the medium porosity, and eventually a clogging of the pores in the vicinity of the inlet ( Figs. 11 c and d).Ultimately, the bacterial colonies far from the entrance will receive less substrate with consequences on their development and finally, causing their extinction.We should recover at large times the observed behavior for low Pe , high Da (see Section 4.3.2 ).

Low Pe, high Da, low τ
We now consider the case of a growing biofilm under conditions of large Damköhler number, i.e., Da = 100 but still at low Péclet number ( Pe = 0 . 1 ).The characteristic reaction time is now significantly lower than the characteristic time associated with nutrient transport.In other words, the electron donor species is consumed before completely passing through the medium.As we can see in Fig. 12 a, we have a very large concentration gradient on the first cells of domain and a concentration that drops down close to zero just after the entrance.This significant reduction of the concentration leads at the inlet, where the substrate concentration is maximum, to a relatively large increase in bacterial concentration ( Fig. 12 b and c).On the contrary, near the outlet where the amount of available substrate is reduced or insufficient, biomass growth is very limited or even non-existent ( Fig. 12 a).Bacteria extinction begins to appear far from the injection zone ( Fig. 12 b).
It is for such conditions that the plugging of the medium may be fast.In our case, a full bio-clogging near the entrance of the porous media occurs after t = 2 ( Fig. 12 c and d).The hydrodynamic properties of the medium are strongly impacted.The porosity of the medium decreases slightly over the entire domain (although it decreases strongly locally) while the permeability falls of several orders of magnitude.

High Pe, low Da, high τ
We end this analysis by the case at high Péclet number, Pe = 500 and low Damköhler, Da = 0 . 1 .As in the case considered in section 4.3.1 ( Pe = 1 -Da = 0 .01 ), the Péclet number is much larger than the Damköhler number, which implies that the transport mechanisms are predominant relative to the reaction rate.At short times, the behaviour is identical with that described in the above-mentioned paragraph characterized by a rather uniform biofilm growth around the grains and uniform distribution of nutrient concentration within the medium ( Fig. 13 a).
However, with the development of the biofilm, the pores gradually shrink.We observe firstly a clogging of the finer pores, whereas only a few channels remain open.These channels rapidly create preferential flow paths which lead to an increase of local flow velocities and hence, of shear stresses exerted on the biofilm surface ( Fig. 13 e).As soon as the shear stresses exceed the imposed critical value τ max (i.e., τ > 1), the biomass is subjected to erosion ( Fig. 13 c).In stagnant flow regions, located between several solid grains, the clogging becomes almost complete ( Fig. 13 d).After some time, t = 140 , the system seems to reach an equilibrium between bacterial growth and detachment processes.The substrate is only available within the channels where the growth is no longer possible due to the intensity of the shear stress ( Fig. 13 b).The preferential flow mechanism encountered in this simulation is very similar to the one observed during dissolution in porous media with the development of wormholes ( Golfier et al., 2002 ).

Growth regime diagram Pe vs Da vs τ
Based on the simulations carried out, three characteristics regimes of biofilm growth in porous media have been observed.These regimes can be related to the dimensionless numbers values of Pe, Da and τ : Reaction rate-limited growth: biofilm develops mainly close to the injection zone of the nutrient, leading rapidly to a total bio-plugging of the medium ( Fig. 12 ).This regime occurs when Pe Da and τ < 1.
Mass transfer-limited growth: biofilm grows uniformly in the medium resulting in a progressive but uniform decrease in porosity ( Fig. 11 ).This regime occurs when Pe Da and τ < 1. Shear rate-limited growth: preferential flow paths develop within the colonised medium, under the effect of shear forces.The obstruction of the medium is never complete ( Fig. 13 ).This regime occurs when τ > 1.

Calculation of macroscopic properties
At this point, we can take advantage of these pore-scale simulations to investigate the impact of physical changes in the microstructure on the effective properties of the porous medium.
To illustrate these biofilm-induced changes on the macroscopic  behaviour of the medium, we will focus on the relation existing between permeability and porosity.Indeed, as a consequence of biofilm growth, the relationship between the macroscopic quantities may be historical.It must be emphasized that this question is of a general interest, and similar problems were found in geochemistry, or in dissolution in porous media ( Golfier et al., 2002 ).An approach, classically made, consists in assuming a direct relationship between the permeability and the porosity, i.e.K ( ) instead of K ( t ) and ( t ).Here, we will adopt this point of view.
In order to test this assumption, the biofilm growth patterns obtained at the pore-scale in the previous Section 4.3 were used to determine the macroscopic parameters.From the series of simulations performed for different values of Pe and Da on the representative elementary volume (REV) of porous medium, the biomass volume fraction is calculated at the different time steps by spatial integration.In a similar way, the permeability is calculated over the REV by Darcy's law: where K is the permeability, L the length of the domain (here equal to l x ) and dP the pressure difference.Fig. 14 shows the temporal evolution of both the biofilm volume fraction b and the dimensionless permeability K = K K 0 , where K 0 is the initial permeability.For different values of Da and Pe , a quick decay of the permeability at the beginning of the simulation is first observed.This sharp reduction of K is due to the choice of initial conditions.We start the simulations with C = 1 everywhere and a high value of the initial bacterial concentration ( ρ bio = 1 ) that enhance a fast spreading of biomass in the pores before reaching a stabilisation of the concentration profile.Then, permeability gradually decreases more or less drastically.For reaction rate-limited regime, permeability falls down to zero due to the final bio-clogging of the domain ( Fig. 13 b and d) and drops of sev- eral orders of magnitude for the mass transfer limited regime.Such a significant permeability decrease is in agreement with the experimental observations of Cunningham et al. (1991) and Kim et al. (2006) who have found a reduction in permeability up to three orders of magnitude due to biofilm formation.Conversely, for the shear rate-limited growth regime, the permeability converges to a plateau ( Fig. 13 f) corresponding to a decrease of about one order of magnitude.This phenomenon is due to the formation of preferential flow channels that preserve partially the permeability of the porous medium.In the same time, the biofilm volume fraction increases quasi linearly for all the simulations after an abrupt increase at short times (related to the sharp decrease of permeability discussed above).For both reaction rate-limited and mass transfer-limited growth regimes, the increase of b is limited by the bioclogging ( Fig. 13 a and b) while for the shear rate-limited growth regime, it stabilizes around a constant value (about 0.4 for Da = 0 . 1 and Pe = 500 ) due to the quasi-equilibrium state between biomass growth and bacterial detachment ( Fig. 13 e).
Finally, we have represented in Fig. 15 the correlations relating the change in permeability with the increase in biofilm volume fraction for all the growth regimes.We observe significant discrepancies between these different relationships due to the differences in biomass distribution for the biofilm growth regimes.For comparison purpose, we have also represented in the same figure the permeability changes predicted by the semi-empirical Kozeny-Carman law classically used in macroscopic models of biofilm growth, (44) where C S represents a geometrical shape factor, 0 the initial porosity of the domain and b the biomass volume fraction.
If this formulation predicts correctly the permeability changes for the mass transfer limited regime since the biofilm is coating gradually the solid grains, it fails to reproduce the behavior observed when the shear stress or the reaction rate is the limiting step for biofilm growth.This result emphasizes the need of further research for predicting bio-clogging and shows that an universal relationship between hydraulic conductivity and effective porosity cannot be formulated.This important point should be kept in mind for large-scale simulations in biofilm-related applications.

Conclusions
In this work, a two-dimensional pore-scale numerical model of biofilm growth was developed.We have adopted a cellular automata model for describing the spatial and temporal evolution of the biomass under the influence of bacterial growth and the mechanisms induced by fluid flow, such as detachment or attachment of bacterial floc.The use of immersed boundary methods allowed us to overcome the main drawback encountered in this type of approach because of the size of the elementary cell that constrains the volume of biomass produced at each time step.Fluid flow is simulated with an immersed boundary-lattice Boltzmann model while solute transport is described with a volume-of-fluid-type approach.First, this model was validated from qualitative comparison with benchmark cases issued from the literature.Under hydrostatic conditions, the numerical solutions are strongly dependent of the Damköhler number.At low Da , the influence of the mass transfer rate within the biofilm is predominant.The large penetration of nutrient concentration inside the biofilm leads to an enhanced microbial growth and a bacterial concentration relatively uniform across the biofilm.At high Da , the biofilm growth is limited by substrate diffusion, generating an unstable development in the form of fingers for bacterial colonies.In hydrodynamic conditions, our interest is focused on the influence, through the Péclet number, of convection on substrate availability and the changes in bacterial population.When the detachment mechanisms are not considered, the bacterial growth is heterogeneous at low Péclet number and tends towards a homogeneous distribution at high Péclet number.In a second step, the model was applied to a 2D complex porous medium, in order to analyze the influence of environmental conditions (reaction kinetics, flow rate) on biofilm patterns.Based on the numerical results obtained, a regime diagram representing the different modes of biofilm growth was proposed.
It depends on the Damköhler Da and Péclet Pe numbers and also on the dimensionless shear stress τ .Thus, for Pe Da and τ < 1, the reaction rate-limited growth regime takes place.The mass transfer-limited growth regime occurs when Pe Da and τ < 1.Finally, if τ > 1, the shear rate-limited growth regime is observed.Moreover, the influence of biofilm growth on the macroscopic properties of the medium was discussed by analyzing the evolution of correlations between permeability and porosity (or biomass volume fraction) as a function of biofilm pattern.Results suggest the non-unicity of these permeability -porosity relationships in the presence of bacterial development characterised by a strong path-dependency.
If the numerical models proposed in this work helped to shed additional light on the mechanisms governing the development of bacteria in porous media, unfortunately we could not validate our numerical model against experimental data.However, recent advances in non-invasive imaging methods for characterizing the distribution of biomass within a porous medium ( Davit et al., 2011 ) should open up new possibilities to better understand the complexity of mechanisms involved and to assess the predictive capabilities of biofilm growth models.
Several ways of improvement for the model have been identified and are under investigation.First, the model is limited to mono-species biofilm and considers only a single nutrient source.Extending the model to several bacterial species using different electron acceptors (aerobic /anaerobic bacteria for example) would be interesting to investigate the effects of competition and stratification of different bacterial colonies within the biofilm.Moreover, the production of EPS and bacterial cells is not differentiated.It would be wise to make evolve this EPS production according to environmental stress conditions and to assess its impact on the hydrodynamic properties of the medium (e.g., resistance to detachment depending on the quantity of EPS produced).

Fig. 1 .
Fig. 1.Illustration of the interpolation procedure in the x direction with the immersed boundary method (IBM).

Fig. 2 .
Fig. 2. Cell ( i, j ) with its four neighbors and their velocity vectors for calculating the shear stress in 2 D .

Fig. 5 .
Fig. 5. Schematic representation of the computational domain.Bacterial colonies are green, the solid substrate is shown in red.(a) biofilm growth in hydrostatic conditions.(b) biofilm growth in hydrodynamic conditions.The thick arrow indicates the direction of nutrient transport.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6 .
Fig. 6.Temporal evolution of the biofilm geometry for Da = 10 .Biofilm is shown in green, fluid in blue and solid substratum in red.The black lines indicate the lines of equal concentration.These lines show the decrease of substrate concentration from the upper fluid surface (blue region) to the biofilm (green), with a variation of 10% between the lines.a) t = 10 , b) t = 50 , c) t = 100 , d) t = 160 .(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 7 .
Fig. 7. Right: Biofilm patterns observed for different Da .The black lines represent the isoconcentrations with a variation of 10% between each one.Left: profiles of the average substrate concentration C in red, and the average biomass concentration ρ in green.(a) Da = 0 .05 , t = 600 , (b) Da = 10 , t = 160 , (c) Da = 10 0 0 , t = 46 .(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 10 .
Fig. 10.Initial geometry (representative region of a porous medium extracted from Wood et al., 2007 ), the biofilm phase is illustrated in green, solid in red and fluid in blue.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Pe − Da values have been chosen as particularly representative of the different growth regimes obtained: Pe = 1 -Da = 10 −2 (low Pe , low Da , low τ ), Pe = 10 −1 -Da = 10 2 (low Pe , high Da , low τ ) and Pe = 500 -Da = 10 −1 (high Pe , low Da , high τ ).The results are discussed below.

Table 1
Simulation data values.

Table 2
Physical and biological parameter values used in the simulations.