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

M. Benioug, F. Golfier, C. Oltéan, M.A. Buès, T. Bahar, J. Cuny

► To cite this version:

M. Benioug, F. Golfier, C. Oltéan, M.A. Buès, T. Bahar, et al. An immersed boundary-lattice Boltzmann model for biofilm growth in porous media. Advances in Water Resources, 2017, 107, pp.65-82. 10.1016/j.advwatres.2017.06.009. hal-02968533

HAL Id: hal-02968533

https://hal.univ-lorraine.fr/hal-02968533

Submitted on 16 Oct 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An immersed boundary-lattice Boltzmann model for biofilm growth in porous media

M. Benioug*, F. Golfier, C. Oltéan, M.A. Buès, T. Bahar, J. Cuny

University of Lorraine / CNRS / CREGU, GeoRessources Laboratory, BP 40, F54501 Vandoeuvre-les-Nancy, France

## ARTICLE INFO

**Article history:**

Received 30 July 2016

Revised 18 May 2017

Accepted 8 June 2017

Available online 15 June 2017

**Keywords:**

Lattice Botzmann method

Immersed boundary method

Biofilm growth

Porous media

Cellular automata

## ABSTRACT

In this paper, we present a two-dimensional pore-scale numerical model to investigate the main mechanisms governing biofilm growth in porous media. The fluid flow and solute transport equations are coupled with a biofilm evolution model. Fluid flow is simulated with an immersed boundary-lattice Boltzmann model while solute transport is described with a volume-of-fluid-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 fluid circulation, nutrient transport and bacterial growth is tested under different hydrostatic and hydrodynamic conditions (i) on a flat medium and (ii) for a complex porous medium. For the second case, different regimes of biofilm growth are identified 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 biofilm growth on the macroscopic properties of the porous medium is investigated and we discuss the unicity of the relationships between hydraulic conductivity and biofilm volume fraction.

## 1. 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., 2001; Prakash et al., 2003; Smith, 2005). 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, 2000; 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, 2000), 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 quantum sensing (Frederick et al., 2011), physiological state of the cells

