\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}
\AtBeginDocument{{\noindent\small
Seventh Mississippi State - UAB Conference on Differential Equations and
Computational Simulations,
{\em Electronic Journal of Differential Equations},
Conf. 17 (2009), pp. 51--69.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2009 Texas State University - San Marcos.}
\vspace{9mm}}
\begin{document} \setcounter{page}{51}
\title[\hfilneg EJDE-2009/Conf/17\hfil A simple bioclogging model]
{A simple bioclogging model that accounts for spatial
spreading of bacteria}
\author[L. Demaret, H. J. Eberl, M. A. Efendiev, P. Maloszewski\hfil EJDE/Conf/17 \hfilneg]
{Laurent Demaret, Hermann J. Eberl, Messoud A. Efendiev, Piotr Maloszewski} % in alphabetical order
\address{Laurent Demaret, Messoud A. Efendiev \newline
Institute of Biomathematics and Biometry\\
Helmholtz Zentrum M\"unchen \\
German Research Center for Environmental Health\\
Ingolst\"adter Landstr 1, 85764 Neuherberg, Germany}
\email{laurent.demaret@helmholtz-muenchen.de}
\email{messoud.efendiyev@helmholtz-muenchen.de}
\address{Hermann J. Eberl\newline
Department of Mathematics and Statistics,
University of Guelph,
Guelph, ON, N1G 2W1, Canada}
\email{heberl@uoguelph.ca}
\address{Piotr Maloszewski \newline
Institute of Groundwater Ecology \\
Helmholtz Zentrum M\"unchen \\
German Research Center for Environmental Health\\
Ingolst\"adter Landstr 1, 85764 Neuherberg, Germany}
\email{maloszewski@helmholtz-muenchen.de}
\thanks{Published April 15, 2009.}
\subjclass[2000]{35K65, 35M10, 68U20, 76S05, 92D25}
\keywords{Bioclogging; biofilm; hydrodynamics; porous medium;
\hfill\break\indent
mathematical model; nonlinear-diffusion; simulation}
\begin{abstract}
An extension of biobarrier formation and bioclogging models is presented
that accounts for spatial expansion of the bacterial population in the soil.
The bacteria move into neighboring sites if locally almost all of the
available pore space is occupied and the environmental conditions are
such that further growth of the bacterial population is sustained.
This is described by a density-dependent, double degenerate
diffusion-equation that is coupled with the Darcy equations and a
transport-reaction equation for growth limiting substrates.
We conduct computational simulations of the governing differential
equation system.
\end{abstract}
\maketitle
\numberwithin{equation}{section}
\section{Introduction}
In soils, aquifers, and many other porous and fractured media
bacteria colonize in abundance. Typically they are sessile, i.e.
attached to the porous material and embedded in a self-produced
extracellular polymeric matrix (EPS) that protects them against
harmful environmental factors and mechanical washout. Naturally
occurring bacteria are a major contributor to the soil's
self-purification capacity. Environmental engineers use these
properties to devise microbially based technologies for soil
remediation or groundwater protection.
A growing bacterial population occupies more and more pore space
and thus alters the local hydraulic conductivity of the porous
medium (bioclogging) \cite{Baveye:1998}. Already moderate
bioclogging can change the local flow velocity and lead to the
development of preferred flow paths, mainly due to the EPS matrix
\cite{STHM:2006}. Changing flow velocities affect the convective
transport of dissolved substrates and thus the local availability
of nutrients and biocides, i.e. the growth conditions for the
bacteria \cite{Thullner:2002}. Thus, soil hydrodynamics, substrate
transport and bacterial population dynamics are strongly coupled
processes. This nonlinear phenomenon not only has implications for
the natural microbial ecology in soils but also has possible use
in engineering applications, cf \cite{Baveye:1998, STHM:2006} and
the references there. For example, biofilm forming bacteria can be
injected to soils and form biobarriers that are impermeable to
water flow and thus prevent contaminants to reach the groundwater.
In many instances these biobarriers are reactive, i.e. they
actively degrade the pollutant. In other engineering applications,
however, such as waste treatment, bioclogging can be detrimental
to process performance \cite{Baveye:1998}.
Mathematical modeling of hydrodynamics and pollutant transport in
porous media has a long history. Traditionally the focus has been
on physical processes, namely hydrodynamics and substrate
transport, while the microbiological aspects often are only
considered in a simplified manner. If the restricting assumption
of a bacterial population in equilibrium is given up, then the
bacterial populations are typically treated like planktonic
populations. On one end of the spectrum, they are assumed to be
suspended and subjected to the same transport mechanisms
(convection and Fickian diffusion) as the chemical species. On the
other hand, when the focus is on the bioclogging properties, it is
often assumed that bacteria are stationary, in the sense that
bacterial biomass may locally aggregate due to reactions but it
does not move into neighboring regions. In such models the
available pore-space is typically not explicitly considered as a
factor that limits the local population size; instead biomass
production is limited only by the usual growth kinetics, i.e.
substrate availability, cf. \cite{Schaefer:1998,Thullner:2004}.
With these models simulations are typically terminated before the
local pore space becomes completely clogged. Alternatively, ad hoc
growth limiters are sometimes used to enforce a shut-down of
biomass production if the available pore-space becomes limited,
e.g. \cite{Chen:2003}.
It seems more appropriate to consider bioclogging as a volume
filling problem. The microbes, immobilized in the EPS matrix,
fill up the space that is locally available. If the pore is close
to be completely filled but the environmental conditions are such
that microbial growth continues, the bacterial communities need to
expand and spill over into neighboring sites. In this paper we
propose a model to describe this spatial movement of bacteria in
porous media. Such volume filling problems can be described by
density-dependent diffusion processes \cite{KhHE:2009,PH:2002}. This leads
to a quasi-linear degenerate transport-reaction equation for the
bacterial population, which needs to be coupled with models for
flow field and substrate transport. The standard model for
macroscopic flow velocities in porous media are the Darcy
equations \cite{Bras:1990} and this is what we use. In a nut-shell
they say that the flow velocity is proportional to the pressure
gradient, where the proportionality factor is the hydraulic
conductivity, which in our case varies in space and time due to
biofilm growth. More formally, the Darcy equations can be derived
from the Stokes equation by homogenization \cite{Hornung:1997}.
Transport equations for dissolved substrates can be derived from
microscopic equations by homogenization or volume averaging
techniques. The particular difficulty for the system at hand is to
correctly account for the biofilm properties of the microbial
population. The development of such transport equations is
currently an active area of research. According to \cite{S:2003}
the most striking difference between biofilms and planktonic
bacterial populations is the diffusive limitation of substrates in
the EPS matrix. Starting from two separate transport equations for
biofilm and aqueous phase in the pore space, it was shown in
\cite{GWOQB:sub} that under a straightforward but simplifying
equilibrium assumptions this can be described by a classical
semi-linear convection-diffusion reaction model, if the parameters
are accordingly interpreted. Models of this type have been used in
bio-barrier studies previously, e.g. \cite{Chen:2003}, and we use
them here as well.
The resulting model is a system of nonlinear, degenerate partial
differential equations of mixed elliptic/parabolic type, which is
studied in computer simulations.
\section{Governing equations}
We propose the following model for bioclogging and biobarrier
formation in the computational domain $\Omega\in\mathbb{R}^d$
($d\in\{2,3\}$)
\begin{equation}\label{model}
\begin{gathered}
V = \frac{1}{\mu \rho}A(M)(F - \nabla P), \\
0 = \nabla V, \\
\mathfrak{P} \partial_t C = \nabla(D_C \nabla C - V C) + \mathfrak{P} \frac{\kappa_1 CM}{\kappa_2+C},\\
\mathfrak{P} \partial_t M = \nabla(D_M(M) \nabla M) +
\mathfrak{P} \big[\frac{\kappa_3 CM}{\kappa_2 +C}-\kappa_4 M\big]
\end{gathered}
\end{equation}
for the dependent variables
\begin{trivlist}
\item $V=V(t,x)$: Darcy velocity vector,
\item $P=P(t,x)$: pressure,
\item $C=C(t,x)$: concentration of growth limiting substrate,
\item $M=M(t,x)$: biomass density.
\end{trivlist}
The independent variables are
\begin{itemize}
\item $t\geq 0$: time,
\item $x\in \Omega$: space.
\end{itemize}
Model \eqref{model} is formulated such that all parameters are
non-negative. By $\mathfrak{P}$ we denote the fraction of the
pores per unit volume of the "empty" (i.e. absence of bacterial
biomass) porous medium, which we assume to be a positive constant.
Note that we distinguish here between three phases: the solid
phase, which occupies the volume fraction $1-\mathfrak{P}$, the
biofilm phase, which occupies the volume fraction $\mathfrak{P}
M/M_{\rm max}$ and the liquid (void) phase, occupying the volume
fraction $\mathfrak{P} - \mathfrak{P} M/M_{\rm
max}=\mathfrak{P}(1-M/M_{\rm max})$. Here $M_{\rm max}$ denotes
the maximum biomass density, i.e. is a measure for the maximum
number of cells that can fit into a unit volume.
{\it Hydrodynamics.} The first equation of \eqref{model} is the
Darcy equation, the standard model for flow in porous media
\cite{Bras:1990,Hornung:1997}. Constants $\mu$ and $\rho$ are the
dynamic viscosity and the density of the fluid, respectively.
Vector $F$ denotes body forces. We will use $F\equiv 0$ for
simplicity and in accordance with \cite{Chen:2003}. The second
equation of \eqref{model} describes conservation of water in the
porous medium. We can assume that this equation holds because
biofilms are largely composed of water. Thus, a growing biofilm
does not push water out of the system but assimilates it. This
argumentation was also used in \cite{Thullner:2004}.
In order to take the effect of bioclogging on the flow field into
account, the hydraulic conductivity $A$ must change locally with
the amount of bacterial biomass. In the absence of biomass, $M=0$,
it is the hydraulic conductivity of the "empty" porous medium,
$A(0)=A_0$. It decreases as $M$ increases, i.e. $A'(M)\leq 0$.
Several expressions have been obtained for $A$ in various
experiments. Some of them are summarized and compared in
\cite{Thullner:2004}; a more theoretical approach is taken in
\cite{SM:2001}. These models, along with others that are used in
the bioclogging literature, e.g. \cite{Chen:2003,HKP:2007}, agree
in that $A$ is in good approximation a cubic polynomial in a
three-dimensional setting (some authors use the exponent
$19/6\approx 3$, cf. \cite{HKP:2007,Thullner:2004}), and a
quadratic polynomial in a two-dimensional setting. The {\it
ansatz}
\begin{equation}\label{permeability}
A(M)=\begin{cases}
A_0\cdot \frac{\big(1-\frac{M}{M_{\rm clog}} \big)^b + a}{1+a},
& \text{if } M\leq M_{\rm clog}\\
A_0 \cdot \frac{a}{1+a}, & \text{if } M> M_{\rm clog}
\end{cases}
\end{equation}
is proposed in \cite{Thullner:2004} for bioclogging induced by
homogeneous biofilms. The main difference between this choice of
$A(M)$ and other proposed models is that $A$ is not entirely
reduced to 0 as the pores clog but still allows for a minimum flow
through the soil. This is described by the two parameters $M_{\rm
clog}$ (the biomass density beyond which no further reduction of
conductivity is observed; note that $M_{\rm clog} < \mathfrak{P}
M_{\rm max}$, we use $M_{\rm clog}= 0.9 \mathfrak{P} M_{\rm
max})$, and $a$, an adjustable parameter that relates the minimum
conductivity to the maximum conductivity (we use $a=0.05$).
{\it Substrate transport and degradation.} The third equation of
\eqref{model} describes convective and diffusive transport of the
limiting dissolved substrate (e.g. a nutrient) as well as its
degradation in biochemical reactions. The pore space can be viewed
as subdivided into an aqueous phase, with convective and diffusive
substrate transport, and a biofilm phase, with diffusive transport
and biochemical reactions. We model it here by a single equation,
as has been suggested by other authors before \cite{ADQ:2007,
Chen:2003, GWOQB:sub}, thereby implicitly assuming an equilibrium
between both phases. $D_C$ is the diffusion coefficient of the
substrate in the porous medium, where we assume that the substrate
does not diffuse through the solids in the porous medium. For
small molecules like oxygen or carbon, the difference of the
diffusion coefficient in water and in the biofilm matrix is small
\cite{BD:1998} and can often be neglected in biofilm modeling
\cite{TGR:2006}. We adopt this strategy here for the sake of
simplicity. Thus, no distinction is made between substrate
diffusion in the biofilm and in the water.
The coefficient $\kappa_1$ in \eqref{model} is the maximum
substrate consumption rate and $\kappa_3$ is the maximum specific
growth rate. The ratio $Y:=\kappa_3/\kappa_1$ is the yield
coefficient that indicates how much mass of substrate is required
to produce a unit mass of biomass. Typically we have $Y<1$, {\it
e.g.} for organic carbon. This means it takes more substrate than
new biomass is produced; the difference in mass is, for example,
oxidized to gain energy. The Monod half-saturation constant
$\kappa_2$ is the value for which $C=\kappa_2$ implies bacterial
growth at half the maximum growth rate. Note that the Monod
kinetics used here implies a saturation effect, in the sense that
even if substrate is available in abundance the specified maximum
growth cannot be exceeded.
{\it Biomass growth and spreading.} The last equation in
\eqref{model} finally describes formation and spreading of the
bacterial population in the porous medium. Production of new
biomass is closely related to substrate degradation as described
in the previous paragraph. A first order decay model is assumed
for cell loss; the parameter $\kappa_4$ is the cell death rate.
The spatial transport terms for biomass were introduced in the
last equation of \eqref{model} to allow for bacteria to spread
into neighboring regions if locally the pores become clogged but
biomass production continues. Similar volume filling problems have
been described by nonlinear diffusion equations \cite{PH:2002},
where the diffusion coefficient $D_M(M)$ depends on the local
population density. It increases with $M$, $D_M'(M)\geq 0$.
Moreover, the biomass should not move into neighboring regions
while space is available locally for more cells. We adopt the
following form for $D_M(M)$ that was first introduced in
\cite{Eberl:2001} for a mesoscopic biofilm model,
\begin{equation}\label{DMcoeff}
D_M(M)=d \frac{ \big( \frac{M}{\mathfrak{P} M_{\rm
max}}\big)^\alpha } {\big(1-\frac{M}{\mathfrak{P} M_{\rm
max}}\big)^\beta}, \quad 1\gg d>0, \; \alpha >1, \; \beta > 1
\end{equation}
where $\mathfrak{P} M_{\rm max}$ is the maximum biomass density
that can be attained in the pore. Accordingly, $M/(\mathfrak{P}
M_{\rm max})$ is the fraction of the unit volume that is occupied
by biomass. It should be pointed out that the behavior of such a
density-dependent diffusion is entirely different than the
behavior of the more familiar linear Fickian diffusion. For one,
for small biomass densities $M\ll \mathfrak{P} M_{\rm max}$, the
biomass equation in \eqref{model} behaves essentially like the
porous medium equation that describes flow in unsaturated soils.
In particular, unlike Fickian diffusion, this effect implies that
initial data with compact support lead to solutions with compact
support. More specifically it describes the formation of a sharp
interface between the region that is at least partly occupied with
biomass, $\Omega_2(t):\{ x\in \Omega: M(t,x)>0\}$, and the region
$\Omega_1(t)=\{x \in \Omega : M(t,x)=0\}$ that is free of biomass.
This interfaces moves at finite speed. On the other hand, for $M
\approx \mathfrak{P} M_{\rm max}$ the density-dependent diffusion
equation behaves like fast diffusion. This guarantees that the
maximum attainable biomass density $\mathfrak{P} M_{\rm max}$ is
not exceeded. It can be shown that $M$ remains separated from this
maximum value \cite{EEZ:2002,EEZ:2008}, i.e. there is a small constant
$\delta>0$, such that $\mathfrak{P} M_{\rm max}-\delta$. Thus, the
fast diffusion mechanism keeps the solution away from the
singularity. The interplay of both nonlinear diffusion effects can
be summarized as follows: As long as $M$ is small, the biomass
will compress locally but does not spread notably. Only as $M\to
\mathfrak{P} M_{\rm max}$, biomass starts to spread out locally
into neighboring regions in order to guarantee that the maximum
biomass density is not exceeded.
Note that both non-linear diffusion effects are needed to safely
guarantee this behavior. The porous medium power law $D(M)=d M^a$
alone is not able to ensure boundedness of the solution by the
physically maximal possible value $\mathfrak{P} M_{\rm max}$,
while the power law $D(M)=d(1-M/\mathfrak{P} M_{\rm max})^{-b}$
leads to an immediate dilution if $M\approx \mathfrak{P} M_{\rm
max}$ instead of an interface that propagates at finite speed.
The solution $M$ is continuous at the interface but not
necessarily differentiable. Thus, $M$ is to be understood as a
weak solution. The biomass equation with nonlinearity
(\ref{DMcoeff}) is a double-degenerate parabolic equation, which
is not yet well understood analytically. We refer to
\cite{DE:2006,ED:sub,EEZ:2008,EM:sub} for first rigorous results.
Including the spatial spreading terms for bacteria is a major
difference to other bioclogging models. In some existing
approaches, the biomass density is allowed to grow unbounded. Thus
eventually it attains locally unphysical values $M> \mathfrak{P}
M_{\rm max}$, which marks the breakdown of the model. One model of
this type is used in \cite{Thullner:2004}. However, the
simulations carried out in that study are terminated before this
breakdown situation occurs, in accordance with the experiments to
which they are compared. In the biobarrier formation model in
\cite{Chen:2003}, $M\to P M_{\rm max}$ somewhere in $\Omega$
triggers an artificial local growth limiter in order to avoid this
breakdown situation. In fact, comparing the equations for $C$ and
$M$ in our model \eqref{model} with the equations for carbon and
the biofilm former {\it Klebsiella oxytoca} in the single species
biobarrier model of \cite[cf. the 2nd and 3rd equations of model
(9)]{Chen:2003}, we note that \eqref{model} is a direct model
extension of this earlier biobarrier model. The general approach
to include spatial spreading of bacterial biomass as
density-dependent diffusion, however, can be included in other
bioclogging models as well, e.g. in \cite{Thullner:2004}. For
moderate biocloging, as long as $M$ remains clearly below
$\mathfrak{P} M_{\rm max}$, the three models \eqref{model},
\cite{Chen:2003}, \cite{Thullner:2004} approximate each other
quite well if the same hydraulic conductivity function $A(M)$ is
chosen.
Model \eqref{model} is a fully transient model,
despite the absence of time-derivatives in the flow equations. We
consider a creeping flow regime in the soil, where inertia effects
do not matter. The Darcy equations in the given form then can be
derived from the Stokes equation by homogenization techniques, as
carried out in \cite{Hornung:1997}. As the soil texture changes,
due to compression and expansion of biomass, the flow field
quickly adapts to the new conditions. Thus, both the flow velocity
$V$ and the pressure $P$ are indeed functions of $t$ and of $x$.
Following the argumentation given in \cite{Thullner:2004}, based
on \cite{Rittmann:1982} and the creeping flow conditions that we
impose, we can neglect detachment and re-attachment of biofilm in
our model in this first study.
In the form \eqref{model}, the model is formulated for general
two- or three-dimensional domains. Lateron in the numerical
simulations we will consider a two-dimensional rectangular domain
$\Omega = [0,L] \times [0,W]$ and drive the flow in the system by
applying a pressure difference at the inflow and outflow
boundaries. Substrate is added at a prescribed concentration at
inflow. This leads to the following boundary conditions to be
imposed on \eqref{model}
\[
\begin{gathered}
P \big|_{x_1=0}=P_\infty,\quad P\big|_{x_1=L}=P_0,\quad
\partial_n P \big|_{\partial \Omega \cap \{ 0 < x_1 < L\}}=0\\
V^Tn\big|_{\partial \Omega \cap \{0P_0$, and $C_\infty=C_\infty(x_2) \geq 0$.
Furthermore we provide the following initial data
\[
C(0,x)=C_0(x),\quad M(0,x)=M_0(x),\quad x\in \Omega,
\]
where we require $C_0 \geq 0$ and $0\leq M_0 < \mathfrak{P} M_{\rm
max}$. Note that no initial data are required for the hydrodynamic
quantities as these always follow the biomass distribution tightly
(quasi steady state assumption).
We re-scale \eqref{model} by introducing the dimensionless
variables
\begin{equation}
\begin{gathered}
\tilde x = \frac{x}{L},\quad \tilde t= \frac{t}{T},\\
p=\frac{P}{P_\infty-P_0},\quad v= \frac{TV}{L \mathfrak{P}}, \quad
c=\frac{C}{\| C_\infty\|_\infty}, \quad m=\frac{M}{\mathfrak{P}
M_{\rm max}}
\end{gathered}
\end{equation}
and, thus, $\tilde \nabla = L \nabla$ and $\partial_{\tilde t}=T
\partial_t$. Furthermore, we make the assumption $F=0$ and obtain
from \eqref{model} the dimensionless form
\begin{equation}\label{dlmodel}
\begin{gathered}
v = -k_0 a(m) \tilde \nabla p \\
0 = \tilde \nabla v \\
\partial_{\tilde t} c = \tilde \nabla \left(d_c \tilde \nabla c - vc\right)
- k_1 \frac{cm}{k_2+c} \\
\partial_{\tilde t} m = \tilde \nabla \left(d_M(m) \tilde \nabla m\right)
+ k_3 \frac{cm}{k_2+c} - k_4 m
\end{gathered}
\end{equation}
where the new dimensionless parameters and coefficients are
\begin{gather*}
d_C= \frac{D_C T}{\mathfrak{P} L^2},\quad d_m(m)= \frac{d T}{L^2}
\frac{m^\alpha}{(1-m)^\beta},\\
a(m)= \frac{A(m\mathfrak{P} M_{\rm max})}{A_0}, \quad
k_0 = \frac{(P_\infty -P_0) A_0 T}{\mu \rho\mathfrak{P}}, \\
k_1 = T \kappa_1 \frac{\mathfrak{P} M_{\rm max}}{\| C_\infty \|_\infty},
\quad k_2=\frac{\kappa_2}{\| C_\infty \|_\infty},
\quad k_3 = T \kappa_3,
\quad k_4 = T \kappa_4.
\end{gather*}
From now on we will refer always to the dimensionless equation
\eqref{dlmodel}. However, for the convenience of the notation we
will drop the tildes.
\section{Computer simulations}
\subsection{General setup}
We conduct a series of computer simulations to illustrate the
behavior of the proposed bioclogging model. The numerical solution
is based on the following standard re-formulation of
\eqref{dlmodel},
\begin{equation}\label{dlmodel2}
\begin{gathered}
v = -k_0 a(m) \nabla p \\
0 = \nabla \left(a(m) \nabla p \right)\\
\partial_{t} c = \nabla \left(d_c \nabla c - vc\right) - k_1 \frac{cm}{k_2+c} \\
\partial_{t} m = \nabla \left(d_M(m) \nabla m\right) + k_3 \frac{cm}{k_2+c} - k_4 m,
\end{gathered}
\end{equation}
which is obtained by taking the divergence of the first equation
of \eqref{dlmodel} and using the incompressibility condition
$\nabla v=0$, i.e. the second equation of \eqref{dlmodel}.
The integration of the equations for $c$ and $m$ follows the
approach that was introduced and studied in \cite{ED:2007} for
the underlying diffusion-reaction biofilm model. This is a
semi-implicit, finite difference based finite volume scheme,
formulated for the concentrations in the centers of the grid
cells, and implemented on a uniform rectangular grid. The main
feature is a non-local time discretisation. This reduces the
computational effort in every time-step to a linear system for
substrate concentration and one for biomass density. The time-step
size is variable and chosen such that $0\leq m < 1$ is enforced
for the numerical approximation.
For our current purpose this numerical scheme must be coupled with
a solver for the Darcy equation. For a given biomass
distribution, the second equation of (\ref{dlmodel2}) is a linear
(non-constant coefficient) elliptic problem for $p$ which is
solved in the cell centers as well, using standard second order
discretization \cite{MM:1994}. The velocity field is obtained from
$p$ by second order finite differences on the cell edges, i.e. on
a staggered grid.
Thus, in every time step, first $p$ is computed from the current
biomass distribution, then $v$. Finally the equations for $c$ and
$m $ are moved to the next time level, and the routine is
repeated.
In total, this procedure requires per time-step the solution of
three sparse (banded diagonal) linear algebraic systems, one for
each of $p, c, m$. The stabilized bi-conjugated gradient method is
used for this purpose.
In our simulations we choose the computational domain to be
rectangular, of size $\Omega=[0,1]\times[0,\frac{1}{2}]$. It is
discretized by a uniform grid of $400 \times 200$ cells. Order of
magnitude estimates for the physical and biological model
parameters were obtained from the hydrological and biofilm
modeling literature \cite{Bras:1990, Thullner:2004, TGR:2006},
leading to the following dimensionless parameters that were used
in the simulations
\begin{gather*}
k_0=2.94, \quad k_1=22.2, \quad k_2=0.1667, \quad k_3=1.0, \\\
k_4= 0.02, \quad d_C=0.002, \quad \alpha=\beta=4, \quad d=10^{-8}.
\end{gather*}
\subsection{Bioclogging initiated by a large single patch of bacterial
biomass}\label{x1_sec}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.42\textwidth]{fig1a}\; % x1001_1
\includegraphics[width=0.42\textwidth]{fig1b} \\[-5pt] %x1001_5
$t=0.01$ \hspace{5cm} $t=2$
\\[3pt]
\includegraphics[width=0.42\textwidth]{fig1c}\; % x1001_11
\includegraphics[width=0.42\textwidth]{fig1d} \\[-5pt] %x1001_15
$t=5$ \hspace{5cm} $t=7$
\\[3pt]
\includegraphics[width=0.42\textwidth]{fig1e}\; %x1001_20
\includegraphics[width=0.42\textwidth]{fig1f} \\[-5pt] % x1001_25
$t=9.5$ \hspace{5cm} $t=12 $
\end{center}
\caption{Bioclogging, induced by a spherical region with initially
low biomass density: flow is from bottom to top.
Shown are: $c$ (background) [left color map],
$p$ (equidistant white iso-lines);
$m$ (colored iso-lines) [right color map], for selected $t$.
}
\label{sphere_illustration_fig}
\end{figure}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.43\textwidth]{fig2a}\; % x1001_30
\includegraphics[width=0.43\textwidth]{fig2b}\\[-5pt] % x1001_35
$t=14.5 $ \hspace{5cm} $t=17$ \\[5pt]
\includegraphics[width=0.43\textwidth]{fig2c}; %x1001_37
\includegraphics[width=0.43\textwidth]{fig2d} \\[-5pt] %x1001_40
$t=18$ \hspace{5cm} $t=19.5$\\[5pt]
\includegraphics[width=0.43\textwidth]{fig2e}\; %x1001_43
\includegraphics[width=0.43\textwidth]{fig2f}\\[-5pt] %x1001_47
$t=21$ \hspace{5cm} $t=23$
\end{center}
\caption{Bioclogging, Fig. \ref{sphere_illustration_fig}
continued.}\label{sphere_illustration_fig2}
\end{figure}
In a first illustration of the bioclogging model, a spherical
region of radius $0.1$ in the center of the domain is inoculated
with bacterial biomass at the low density $m=0.01$. Thus,
initially only 6.3\% of the domain is occupied by a low density
biomass core. In Figure \ref{sphere_illustration_fig}, the biomass
density, substrate concentration and pressure are shown for
selected subsequent time instances. The concentration field is
shown color-coded in the background. For the biomass we plot
color-coded equidistant iso-lines, while the equidistant iso-lines
for the pressure are plotted in white. Note that the scale of the
biomass iso-lines changes, as the bacterial population increases,
while the pressure field always is confined to the interval
$[0,1]$. The curvature of the isobars indicates the direction of
the local flow velocity, due to the first of (\ref{dlmodel2}). In
the simulation two distinctive phases can be distinguished.
{\it Compression phase.} Initially ($t<5$) the biomass density
increases homogeneously in the spherical core and no spatial
spreading of biomass takes place. Even during this compression
phase the increasing biomass density increases the flow resistance
and the flow field seeks preferred flow directions around the
spherical core, as indicated by the curvature of the pressure
lines. Production of biomass is at the expense of nutrient
consumption. In the biomass core, substrate is transported by
diffusion and by convection (the convective contribution decreases
with increasing biomass density), and consumed by the bacteria.
Accordingly, substrate gradients are observed into the core and in
main flow direction. Downstream of the biomass core, the substrate
concentration is lower than in the biomass but nutrient is
initially replenished, both by convection (flow around the biomass
patch) and diffusion. At $t=5$ substrate in the downstream end of
the biomass core is decreased to about 50\% of the inflow
concentration value.
{\it Expansion phase.} Eventually the entire pore space in the
initially spherical region is almost completely occupied by
biomass and the biomass region starts expanding. Due to substrate
availability, the living conditions are best in the upstream rim
of this pocket and the biomass develops and expands faster there
than in the downstream region. In the growing biomass region
substrate demand increases and the concentration drops to low
values inside the domain quickly. As a consequence of the
drastically different and heterogeneous growth conditions inside
the biomass region, the spatial homogeneity is perturbed. In
particular in the downstream trail end the biomass density is
smaller than in the upstream rim. The expanding biomass region
leads also to an acceleration of the flow field around the pocket,
as implied by the stronger curvature of and smaller distance
between the isobars. We observe a flow induced substrate
concentration boundary layer that forms around the biomass region
and extends in the downstream region of the domain along the
boundaries all the way down to the outflow region, as long as free
flow around the biomass core is possible. In this concentration
boundary layer microbial activity is highest and biomass density
largest.
Eventually, at $t\approx 15$ the biomass core spans over the width
of the entire domain. The concentration boundary layer is now
found completely inside the biomass region, in the upstream rim.
The biomass front propagates upstream, toward the nutrient source.
The complete clogging of the flow channel implies that the flow
rate drops drastically, since the pressure difference between
inflow and outflow remains constant. Thus, substrate supply becomes
more and more diffusion dominated. Substrate becomes limited
inside the bacterial core over a short penetration distance,
leading to an initially almost stationary downstream interface.
Eventually the biomass front reaches the upstream boundary. Due to
favorable conditions ($c\approx 1$), more new biomass is produced
here than decays in the substrate depleted downstream zones. As a
consequence of space limitations, active biomass is forced
downstream. Eventually a biomass front travels in flow direction.
In order to quantify and summarize the effect of bioclogging on
the reactor scale, we plot in Figure \ref{sphere_lumped_fig} the
following lumped parameters as functions of time:
\begin{itemize}
\item {\it occupancy of the domain}, i.e. the region of the domain in
which bacteria can be found, regardless of the biomass density
\[
\omega(t):=\frac{1}{|\Omega|} \int_{\{x\in \Omega:
\;m(t,x)>0\}}\!\!\! dx
\]
\item the {\it total bacterial biomass} in the system
\[
M_{\rm tot}(t):=\int_\Omega m(t,x) dx
\]
\item the hydrodynamic {\it flow rate} on inflow and outflow
\[
Q_{in}(t):=\int_{0}^{1/2} v(t,0,y) dy,\quad
Q_{out}(t):=\int_{0}^{1/2} v(t,1,y) dy.
\]
Our flow model is based on a quasi-steady state assumption, due to
a time-scale separation as outlined above. Therefore, in every
time-step and for every biomass distribution the flow field $v$ is
stationary. Thus, we must have, $Q_{in}(t)\equiv Q_{out}(t)$. This
can be used as a criterion to test mass conservation properties of
the flow solver.
\item the {\it substrate flux} on inflow
\[
J_{in}(t):=\int_{0}^{L/2} v(t,0,y)c(t,0,y) dy - d_C
\partial_x c\big|_{x=0}
\]
and on outflow, which due to the Neuman boundary condition there reduces
to the convective flux
\[
J_{out}(t):=\int_{0}^{L/2} v(t,1,y)c(t,1,y) dy.
\]
The difference $J_{in}-J_{out}$ indicates how much substrate is degraded
by the bacteria.
\end{itemize}
\begin{figure}[ht] % figure 3
\begin{center}
\includegraphics[width=0.5\textwidth]{fig3a}\\ %x1MtotW
\includegraphics[width=0.5\textwidth]{fig3b}\\ %x1Flow
\includegraphics[width=0.5\textwidth]{fig3c} %x1Fluxes
\end{center}
\caption{Some lumped parameters for the simulation shown in Figure
\ref{sphere_illustration_fig}}\label{sphere_lumped_fig}
\end{figure}
Integrating the last equation of \eqref{dlmodel} over the
computational domain $\Omega$ yields the differential inequality
for $M_{\rm tot}(t)$
\[
\frac{ d M_{\rm tot}}{dt} < \Big(\frac{k_3}{k_2+1}-k_4\Big)
M_{\rm tot}
\]
and thus in particular the upper bound
\begin{equation}\label{MUBeq}
M_{\rm tot}(t)< M_{\rm tot}(0) e^{(\frac{k_3}{k_2+1}-k_4)t}
\end{equation}
where $\frac{k_3}{k_2+1}-k_4$ is the maximum net growth rate. The
growth curve for $M_{\rm tot}$ will stay below this maximum growth
rate when substrate $c$ becomes limited in the biomass pocket.
A similar simple {\it a priori} estimate can be derived for $\omega(t)$,
using that due to the fast diffusion nonlinearity $m\leq 1$. We have
\[
\omega(t) < \omega(0) e^{(\frac{k_3}{k_2+1}-k_4)t}
\]
Note that initially $m\ll 1$, therefore as seen above, the increase in
$\omega$ is delayed and starts only after the end of the initial
compression period at $t=t_{compression}$. We can derive a quantitatively
probably better but less rigorous estimate
\[
\omega(t) \approx \omega(0) e^{(\frac{k_3}{k_2+1}-k_4)(t-t{\rm
compression})}
\]
which is valid for some time, until biomass production becomes
severely limited inside the biomass pocket.
According to Figure \ref{sphere_lumped_fig}.a, biomass $M_{\rm
tot}(t)$ initially grows exponentially (up to $t\approx 5$), at
maximum growth rate, indicating that substrate availability is not
limited during this period. During this time the biomass occupied
region does not spread, but the biomass in the core region
solidifies and compresses. After that the biomass core starts
spreading and the biomass growth rate remains sub-linear,
indicating growth limitations due to substrate limitation.
The flow in the system is driven by a constant pressure difference
that is applied at the inflow and outflow boundaries. Therefore,
reductions in the flow rate $Q_{in}(t)$, as plotted in Figure
\ref{sphere_lumped_fig}.b, are entirely due to an increased flow
resistance, i.e. due to an increase in bacterial biomass which
clogs flow pathways in the porous medium. This is expressed in
terms of a decreasing permeability $a(m)$. While the decline in
the flow rate is first gradually, an almost instantaneous drop is
observed at around, $t\approx 15$, when the biomass front reaches
the lateral boundaries (cf. also Fig.
\ref{sphere_illustration_fig}). From then on the biomass region
spans the total width of the domain. While before that water could
bypass the biomass region at increased flow velocities, this is
not possible anymore from then on. Thus all flow after that time
is flow through the biomass region.
In Figure \ref{sphere_lumped_fig}.c the substrate fluxes
$J_{in}(t)$ on inflow and $J_{out}(t)$ on outflow are plotted, as
well as the difference between them. Due to the Dirichlet
condition $c=1$ on inflow, $J_{in}(t)$ is closely aligned with
$Q_{in}(t)$ (plotted in Fig. \ref{sphere_lumped_fig}.b).
Differences between both functions are due to diffusive fluxes on
inflow. As already observed in Figure
\ref{sphere_illustration_fig}, those are close to negligible. Note
that this behavior is quite different from systems in which
transport is only diffusive, i.e. the flow induced convective
transport contribution matters. The inflow flux $J_{in}(t)$
remains positive over the entire simulation interval. Naturally,
$J_{out}(t)$ lies always below $J_{in}(t)$ due to substrate
consumption in the biomass core. As already indicated in Fig.
\ref{sphere_illustration_fig}, substrate is eventually completely
depleted in the downstream region, leading to $J_{out}\approx 0$
eventually ($t\approx 16$ in our simulation). From then on, all
substrate that is supplied is completely consumed by the bacteria.
Before the biomass core spans across the whole domain ($t<15$),
the flow field is able to transport substrate around the biomass
region, leading to $J_{out}>0$ even if some bacteria experience
already substrate limitation in the biomass core, as observed in
Fig. \ref{sphere_illustration_fig}. Note that the declining
consumption rate $J_{in}(t)-J_{out}(t)$ for large $t$ indicates a
slow-down in microbial activity, i.e. decreased growth rates, as
evident also from $M_{\rm tot}(t)$ in
Fig.\ref{sphere_lumped_fig}.a.
Note that other bioclogging models that do not account for spatial
spreading of biomass, e.g. \cite{Chen:2003, Thullner:2004}, are
only able to simulate the process over the initial period, i.e. up
to $t\approx 5$.
\subsection{Non-uniform substrate supply and irregular biomass distribution}
\label{P1_sec}
In a small numerical experiment we study the spatio-temporal
development of biomass and substrate patterns, as well as and
preferential flow direction, under less regular conditions than
in the previous Section \ref{x1_sec}. To this end we perturb the
symmetry in the initial distribution of biomass and replace the
constant (across the flow channel width) inflow concentration by a
variable one. The simulation experiment mimics some characteristic
features of experimental setups that are commonly used in
laboratory experiments \cite{STHM:2006,Thullner:2004}.
The hydrodynamic conditions remain unchanged but the upstream
substrate concentration profile is chosen as the Gaussian
distribution
\[
c(t,0,y)=c_\infty(y)= e^{-40(y-1/2)^2},
\]
mimicking a substrate plume that develops after a point injection
into a well developed flow field upstream.
The biomass is initially distributed in four spherical pockets,
the centers of which are chosen randomly in the subdomain
$[0.3,1]\times [0,0.5]$. The upstream offset was chosen to allow
for the development of a plume-like concentration field, i.e. to
keep upstream boundary effects small. In these pockets the biomass
density is set to $m=0.01$. The spherical regions are allowed to
overlap and are not necessarily fully contained in $\Omega$. Thus,
the initial occupancy $\omega(0)$ and total biomass $M_{\rm
tot}(0)$ can vary between different simulations.
In Figures \ref{plume1_fig} and \ref{plume2_fig} we show the
simulations for two different initial distributions of biomass.
The {\it compression phase} and the {\it expansion phase} as
described above are observed in these simulations as well.
Initially, the substrate concentration appears unaffected by the
bacterial biomass but soon we observe concentration boundary
layers at the rim of the biomass pockets as described in the
previous example. Depending on their locations relative to the
flow and bulk concentration field, the individual biomass pockets
are affected differently. In particular in the simulation of
Figure \ref{plume2_fig} one colony lies outside the plume and,
therefore, develops slower than the other three. Nevertheless, its
contribution to the development of an irregular flow field (as
indicated by the pressure iso-lines) is notable. During the {\it
expansion phase}, initially disjoint neighboring biomass regions
merge. In Figure \ref{plume1_fig} substrate is carried past the
biomass core by convection, although at a largely reduced
(compared to inflow) concentration level.
\begin{figure}[ht] % figure 4
\begin{center}
\includegraphics[width=0.42\textwidth]{fig4a}\; % P1003_02
\includegraphics[width=0.42\textwidth]{fig4b}\\[-5pt] %P1003_12.eps
$t=0.5$ \hspace{4cm} $t=5.5$\\[5pt]
\includegraphics[width=0.42\textwidth]{fig4c}\; %P1003_18
\includegraphics[width=0.42\textwidth]{fig4d}\\[-5pt] %P1003_24
$t=8.5$ \hspace{4cm} $t=11.5$\\[5pt]
\includegraphics[width=0.42\textwidth]{fig4e} \; %P1003_32
\includegraphics[width=0.42\textwidth]{fig4f}\\[-5pt] %P1003bf
$t=15.6$ \hspace{5cm} $\quad$
\end{center}
\caption{Bioclogging, initiated by randomly distributed patches of
biomass: flow is from bottom to top.
Shown are: $c$ (background) [left color map],
$p$ (equidistant white iso-lines),
$m$(colored iso-lines) [right color map],
for selected $t$; bottom right: $\omega(t)$ and $M_{\rm tot}(t)$.
}\label{plume1_fig}
\end{figure}
\begin{figure}[ht] % figure 5
\begin{center}
\includegraphics[width=0.43\textwidth]{fig5a}\; %P1002_02
\includegraphics[width=0.43\textwidth]{fig5b}\\[-5pt] %P1002_12
$t=0.5$ \hspace{4cm} $t=5.5$\\[5pt]
\includegraphics[width=0.43\textwidth]{fig5c}\; %P1002_18
\includegraphics[width=0.43\textwidth]{fig5d} \\[-5pt] %P1002_24
$t=8.5$ \hspace{4cm} $t=11.5$\\[5pt]
\includegraphics[width=0.43\textwidth]{fig5e}\; %P1002_32
\includegraphics[width=0.43\textwidth]{fig5f} \\[-5pt] %P1002bf
$t=15.6$ \hspace{5cm} $\quad$
\end{center}
\caption{Bioclogging, as in Fig. \ref{plume1_fig}, but with
different initial distribution of biomass.}\label{plume2_fig}
\end{figure}
It is naturally expected that different initial biomass
distributions lead to different spatio-temporal patters of biomass
and substrate, as well as different preferential flow paths, as
shown in Figures \ref{plume1_fig} and \ref{plume2_fig}. Whether
this carries over to the lumped parameters $\omega(t)$ and $M_{\rm
tot}(t)$ is investigated by comparing several (in this case 14)
simulations, differing by initial distribution of biomass, cf.
Figure \ref{P1_lumped_fig}. While there is some variation in the
time evolution, in all cases the total biomass in the system
$M_{\rm tot}(t)$ levels off at about the same value, determined by
availability of substrate, as a consequence of flow conditions
(pressure difference between inflow and outflow) and substrate
concentration on inflow. The growth period before the plateau
phase is exponential in most cases, i.e. in good agreement with
the lumped estimate (\ref{MUBeq}), although for few simulations
sub-exponential growth curves evolve from the start, indicating a
local substrate limitation somewhere in the domain already in the
initial phase of the experiment. The spatial distribution of
biomass in the system, expressed in terms of the occupancy
function $\omega(t)$, appears to be more sensitive to initial
distribution of biomass than the amount of bacterial biomass
$M_{\rm tot}(t)$.
\begin{figure}[ht] % figure 6
\includegraphics[width=0.45\textwidth]{fig6a}\; %P1bf_oc
\includegraphics[width=0.45\textwidth]{fig6b} %P1bf_Mtot
\caption{Occupancy $\omega(t)$ [left] and total biomass
$M_{\rm tot}(t)$ [right] for
14 simulations of Section \ref{P1_sec} }\label{P1_lumped_fig}
\end{figure}
\section{Summary and conclusion}
We proposed a bioclogging model, that is based on a volume-filling argument.
Unlike other models of bioclogging, it does not require an artificial
biomass production limiter when the available pore space is filled by
biomass but allows for continued local production of new biomass as
long as growth conditions are favorable. Instead, the bacteria move
into neighboring void regions when the space becomes locally limited.
This is described by a double-degenerate diffusion-mechanism for biomass
that encompasses two nonlinear diffusion effects, namely porous medium
degeneracy and fast diffusion singularity. This spatial model extension
could also be incorporated into other bioclogging models such
as \cite{Chen:2003, Thullner:2004}.
The model for spatio-temporal population dynamics needs to be
coupled with a model for hydrodynamics and for substrate
transport. On the macro-scale, the Darcy equations are the
standard model to describe the flow field in a porous medium and
this is what was used here, taking into account that the hydraulic
conductivity of the soil changes locally when the amount of
biomass in the pore space changes. Developing macroscopic
substrate transport equations for biofilm formation in porous
media is not as straightforward and in fact currently an active
area of research. The difficulty here lies in taking into account
that the pore space is divided into a liquid and a biofilm phase,
which change as the biomass content increases. We use here a
simple one-equation model that incorporates both, biofilm and
aqueous phase in the pores, based on an equilibrium assumption
\cite{GWOQB:sub}. It is to be seen which modifications, if any at
all, need to be made in the bioclogging model once less
restrictive substrate transport models become available.
The bioclogging model is formulated here in a first version for a
rather simple prototype system with only one bacterial species
and one growth limiting substrate, but it generalizes easily to
more involved systems that involve several biomass fractions or
several dissolved substrates (nutrient, pollutants, biocides).
In computer simulations we could study the dependence between
hydrodynamics, substrate transport and population dynamics, and
illustrate some qualitative properties of the mathematical model.
A rigorous solution theory for this class of models, however, is
yet to be established. First results for such double degenerate
diffusion-reaction equations exist in the biofilm modeling
literature, upon which future efforts in this direction will
build.
\begin{thebibliography}{99}
\bibitem{ADQ:2007}
Aspa Y, Debenest G, Quintard M.;
Effect transport properties of
porous biofilms,{\it Eurotherm Seminar No 81: Reactive Heat
Transfer in Porous Media, Ecole des Mines d'Albi}, 8pp, 2007
\bibitem{Baveye:1998}
Baveye P, Vandevivere P, Hoyle BL, DeLeo PC, Sanchez de Lozada D.;
Environmental impact and mechanisms of the biological clogging of
saturated soils and aquifer materials, {\it Crit. Rev. Env. Sci.
\& Techn.}, 123(2):123-191, 1998
\bibitem{Bras:1990}
Bras RL. {\it Hydrology}, 643pp, Addison Wesley, Reading, Mass., 1990
\bibitem{BD:1998}
Bryers J. D., Drummond F.; Local macromolecule diffusion
coefficients in structurally non-uniform bacterial biofilms using
fluorescence recovery after photobleaching (FRAP). {\it
Biotechnol. \& Bioeng.} 60(4):462-473, 1998
\bibitem{Chen:2003}
Chen-Carpentier B. M., Kojouharov H.V.; Numerical simulation of
dual-species biofilm in porous media, {\it Appl. Num. Math},
47:377-389, 2003
\bibitem{DE:2006}
Duvnjak A., Eberl H. J.; Time-discretisation of a degenerate
reaction-diffusion equation arising in biofilm modeling, {\it El.
Trans. Num. Analysis}, 23:15-38, 2006
\bibitem{Eberl:2001}
Eberl H. J., Parker D. F., van Loosdrecht M. C. M.;
A new Deterministic Spatio-Temporal Continuum Model For Biofilm
Development. {\it J. Theor. Medicine}, 3:161-175, 2001
\bibitem{ED:2007} Eberl H. J., Demaret L.;
A finite difference scheme for a doubly
degenerate diffusion-reaction equation arising in microbial
ecology. Electron. J. Diff. Equs. Conf. 15: 77-95, 2007
\bibitem{ED:sub}
Efendiev M. A., Demaret L.; On the structure of attractors for a
class of degenerate reaction-diffusion systems, {\it Adv. Math.
Sci. Appls.}, in press.
\bibitem{EEZ:2002}
Efendiev MA, Eberl H. J., Zelik S. V.; Existence and longtime
behavior of solutions of a nonlinear reaction-diffusion system
arising in the modeling of biofilms. {\it RIMS Kokyuroko},
1258:49-71, 2002
\bibitem{EEZ:2008}
Efendiev M.A., Zelik S., Eberl H. J.; Existence and longtime
behavior of a biofilm model, {\it Comm. Pure Appl. Analysis}, 8(2):509-531, 2009
press
\bibitem{EM:sub}
Efendiev M.A., M\"uller J.; Classification of existence and
non-existence of running fronts in the case of fast diffusion,
{\it Adv. Math. Sci. Appls.}, in press
\bibitem{GWOQB:sub}
Golfier F, Wood B.D., Orgogozo L., Quintard M., Bues M.;
Biofilms in porous media: development of macroscopic transport equations
via volume averaging with closure for local mass equilibrium
conditions, {\it preprint}
\bibitem{HKP:2007}
Ham Y. J., Kim S. B., Park S. J.; Numerical experiments for
bioclogging in porous media, {\it Env. Technology}, 28:1079-1089,
2007
\bibitem{Hornung:1997}
Hornung U.; {\it Homogenization and porous media}, 278pp,
Springer, New York, 1997
\bibitem{KhHE:2009}
Khassehkhan H., Hillen T., Eberl H.J.; A Nonlinear Master Equation for a Degenerate Diffusion Model of Biofilm Growth,
{\it Lecture Notes in Computer Science}, in press
\bibitem{MM:1994}
Morton K. W., Mayers D. F.; {\it Numerical solution of partial
differential equations}, 228pp, Cambridge Univ. Press, Cambridge,
1994
\bibitem{PH:2002}
Painter K.J., Hillen T.; Volume-filling and quorum-sensing in
models for chemosensitive movements, {\it Can. Appl. Math. Quart.}
10(4):501-543, 2002
\bibitem{Rittmann:1982}
Rittmann B. E.; The effect of shear stress on biofilm loss rate,
{\it Biotech. \& Bioeng.} 24:501-506, 1982
\bibitem{Schaefer:1998}
Sch\"afer D., Sch\"afer W., Kinzelbach W.; Simulation of reactive
processes related to biodegradation in aquifers: 1. Structure of
the three-dimensional reactive transport model. {\it J.Cont.
Hydr.} 31:167-186, 1998
\bibitem{SM:2001}
Seki K., Miyazaki T.; A mathematical model for biological clogging
of uniform porous media, {\it Wat. Res. Res.} 37(12):2995-2999,
2001
\bibitem{STHM:2006}
Seki K., Thullner M., Hanada J., Miyazaki T.; Moderate bioclogging
leading to preferential flow paths in biobarriers, {\it Ground
Water Monitoring \& Remediation}, 26(3):68-76, 2006
\bibitem{S:2003}
Stewart P.S.; Diffusion in biofilms, {\it J. Bacteriol.}
185:1485-1491, 2003
\bibitem{Thullner:2002}
Thullner M., Mauclaire L., Schroth M. H., Kinzelbach W., Zeyer J.;
Interaction between flow and spatial distribution of microbial
growth in a two-dimensional flow field in saturated porous media,
{\it J. Cont. Hydr.}, 58:169-189, 2002
\bibitem{Thullner:2004}
Thullner M, Schroth M. H., Zeyer J., Kinzelbach W.; Modeling of a
microbial growth experiment with bioclogging in a two-dimensional
saturated porous media flow field. {\it J. Cont. Hydr.} 70:37-62,
2004
\bibitem{TGR:2006}
Wanner O., Eberl H., Morgenroth E., Noguera D. R., Picioreanu C.,
Rittmann B., van Loosdrecht M. C. M.; {\it Mathematical modeling
of biofilms}, 178pp, IWA Puplishing, London, 2006
\end{thebibliography}
\end{document}