\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}