<footnote>* Corresponding author.
E-mail address: benioug@gmail.com (M. Benioug).</footnote>
(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 moti-
vated 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 descrip-
tions have been achieved (Ebigbo et al., 2013; 2010; Kapellos et al.,
2007; Knutson et al., 2005; Peszynska et al., 2016; von der Schu-
lenburg et al., 2009) and have given additional insight on the
bio-chemical and physical mechanisms that drive the biofilm dy-
namics, important questions on how the morphology and distribu-
tion 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 cap-
pable 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 kinet-
ics 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 hydrody-
namic properties of the medium. In particular, the relationship be-
tween porosity and permeability has been discussed as a function
of the biofilm growth regime.

## 2. 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 dis-
crete (or individual-based) approach (including multi-agent and
cellular automata models) and continuous approach. In the multi-
agent 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 re-
mains 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., 2006; 2006; 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 multi-
species simulation by Kreft et al. (2001). Multi-agent models have
been successfully applied to describe and optimize various ap-
plications 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 rel-
evant if the interest is focused on interactions between individ-
uals (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 com-
putation time particularly in a 3D configuration (von der Schu-
lenburg et al., 2009). Conversely, cellular automata model divides
the space into grid cells, each cell being considered as individ-
ual entities that can take several possible states (e.g., solid, fluid
or biofilm). The state of a grid cell is determined by the previ-
ous 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 aggre-
gation (DLA) model used by Schindler and Rataj (1992); Schindler
and Rovensky (1994) and Fujikawa (1994). Ever since, several mod-
els have been developed using this approach (Hermanowicz, 1998;
2001; Khassehkhan et al., 2009; Laspidou et al., 2012; Laspidou
and Rittmann, 2004; Laspidou et al., 2005; Liao et al., 2012; Pi-
cioreanu 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 2D configuration) until reach-
ing 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 ap-
proach 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 tak-
ing 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 bacte-
rial growth causes the development of a pressure field inside the
biofilm (Eberl et al., 2001). This pressure field constitutes the driv-
ing 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 ve-
locity of the fluid/biofilm interface due to biofilm growth in a one-
dimensional model. The average velocity of biomass is thus deter-
mined by mechanical laws, and more precisely by the Euler equa-
tions 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 meth-
ods 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 automa-
ton. This cell size constrains the minimal change in biofilm thick-
ness 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 ac-
curacy 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 cel-
lular automaton model will take advantage of these numerical ef-
forts to work on entities smaller than the mesh size.
### 3. 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.

#### 3.1. 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:

$$ \nabla \cdot \mathbf{u} = 0, \qquad (1) $$

$$ \rho_f \left[ \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right] = -\nabla p + \mu_f \nabla^2 \mathbf{u}, \qquad (2) $$

with non-slip boundary conditions for the fluid-solid and fluid-biofilm interfaces, respectively $A_{fs}$ and $A_{fb}$,

$$ C.L.1 \qquad \mathbf{u} = 0, \qquad \text{at } A_{fb} \qquad (3) $$

$$ C.L.2 \qquad \mathbf{u} = 0, \qquad \text{at } A_{fs} \qquad (4) $$

where $\mathbf{u}$, $p$, $\mu_f$, $\rho_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 $\mathbf{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 $\Omega_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:

$$ f_k(\mathbf{x} + \mathbf{e}_k \Delta t, t + \Delta t) = f_k(\mathbf{x}, t) - \frac{\Delta t}{\tau} \left[ f_k(\mathbf{x}, t) - f_k^{(\text{eq})}(\mathbf{x}, t) \right] \qquad (5) $$

where $\tau$ represents the relaxation time and $f_k^{(\text{eq})}$ the equilibrium distribution function associated with the velocity vector $\mathbf{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:

$$ \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} = - \frac{\nabla p}{\rho_f} + \nu \nabla^2 \mathbf{u} + \frac{\mathbf{g}}{\rho_f}, \qquad (6) $$

where $\nu$ is the kinematic viscosity and $\mathbf{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 $\mathbf{e}_k$ is given by Kang and Hassan (2011):

$$ g_k(\mathbf{x}, t) = \left(1 - \frac{1}{2\tau}\right) \omega_k \left[ 3 \frac{\mathbf{e}_k - \mathbf{u}(\mathbf{x}, t)}{c^2} + 9 \frac{\mathbf{e}_k \cdot \mathbf{u}(\mathbf{x}, t)}{c^4} \mathbf{e}_k \right] \cdot g_2(\mathbf{x}, t) \qquad (7) $$

with:

$$ g_2(\mathbf{x}, t) = 2\rho(\mathbf{x}, t) \frac{\bar{\mathbf{U}}(\mathbf{x}, t) - \mathbf{u}^*(\mathbf{x}, t)}{\Delta t} \qquad (8) $$

and $\mathbf{u}^*$ the velocity in the bulk fluid and $\bar{\mathbf{U}}$ the interpolated velocity at the boundary of the immersed cell. As illustrated in Fig. 1, the linear interpolation of the component $\bar{\mathbf{U}}_j$ along the direction normal to the immersed boundary and located at a distance $(1 - \phi(\mathbf{x}_f, t)) \Delta x$ from the interface, leads to:

$$ \bar{U}_j = \frac{(1 - \phi(\mathbf{x}_f, t)) \Delta x u_j^*}{\Delta x + (1 - \phi(\mathbf{x}_f, t)) \Delta x} \qquad (9) $$
<figure><img src="image_2.png" /><figcaption>Fig. 1. Illustration of the interpolation procedure in the x direction with the immersed boundary method (IBM).</figcaption></figure>

where $\phi(\mathbf{x}, t)$ describes the biofilm volume fraction for the grid cell $\mathbf{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.

## 3.2. Solute transport

The transport of chemical species both in the fluid and biofilm phases is governed by the following conservation equations:

$$ \nabla \cdot (\mathbf{u}\mathbf{c}_f) = \nabla \cdot (D_f \nabla c_f), \quad \text{in the fluid phase} \qquad (10) $$

$$ 0 = \nabla \cdot (D_b \nabla c_b) - R_b, \quad \text{in the biofilm phase} \qquad (11) $$

$$ D_f \nabla c_f \cdot \mathbf{n}_{fb} = D_b \nabla c_b \cdot \mathbf{n}_{fb}, \quad \text{at } A_{fb} \qquad (12) $$

$$ c_f = K_{bf} c_b, \quad \text{at } A_{fb} \qquad (13) $$

$$ 0 = D_f \nabla c_f \cdot \mathbf{n}_{fs}, \quad \text{at } A_{fs} \qquad (14) $$

$$ 0 = D_b \nabla c_b \cdot \mathbf{n}_{bs}, \quad \text{at } A_{bs} \qquad (15) $$

where $c_f$ (resp. $c_b$) is the concentration of the solute within $\Omega_f$ (resp. $\Omega_b$), $\mathbf{u}$ is the velocity vector field as predicted from Eqs. (2)-(4) and $D_f$ and $D_b$ are the molecular diffusion coefficients in $\Omega_f$ and $\Omega_b$. $K_{bf}$ is the partition coefficient between the two phases $\Omega_f$ and $\Omega_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:

$$ R_b = \rho_{bio} \mu_{max} \frac{c_b}{c_b + K} \qquad (16) $$

where $\rho_{bio}$, $\mu_{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 $\Omega_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 $\Omega_f$ and $\Omega_b$ can be defined. This

approximation leads to a one-equation transport model valid everywhere expressed as:

$$ \nabla \cdot (\mathbf{u}_{\phi}C) = \nabla \cdot (D_{\phi} \nabla C) - R_{\phi}, \quad \text{in } \Omega = \Omega_f \cup \Omega_b \qquad (17) $$

where $\phi(\mathbf{x}, t)$ describes the volume fraction of $\Omega_b$ over $\Omega$. $D_{\phi}$, $R_{\phi}$ and $\mathbf{u}_{\phi}$ represent respectively the molecular diffusion coefficient, the reaction term and the velocity vector within $\Omega$. They are expressed as a function of $\phi$ so that, for the two limit cases, i.e., $\phi = 0$ (fluid cell) and $\phi = 1$ (biofilm cell), Eqs. (10) and (11) are recovered. Thus, we have:

$$ D_{\phi}(\mathbf{x}, t) = \phi(\mathbf{x}, t)D_r + (1 - \phi(\mathbf{x}, t))D_f \qquad (18) $$

$$ R_{\phi}(\mathbf{x}, t) = \phi(\mathbf{x}, t) R_b(\mathbf{x}, t) \qquad (19) $$

$$ \mathbf{u}_{\phi}(\mathbf{x}, t) = (1 - \phi(\mathbf{x}, t)) \mathbf{u}(\mathbf{x}, t) \qquad (20) $$

## 3.3. Biomass growth and decay

The time evolution of the biomass concentration $\rho_{bio}$ is governed by the following mass balance equation:

$$ \frac{\partial \rho_{bio}}{\partial t} = F_{\lambda} \rho_{bio} \mu_{max} \frac{c_b}{c_b + K} - k_d \rho_{bio} \qquad (21) $$

$F_{\lambda}$ 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 $\rho_{bio}$ calculated at each time step from Eq. (21) is then used to determine the value of the volume fraction $\phi$ 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 $\phi$ of 0.01 based on the maximum reaction rate $\mu_{max}$. As soon as the $\rho_{bio}$ value exceeds locally, within a biofilm cell, an arbitrary value denoted by $\rho_{biomax}$ (corresponding to a value of $\phi = 1$), the biomass density is distributed spatially by using the cellular automata model (described in more details in the next section).

## 3.4. 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 $\phi(x, y) > 1$ (i.e., $\rho_{bio}(x, y) > \rho_{biomax}$), the excess of biomass $\Delta\rho_{bio} = \rho_{bio}(x, y) - \rho_{biomax}$ is moved towards an adjacent cell. This algorithm is
repeated for each cell in turn until all the excess biomass is redis- tributed.

During biomass spreading, depending on the nature of neigh- boring cells, three different cases can be identified:

1. If there is one fluid cell or more among the neighboring cells, Δρ<sub>bio</sub> is placed in one of them, randomly chosen with equal probability, and this one becomes a partially-filled cell with a volume fraction φ = Δρ<sub>bio</sub>/ρ<sub>biomax</sub>.

2. If all the neighboring cells are biofilms but, at least one of them contains a biofilm volume fraction less than 1, Δρ<sub>bio</sub> 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 to- wards adjacent cells.

3. If all the neighboring cells are biofilms and have a biofilm vol- ume fraction equal or greater than 1, Δρ<sub>bio</sub> 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 ap- proach employed for solute transport. However, we must note that this algorithm does not allow to specify the shape of the fluid- biofilm 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 follow- ing that these fictitious interfaces are parallel to the mesh facets and perpendicular to the direction of growth.

3.5. 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 as- sume 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 τ<sub>max</sub> (Choi and Morgenroth, 2003; Stood- ley 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 φ<sub>det</sub>) 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 2D for instance:

$$
\tau = \max(\tau_{\text{north}}, \tau_{\text{south}}, \tau_{\text{east}}, \tau_{\text{west}}) \qquad (22)
$$

$$
\tau = \mu_f \times \max \left( \frac{\partial u_y}{\partial x}, \frac{\partial u_x}{\partial y} \right) \quad (23)
$$

<figure><img src="image_3.png" /><figcaption>Fig. 2. Cell (i, j) with its four neighbors and their velocity vectors for calculating the shear stress in 2D.</figcaption></figure>

where the cells east, west, south and north represent the neigh- boring 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<sub>stress</sub>. Thus, for the shear stress on the upper face, i.e., the north face, we obtain:

$$
\tau_{\text{north}} = \mu_f \times \frac{u_x(i, j+1) - 0}{l_{\text{stress}}} \qquad (24)
$$

In the particular instance illustrated in Fig. 2, l<sub>stress</sub> is given by:

$$
l_{\text{stress}} = \left( \frac{3}{2} - \phi(i, j) \right) \Delta y \qquad (25)
$$

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 re- cover l<sub>stress</sub> = Δy/2.

To model the process of detachment, we are based on the ex- perimental 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 em- pirical formulation relating the detached biofilm volume fraction φ<sub>det</sub> to the shear stress τ:

$$
\phi_{\text{det}} = 0 \quad \text{for } \tau < \tau_{\text{max}} \qquad (26)
$$

$$
\phi_{\text{det}} = e^{C_{\text{det}} \frac{(\tau - 2 \times \tau_{\text{max}})}{\tau_{\text{max}}}} \quad \text{for } \tau_{\text{max}} \le \tau \le 2 \times \tau_{\text{max}} \qquad (27)
$$

$$
\phi_{\text{det}} = 1 \quad \text{for } \tau \ge 2 \times \tau_{\text{max}} \qquad (28)
$$

with C<sub>det</sub> an empirical constant. The remaining volume fraction of biomass after detachment can be derived easily as follows (we as- sume that the biomass keeps a constant density):

$$
\phi^{t+1} = (\phi^t - \phi_{\text{det}}) \qquad (29)
$$

such as φ<sup>t</sup> 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.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Schematic representation of stages of bacterial attachment.</figcaption></figure>

## 3.6. 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., $\tau \le \tau_{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. 3a). 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 ($\tau \le \tau_{attach}$), then the virtual biofilm particle is attached, i.e., the fluid cell is converted into partially-filled cell (Fig. 3b). 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.

## 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; 2000; 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-

<figure><img src="image_5.png" /><figcaption>Fig. 4. General algorithm of model.</figcaption></figure>

ditions and verify that the model captures the dynamics of the developing biofilm pattern.

Finally, we apply our model to a 2D 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:

$$x' = \frac{x}{l}, \quad z' = \frac{z}{l}, \quad t' = t \frac{D_f}{l^2}, \quad c_i' = \frac{c_i}{C_0}, \quad (i = f, b) \quad (30)$$

$$p' = \frac{pl}{\mu_f u_{moy}}, \quad \mathbf{u}' = \frac{\mathbf{u}}{u_{moy}}, \quad \rho'_{bio} = \frac{\rho_{bio}}{\rho_{biomax}}, \quad \tau' = \frac{\tau l}{\mu_f u_{moy}} \quad (31)$$

$$K' = \frac{K}{C_0}, \quad \mu'_{max} = \mu_{max} \frac{l^2}{D_f}, \quad k'_{d} = k_d \frac{l^2}{D_f}, \quad D'_{i} = \frac{D_i}{D_f} \quad (32)$$
<figure><img src="image_6.png" /><figcaption>**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.)</figcaption></figure>

$\rho_{biomax}$, $C_0$, $u_{moy}$ and $l$ respectively represent the maximum biomass concentration, the substrate concentration imposed at the inlet, the fluid averaged velocity and the characteristic length of the medium. We also introduce the following dimensionless numbers, Pe and Da that characterize mass transport and bacterial reactivity and the dimensionless relative shear stress, denoted by $\tau^*$:

$$ \text{Péclet number: } Pe = \frac{u_{moy} \cdot l}{D_f} \qquad (33) $$

$$ \text{Damköhler number: } Da = \frac{\mu_{\max} \cdot l^2 \cdot \rho_{biomax}}{D_f K} = \frac{\mu'_{\max}}{K'} \Gamma \qquad (34) $$

$$ \text{Dimensionless shear stress: } \tau^* = \frac{\tau_\infty}{\tau_{\max}} \qquad (35) $$

where $\tau_\infty$ represents the maximum shear stress value encountered in the whole domain.

### 4.1. Case 1. Biofilm growth under hydrostatic conditions

The computational domain $\Omega = [0, lx] \times [0, lz]$ is discretized with a uniform mesh size $(Nx - 1) \times (Nz - 1)$ and a resolution $\Delta x = \frac{l_x}{Nx-1} = \Delta z = \frac{l_z}{Nz-1}$ 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. 5a). 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 $\rho'_{bio}(t' = 0) = 1$. A constant concentration $C_0$ in nutrient is imposed on the upper limit of the domain,

<table><caption>**Table 1**
Simulation data values.</caption><thead><tr><th>Parameter</th><th>value</th></tr></thead><tbody><tr><td>lx'</td><td>3</td></tr><tr><td>lz'</td><td>1</td></tr><tr><td>Nx</td><td>151</td></tr><tr><td>Nz</td><td>51</td></tr><tr><td>n<sub>0</sub></td><td>70</td></tr><tr><td>ρ<sub>bio</sub></td><td>1</td></tr><tr><td>D<sub>f</sub></td><td>1</td></tr><tr><td>D<sub>b</sub></td><td>0.25</td></tr><tr><td>F<sub>λ</sub></td><td>0.7</td></tr><tr><td>k'<sub>d</sub></td><td>5 × 10<sup>-3</sup></td></tr><tr><td>Γ</td><td>300</td></tr></tbody></table>

i.e., at $z' = 1$. Periodic boundary conditions are imposed in the $x$-direction 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 $\mu_{\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 isocentric concentrations are practically horizontal along the $x$-direction. We observe a homogeneous development of biofilm which occupies practically the entire solid surface (Fig. 6a). The biofilm growth is progressively causing a change in the availability of nutrient (Fig. 6b).
<figure><img src="image_7.png" /><figcaption>**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.)</figcaption></figure>

This disturbance of the concentration field intensifies progressively and destabilizes the biomass development (Fig. 6b 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. 6d), leads to the formation of finger-like structures (Fig. 6c 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 \cdot 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. 7a-right). This observation is corroborated by Fig. 7a-left, where the average concentration $\bar{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. 7a), 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. 7b 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 $\bar{\rho}$ is also more spread out over the z-direction and decreases when approaching the solid surface where the most mature bacteria layers are affected by the cellular extinction due to nutrient depletion.

## 4.2. Case 2. Biofilm growth under hydrodynamic conditions

### 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. 5b). These conditions are expressed as follows:

$$u_x'(0, z') = (6z/lz - 6z^2/lz^2), \quad z' \in [0, 1] \qquad (36)$$

$$u_z'(0, z') = 0, \quad z' \in [0, 1] \qquad (37)$$

