\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}
\AtBeginDocument{{\noindent\small Sixth Mississippi State
Conference on Differential Equations and Computational
Simulations, {\em Electronic Journal of Differential Equations},
Conference 15 (2007), pp. 77--95.\newline ISSN: 1072-6691. URL:
http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu (login: ftp)}
\thanks{\copyright 2007 Texas State University - San Marcos.}
\vspace{9mm}}
\begin{document} \setcounter{page}{77}
\title[\hfilneg EJDE-2006/Conf/15\hfil A finite difference scheme]
{A finite difference scheme for a degenerated diffusion equation
arising in microbial ecology}
\author[H. J. Eberl, L. Demaret \hfil EJDE/Conf/15 \hfilneg]
{Hermann J. Eberl, Laurent Demaret} % in alphabetical order
\address{Hermann J. Eberl \newline
Dept. Mathematics and Statistics \\
University of Guelph, Guelph, On, Canada}
\email{heberl@uoguelph.ca}
\address{Laurent Demaret \newline
Inst. Biomathematics and Biometry \\
GSF - National Research Centre for Environment and Health \\
Neuherberg, Germany}
\email{laurent.demaret@gsf.de}
\thanks{Published February 28, 2007.}
\subjclass[2000]{35K65, 65M06, 92C17}
\keywords{Finite differences; nonlinear diffusion;
non-local representation; \hfill\break\indent
non-standard discretisation; numerical simulation; biofilm}
\begin{abstract}
A finite difference scheme is presented for a density-dependent
diffusion equation that arises in the mathematical modelling
of bacterial biofilms. The peculiarity of the underlying model
is that it shows degeneracy as the dependent variable vanishes,
as well as a singularity as the dependent variable approaches
its a priori known upper bound. The first property leads to
a finite speed of interface propagation if the initial data have
compact support, while the second one introduces counter-acting
super diffusion. This squeezing property of this model leads to
steep gradients at the interface.
Moving interface problems of this kind are known to be problematic
for classical numerical methods and introduce non-physical and
non-mathematical solutions. The proposed method is developed to
address this observation. The central idea is a non-local (in time)
representation of the diffusion operator.
It can be shown that the proposed method is free of oscillations
at the interface, that the discrete interface satisfies a discrete
version of the continuous interface condition and that the effect
of interface smearing is quantitatively small.
\end{abstract}
\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Corollary}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{remark}[theorem]{Remark}
\section{Introduction}
\subsection{Biological background}
The mathematical model in the focus of this study describes the
formation and growth of bacterial biofilms. These are communities
of migroorganisms that live embedded in a self-produced layer of
extracellular polymeric substances ({\it EPS}). The EPS offers
protection to the bacteria in the biofilm, which develop
behavioral patterns that differ from those of planktonic bacteria
\cite{Costerton}, \cite{Flemming}. As a consequence, harmful
biofilm bacteria are more difficult to remove, mechanically or by
antimicrobial agents, than suspended bacteria \cite{Stewart},
\cite{Watnick}. Biofilms occur naturally on surfaces and
interfaces in aquatic systems wherever enough nutrients are
available to sustain microbial growth \cite{Flemming}. Harmful
occurrences of biofilms include bacterial infections, dental
plaque and associated diseases, biocorrosion of drinking water
pipes or industrial facilities, and food spoilage. On the other
hand, the sorption properties and enhanced mechanical stability of
biofilms are beneficially used by environmental engineers, {\it
e.g.} in wastewater treatment, soil remediation, and groundwater
protection \cite{Wanner:2006}.
While the term biofilm indicates a homogeneous film-like layer,
biofilms on the meso-scale ($10 \mu m \sim 1mm $, the actual
biofilm scale) in reality can develop in structurally highly
complicated architectures. The morphology of a biofilm is
triggered by environmental conditions, such as the availability of
required substrates like nutrients and, in the case of aerobic
biofilms, oxygen \cite{vanL:1995}, \cite{vanL:1997},
\cite{Wanner:2006}, \cite{Xu:1998}. Well-known are the so-called
mushroom or pillar shaped biofilm morphologies. These are clusters
of biofilm colonies, trenched by water-filled pores and channels,
cf. the confocal laser scanning micrograph in Figure \ref{CLSM}
and the schematic in Figure \ref{schematic}. Biofilm morphologies
of this type are observed in nutrient limited regimes. Initially,
at low population density, the bacterial population develops
unhindered. As the population grows, demand of food increases and
inner-species competition for food sets in when nutrients become
depleted. Large colonies have an advantage allowing them to grow
faster towards the nutrient source and thus out-compete smaller
colonies that stay behind in their development. As these biofilms
grow under nutrient limitation they typically develop slowly. On
the other hand, biofilms developing in a nutrient rich, un-limited
environment tend to grow quickly and form homogeneous, compact,
flat yet thick structures. In order to capture these features,
spatially resolving techniques are required for biofilm research.
In addition to modern non-invasive microscopy methods in
experimental studies, mathematical models and computer simulations
for theoretical studies become more and more accepted
\cite{Wanner:2006}, \cite{Wilderer:2002}.
\begin{figure}[ht]
\begin{center}
%\epsfxsize=10truecm\epsffile{
\includegraphics[width=0.8\textwidth]{figures/frame3gs}
\end{center}
\caption{Microscopy snapshot of a very young spatially heterogeneous
biofilm of Listeria monocytogenes and Pseudomonas putida developing in
a cluster-and-channel morphology. The CLSM figure was made available
by Heidi Schraft, Lakehead University, Thunder Bay, On, Canada.
The original colors were modified for better greyscale printing.}
\label{CLSM}
\end{figure}
\subsection{Mathematical models and numerical simulation}
Traditional biofilm models were formulated as one-dimensional
free-boundary value problems, based on the assumption that newly
produced biomass is converted into new biofilm volume. This leads
to an increasing thickness of the biofilm and the speed of
propagation of the biofilm/liquid interface normal to the
substratum can be calculated from the production terms by
integration over the biofilm thickness \cite{Wanner:2006}. As this
concept is inherently one-dimensional, it cannot predict spatially
heterogeneous biofilm morphologies. To fill this gap, a variety of
models for spatially heterogeneous biofilms has been proposed in
recent years, ranging from stochastic individual based models to
cellular automata models to deterministic continuum models
\cite{Pici:2003}. The underlying mathematical principles of these
models are quite different and, therefore, the mathematical and
computational challenges vary from model to model. Nevertheless,
they all show the same qualitative behavior in predicting the
development of biofilm morphologies in response to substrate
limitation, cf. \cite{Pici:2003} and the references therein for a
selection and overview.
The model for biofilm formation that we consider was introduced
originally in \cite{Eberl:2001} and studied in a small series of
papers analytically, numerically and in applications
\cite{Duvnjak:2004,Eberl:2003,Eberl:2004,
Eberl:2006,Efendiev:2005,Khassehkhan:2006}. It is a quasi-linear
diffusion equation for the dependent variable biomass density $u$,
with two interacting causes of degeneracy. One is degeneracy as in
the porous medium equation, i.e. the density dependent diffusion
coefficient $D(u)$ vanishes as the dependent variable vanishes,
$D(0)=0$; the other one is a singularity in the density-dependent
diffusion coefficient as the dependent variable approaches its a
priori known maximum density, $u\to 1$. In the current study we
propose a numerical discretisation scheme that is able to handle
these two simultaneously occurring effects and to render important
qualitative features of the continuous model. The actual biofilm
is the region in which the biomass density is positive. The
transition from the surrounding aqueous phase to the biofilm is
sharp, that is the biomass gradients are steep at the
biofilm/liquid interface. Since the biofilm grows, this interface
is not stationary but its location changes over time. The speed
and direction of the interface propagation is not a priori known
but an implicit consequence of the solution of the equation;
furthermore, as two or more biofilm colonies merge into a bigger
one the interfaces collide and become dissolved. Interface
problems of this type are known to pose difficulties for many
numerical methods, often leading to either un-physical and
un-mathematical oscillations in the solution or smearing of the
sharp interface. The underlying reason is that discretisation
schemes often assume more regularity than the initial-boundary
value problem offers. This effect can be observed already in very
simple examples that possess classical smooth solutions, e.g. the
heat equation with a discontinuous interface in the initial data
\cite{Hundsdorfer}.
A variety of numerical methods were used in previous studies of
the biofilm model. In \cite{Eberl:2001,Eberl:2003,Eberl:2004},
based on a time-scale argument, the faster (in
terms of characteristic time scales) substrate equations were
solved implicitly while the slower biomass equation is solved with
explicit schemes, e.g. the adaptive 4th/5th order
Runge-Kutta-Fehlberg method, cf\cite{Epperson}. The maximum
feasible time-step in this approach depends on the solution of the
equation and becomes very small as the biofilm grows and maximum
biomass density $u\to 1$ is approached somewhere. A very different
idea was followed in \cite{Duvnjak:2004}, where the degenerated
equation was transformed into new dependent variables such that
the spatial operator becomes the Laplacian and all non-linear
effects appear in the time-derivative. A fully implicit method for
this equation was devised that has only a very mild time-step
constraint. This method handles the diffusion singularity effect
very well, but it was found that in some cases it emphasises
interface oscillations, i.e. is not able to correctly deal with
the porous medium degeneracy. In \cite{Khassehkhan:2006} the model
equation was re-formulated in a weak form in the moving frame and
subsequently discretised by a moving grid finite element method
that explicitly tracks the biofilm/liquid interface, following an
idea that was originally developed in \cite{Baines} for degenerate
parabolic equations. This method describes the porous medium
degeneracy very well. It is fully adaptive in time and space (with
explicit time integration) and builds its biofilm grid adaptively.
However, re-meshing and grid-refinement is required as new biomass
is produced and the biofilm structure grows, in order to avoid
finite elements becoming too large. This requires interpolation
and becomes inefficient in higher dimensions.
In this paper we will present a numerical method for the biofilm
model that is based on a non-local (in time) representation of
the nonlinear density dependent diffusive biomass flux. This leads
to a semi-linear method that inherits the mild time-step
constraint of the implicit method in \cite{Duvnjak:2004} and the
advantages of the explicit method. In every time step the solution
of only one linear algebraic system is required. In particular the
method deals well with the moving interface, preserves the
positivity of the solutions and is free of oscillations.
\section{The continuous model}\label{sec_modeldesc}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.85\textwidth]{figures/schematic}
\end{center}
\caption{The domain $\Omega \subset \mathbb{R}^2$ with liquid region
$\Omega_1(t)=\{x \in \Omega: u(t,x)=0\}$ and biofilm region
$\Omega_2(t)=\{ x \in \Omega: u(t,x)>0\}$. As $u$ changes with time,
the interface between $\Omega_1$ and $\Omega_2$ starts moving.}
\label{schematic}
\end{figure}
We consider the biofilm as a continuum and study a
density-dependent diffusion-reaction equation that describes the
development of a spatially structured bacterial biofilm. The
model is formulated as an evolution equation over the domain
$\Omega \subset R^d$, $d=1,2,3$, in terms of the biomass density
$u$, normalised with respect to the maximum biomass density. Under
this {\it ansatz}, bacteria and EPS are considered together in one
fraction. This is the classical approach in biofilm modelling
\cite{Wanner:2006}.
The region $\Omega_1(t):=\{x \in \Omega: u(t,x)=0\}$ describes the
surrounding aquatic environment (bulk liquid, channels and pores
of a biofilm) without biomass, while $\Omega_2(t)=\{ x \in \Omega:
u(t,x)>0\}$ is the actual biofilm (cf. Figure \ref{schematic}).
The model under consideration was introduced in \cite{Eberl:2001}.
In the notation used here it reads
\begin{equation}\label{model}
u_t = \nabla_x (D(u) \cdot \nabla_x u) + ku
\end{equation}
where
\begin{equation}\label{diffcoeff}
D(u) = \delta {u^b \over (1-u)^a},\quad a,b \geq 1\gg \delta >0
\end{equation}
together with Dirichlet, Neumann, or mixed boundary conditions. In
(\ref{model}) the quantity $k$ describes the production rate. For
our purpose it is sufficient to assume that it is bounded. If $k$
is a positive constant (\ref{model}) corresponds to a biofilm
system in which nutrients are not limited. In this case we expect
a homogeneous biofilm morphology for large $t$. In the case of
nutrient limitations, $k(t,x)$ is a non-constant function. We
shall present simulations for both scenarios below.
Some analytical results about this model can be found in
\cite{Duvnjak:2004}, \cite{Efendiev:2005}; biofilm applications
of this type of equation and its generalisations were discussed in
\cite{Eberl:2001}, \cite{Eberl:2003}, \cite{Eberl:2004}.
The long term behavior of the initial-boundary value problem
associated with (\ref{model}) was studied in
\cite{Efendiev:2005}. It was shown that for initial data $u_0$ with
\[
u_0\in L^\infty (\Omega), \quad F(u_0) \in H_0 ^1
(\Omega),\quad u_0\ge 0,\quad\Vert u_0 \Vert _{L^\infty}
(\Omega)<1
\]
where
\[
F(u)=\int_0 ^u \frac{v^\alpha}{(1-v)^\beta} dv,\quad
0\le u<1
\]
there exists a unique solution $u$ of (\ref{model}) in the class
of functions
\begin{gather*}
u\in L^\infty (R_+ \times \Omega)\cap \, C([0,\infty),L^2
(\Omega)))
\\
F(u)\in L^\infty(R_+,H^1(\Omega)\cap C([0,\infty),L^2(\Omega))
\\
0\le u(t,x)\le 1 \quad \Vert u \Vert_{L^\infty (R_+ \times \Omega)}<1
\end{gather*}
If the boundary conditions are homogeneous Neumann conditions
everywhere, then the solution $u$ reaches its maximum density $1$
almost everywhere in finite time. Physically this corresponds to
the situation where the biofilm growths unlimited in a closed
vessel from where no biomass can escape. If Dirichlet conditions
$u=0$ are specified somewhere on the boundary of $\Omega$, then the
solution will remain below the maximum density for all $t>0$
almost everywhere. Further results include the existence of a global attractor \cite{Efendiev:2005} and the construction of a
Lyapunov functional \cite{Duvnjak:2004}.
An essential property of the solutions of (\ref{model}) is that initial data with compact support imply solutions with compact support and that the solution is indeed bounded by the maximum biomass density. The first effect is due to the porous medium degeneracy $u^b$ in (\ref{diffcoeff}), while the latter is due to the singularity $(1-u)^{-a}$ in (\ref{diffcoeff}).
It is the production term $ku$ together with this second effect that drives the spatial spreading of $\Omega_2(t)$. This is counteracted by the degeneracy as $u=0$ at the interface. As a consequence, $u$ squeezes in $\Omega_2(t)$ and, in the absence of loss terms, approaches its maximum value 1.
Hence, the interaction of both non-linear diffusion effects with the growth term in (\ref{model}) is needed to describe spatial biomass spreading. It leads to very steep gradients of $u$ at the moving interface between $\Omega_1(t)$ and $\Omega_2(t)$. If $\Omega_2(t)$ is initially not connected but consists of several ``sub-domains'' (each of which connected), these will eventually join if $k$ is positive and large enough, due to the expansion of each individual sub-domain. In biofilm applications this describes the merging of two initially isolated bacterial colonies. Thus, if two $\Omega_1$-$\Omega_2$ interfaces collide, they are dissolved and the sub-domains become connected.
In \cite{Duvnjak:2004} a semi-discretisation in time was discussed
that aimed at overcoming the obvious numerical difficulties
arising when $u \to 1$. After transformation of the dependent
variable $u \mapsto v:=F(u) $, the nonlinearities appear in the
time-derivative, while the spatial spreading effect is described
by the linear diffusion operator. Backward Euler discretisation
leads to an elliptic boundary value problem in every time-step.
Convergence results showed that the time-step restriction of this
model is not critical, {\it i.e.} the maximum time-step size
permitted for convergence is larger than the characteristic
time-scale of the biofilm formation process, i.e. for given
spatial discretisation it behaves like $\mathcal{O}(1/k)$.
However, it was found that the numerical discretisation of the
spatial terms with standard Finite Element or Finite Difference
methods introduces oscillations with negative values at the
biofilm/liquid interface. Local grid refinement can postpone and
dampen this effect but not avoid it. In the current work we aim at
a generally applicable numerical scheme for (\ref{model}) and
model extensions that is free of this non-physical effect.
\section{Numerical Method}
\subsection{Finite Difference Scheme For The Continuous Problem}
We introduce a finite difference scheme for (\ref{model}) on a
regular grid based on a first order difference approximation for
the time derivative and the usual second order spatial difference
for self-adjoint diffusion operators, {\it cf} \cite{Morton:1994}.
The key idea of the method is to represent the nonlinear diffusive
flux $D(u)\cdot \nabla_x u$ in (\ref{model}) and the reaction
term $k \cdot u$ non-local in time: The diffusion coefficient is
evaluated in time level $t$ while the gradient is evaluated at the
new time level $t+\Delta t$, {\it i.e.} one has
\[
D(u)\cdot \nabla_x u \approx D\bigl(u(t,\cdot)\bigr)\cdot
\nabla_x u(t+\Delta t, \cdot)
\]
and similarly for the reaction term
\[
k\cdot u \approx k(t,\cdot) \cdot u(t +\Delta t,\cdot)
\]
Finite difference methods for (ordinary or partial) differential
equations utilising such a non-local representation of non-linear
terms are called {\it non-standard finite difference schemes}
\cite{Lubuma:2001}. Despite the low order convergence of the
chosen derivative approximations, such a discretisation is known
to be optimal among linear finite difference schemes for certain
problems under the aspect of positivity and monotonicity (cf.
\cite{Hundsdorfer} for details). In the context of non-standard
finite difference schemes this strategy of low order
discretisation of derivatives is routinely used for the
construction of discretisation schemes and referred to as {\it
Rule 1} in \cite{Mickens:2000}.
In the sequel, we denote by $u_j^n$ the numerical approximation of
the exact solution $u(t^n, x_j)$ in the grid points $x_j \in R^d$.
Here $j \in J$ can be a multi-index or a grid ordering in the
multi-dimensional case $d>1$; then $J$ is an appropriate index
set. The vector $U^n=(u_j^n)_{j\in J} \in R^{| J|}$ is the
(ordered) vector, the coefficients of which are the numerical
approximations of the solution at time level $n$. Employing a
regular grid with spatial step size $\Delta x$ one obtains for the
discretised diffusion operator
\begin{equation}\label{diffdisc}
\nabla_x \left(D(u) \nabla_x u\right) \dot= {1 \over 2 \Delta
x^2}\sum_{k\in N_j} (D(u_k^n) +D(u_j^n))(u_k^{n+1} - u_j^{n+1})
\end{equation}
where $N_j \subset J$ is the index set pointing at the direct
neighbors of $x_j$ on the grid. For $d=1$ it has two elements, for
$d=2$ there are four and for $d=3$ there are six. The non-standard
finite difference scheme that is derived under the above
assumptions reads
\begin{equation}\label{imnsfd}
U^{n+1} - U^n = \Delta t \mathcal{D}(U^n) \cdot U^{n+1} + \Delta
t \mathcal{K} U^{n+1}
\end{equation}
where $\mathcal{K}=diag(k_j^n)$ is the diagonal matrix the
coefficients of which are the net growth rates $k(t^n,x_j)$. The
matrix $\mathcal{D}(U^n)$ is the banded matrix obtained by second
order standard finite difference discretisation of the diffusion
operator evaluated at the previous time-step $t^n$, according to
(\ref{diffdisc}).
System (\ref{imnsfd}) can be re-written in matrix-vector form as
\begin{equation}\label{nsfd}
U^{n+1} = \left[\mathcal{I} - \Delta t \mathcal{D}(U^n) - \Delta
t \mathcal{K}\right]^{-1} U^n
\end{equation}
if the inverse matrix exists. This is the case if the time-step
$\Delta t$ is chosen such that $1/\Delta t$ is not an eigenvalue
of the matrix $\mathcal{D}(U^n) + \mathcal{K}$. Then, the
numerical solution $U^{n+1}$ is unique. Hence, a sufficient
condition for the existence of a unique solutions is the
restriction of the time-step size $t^n \to t^n+\Delta t =:
t^{n+1}$ by
\begin{equation}\label{dtrestrict}
\Delta t \cdot \max_i k_i^n <1
\end{equation}
This ensures that $I-\Delta t \mathcal{K}$ is diagonal with
strictly positive entries. Then in (\ref{nsfd}) the matrix
$\mathcal{ I} - \Delta t \mathcal{D}(U^n) - \Delta t \mathcal{K}$
is diagonally dominant. Inequality (\ref{dtrestrict}) poses a
time-step restriction for the numerical method. From the
application point of view this is a weak restriction, as $1/k$ is
the characteristic time scale of biomass production. Accordingly,
time-steps larger than $1/k$ would be at the expense of inaccurate
description of the growth process. The same time stepping
condition (\ref{dtrestrict}) was obtained for the transformation
method in \cite{Duvnjak:2004}. If $k$ is not a constant,
(\ref{dtrestrict}) leads to a variable time-stepping strategy.
For a given $U^n$ with $0\leq U^n <1$ (where all vector
inequalities are understood coefficient-wise), $U^{n+1}$
according to (\ref{nsfd}) is a continuous function of $\Delta t$.
Hence, for $\Delta t$ small enough $U^n <1$ implies $U^{n+1}<1$. A
more quantitative statement will be presented in the next section.
The non-standard finite difference scheme (\ref{nsfd}) is a system
of explicit non-linear difference equations, {\it i.e.} it can be
understood as a discrete dynamic system. In the numerical
simulation the solution of one coupled linear system is required
in every time-step. The matrix is sparse, banded and diagonally
dominant. In two- and three-dimensional simulations, the system
can be solved iteratively, e.g. by the stabilised bi-conjugate
gradient method, which is described in \cite{Saad}; due to
diagonal dominance, pre-conditioning is not required for
convergence. In the one-dimensional case the system is tridiagonal
and the Thomas algorithm (e.g. \cite{Epperson}) can be used. These
are the methods that we employ in the simulations below.
In the next section we will show that the numerical solution
according to (\ref{nsfd}) inherits the following important
qualitative properties of the continuous model:
\begin{itemize}
\item positivity
\item boundedness $U^n < 1$ of the solution
\item finite speed of propagation of the discrete liquid/biofilm interface
\item the solution of (\ref{nsfd}) satisfies a discrete interface
condition that corresponds to the interface condition for the
continuous model (\ref{model}).
\item the numerical interface is sharp, i.e. only weak interface
smearing takes place
\item the solution is monotonous at the interface and merging of two
colonies is well-posed
\end{itemize}
\subsection{Properties Of The Scheme}
\begin{lemma} [Positivity] \label{posLemma}
Let (\ref{dtrestrict}) be always satisfied. If $0\leq U^0 <1$
then $U^{n}\geq 0$
\end{lemma}
\begin{proof}
We use an argument from the theory of $M$-matrices (cf.
\cite{Hackbusch}) and mathematical induction: If $\Delta t$ is
chosen according to (\ref{dtrestrict}) then the matrix in
(\ref{nsfd}) is diagonally dominant with positive diagonal
elements and non-positive off-diagonals. Hence, its inverse is
positive and, thus, positivity of $U^n$ implies positivity of
$U^{n+1}$ according to iteration (\ref{nsfd}).
\end{proof}
\begin{remark} \label{rmk3.2} \rm
An alternate proof can be derived based on the following
observation:
In the time step $t^n \to t^{n+1}$ method (\ref{imnsfd}) can be
understood as an implicit Euler method for the linear system of
ordinary differential equations
\begin{equation}\label{semidiscrete}
\dot v = \left(\mathcal{D}(U^n)+\mathcal{K}\right) v, \quad
v(t^n)=U^n, \quad t\in [t^n,t^{n+1}]
\end{equation}
where $v$ is a time-continuous vector valued function. Thus,
method (\ref{imnsfd}) is a piecewise linear problem and in every
time-step linear theory applies, which can be found {\it e.g.} in
\cite{Hundsdorfer}. Thus in every time-step positivity is ensured
if the matrix in (\ref{semidiscrete}) has only non-negative
off-diagonal entries and if the diagonal is bounded from below,
{\it i.e.} if there is an $\alpha_n>0$ such that
$\left(\mathcal{D}(U^n)+\mathcal{K}\right)_{jj} > -\alpha_n$. That
this is the case follows from a recursive argument: It is easily
verified that these conditions are satisfied for $n=0$. If they
hold for one $n$ then it follows from a simple calculation that it
also holds for $n+1$, due to the assumption that $\Delta t$ is
small enough to warrant $U^n<1$ and the definition of $\mathcal{
D}(U)$. Note that both methods of proof for Lemma \ref{posLemma}
use the same property of the system matrix in (\ref{nsfd}), and
thus are essentially equivalent.
\end{remark}
The M-matrix property is also used to derive the following
sufficient (but not necessary) time-step condition for
boundedness of the solution.
\begin{lemma} [Boundedness] \label{bounded}
The choice of a time-step $\Delta t$ such that
(\ref{dtrestrict}) holds and
\[
\Delta t < \min_i {1 -u_i^n \over \sum_j[\mathcal{D}(U^n) +
\mathcal{ K}]_{ij}}
\]
guarantees that the numerical solution obeys the upper bound, i.e.
$U^{n+1} \leq 1$ if $U^n <1$.
\end{lemma}
\begin{proof} We introduce $W^n:=e-U^n$ where $e:=(1,\dots ,1)^T$.
Hence, $U^n\leq 1$ is equivalent to $W^n\geq 0$. Then (\ref{nsfd})
gives
\begin{equation}\label{wnsfd}
\left[\mathcal{I} - \Delta t \mathcal{D}(U^n) - \Delta t \mathcal{
K}\right](e-W^{n+1})=(e-W^n)
\end{equation}
and, therefore,
\begin{equation}
\left[\mathcal{I} - \Delta t \mathcal{D}(U^n) - \Delta t \mathcal{
K}\right] W^{n+1} =W^n- \Delta t \left[\mathcal{D}(U^n) +
\mathcal{ K}\right]e
\end{equation}
If $W^n>0$, positivity of $W^{n+1}$ is ensured if $\Delta t$ is
chosen such that the right hand side of (\ref{wnsfd}) is
non-negative. This is certainly the case if
\begin{equation}\label{dtrestrict2}
\Delta t < \min_i {1 -u_i^n \over \sum_j[\mathcal{D}(U^n) +
\mathcal{ K}]_{ij}}
\end{equation}
\end{proof}
It is easy to verify that criterion (\ref{dtrestrict2}) can be
more strict than (\ref{dtrestrict}). It can be relaxed by the
observation that it is sufficient but not necessary. Indeed, in
the numerical simulations conducted below, this time-step
criterion was implemented as an emergency strategy to guarantee
the upper bound if other time-step control mechanisms fail, but
never activated.
Note that if in one coefficient $U^n$ reaches 1, (\ref{dtrestrict})
implies the termination of the simulation. This is in accordance
with the continuous problem (\ref{model}), since it was shown in
\cite{Efendiev:2005} that the solution of (\ref{model}) can reach
1 in finite time in dependence of boundary conditions and reaction
rates.
To formulate the next result we introduce the following concept.
\noindent{\bf Definition.} We denote by $\Omega_1^n$ a discrete
grid analogy of $\Omega_1(t^n)$ as
\[
\Omega_1^n = \{ x_j : u_j^n=0 \}
\]
An interior point of $\Omega_1^n$ is a point in $\Omega_1^n$ whose
direct neighbors are in $\Omega_1^n$ as well, i.e. a point $x_j
\in \Omega_1^n$ such that $x_i \in \Omega_1^n$ for all $i \in
N_j$. Similar definitions can be made for the biofilm region
$\Omega_2$
Note that this definition is based on the numerical solution
$u_j^n$. Thus it does not imply that $\Omega_1^n$ is contained in
$\Omega_1(t^n)$ or {\it vice versa}. For now, the discrete
interface at time level $t^n$ is defined by grid points of
$\Omega_1^n$ with a direct neighbor in $\Omega_2^n$. The following
result states that the discrete interface travels at most with
finite speed. In particular it states that in one time-step it
moves across at most one grid cell.
\begin{lemma} [Finite speed of propagation] \label{finspeedLemma}
Let \eqref{dtrestrict} hold. If $x_j$ is an interior point of
$\Omega_1^n$ then $x_j \in \Omega_1^{n+1}$.
\end{lemma}
\begin{proof}
If $x_j\in \Omega_1^n$ is interior then the corresponding $j$th
row of the matrix $\mathcal{D}(U^n)$ contains only zeros. Thus one
has from (\ref{nsfd})
\[
(1-\Delta t k_j^{n}) u_j^{n+1} = u_j^n= 0
\]
with (\ref{dtrestrict}) it follows that $u_j^{n+1}=0$ and thus
$x_j \in \Omega_1^{n+1}$.
\end{proof}
In the one-dimensional case $d=1$ an explicit expression can be
derived for the speed of interface propagation in the continuous
model (\ref{model}). An interface between $\Omega_1(t)$ and
$\Omega_2(t)$ can be parameterised locally by a curve $x^\ast(t)$.
The speed of the interface is $\dot x^\ast(t)$. Due to continuity
and $u(x^{\ast}(t),t)=0$ we obtain from (\ref{model})
\begin{equation}\label{intcond}
\dot x^\ast(t) = -[u_t/u_x]_{x^\ast(t)}
\end{equation}
where the right hand side is to be understood in the sense of a limit of
the quotient of the one-sided derivatives as one approaches the
interface. The theory developed in \cite{Efendiev:2005} implies
that this is finite. We show that the numerical solution satisfies
a corresponding discrete interface condition if the derivatives in
(\ref{intcond}) are replaced by the usual difference quotients.
Without loss of generality, we restrict ourselves to the case of
an expanding biofilm (i.e. $\Omega_2$ increases) with an interface
with positive speed. We define the location $x_{j^{\ast,n}}$ of a
discrete interface at time $t^n$ as follows: $x_{j^{\ast,n}-1} \in
\Omega_2^n$, $x_{j^{\ast,n}} \in \Omega_2^n$, $x_{j^{\ast,n}+1}
\in \Omega_1^n$ and $x_{j^{\ast,n}+2}$ an interior point of
$\Omega_1^n$. That is, $j^{\ast,n}$ is the index of the last
biofilm point before the liquid/biofilm interface at time $t^n$.
\begin{lemma} [Moving interface condition]
The following discrete interface condition holds
\begin{equation}\label{dintcond}
{ x_{j^{\ast,n+1}} - x_{j^{\ast,n}} \over \Delta x}
= - { u_{j^{\ast,n+1}}^{n+1}
- u_{j^{\ast,n+1}}^{n} \over u_{j^{\ast,n+1}+1}^{n+1}
- u_{j^{\ast,n+1}}^{n+1}}
\end{equation}
\end{lemma}
\begin{proof} For brevity of the notation let $i:=j^{\ast,n}$.
We obtain from the $(i+1)$th coefficient of (\ref{nsfd}) and the
hypothesis on the position of the interface
\begin{equation}\label{interface}
-{\Delta t \over 2} D(u_i^n) u_i^{n+1}
+\Big(1 +{\Delta t \over 2} D(u_i^n) - \Delta t
k_{i+1}^{n}\Big) u_{i+1}^{n+1}=0
\end{equation}
and thus $u_{i+1}^{n+1} >0$. With the previous Lemma we have
that $j^{\ast,n+1}=i+1=j^{\ast,n}+1$. Thus $u_{j^{\ast,n+1}}^{n}=0$.
Furthermore we have by hypothesis $u_{j^{\ast,n+1}+1}^{n+1}=0$.
The assertion follows by substituting these observations
into (\ref{dintcond}).
\end{proof}
Note that (\ref{dintcond}) is a discrete version of (\ref{intcond})
after multiplying with $\Delta x/\Delta t$.
The proof showed that it could have been formulated in a stronger version
that essentially states that the speed of propagation of the discrete
interface is $\Delta x/ \Delta t$. This is a consequence of the
smearing around the interface that is introduced by the
spatial discretisation of the diffusion operator.
However, from (\ref{interface}) it follows
\[
{\Delta t \over 2} D(u_i^n)\left[u_{i+1}^{n+1} - u_i^{n+1}\right]
+ \left[1 - \Delta t k_{i+1}^{n}\right] u_{i+1}^{n+1} =0
\]
and thus with (\ref{dtrestrict}) and the positivity of (\ref{nsfd})
we obtain $0< u_{i+1}^{n+1} < u_i^{n+1}$. Moreover, after re-arranging
terms we have
\[
u_{i+1}^{n+1} = \mathcal{O}\left((u_i^n)^b u_i^{n+1}\right)
\]
with the definition (\ref{diffcoeff}), since $u_i^n$ small at the
interface implies $(1-u_i^n) \gg 0$. In biofilm applications one
has $b$ in the range 2 through 6. Thus we may expect the interface
smearing effect to be quantitatively small. This will be
demonstrated in the numerical simulations in the next section.
The last theoretical result in this section shows that the
numerical solution at time $t^{n+1}$ in a grid point in
$\Omega_1^n$ is bounded from above by the values of the solution
in neighboring grid points. This is of particular interest in the
case where the grid point changes its membership from $\Omega_1^n$
to $\Omega_2^{n+1}$ as $t^n \to t^{n+1}$, i.e. is passed by a
biofilm/water interface.
\begin{lemma} [Merging of two colonies] \label{mergeLemma}
Let $x_i \in \Omega_1^n$ and assume that (\ref{dtrestrict}) holds.
Then $u_i^{n+1}$ is bounded by the values of $u$ in the neighboring nodes,
i.e.
\[
0\leq u_i^{n+1} \leq \max_{j\in N_i}\{u_j^{n+1} \}
\]
\end{lemma}
\begin{proof} We have $x_i^n \in \Omega_1^n \Longrightarrow u_i^n =0
\Longrightarrow D(u_i^n)=0$.
Thus (\ref{imnsfd}) and (\ref{diffdisc}) imply
\[
u_i^{n+1} = {\Delta t \over \Delta x^2}\sum_{j\in N_i}
D(u_j^n)(u_j^{n+1} - u_i^{n+1}) + \Delta t k_i^{n} u_i^{n+1}
\]
and, hence,
\begin{align*}
u_i^{n+1} &={{\Delta t \over \Delta x^2}\sum_{j\in N_i}
D(u_j^n)u_j^{n+1} \over 1+ {\Delta t \over \Delta x^2} \sum_{j\in
N_i} D(u_j^n) - \Delta t k_i^{n} } \\
&\leq {{\Delta t \over \Delta x^2}\sum_{j\in N_i} D(u_j^n) \cdot
\max_{j\in N_i} u_j^{n+1} \over 1+ {\Delta t \over \Delta x^2}
\sum_{j\in N_i} D(u_j^n) - \Delta t k_i^{n} } \\
&\leq \max_{j\in N_i} u_j^{n+1}
\end{align*}
The last inequality follows from (\ref{dtrestrict}) with $1-\Delta t k_j^n>0$.
From Lemma \ref{posLemma}, we have $0\leq u_i^{n+1}$.
\end{proof}
Thus, if $x_i \in \Omega_1^n$ separates two approaching interfaces
then Lemma \ref{mergeLemma} implies that the colonies merge
smoothly without forming bulges or inducing oscillations. In the
case of a single moving interface, Lemma \ref{mergeLemma} implies
monotonicity, i.e. the absence of non physical oscillations. Note
that Lemma \ref{mergeLemma} includes the previous Lemma
\ref{finspeedLemma} for the trivial case where $x_i$ is interior
to $\Omega_1^n$.
\section{Numerical Simulations}\label{simulation_sec}
We present three sets of simulations of (\ref{model}) with the
nonstandard finite difference scheme (\ref{nsfd}). The first one
is one-dimensional in space. While this does not fully capture the
features of biofilm formation and growth it allows to investigate
the behavior of the numerical scheme. In the second part of this
section we carry out some two-dimensional numerical illustrations
of biofilm growth for constant net growth rate $k$. In both cases
we had for the parameters of the dimensionless equation $a=b=4$,
$k=0.1$, $\delta=10^{-8}$. In the third simulation study we apply
the discretisation scheme to a more complicated system, where
biofilm formation is controlled by two dissolved substrates. All
simulations were carried out using double accuracy arithmetics.
Codes were implemented in Fortran 95 on Linux based personal
computers. In all simulations the time-step was variable, as
\begin{equation}\label{dtchoice}
\Delta t = \min\Big\{1/k, \min_j {M \Delta x^2 \over D(u_j^n)} ,
0.1 \Big\}
\end{equation}
The first condition is according to (\ref{dtrestrict}), the last
one to keep the time-step bounded for the sake of accuracy. This
was also the reason to introduce the second condition, where $M$
is a user-defined constant. It renders the typical $\Delta t /
\Delta x^2$ dependency for parabolic problems and allows
comparability of results for different choices of $\Delta x$. In
our simulations, we use $M=40$; this choice of time-step is well
above the limitations of the explicit method. That said, this
time-step constraint was introduced as a trade-off between between
fast and accurate computation. The condition (\ref{dtrestrict2})
was implemented as an emergency brake, such that it becomes active
only if (\ref{dtchoice}) would lead to $\quad \max U^{n+1}>1$. It
never became activated.
\subsection{One-dimensional simulations}
\begin{figure}
\begin{center}
\includegraphics[width=0.488\textwidth]{figures/nsfd1}
\includegraphics[width=0.488\textwidth]{figures/nsfd2}
\end{center}
\caption{Simulation of biofilm formation in time: the inoculum at $t=0$
consist of (right) three sinusoidal colonies of different height and
(left) two colonies with constant densities of different height.}
\label{nsfdfig}
\end{figure}
The results of two numerical simulations of (\ref{model}) with the
finite difference method presented here are shown in Figure
\ref{nsfdfig}. In case the (I), the initial data were specified by the
function
\[
u(x,0)=\max( 0, -0.8\sin(7\pi x)(1-x^4)).
\]
In the case (II), the initial data were specified by
\[
u(x,0)=\begin{cases}
0.8, & 0.3\leq x \leq 0.4 \\
0.9, & 0.5 \leq x\leq 0.6, \\
0, & \mbox{elsewhere}.
\end{cases}
\]
In the case (III), finally one sinusoidal colony was
placed around the centre of the interval $(0,1)$.
This scenario is included to study the dependence of the interface
on the spatial resolution.
That is, in the first and third case the initial data are continuous
but not differentiable at the interface, while in the second
case they are only piecewise continuous. Both scenarios are known to be problematic if treated with standard techniques.
It is easily verified that in the first case the biomass is organised in three colonies originally.
Figure \ref{nsfdfig} shows the simulations of cases I and II.
They were carried out with a spatial step-size $\Delta x = 1/200$.
In the third scenario the dependency of the simulation results on
the spatial resolution was tested. To this end, the location of
the actual interface between $\Omega_1(t)$ and $\Omega_2(t)$ for
different step-sizes $\Delta x=1/N,\quad N=50, 100, 200, 400, 800$
is plotted in the $x$-$t$ plane. The second panel in Fig.
\ref{interfacefig} shows the solution surface for various $\Delta
x$. In Figures \ref{nsfdfig} and \ref{interfacefig} (right panel),
the solution surfaces $u(t,x)$ are plotted where $| u | >
10^{-16}$ for the sake of readability, assuming that $10^{-16}$ is
a reasonably good approximation of $0$ for our purpose.
\begin{figure}[t]
\includegraphics[width=0.48\textwidth]{figures/interfc}
\includegraphics[width=0.48\textwidth]{figures/nsfdsc3}
\caption{left: Location of the left interface between liquid and biofilm
in the third scenario for different choices of
$\Delta x=1/N, \quad N=50, 100, 200, 400, 800$; right:
solution of scenario 3 for various choices of $N$.} \label{interfacefig}
\end{figure}
As predicted by the theory outlined above, the new numerical
method gives non-negative and oscillation-free solutions. The
simulations also confirm the above statement that the diffusive
smearing around the interface is quantitatively small, {\it i.e.}
negligible. In both first scenarios the colonies of the inoculum
compress initially, that is the biomass density $u$ increases
while only very little spatial expansion is observed. As $u$
approaches 1, $\Omega_2$ expands exponentially due to the first
order growth term in (\ref{model}) and the squeezing property. In
both cases the colonies eventually merge. In the first scenario,
the first initial colony is bigger than the second one, which is
bigger than the third one. Accordingly the first and second colony
merge first. The actual interface was defined somewhat arbitrary
as the location at which $u$ becomes bigger than $10^{-5}$ (as an
approximation of $0$). Both figures show a very good agreement. It
can be seen that the computed interface converges as $\Delta x$
becomes smaller. Even for the coarsest discretisation the
deviations are accurate within $\Delta x$. This again confirms the
above comment that the smearing effect at the interface is small.
\subsection{A two-dimensional illustration: Formation of a homogeneous biofilm}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.48\textwidth]{figures/2d1}
\includegraphics[width=0.48\textwidth]{figures/2d12}
\\
\includegraphics[width=0.48\textwidth]{figures/2d20}
\includegraphics[width=0.48\textwidth]{figures/2d28}
\\
\includegraphics[width=0.48\textwidth]{figures/2d36}
\includegraphics[width=0.48\textwidth]{figures/2d42}
\\
\includegraphics[width=0.48\textwidth]{figures/2d48}
\includegraphics[width=0.48\textwidth]{figures/2d59}
\end{center}
\caption{Formation of a homogeneous, thick biofilm from originally heterogeneously distributed biomass.
The substratum $y=0$ is at the top of each picture.
Shown are snap shots at $t=0.01, 1.22, 2.1, 2.98, 3.85, 4.45, 5.06, 6.16$
(top left to bottom right). Biomass density $u$ is coded in a linear
greyscale: black corresponds to $u=0$, white to $u=1$.} \label{Fig2D}
\end{figure}
The simplified biofilm model (\ref{model}) with positive $k={\rm const}$
describes biofilm formation under non-limiting conditions. This
can be the case for nutrient rich environments, e.g. bacterial
growth in foods or relatively thin biofilms (in the initial
stages) in wastewater treatment plants. For many other biofilms
this is not a realistic assumption because often either nutrients
or oxygen or ph-value or other controlling substances are not
constant inside the biofilm. Typically, the thicker a biofilm is,
the more limited become nutrients in the deeper layers. Often,
decay due to lysis will even become dominant, leading to negative
local net growth rates. This spatially heterogeneous situation
will be addressed in the next section. Since we want to focus on
the spatial discretisation of models of this kind first, we
nevertheless accept the case $k={\rm const}$ for an illustrative
example. Under such non-limiting conditions it is expected (cf.
\cite{Eberl:2001}, \cite{Wanner:2006}) that even starting from a
heterogeneous initial distribution of biomass, a homogeneous
biofilm will form eventually and that the biomass distribution
within $\Omega_2$ will be more or less homogeneous. This will be
verified in the following simulations.
Initially biomass is randomly distributed in 15 different
locations across the substratum, {\it i.e.} the surface on which
the biofilm forms. That is, the biomass is placed initially in the
first layer of grid cells. The biomass density in each of these
locations is chosen randomly between 0 and 1. The computational
domain is the interval $[0,1] \times [0,0.3]$. We specify Neumann
boundary conditions at the substratum $y=0.3$ and on the
boundaries $x=0$, $x=1$. The first assumptions ensures that no
biomass leaves the system across the solid surface, while the
second assumption allows us to consider the computational domain
as part of a larger system (periodic conditions would do this as
well). At $y=0$ we specify homogeneous Dirichlet conditions. The
results of a simulation are shown in Figure \ref{Fig2D}. The rate
at which the biofilm colonies grow depends initially on the
biomass density at $t=0$. Colonies with larger $u(x,0)$ grow
faster and merge earlier with their neighbors. Eventually all
colonies merge and form a compact film as expected. As it was
already observed in the one-dimensional simulations, the biomass
density inside the biofilm approaches the maximum density rapidly
everywhere, due the homogeneous growth rate.
\subsection{A two-dimensional illustration: Formation of a
cluster-and-channel biofilm}\label{mushroom_sec}
We investigate now the more general case where the local growth
of biomass is controlled by the availability of required
substances. As an example we choose a heterotrophic biofilm that
is limited by two dissolved substrates, nutrients and oxygen.
This is a standard example of biofilm systems in wastewater
engineering and was defined by the International Water
Association's taskgroup on biofilm modeling as a first benchmark
system (BM1) for biofilm models, cf. \cite{Morgenroth:2004},
\cite{Wanner:2006}. The local net growth rate $k(t,x)$ in this
example is given by
\begin{equation}\label{reacs}
k(t,x)= k_1 {c_1(t,x) \over k_{s,1} + c_1(t,x)}
{c_2(t,x) \over k_{s,2} + c_2(t,x)} - k_3
\end{equation}
where the positive constant $k_1$ is the maximum specific growth
rate and the positive constant $k_3$ is the lysis rate describing
biomass deactivation. The positive parameters $k_{s,1}$ and
$k_{s,2}$ are the Monod half saturation constants. By $c_1$ and
$c_2$ we denote the concentrations of the dissolved substrates.
They are governed by the semi-linear diffusion-reaction equations
\begin{equation}\label{c1eq}
c_{1,t} = D_1 \Delta_x c_1 - k_4 u {c_1 \over k_{s,1} + c_1}
{c_2 \over k_{s,2} + c_2}
\end{equation}
and
\begin{equation}\label{c2eq}
c_{2,t} = D_2 \Delta_x c_2 - k_5 u {c_1 \over k_{s,1} + c_1} {c_2
\over k_{s,2} + c_2}
\end{equation}
where the consumption rates $k_4$ and $k_5$ are essentially the maximum
specific growth rate divided by the respective yield factors,
which quantify how much substrate is needed to produce one unit mass
of microbial biomass. The positive parameters $D_1$ and $D_2$ are the
diffusion coefficients of both substrates. For small molecules these
are constants, attaining the same values in the biofilm and in the
surrounding bulk phase.
For our simulation we choose the same reaction and diffusion parameters
as in the benchmark problem specification in \cite{Morgenroth:2004}.
The model is completed by initial and boundary conditions for
$c_1$ and $c_2$. We choose the computational domain to be
rectangular of size $L_1 \times L_2$ and specify the mixed
boundary conditions for $j=1,2$
\[
c_{j}=c_{\infty,j} \quad\mbox{for}\quad x_2=L_2
\]
and
\[
{\partial c_j \over \partial n} =0 \quad\mbox{elsewhere on}\quad \partial \Omega
\]
The initial conditions for $c_{1,2}$ are
\[
c_j(x,0)= c_{\infty,j}
\]
For the biomass density $u$ we specify the same initial and
boundary conditions as in the previous example.
The purpose of this example is to show that the numerical scheme
is able to simulate the formation of spatially heterogeneous of
mushroom-shaped cluster-and-channel biofilm architectures. It is
well understood in the biofilm literature, cf \cite{Wanner:2006},
that a good predictor for the surface heterogeneity of a biofilm
is how biomass growth terms relate to substrate availability. If
substrate does not become limited, a flat, spatially homogeneous
biofilm is eventually obtained, as in the previous example. If
substrate in the biofilm becomes limited, growth of biomass slows
down. Microbial inner-species competition for nutrients leads to
irregular biofilm morphologies. Eventually larger biofilm colonies
closer to the food source become dominating and smaller colonies
stay back in their development. With our choice of parameters from
\cite{Morgenroth:2004}, the growth terms are fixed. Hence, we use
substrate availability to control the biofilm structure. It is
quantified in terms of the diffusion coefficients (also fixed from
\cite{Morgenroth:2004}) and the environmental conditions, in
particular the Dirichlet values $c_{\infty,j}$ and the diffusion
length $L_2$. We choose the bulk concentrations in the same range
as the Monod half saturation concentrations,
$c_{\infty,1}=k_{s,1}$ and $c_{\infty,2}=2k_{s,2}$. The maximum
principle for parabolic equations implies that these values are
upper estimates for the concentrations $c_1$ and $c_2$ and that
inside the domain and, hence, in the biofilm lower concentration
values will be obtained. Thus, our choice of environmental
concentrations ensures that substrates are not available in
abundance. The system height $L_2$ is the distance over which
substrates need to be transported by diffusion to reach the
bacteria at the substratum. The smaller this value is, the more
compact a biofilm develops \cite{Wanner:2006}. We choose
$L_2=L_1=500\mu m$, i.e. carry out our simulations in a square
domain.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.47\textwidth]{figures/BM1_001}
\includegraphics[width=0.47\textwidth]{figures/BM1_005}
\\
\includegraphics[width=0.47\textwidth]{figures/BM1_007}
\includegraphics[width=0.47\textwidth]{figures/BM1_010}
\\
\includegraphics[width=0.47\textwidth]{figures/BM1_020}
\includegraphics[width=0.47\textwidth]{figures/BM1_040}
\end{center}
\caption{Formation of a cluster-and-channel biofilm morphology
(top left to bottom right): Shown are the biofilm/liquid
interface and the limiting oxygen concentration $c_2(t,x)$
in time $t= 0.02, 4.03, 6.05, 9.08, 19.12, 39.13 d$. The food source
is located at the right boundary $x_2=L_2$.} \label{BM1_sim_fig}
\end{figure}
\begin{figure*}
\begin{center}
\includegraphics[width=0.48\textwidth]{figures/BM1_080}
\includegraphics[width=0.48\textwidth]{figures/BM1_120}
\\
\includegraphics[width=0.48\textwidth]{figures/BM1_160}
\includegraphics[width=0.48\textwidth]{figures/BM1_200}
\\
\includegraphics[width=0.48\textwidth]{figures/BM1_230}
\includegraphics[width=0.48\textwidth]{figures/BM1_250}
\end{center}
\caption{Figure \ref{BM1_sim_fig} continued for
$t= 79.13, 119.13, 159.13, 199.13, 229.20, 249.35 d$}
\end{figure*}
For the numerical solution of (\ref{c1eq}) and (\ref{c2eq}) we use
a standard finite volume scheme on the same grid on which
(\ref{model}) is discretised, 2nd order in space and 1st order in
time. The simulation results of this example are shown in Figure
\ref{BM1_sim_fig}; plotted are the biofilm/liquid interface and
the concentration field of oxygen, $c_2$. The nutrient
concentration field $c_1$ is qualitative similar, albeit on a
higher quantitative level. Initially the substratum is inoculated
by few small bacterial colonies. As time increases, first a
smooth, almost homogeneous, smooth biofilm layer develops which
than changes in a cluster-and-channel morphology as substrates
become limited. It should be noted that due to the lysis term in
(\ref{reacs}) and due to substrate limitations the active biomass
density remains distinctively separated from the maximum biomass
density, i.e. $u<1$ clearly. Thus, the simulation results agree
with our a priori expectations on model behavior.
\section{Conclusion}
A finite difference scheme was developed for a parabolic evolution
equation with two distinct nonlinear effects in the
density-dependent diffusion coefficient. One is degeneracy as in
the porous medium equation, while the other one is a singularity
as the dependent variable approaches its {\it a priori} known
maximum value. In the numerical method, the nonlinearity of the
diffusion operator was handled in a non-local representation in
time. It was shown that the numerical method renders the essential
qualitative features of the solution of the continuous problem. In
particular it is free of non physical oscillations that often are
observed in problems of this kind, it shows the finite speed of
propagation of interfaces between the regions $\Omega_1(t)$ where
$u=0$ and $\Omega_2(t)$ where $u>0$, and it guarantees that the
upper bound of the solution.
The nonstandard finite difference scheme for the proto-type
biofilm model (\ref{model}) is based on a non-local in time
representation of the nonlinear diffusive flux. The method was
constructed in a manner that allows a straightforward application
to more complicated biofilm systems with several particulate
substances as in section \ref{simulation_sec}\ref{mushroom_sec}
and dissolved substrates, {\it i.e.} for models as studied in
{\it e.g.} \cite{Eberl:2003}, \cite{Eberl:2004}. Although only
one- and two-dimensional simulations have been carried out here
for the illustration of the new method, an application to more
realistic three-dimensional biofilm descriptions is possible. In
fact, the method (\ref{nsfd}) itself is formulated in the general
setup and the analytical results were derived for the
three-dimensional case. This makes the new method attractive for
further studies in mathematical modeling of biofilms.
\subsection*{Acknowledgments}
This study was supported by the
Volkswagen Stiftung, Hannover, Germany with a grant in {\it
Interdisciplinary Environmental Research}, by the NSERC Canada with a
{\it Discovery Grant} and by the {\it Advanced Foods and
Materials Network of Centres of Excellence (AFMnet)} with a
project grant. The original draft version of this manuscript was
written while the first author visited the GSF Research Centre in
Neuherberg, Germany. The authors also wish to thank Messoud
Efendiev for many stimulating discussions. The CLSM micrograph was
made available by Heidi Schraft, Dept. Biology, Lakehead
University, Thunder Bay, ON.
\thebibliography{00}
\bibitem{Lubuma:2001}
Anguelov R, Lubuma J. M. S.
Contributions to the mathematics of the nonstandard finite difference
method and its applications, {\it Num. Meth. PDE} 17:518-543, 2001
\bibitem{Baines}
Baines M. J., Hubbard M. E., Jimack P. K., A moving mesh finite
element algorithm for the adaptive solution of time dependent PDEs
with moving boundaries, {\it Appl. Num. Math.}, 54:450-469, 2005
\bibitem{Costerton}
Costerton JW, Lewandowsky Z, Caldwell DE, Korber DR, Lappin-Scott HM,
Microbial Biofilms, {\it Annu. Rev. Microbiol.} 49:711-745, 1995
\bibitem{Duvnjak:2004}
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(3), 2001
\bibitem{Eberl:2003}
Eberl H. J., Efendiev M. A. A Transient Density Dependent
Diffusion-Reaction Model for the Limitation of Antibiotic
Penetration in Biofilms. {\it El. J. Diff Eq}. CS10:123-142, 2003
\bibitem{Eberl:2004}
Eberl H. J., A deterministic continuum model for the formation of EPS
in heterogeneous biofilm architectures, {\it Proc. Biofilms
2004}, Las Vegas, 2004
\bibitem{Eberl:2006}
Eberl H. J., Schraft H., A diffusion-reaction model of a mixed culture
biofilm arising in food safety studies, {\it accepted}
\bibitem{Efendiev:2005}
Efendiev M. A., 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 Kyoto}, 1258:49-71, 2002
\bibitem{Epperson}
Epperson J. F., {\it An Introduction to Numerical Methods and Analysis},
Wiley \& Sons, 2002
\bibitem{Flemming}
Flemming H. C., Biofilme -- das Leben am Rande der Wasserphase, {\it
Nachr. Chemie}, 48:442-447, 2000
\bibitem{Hackbusch}
Hackbusch W., {\it Theorie und Numerik elliptischer
Differentialgleichungen}, Teubner, Stuttgart, 1986
\bibitem{Hundsdorfer}
Hundsdorfer W., Verweer J. G., {\it Numerical Solution of
Time-Dependent Advection-Diffusion-Reaction Equations}, Springer,
2003
\bibitem{Khassehkhan:2006}
Khassekhan H., Eberl H. J., Interface tracking for a non-linear
degenerated diffusion-reaction equation describing biofilm
formation, {\it accepted}
\bibitem{vanL:1995}
Van Loosdrecht M. C. M., Eikelboom D., Gjaltema A., Mulder A., Tijhuis L.,
Heijnen J. J., Biofilm Structures, {\it Wat. Sci. Tech.} 32(8):35-43, 1995
\bibitem{vanL:1997}
Van Loosdrecht M. C. M. , Picioreanu C., Heijnen J. J., A more unifying
hypothesis for the structure of microbial biofilms, {\it FEMS
Microb. Ecol.}, 24:181–183, 1997
\bibitem{Mickens:2000}
Mickens R. E. Nonstandard finite difference schemes, {\it in}
Mickens R. E. (ed), {\it Applications of nonstandard finite
difference schemes}, World Scientific, Singapore, 2000
\bibitem{Morgenroth:2004}
Morgenroth E., Eberl H. J., Van Loosdrecht M.C. M., Noguera D. R., Picioreanu
C, Rittmann BE, Schwarz AO, Wanner O. Results from the single
species benchmark problem (BM1), {\it Wat. Sci. Tech.},
49(11/12):145-154, 2004
\bibitem{Morton:1994}
Morton K. W., Mayers D. F., {\it Numerical solution of partial
differential equations}, Cambridge Univ. Press, 1994
\bibitem{Pici:2003}
Picioreanu C., Van Loosdrecht M. C. M., Use of mathematical modelling to
study biofilm development and morphology. In: Lens P et al (eds).
{\it Biofilms in Medicine, Industry and Environmental
Biotechnology – Characteristics, Analysis and Control}, IWA
Publishing. pp 413-437, 2003
\bibitem{Saad}
Saad Y., {\it Iterative Methods for Sparse Linear Systems}, 2nd edition, 2000
\bibitem{Stewart}
Stewart P. S., Multicellular nature of biofilm protection from
antimicrobial agents, in: McBain A et al (eds), {\it Biofilm
Communities: Order from Chaos}, BioLine, 2003
\bibitem{Wanner:2006}
Wanner O., Eberl H. J., Van Loosdrecht M. C. M., Morgenroth E., Noguera D. R.,
Picioreanu C, Rittmann BE. {\it Mathematical Modeling of
Biofilms}, IWA Publishing, 2006
\bibitem{Watnick}
Watnick P., Kolter R., Biofilm -- City of Microbes (Minireview),
{\it J. Bacteriol.}, 182(10):2675-2679, 2000
\bibitem{Wilderer:2002}
Wilderer P. A., Bungartz H. J., Lemmer H., Wagner M., Keller J., Wuertz S.,
Modern Scientific Methods and their Potential in Wastewater
Science and Technology, {\it Water Res.} 36:370-393, 2002
\bibitem{Xu:1998}
Xu K. D., Stewart P. S., Xia F., Huang C.-T., McFeters G. A., Spatial
physiological heterogeneity in Pseudomonas aeruginosa biofilm is
determined by oxygen availability, {\it Appl. and Env.
Microbiol.}, 64(10):4035-4039, 1998
\endthebibliography
\end{document}