$$\frac{\partial u_x'}{\partial z'}(x', 1) = 0, \quad x' \in [0, 1] \qquad (38)$$

$$u_z'(x', 1) = 0, \quad x' \in [0, 1] \qquad (39)$$

$$\frac{\partial u_x'}{\partial x'}(1, z') = 0, \quad z' \in [0, 1] \qquad (40)$$

$$u_z'(1, z') = 0, \quad x' \in [0, 1] \qquad (41)$$

$$\left. \frac{\partial C}{\partial n} \right|_{z'=1} = 0 \qquad (42)$$

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. 8c), 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. 8a 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
<figure><img src="image_8.png" /><figcaption>**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 $\bar{C}$ in red, and the average biomass concentration $\bar{\delta}$ in green. (a) $Da = 0.05$, $t' = 600$. (b) $Da = 10$, $t' = 160$. (c) $Da = 1000$, $t' = 46$. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

growth of colonies closest to the source. These colonies create an obstacle to the flow and change the streamline distribution. These numerical results are in perfect agreement with those obtained by Picioreanu et al. (2000).

### 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 \times Nz = 100 \times 33$ and $80 \times 26$. The Péclet number is $Pe = 50$ with a critical shear stress $\tau'_{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. 9a and b). The analysis of the increase of biomass concentration versus time

(Fig. 9d) 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. 9c and d) highlights the role of this regulatory mechanism on biomass development. Impact of detachment becomes clearly visible after $t' = 22$ (Fig. 9d). After $t' = 32$, the biofilm height reaches a plateau and stabilizes around 0.4 (Figs. 9a and b) while the biomass continues to grow infinitely in the absence of bacterial detachment, at least if the substrate availability is sufficient (Fig. 9c). 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 $\tau_{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.
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Biofilm growth patterns for $Da = 1$ and different $Pe$ values. The black lines represent the isoconcentrations with a variation of 10% between each one. The velocity vectors are also indicated. (a) $Pe = 0.1$, $t' = 1400$, (b) $Pe = 1$, $t' = 1100$, (c) $Pe = 100$, $t' = 450$.</figcaption></figure>

### 4.3. 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 \times lz = 12 \times 12\text{mm}^2$ with an initial porosity $\epsilon_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 $\Gamma$ and $F_\lambda$ 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 \times 10^{-8} \text{s}$ and the half-saturation constant $K$ varies between 0.0032 mg/l and 0.32 mg/l. The maximum growth rate $\mu_{\text{max}}$ was imposed via the Damköhler number.

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Biofilm growth for $Pe = 50$ and $Da = 1$. The two black lines represent the biomass concentration contour, solid line at $\rho'_{\text{bio}} = 0.1$, dashed line at $\rho'_{\text{bio}} = 0.9$ (a) Biofilm growth with bacteria detachment at $t' = 35$ (mesh 100 $\times$ 33), (b) Biofilm growth with bacteria detachment at $t' = 35$ (mesh 80 $\times$ 26), (c) Biofilm growth without bacteria detachment at $t' = 35$. (d) temporal variation of biofilm volume fraction.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**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.)</figcaption></figure>

<caption>**Table 2**</caption>

Physical and biological parameter values used in the simulations.

<table><thead><tr><th>Parameter</th><th>physical value</th><th>dimensionless value</th></tr></thead><tbody><tr><td>lx</td><td>12 mm</td><td>lx' = 1</td></tr><tr><td>lz</td><td>12 mm</td><td>lz' = 1</td></tr><tr><td>Nx</td><td>150</td><td></td></tr><tr><td>Nz</td><td>150</td><td></td></tr><tr><td>ρ<sub>bio</sub> (t = 0)</td><td>10<sup>4</sup> mg/l</td><td>ρ'<sub>bio</sub> (t' = 0) = 1</td></tr><tr><td>C<sub>0</sub></td><td>32 mg/l</td><td>C<sub>0</sub> = 1</td></tr><tr><td>D<sub>f</sub></td><td>10<sup>-9</sup> m<sup>2</sup>/s</td><td>D<sub>f</sub> = 1</td></tr><tr><td>D<sub>b</sub></td><td>0.25 × 10<sup>-9</sup> m<sup>2</sup>/s</td><td>D'<sub>b</sub> = 0.25</td></tr><tr><td>F<sub>λ</sub></td><td>0.7</td><td></td></tr><tr><td>k<sub>d</sub></td><td>3.4 × 10<sup>-8</sup> s</td><td>k'<sub>d</sub> = 5 × 10<sup>-3</sup></td></tr><tr><td>K</td><td>0.0032 – 0.32 mg/l</td><td>K' = 10<sup>-4</sup> – 10<sup>-2</sup></td></tr><tr><td>Γ</td><td>300</td><td></td></tr><tr><td>τ<sub>max</sub></td><td>1.39 Pa</td><td></td></tr><tr><td>C<sub>det</sub></td><td>5</td><td></td></tr></tbody></table>

Note that the dimensionless shear stress value, $\tau^*$, implicitly related to the Péclet number (through the pore velocity magnitude) vary also with the numerical computations performed ($\tau_{max}$ is kept constant). From all the simulations carried out, three pairs of Pe-Da values have been chosen as particularly representative of the different growth regimes obtained: Pe = 1 - Da = 10⁻² (low Pe, low Da, low $\tau^*$), Pe = 10⁻¹ - Da = 10² (low Pe, high Da, low $\tau^*$) and Pe = 500 - Da = 10⁻¹ (high Pe, low Da, high $\tau^*$). The results are discussed below.

### 4.3.1. Low Pe, low Da, low $\tau^*$

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. 11a) and the biomass grows in a relative uniform manner around the solid grains (Fig. 11c).

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. 11b) 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. 11c 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).

### 4.3.2. Low Pe, high Da, low $\tau^*$

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. 12a, we have a very large concentration gradient on the first cells of domain and a concentration that drops
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Simulation results for Pe = 1, Da = 0.01. (a) substrate concentration field at t' = 100. (b) biomass concentration field at t' = 100. (c) geometry of the medium at t' = 100. (d) geometry of the medium at t' = 9800.</figcaption></figure>

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. 12b 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. 12a). Bacteria extinction begins to appear far from the injection zone (Fig. 12b).

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. 12c 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.

#### 4.3.3. *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. 13a).

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
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 12.** Simulation results for Pe = 0.001, Da = 100. (a) substrate concentration field at *t*' = 1, (b) biomass concentration field at *t*' = 1, (c) medium geometry at *t*' = 1, (d) medium geometry at *t*' = 2.</figcaption></figure>

surface (Fig. 13e). As soon as the shear stresses exceed the imposed critical value τmax (i.e., τ* > 1), the biomass is subjected to ero-sion (Fig. 13c). In stagnant flow regions, located between several solid grains, the clogging becomes almost complete (Fig. 13d). Af-ter 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. 13b). The pref-erential 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).

4.4. 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 to-tal 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.

4.5. 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
<figure><img src="image_18.png" /><figcaption>**Fig. 13.** Simulation results for Pe = 500 and Da = 0.1. (a) Substrate concentration at *t*' = 70. (b) substrate concentration at *t*' = 140. (c) medium geometry at *t*' = 70. (d) medium geometry at *t*' = 140. (e) shear stress at *t*' = 700. (f) shear stress at *t*' = 140.</figcaption></figure>
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 14. Temporal variation of biofilm volume fraction and permeability. (a)-(b) Pe = 0.001, Da = 100. (c)-(d) Pe = 1, Da = 0.01. (e)-(f) Pe = 500, Da = 0.1.</figcaption></figure>

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

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

the REV by Darcy's law:

$$K = u_{moy} \frac{\mu_f \cdot L}{dP}, \qquad (43)$$

where K is the permeability, L the length of the domain (here equal to l<sub>x</sub>) and dP the pressure difference.

Fig. 14 shows the temporal evolution of both the biofilm volume fraction ε<sub>b</sub> and the dimensionless permeability K' = K/K<sub>0</sub>, where K<sub>0</sub> 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 (ρ'<sub>bio</sub> = 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. 13b and d) and drops of sev-
<figure><img src="image_25.png" /><figcaption>Fig. 15. Relationships between permeability and biomass volume fraction.</figcaption></figure>

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. 13f) 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 $\epsilon_b$ is limited by the bioclogging (Fig. 13a 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. 13e).

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,

$$ 
\begin{aligned}
K' &= C_s \frac{\epsilon^3}{(1 - \epsilon)^2} \\
    &= C_s \frac{(\epsilon_0 - \epsilon_b)^3}{(1 - \epsilon_0 + \epsilon_b)^2},
\end{aligned} 
\qquad (44) 
$$

where $C_s$ represents a geometrical shape factor, $\epsilon_0$ the initial porosity of the domain and $\epsilon_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.

## 5. 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 $\tau^*$. Thus, for *Pe* << *Da* and $\tau^* < 1$, the reaction rate-limited growth regime takes place. The mass transfer-limited growth regime occurs when *Pe* >> *Da* and $\tau^* < 1$. Finally, if $\tau^* > 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).

## Acknowledgments

This work was supported by the French National Research Agency (ANR) through the MOBIOPOR project, with the reference ANR-10-BLAN-0908 and was performed as part of the French Scientific Interest Group-Industrial Wasteland (GISFI) program.

## References

Al-Bader, D., Kansour, M., Rayan, R., Radwan, S., 2013. Biofilm comprising phototrophic, diazotrophic, and hydrocarbon-utilizing bacteria: a promising consortium in the bioremediation of aquatic hydrocarbon pollutants. Environ. Sci. Pollut. Res. 20 (5), 3252-3262.

Alpkvist, E., Klapper, I., 2007. A multidimensional multispecies continuum model for heterogeneous biofilm development. Bull. Math. Biol. 69 (2), 765-789.

Alpkvist, E., Klapper, I., 2007. Description of mechanical response including detachment using a novel particle model of biofilm/flow interaction. Water Sci. Technol. 55, 265-273.

Alpkvist, E., Picioreanu, C., van Loosdrecht, M., Heyden, A., 2006. Three-dimensional biofilm model with individual cells and continuum EPS matrix. Biotechnol. Bioeng. 94, 961-979.

Batstone, D., Picioreanu, C., van Loosdrecht, M., 2006. Multidimensional modelling to investigate interspecies hydrogen transfer in anaerobic biofilms. Water Resour. 40, 3099-3108.

Benioug, M., Golfier, F., Tinet, A.-J., M.A. B., Oltéan, C., 2015. Numerical efficiency assessment of IB-LB method for 3d pore-scale modeling of flow and transport. Transp. Porous Media 109, 1-23.

Bottero, S., Storck, T., Heimovaara, T.J., van Loosdrecht, M.C.M., Enzien, M.V., Picioreanu, C., 2013. Biofilm development and the dynamics of preferential flow paths in porous media. Biofouling 29 (9), 1069-1086.

Brackman, G., Hillaert, U., Calenbergh, S.V., Nelis, H., Coenye, T., 2009. Use of quorum sensing inhibitors to interfere with biofilm formation and development in burkholderia multivorans and burkholderia cenocepacia. Res. Microbiol. 160 (2), 144-151.

Bunt, C., Jones, D., Tucker, I., 1993. The effects of pH, ionic strength and organic phase on the bacterial adhesion to hydrocarbons (bath) test. Int. J. Pharm. 99 (2-3), 93-98.

Chmielewski, R., Frank, J., 2003. Biofilm formation and control in food processing facilities. Compr. Rev. Food Sci. Food Saf. 2 (1), 22-32.

Choi, Y., Morgenroth, E., 2003. Monitoring biofilm detachment under dynamic changes in shear stress using laser-based particle size analysis and mass fractionation. Water Sci. 47, 69-76.

Colasanti, R., 1992. Cellular automata models of microbial colonies. Binary, Comput.Microbiol. 4, 191-193.

Costerton, J., Caldwell, D., Korber, D., Lappin-Scott, H., 1994. Microbial biofilms. Ann. Rev. Microbiol. 49, 711-745.

Cowan, M., Warren, T., Fletcher, M., 1991. Mixed species colonization of solid surfaces in laboratory biofilms. Biofouling 3, 23-34.

Cunningham, A., Characklis, W., Abedeen, F., D., C., 1991. Influence of biofilm accumulation on porous media hydrodynamics. Environ. Sci. Technol. 25, 1305-1311.

Davies, D., Parsek, M., Pearson, J., Iglewski, B., Costerton, J., Greenberg, E., 1998. The involvement of cell-to-cell signals in the development of a bacterial biofilm. Science 280 (5361), 295-298.

Davit, Y., Iltis, G., Debenest, G., Veran-Tissoires, S., Wildenschild, D., Gerino, M., Quintard, M., 2011. Imaging biofilm in porous media using x-ray computed tomography. J. Microsc. 242 (1), 15-25.

Decho, A., 2000. Microbial biofilms in intertidal systems: an overview. Cont. Shelf Res. 20 (10-11), 1257-1273.

Derlon, N., Coufort-Saudejaud, C., Queinnec, I., Paul, E., 2013. Growth limiting conditions and denitrification govern extent and frequency of volume detachment of biofilms. Chem. Eng. J. 218, 368-375.

Dupin, H.J., Kitanidis, P.K., McCarty, P.L., 2001. Simulation of two-dimensional modeling of biomass aggregate growth in network models. Water Resour. Res. 37 (12), 2981-2994.

Eberl, H., Parker, D., van Loosdrecht, M., 2001. A new deterministic spatio-temporal continuum model for biofilm development. J. Theor. Med. 3, 161-175.

Ebigbo, A., Golfier, F., Quintard, M., 2013. A coupled pore-scale model for methanogenic microbial activity in underground hydrogen storage. Adv. Water Resour. 61, 74-85.

Ebigbo, A., Helmig, R., Cunningham, A.B., Class, H., Gerlach, R., 2010. Modelling biofilm growth in the presence of carbon dioxide and water flow in the subsurface. Adv. Water Resour. 33 (7), 762-781.

Frederick, M.R., Kuttler, C., Hense, B.A., Eberl, H.J., 2011. A mathematical model of quorum sensing regulated EPS production in biofilm communities. Theor. Biol. Med. Modell. 8. 8-8

Fujikawa, H., 1994. Diversity of the growth patterns of bacillus subtilis colonies on agar plates. FEMS Microbiol. Ecol. 13 (3), 159-168.

Golfier, F., Wood, B., Orgogozo, L., Quintard, M., Bus, M., 2009. Biofilms in porous media: development of macroscopic transport equations via volume averaging with closure for local mass equilibrium conditions. Adv. Water Resour. 32 (3), 463-485.

Golfier, F., Zarcone, C., Bazin, B., Lenormand, R., Lasseux, D., Quintard, M., 2002. On the ability of a Darcy-scale model to capture wormhole formation during the dissolution of a porous medium. J. Fluid Mech. 457, 213-254.

Hellweger, F.L., Bucci, V., 2009. A bunch of tiny individuals-individual-based modeling for microbes. Ecol. Modell. 220 (1), 8-22.

Hellweger, F.L., Clegg, R.J., Clark, J.R., Plugge, C.M., Kreft, J.-U., 2016. Advancing microbial sciences by individual-based modelling. Nat. Rev. Microbiol. 14, 461-471.

Hermanowicz, S., 1998. A model of two-dimensional biofilm morphology. Water Sci. Technol. 37, 219-222.

Hermanowicz, S., 2001. A simple 2D biofilm model yields a variety of morphological features. Math. Biosci. 169, 1-14.

Horn, H., Hempel, D.C., 1997. Growth and decay in an auto-/heterotrophic biofilm. Water Res. 31 (9), 2243-2252.

Huang, T., Xu, J., Cai, D., 2011. Efficiency of active barriers attaching biofilm as sediment capping to eliminate the internal nitrogen in eutrophic lake and canal. J. Environ. Sci. 23 (5), 738-743.

Kalmokoff, M., Austin, J., Wan, X.-D., Sanders, G., Banerjee, S., Farber, J., 2001. Adsorption, attachment and biofilm formation among isolates of listeria monocytogenes using model conditions. J. Appl. Microbiol. 91 (4), 725-734. ISSN 1365-2672

Kang, S.K., Hassan, Y.A., 2011. A comparative study of direct-forcing immersed boundary-lattice Boltzmann methods for stationary complex boundaries. Int. J. Numer. Meth. Fl 66 (9), 1132-1158. ISSN 1097-0363

Kapellos, G., Alexiou, T., Payatakes, A., 2007. Hierarchical simulator of biofilm growth and dynamics in granular porous materials. Adv. Water Resour. 30 (6-7), 1648-1667.

Khassekhhan, H., Efendiev, M., Eberl, H., 2009. A degenerate diffusionreaction model of an ammalistic biofilm control system : existence and simulation of solutions. Discrete Continuous Dyn. Syst. Ser. B 12, 371-388.

Kim, G., Lee, S., Kim, Y., 2006. Subsurface biobarrier formation by microorganism injection for contaminant plume control. J. Biosci. Bioeng. 101, 142-148.

Knutson, C.E., Werth, C.J., Valocchi, A.J., 2005. Pore-scale simulation of biomass growth along the transverse mixing zone of a model two-dimensional porous medium. Water Resour. Res. 41 (7), W07007. ISSN 1944-7973

Kone, T., Golfier, F., Orgogozo, L., Oltean, C., Lefevre, E., Block, J., Buës, M., 2014. Impact of biofilm-induced heterogeneities on solute transport in porous media. Water Resour. Res. 50, 9103-9119.

Kreft, J.-U., Booth, G., Wimpenny, J., 1998. BacSim, a simulator for individual-based modelling of bacterial colony growth. Microbiology 144, 3275-3287.

Kreft, J.-U., Picioreanu, C., Wimpenny, J.T., van Loosdrecht, M., 2001. Individual-based modelling. Microbiology 147 (11), 2897-2912.

Lackner, S., Terada, A., Smet, B., 2008. Heterotrophic activity compromises autotrophic nitrogen removal in menbrane aerated biofilm: results of a modeling study. Water Res. 42, 1102-1112.
Laspidou, C., Kungolos, A., Samaras, P., 2010. Cellular-automata and individual-based approaches for the modeling of biofilm structures: pros and cons. Desalination 250 (1), 390-394.

Laspidou, C.S., Liakopoulos, A., Spiliotopoulos, M.G., 2012. A 2D Cellular Automaton Biofilm Detachment Algorithm. In: Sirakoulis, G.C., Bandini, S. (Eds.). Cellular Automata. ACRI 2012. Lecture Notes in Computer Science, vol. 7495. Springer, Berlin, Heidelberg.

Laspidou, C., Rittmann, B., 2004. Modeling the development of biofilm density including active bacteria, inert biomass, and extracellular polymeric substance. Water Resour. 38, 3349-3361.

Laspidou, C., Rittmann, S., Karaman, B., 2005. Finite element modeling to expand the UMCCA model to describe biofilm mechanical behavior. Water Sci. Technol. 52, 161-166.

Liao, Q., Wang, Y., Wang, Y., Chen, R., Zhu, X., Pu, Y., Lee, D., 2012. Twodimension mathematical modeling of photosynthetic bacterial biofilm growth and formation. Int. J. Hydrogen Energy 37, 15607-15615.

Mah, T., O'Toole, G., 2001. Mechanisms of biofilm resistance to antimicrobial agents. Trends Microbiol. 9, 39-41.

Matsumoto, S., Katoku, M., Saeki, G., Terada, A., Aoi, Y., Tsuneda, S., 2010. Microbial community structure in autotrophic nitrifying granules characterized by experimental and simulation analyses. Environ. Microbiol. 12, 192-206.

Matsumoto, S., Terada, A., Aoi, Y., Tsuneda, S., Alpkvist, E., Picioreanu, C., 2007. Experimental and simulation analysis of community structure of nitrifying bacteria in a membrane-aerated biofilm. Water Sci. Technol. 55, 283-290.

MeeGee, R., Kinoshita, S., Fredrickson, A., Tsuchiya, H., 1970. Differentiation and product formation in molds. Biotechnol. Bioeng. 12, 771-801.

Morgenroth, E., Wilderer, P.A., 2000. Influence of detachment mechanisms on competition in biofilms. Water Res. 34 (2), 417-426.

Noguera, D., Pizarro, G., Stahl, D., Rittmann, B., 1999. Simulation of multispecies biofilm development in three dimensions. Water Sci. Technol. 39 (7), 123-130.

Paul, E., Ochoa, J., Pechaud, Y., Liu, Y., Liné, A., 2012. Effect of shear stress and growth conditions on detachment and physical properties of biofilms. Water Res. 46 (17), 5499-5508.

Peszynska, M., Trykozko, A., Iltis, G., Schlueter, S., Wildenschild, D., 2016. Biofilm growth in porous media: experiments, computational modeling at the porescale, and upscaling. Adv. Water Resour. 95, 288-301.

Pfeiffer, T., Schuster, S., Bonhoeffer, S., 2001. Cooperation and competition in the evolution of ATP-producing pathways. Science 292, 504-507.

Picioreanu, C., Batstone, D.J., van Loosdrecht, M., 2005. Multidimensional modelling of anaerobic granules. Water Sci. Technol. 52, 501-507.

Picioreanu, C., Kreft, J.U., van Loosdrecht, M.C.M., 2004. Particle-based multidimensional multispecies biofilm model. Appl. Environ. Microbiol. 70 (5), 3024-3040.

Picioreanu, C., van Loosdrecht, M., Heijnen, J., 1998. Mathematical modeling of biofilm structure with a hybrid differential-discrete cellular automaton approach. Biotechnol. Bioeng. 58, 101-106.

Picioreanu, C., van Loosdrecht, M., Heijnen, J., 1998. A new combined differential-discrete cellular automaton approach for biofilm modeling: application for growth in gel beads. Biotechnol. Bioeng. 57, 718-731.

Picioreanu, C., van Loosdrecht, M., Heijnen, J., 2000. Effect of diffusive and convective substrate transport on biofilm structure formation: a two-dimensional modeling study. Biotechnol. Bioeng. 69, 504-515.

Picioreanu, C., van Loosdrecht, M.C.M., Heijnen, J.J., 1998. Mathematical modelling of biofilm structure with a hybrid differential-discrete cellular automaton approach. Biotechnol. Bioeng. 58, 101-102.

Picioreanu, C., Xavier, J., van Loosdrecht, M., 2004. Advances in mathematical modeling of biofilm structure. Biofilms 1, 1-13.

Potera, C., 1999. Forging a link between biofilms and disease. Science 283, 1837-1839.

Prakash, B., Veeregowda, B.M., Krishnappa, G., 2003. Biofilms: a survival strategy of bacteria. J. Curr. Res. Sci. 85, 9-10.

Rupp, C., Fux, C., Stoodley, P., 2005. Viscoelasticity of staphylococcus aureus biofilms in response to fluid shear allows resistance to detachment and facilitates rolling migration. Appl. Environ. Microbiol. 71 (4), 2175-2178.

Scheuerman, T., Camper, A., Hamilton, M., 1998. Effects of substratum topography on bacterial adhesion. J. Colloid Interface Sci. 208 (1), 23-33.

Schindler, J., Rataj, T., 1992. Fractal geometry and growth models of a bacillus subtilis colony. Binary, Comput.Microbiol. 4, 66-72.

Schindler, J., Rovensky, L., 1994. A model of intrinsic growth of a bacillus colony. Binary, Comput.Microbiol. 6, 105-108.

Seo, Y., Lee, W.-H., Sorial, G., Bishop, P., 2009. The application of a mulch biofilm barrier for surfactant enhanced polycyclic aromatic hydrocarbon bioremediation. Environ. Pollut. 157 (1), 95-101.

Shaw, T., Winston, M., Rupp, C., Klapper, I., Stoodley, P., 1996. Commonality of elastic relaxation times in biofilms. Phys. Rev. Lett. 93 (9), 1-4.

Singh, R., Paul, D., Jain, R., 2006. Biofilms: implications in bioremediation. Trends Microbiol. 14, 389-397.

Skowlund, C., 1990. Effect of biofilm growth on steady-state biofilm models. Biotechnol. Bioeng. 35 (5), 502-510.

Smith, A., 2005. Biofilms and antibiotic therapy: is there a role for combating bacterial resistance by the use of novel drug delivery systems? Adv. Drug Deliv. Rev. 57 (10), 1539-1550.

Stoodley, P., Cargo, R., Rupp, C., Wilson, S., Klapper, I., 2002. Biofilm material properties as related to shear-induced deformation and detachment phenomena. J. Ind. Microbiol. Biotechnol. 29 (6), 361-367. ISSN 1367-5435

Stoodley, P., Wilson, S., Hall-Stoodley, L., Boyle, J., Lappin-Scott, H., Costerton, J., 2001. Growth and detachment of cell clusters from mature mixed-species biofilms. Appl. Environ. Microbiol. 67, 5608-5613.

Suchomel, B.J., Chen, B.M., Allen, M.B., 1998. Network model of flow, transport and biofilm effects in porous media. Transp. Porous Media 30, 1-23.

Sudarsan, R., Milferstedt, K., Morgenroth, E., Eberl, H., 2005. Quantification of detachment forces on rigid biofilm colonies in a roto-torque reactor using computational fluid dynamics tools. Water Sci. Technol. 52, 149-154.

Telgmann, U., Horn, H., Morgenroth, E., 2004. Influence of growth history on sloughing and erosion from biofilms. Water Res. 38 (17), 3671-3684.

Vega, L., Mathieu, J., Yang, Y., Pyle, B., McLean, R., Alvarez, P., 2014. Nickel and cadmium ions inhibit quorum sensing and biofilm formation without affecting viability in burkholderia multivorans. Int. Biodeterior. Biodegrad. 91 (0), 82-87.

von der Schulenburg, D., Pintelon, T., Picioreanu, C., van Loosdrecht, M., Johns, M., 2009. Three-dimensional simulations of biofilm growth in porous media. AIChE J. 55 (2), 494-504.

Walter, M., Safari, A., Ivankovic, A., Casey, E., 2013. Detachment characteristics of a mixed culture biofilm using particle size analysis. Chem. Eng. J. 228 (0), 1140-1147.

Wang, Q., Zhang, T., 2010. Review of mathematical models for biofilms. Solid State Commun. 150 (21 - 22), 1009-1022.

Wanner, O., Ebert, H., Morgenroth, E., Noguera, D., Picioreanu, C., Rittmann, B., van Loosdrecht, M., 2006. Mathematical modeling of biofilms. IWA Scientific and Technical Report No.18.

Wanner, O., Gujer, W., 1986. A multispecies biofilm model. Biotechnol. Bioeng. 28 (3), 314-328.

Wimpenny, J., Colasanti, R., 1997. A unifying hypothesis for the structure of microbial biofilms based on cellular automaton models. FEMS Microbiol. Ecol. 22 (1), 1-16.

Wingender, J., Neu, T., Flemming, H.-C., 1999. What are bacterial extracellular polymeric substances? In: Wingender, J., Neu, T., Flemming, H.-C. (Eds.) Microbial Extracellular Polymeric Substances. Springer Berlin Heidelberg, pp. 1-19.

Wood, B., Radakovich, K., Golfier, F., 2007. Effective reaction at a fluid - solid interface: application to biotransformation in porous media. Adv. Water Resour. 30, 1630-1647.

Wood, B., Whitaker, S., 1998. Diffusion and reaction in biofilms. Chem. Eng. Sci. 53 (3), 397-425.

Xavier, J., de Kreuk, M., Picioreanu, C., van Loosdrecht, M., 2007. Multi-scale individual-based model of microbial and bioconversion dynamics in aerobic granular sludge. Environ. Sci. Technol. 41, 6410-6417.

Xavier, J., Picioreanu, C., van Loosdrecht, M., 2005. A general description of detachment for multidimensional modelling of biofilms. Biotechnol. Bioeng. 91, 651-669.

Xavier, J., Picioreanu, C., Rani, S., van Loosdrecht, M., Stewart, P., 2005. Biofilm-control strategies based on enzymic disruption of the extracellular polymeric substance matrix - a modelling study. Microbiology 151, 3817-3832.

Zhang, K., Ou, M., Wang, W., Ling, J., 2009. Effects of quorum sensing on cell viability in streptococcus mutans biofilm formation. Biochem. Biophys. Res. Commun. 379 (4), 933-938.