# Numerical Analysis of a Parabolic Variational Inequality System Modeling
Biofilm Growth at the Porescale

Azhar Alhammali and Malgorzata Peszynska  Department of Mathematics, Oregon
State University, Corvallis, OR 97330
email<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
## Abstract

In this paper we consider a system of two coupled nonlinear diffusion–reaction
partial differential equations (PDEs) which model the growth of biofilm and
consumption of the nutrient. At the scale of interest the biofilm density is
subject to a pointwise constraint, thus the biofilm PDE is framed as a
parabolic variational inequality. We derive rigorous error estimates for a
finite element (FE) approximation to the coupled nonlinear system and confirm
experimentally that the numerical approximation converges at the predicted
rate. We also show simulations in which we track the free boundary in the
domains which resemble the pore scale geometry and in which we test the
different modeling assumptions.

## Keywords

Parabolic variational inequality, coupled system, biofilm growth, finite
elements, error estimates, semismooth Newton solver.

##  1\. Introduction

Coupled bio-chemical interactions involving microbial cells are very important
in many applications including the food industry and medicine. In this paper
we focus on the applications involving “selective plugging” by biofilm in
porous media such as, e.g., in microbial enhanced oil recovery (MEOR) and in
carbon sequestration. In these applications the microbes coat the grains of
the porous medium and plug the high-permeability zones thereby enhancing
overall recovery [16] or preventing leakage by promoting mineralization [10,
25, 7].

The mathematical model we consider was proposed in [18] to describe some
microorganisms and nutrient suspended in ambient fluid. Their densities (or
concentrations) are denoted by $B(x,t)$ and $N(x,t)$, respectively, and these
suspensions spread by diffusion as in

(1a) $\displaystyle\frac{\partial B}{\partial t}-\nabla\cdot(D_{B}\;\nabla
B)-\Lambda$ $\displaystyle=F(B,N)+f(x,t),$ (1b) $\displaystyle\frac{\partial
N}{\partial t}-\nabla\cdot(D_{N}\;\nabla N)$ $\displaystyle=G(B,N)+g(x,t).$

The nonnegative diffusivities $D_{B},D_{N}$ are discussed in the sequel, and
their choice plays an important role in the dynamics of spreading. For the
growth $F$ and consumption $G$ a common choice are Monod functions

(2) $\displaystyle
F(B,N)=\kappa_{B}P(N)B,\;G(B,N)=-\kappa_{N}P(N)B,\;\;P(N)=\frac{N}{N+N_{0}},$

where $k_{B}\geq k_{N}>0,N_{0}>0$ are known constants. The external source and
sink terms are denoted by $f,g$. The model is completed by boundary and
initial conditions specified later. Here we assume that the fluid is at rest,
so the model (1) does not include transport by advection. See references to
the biofilm–nutrient models in, e.g., [10, 9, 18].

The term $\Lambda$ in (1a) enforces the constraint $B(x,t)\leq B^{*}$ written
as Nonlinear Complementarity Condition (NCC)

(3) $\displaystyle B(x,t)\leq B^{*},\;\;\Lambda\leq 0,\;\;\Lambda(B-B^{*})=0.$

This salient feature of the model in [18] leads to the main challenges
addressed in this paper. We explain the need for (3), and the role and the
relationship between $\Lambda$, $B$, and $B^{*}$ in what follows.

Without the constraint (3) is not active, i.e., with $\Lambda=0$ in (1), the
growth of the solution $B(x,t)$ predicted by (1) under some conditions such as
no-flux boundary conditions and with abundant nutrient can be exponential.
This is only realistic at large spatial scales of m or km, typical for
laboratory containers and reservoirs, or at small concentrations. Other
features of (1), e.g., of dynamics of $B+\tfrac{\kappa_{B}}{\kappa_{N}}N$
under various assumptions on $f$ and $g$ can be analyzed, but are out of our
scope.

We are interested in the growth of microbes at the mm scale relevant, e.g., to
biomedical applications, or in the biofilm phase at the pore-scale of $\mu$m.
At that scale the constraint (3) comes from the fact that the microbial cells
have a finite size, so their density is limited. In addition, when $B\approx
B^{*}$ is reached, the cells form the biofilm, a gel-like substance made of
microbial cells which are surrounded by an extracellular polymer substance
(EPS) matrix. The EPS matrix protects the cells, which tend to aggregate, and
it is not penetrated by larger molecules, e.g., of imaging contrast. The
interface between the biofilm domain denoted by
$\Omega^{*}(t)=\\{x\in\Omega\mathrel{\mathop{\ordinarycolon}}B(x,t)=B^{*}\\}$,
and the surrounding fluid domain
$\Omega^{-}(t)=\\{x\mathrel{\mathop{\ordinarycolon}}B(x,t)<B^{*}\\}$ is a free
boundary $\Gamma(t)=\partial\Omega^{*}(t)\cap\partial\Omega^{-}(t)$ visible to
human eye and easily recognized by x-ray $\mu$-CT tomography or other
microscopy equipment; see [18, 25, 7]. The growth of $\Omega^{*}(t)$ proceeds
through interface only; this explains the tapering of exponential growth
reported, e.g., in [7]. Finding $\Lambda\neq 0$ and the free boundary
$\Gamma(t)$ is part of the evolution problem, and the PDE (1a) is formally a
nonlinear parabolic variational inequality (PVI) coupled to the nutrient
dynamics.

Theoretical and practical challenges and contributions of this paper. First,
the solutions to free-boundary problems and variational inequalities such as
(1) have low regularity. For this reason we consider only a linear Galerkin
Finite Element approximation coupled with fully implicit time-stepping for
which we prove convergence roughly of $O(h)$ in the $L^{2}$ norm if $\Delta
t=O(h)$. The main theoretical challenges are the coupled nonlinear nature of
the system and the constraint. We confirm the rate with numerical experiments
and describe a robust and efficient nonlinear solver. Second, we work in
complicated geometries obtained from imaging, such as the void space between
the grains of a porous medium at the porescale, and study some model variants
which support qualitative behavior in experiments [25, 7, 18]. The
computational challenge is to carry out the simulations in complicated domains
which require careful grid generation and pre- and post- processing.

Literature remarks. Numerical simulation of biofilm growth is of interest
because the physical experiments with biofilm and its imaging are quite
difficult [18, 10, 25, 7]. Various models for biomass growth and biofilm
growth address the growth at the porescale and the permeability of the porous
domains plugged up by biofilm; however, their focus is selective. In
particular, the models [10, 19] do not impose the constraint (3), while the
model in [24] considers a simplified free boundary geometry in a strip. In
turn, the hybrid differential–lattice models in [9] and [21] enforce (3), but
their structure does not permit numerical convergence analysis. A complex
family of phase field models is considered in [29], but its computational
complexity seems prohibitive at the porescale. Further considerations of
complexity and uncertainty in upscaling are in [8]. Finally, the model in [18]
is more general than (1), since it considers advective transport, coupled
Navier-Stokes flow solver, as well as dynamic upscaling of the flow solutions
to obtain permeability. However, there is no numerical analysis in [18].

The literature on the topic of numerical approximation of parabolic systems is
enormous; we recall only a few results which directly guide our work. The
classical results by Wheeler [28] for unconstrained scalar problems (see also
[22], Chapter 13) establish second order $l^{\infty}(L^{2})$ error estimate of
the linear finite approximation in space and backward Euler in time, but
require high regularity of the solution $B$, such as
$\frac{\partial^{2}B}{\partial t^{2}}\in L^{2}(Q)$. The works on elliptic and
parabolic variational inequalities (EVI and PVI), such as on first order
convergence of Galerkin FE [6, 11, 2], account for the lack of regularity of
their solutions but handle only scalar linear problems. The analyses of
Johnson [14] for the subproblem (1a) are the closest to what we consider here
but they require $F=0$ and $D_{B}=$const. The paper [27] considers more
general finite differencing in time, but requires stronger regularity to
obtain first order of convergence in $l^{\infty}(L^{2})$. For systems without
constraints, the analysis in [4] allows degeneracy in the diffusivity in
$D_{B}(B)$ and gives $O(h^{1/2})$ order of convergence in
$l^{\infty}(L^{2})$-norm. Finally there is more work on FE approximation to
free boundary problems [26, 17, 15]; many recent papers focus on adaptivity
which we do not consider here.

Plan of the paper. In Sec. 2 we provide the details on the model. In Sec. 3 we
define the FE approximation and prove our estimate. In Sec. 4 we describe the
solution method framed as a semi-smooth Newton solver for the NCC. In Sec. 5
we present numerical experiments which confirm the predicted convergence rate
and illustrate different qualitative behavior depending on the assumed
nonlinear diffusivity models.

##  2\. Background and formal setting

In this section we provide details of our model which is a parabolic system of
PDEs under a constraint, i.e., of parabolic variational inequalities. We
discuss the rigorous setting including the PDE model in the sense of
distributions and its weak formulation.

Let $\Omega$ be an open bounded domain in $\mathbb{R}^{d}$; $d$=$2$, with a
smooth boundary $\partial\Omega$. (In our numerical experiments cover we
consider $d=1,2,3$.) Throughout the paper, we use standard notation on
$L^{2}(\Omega)$ and Sobolev spaces $H^{2}(\Omega)$ and $H_{0}^{1}(\Omega)$;
see, e.g., [20]. Let $\|\cdot\|_{0}$ denote the norm on $L^{2}(\Omega)$, and
$\|\cdot\|_{s}$ the norm on $H^{s}(\Omega)$ where $s$ is a nonnegative
integer. $J=[0,T]$ denotes the time interval with final time $T>0$. We write
$u\in L^{2}(H_{0}^{1})$ to mean $u\in L^{2}(J;H_{0}^{1}(\Omega))$; similar
shorthand is applied to other notation on functional spaces. In particular,
the notation $L^{2}(L^{2})$ means $L^{2}(J;L^{2}(\Omega))=L^{2}(Q)$ with
$Q=\Omega\times J$.

First we discuss the constraint (3), and denote
$\mathbb{R}^{*}={(-\infty,B^{*}]}$. The indicator function
$I_{\mathbb{R}^{*}}$ takes value $0$ in the set $R^{*}$ and $\infty$ outside:
$I_{\mathbb{R}^{*}}(B)=\left\\{\begin{array}[]{ll}0&B\in\mathbb{R}^{*},\\\
+\infty&\mbox{otherwise}.\end{array}\right.$ Denote by $\partial
I_{\mathbb{R}^{*}}$ the subgradient of $I_{\mathbb{R}^{*}}$, the multivalued
constraint graph

$\partial
I_{\mathbb{R}^{*}}=\left((-\infty,B^{*})\times\\{0\\}\right)\cup\left(\\{B^{*}\\}\times[0,\infty)\right).$

Now we write the model (1) formally in the sense of distributions using the
operator $\partial I_{\mathbb{R}^{*}}(B)$ instead of $-\Lambda$. Formally,
when $\partial I_{\mathbb{R}^{*}}(B)$ appears in an equation, $B$ must satisfy
(3), i.e., be in the domain $K$ of $\partial I_{\mathbb{R}^{*}}(B)$, where the
set $K\mathrel{\mathop{\ordinarycolon}}=\\{B\in
H_{0}^{1}(\Omega);\;B\leq{B^{\ast}}\;\;{a.e}\;\;\mbox{on}\;\;\Omega\\}$ is a
convex subset of $H_{0}^{1}(\Omega)$. The model reads

(4a) $\displaystyle\frac{\partial B}{\partial t}-\nabla\cdot(D_{B}\;\nabla
B)+\partial I_{(-\infty,B^{*}]}(B)$ $\displaystyle\ni F(B,N)+f(x,t),$
$\displaystyle\;\;\mathrm{a.e.\ in\ }\Omega,\;t>0,$ (4b)
$\displaystyle\frac{\partial N}{\partial t}-\nabla\cdot(D_{N}\;\nabla N)$
$\displaystyle=G(B,N)+g(x,t),$ $\displaystyle\mathrm{a.e.\ in\ }\Omega,\;t>0.$
The operator $\partial I_{\mathbb{R}^{*}}(B)$ is multivalued, so we use the
symbol $\ni$. However, there is no ambiguity in the choice of the particular
selection $\Lambda\in\partial I_{\mathbb{R}^{*}}(B)$, since this selection is
always exactly the value that keeps $B(x,t)$ in the convex set $K$.

The model is completed with the initial and boundary conditions as follows

(4c) $\displaystyle B(x,0)$ $\displaystyle=B_{init}(x)$
$\displaystyle\;\;x\in\Omega,$ (4d) $\displaystyle N(x,0)$
$\displaystyle=N_{init}(x)$ $\displaystyle\;\;x\in\Omega,$ (4e) $\displaystyle
B(s,t)$ $\displaystyle=0,$ $\displaystyle\;\;s\in\partial\Omega,\;t>0,$ (4f)
$\displaystyle N(s,t)$ $\displaystyle=0,$
$\displaystyle\;\;s\in\partial\Omega,\;t>0.$

######  Remark 1.

The choice of Dirichlet conditions (4e)–(4f) is merely for the easiness of
analysis. In practice, Neumann conditions are realistic in simulations of a
closed system. Furthermore, while our analysis is restricted to the case where
$D_{B},D_{N}$ do not depend on $B,N$, we consider other variants in numerical
experiments.

The weak form of (4) is the variational inequality

(5a) $\displaystyle B(t)\in
K\mathrel{\mathop{\ordinarycolon}}\left(\frac{\partial B}{\partial
t},\psi-B\right)+(D_{B}\nabla B,\nabla(\psi-B))$ $\displaystyle\geq$
$\displaystyle(F(B,N){+f},\psi-B),\forall\psi\in K,$ (5b) $\displaystyle
N(t)\in H_{0}^{1}(\Omega)\mathrel{\mathop{\ordinarycolon}}\left(\frac{\partial
N}{\partial t},\chi\right)+(D_{N}\nabla N,\nabla\chi)$ $\displaystyle=$
$\displaystyle(G(B,N)+g,\chi),\forall\chi\in H_{0}^{1}(\Omega),$ (5c)
$\displaystyle B(x,0)$ $\displaystyle=$ $\displaystyle{B_{init}(x)},$ (5d)
$\displaystyle N(x,0)$ $\displaystyle=$ $\displaystyle{N_{init}(x)},$

where $(\cdot,\cdot)$ is the duality pairing; (also, the scalar product on
$L^{2}(\Omega)$); see [5, 20, 23]. The symbol $\geq$ in (5a) expresses the
fact that the admissible solutions and test functions are defined by
inequalities, thus they constitute a convex set $K$ rather than the space
$H_{0}^{1}(\Omega)$. In the numerical model the symbol “$\geq$” is replaced by
“$=$” and the presence of a Lagrange multiplier.

For the needs of subsequent numerical analysis model we make some assumptions
on the data and the regularity of the solutions typical for parabolic
variational inequalities.

######  Assumption 1.

  1. (A)

$D_{B},D_{N}\mathrel{\mathop{\ordinarycolon}}\Omega\to\mathbb{R}$ are
Lipschitz continuous with a Lipschitz constant $R$, and $0<\mu_{1}\leq
D_{B}(x)\leq\mu_{2},\;\;0<\nu_{1}\leq D_{N}(x)\leq\nu_{2}\;\;\mbox{for
$x\in\Omega$}$ and some $\mu_{1},\mu_{2},\nu_{1},\nu_{2}$.

  2. (B)

$F(B,N)$ and $G(B,N)$ are smooth functions, Lipschitz with respect to $B$ and
to $N$ with a Lipschitz constant $M$. Further assume $F$ and $G$ are linear in
$B$ and uniformly bounded in $N$ on $\mathbb{R}\times\mathbb{R}^{+}$,

(6) $\displaystyle F(B,0)=0=F(0,N),\;G(B,0)=0=G(0,N)\;\forall
B,N\in\mathbb{R}.$

[We note that Monod growth functions defined by (2) satisfy these conditions.
]

  3. (C)

$f,g\in C(L^{\infty})$, $\frac{\partial f}{\partial t}\in L^{2}(L^{\infty})$.

  4. (D)

$B^{*}>0$ is given.

  5. (E)

$B_{init},N_{init}\in W^{2,\infty}$, and $B_{init}\leq B^{*}$.

  6. (F)

$B,N\in L^{\infty}(W^{2,p})$; $1\leq p<\infty$, and $\frac{\partial
B}{\partial t},\frac{\partial N}{\partial t}\in L^{2}(H_{0}^{1})\cap
L^{\infty}(Q)$. Also, $\frac{\partial^{2}N}{\partial t^{2}}\in L^{2}(Q)$.

######  Remark 2.

The assumption $\frac{\partial^{2}N}{\partial t^{2}}\in L^{2}(Q)$ can be
easily dropped. This is useful if the second PDE of system (4) is constrained.
The assumptions F on regularity of $B$ are realistic; in general
$\tfrac{\partial^{2}B}{\partial t^{2}}\not\in L^{2}(H^{-1})$.

We denote by $\chi_{A}$ the characteristic function of a set $A$. We shall use
the elementary inequality

(7) $ab\leq\frac{\epsilon}{2}a^{2}+\frac{1}{2\epsilon}b^{2}\;\;\forall
a,b\in\mathbb{R},\;\epsilon>0,$

and Poincaré-Friedrichs inequality

(8) $\|\psi\|_{0}\leq C_{PF}\|\nabla\psi\|_{0},\;\;\psi\in
H_{0}^{1}(\Omega);\;\;C_{PF}>0,$

##  3\. The Approximation and Error Estimate

In this section we formulate a fully discrete approximation to (5) using
backward Euler scheme in time and piecewise linear finite element method in
space. The notation is fairly standard; see, e.g., [22].

###  3.1. Discrete Problem

For $h>0$, let $\mathcal{T}_{h}=\\{T_{i}\\}$ be a conforming triangulation of
$\Omega$ such that the length of each side of any element $T_{i}$ is at most
$h$. We assume that no vertex of any triangle lies in the interior of another
triangle, and all angles of the triangles are bounded below by a fixed
positive constant, and that $\Omega=\bigcup_{T\in\mathcal{T}_{h}}$.

Define the space $V_{h}$ for piecewise linear approximations and its convex
subset $K_{h}$

$\displaystyle V_{h}=\\{\psi\in
C(\bar{\Omega})\mathrel{\mathop{\ordinarycolon}}\psi\;\mbox{ is linear on each
}T_{i},\;\psi=0\mbox{ on }\partial\Omega\\},\;\;K_{h}=V_{h}\cap K.$

In what follows we work with the product space
$\mathbf{H}=H_{0}^{1}(\Omega)\times H_{0}^{1}(\Omega)$, with the norm
$\|(B,N)\|_{\mathbf{H}}=\sqrt{\|B\|_{1}^{2}+\|N\|_{1}^{2}}$. Let
$\mathbf{V}_{h}=V_{h}\times V_{h}\subset\mathbf{H}$.

We use uniform time-stepping on $J$ for simplicity of the exposition. One can
also easily use non-uniform or adaptive time-stepping, but we skip details.
Let $\Delta t=N_{T}^{-1}T;\;N_{T}$ a positive integer, $t_{n}=n\Delta t$,
$J_{n}=(t_{n},t_{n+1}]$. We denote $\psi^{n}=\psi(t_{n})$. In this section
$\partial$ refers to a difference
$\partial\psi^{n}=(\psi^{n+1}-\psi^{n})/\Delta t$ for any $\psi$. Let
$\Upsilon=\\{t_{0},\ldots,t_{N_{T}}\\}$ be the set of time steps.

We approximate (5) as follows. We seek
$B_{h}\mathrel{\mathop{\ordinarycolon}}\Upsilon\to K_{h}$ and
$N_{h}\mathrel{\mathop{\ordinarycolon}}\Upsilon\to V_{h}$ which satisfy, for
$n=0,\ldots,N_{T}-1$

(9a) $(\partial B_{h}^{n},\psi-B_{h}^{n+1})+(D_{B}\nabla
B_{h}^{n+1},\nabla\psi-\nabla B_{h}^{n+1})\\\
\geq(F(B_{h}^{n+1},N_{h}^{n+1})+f^{n+1},\psi-B_{h}^{n+1})\;\;\forall\psi\in
K_{h},$ (9b) $\displaystyle(\partial N_{h}^{n},\chi)+(D_{N}\nabla
N_{h}^{n+1},\nabla\chi)$
$\displaystyle=(G(B_{h}^{n+1},N_{h}^{n+1})+g^{n+1},\chi)$
$\displaystyle\forall\chi\in V_{h}.$ (9c)
$\displaystyle\|B_{h}^{0}-B_{init}\|_{0}$ $\displaystyle\leq Ch,$ (9d)
$\displaystyle B_{h}^{0}$ $\displaystyle=I_{h}B_{init}$ (9e) $\displaystyle
N_{h}^{0}$ $\displaystyle=I_{h}N_{init}.$

Here $I_{h}$ is the local smoother defined in the sequel.

###  3.2. Existence and uniqueness of solutions

Here we prove that the fully discrete problem (9) has a unique solution under
mild assumptions on the size of $\Delta t$.

######  Lemma 3.1.

Assume that (i) $\gamma=\min\\{\mu_{1},\nu_{1}\\}>2MC_{PF}^{2}$. If this does
not hold, (ii) assume that $\Delta t$ is sufficiently small, and in particular
that $\Delta t<\frac{C_{PF}^{2}}{2MC_{PF}^{2}-\gamma}$. Then the problem (9)
has a unique solution.

###### Proof.

For $n=0,\ldots,N_{T}-1$, (9) can be written as

$\displaystyle a_{\Delta t}(N_{h}^{n+1};B_{h}^{n+1},\psi-B_{h}^{n+1})$
$\displaystyle\geq(B_{h}^{n}+\Delta tf^{n+1},\psi-
B_{h}^{n+1})\;\;\forall\psi\in K_{h},$ $\displaystyle b_{\Delta
t}(B_{h}^{n+1};N_{h}^{n+1},\chi)$ $\displaystyle=(N_{h}^{n}+\Delta
tg^{n+1},\chi)\;\forall\chi\in V_{h},$

where $a_{\Delta t}\mathrel{\mathop{\ordinarycolon}}V_{h}\times
V_{h}\longrightarrow\mathbb{R}$, $b_{\Delta
t}\mathrel{\mathop{\ordinarycolon}}V_{h}\times V_{h}\longrightarrow\mathbb{R}$
are defined as

$\displaystyle a_{\Delta t}(N;B,\psi)$ $\displaystyle=(B,\psi)-\Delta
t(F(B,N),\psi)+\Delta t(D_{B}\nabla B,\nabla\psi),\;\;\forall(B,\psi)\in
V_{h}\times V_{h},$ $\displaystyle b_{\Delta t}(B;N,\chi)$
$\displaystyle=(N,\chi)-\Delta t(G(B,N),\chi)+\Delta t(D_{N}\nabla
N,\nabla\chi),\;\;\forall(N,\chi)\in V_{h}\times V_{h}.$

In what follows we suppress the indices $h$, and $n$+$1$ in
$N_{h}^{n+1},B_{h}^{n+1}$. Define
$L\mathrel{\mathop{\ordinarycolon}}\mathbf{H}\longrightarrow\mathbf{H}^{\prime}$
as

$(L(B,N),(\psi,\chi))=a_{\Delta t}(N;B,\psi)+b_{\Delta
t}(B;N,\psi)\;\;\forall(B,N),(\psi,\chi)\in\mathbf{H}.$

We will show that $L$ is monotone, continuous, and coercive. From this it
follows that $L$ is demicontinuous, and we can apply ([3] Corollary 2.8, page
73) to prove existence of the solution in $K_{h}\times V_{h}$. If $L$ is also
strictly monotone, then the solution is unique; see, e.g, the existence proof
in ([20], Corollary 7.1, page 84).

To show continuity, we apply Assumption 1, parts A and B, and Cauchy-Schwarz
inequality and get

$|\left(L(B,N),(\psi,\chi)\right)|\leq
C\|(B,N)\|_{\mathbf{H}}\|(\psi,\chi)\|_{\mathbf{H}},\;\mbox{for some constant
$C>0$}.$

Recall $L$ is strictly monotone if

(12)
$\displaystyle\left(L(B_{1},N_{1})-L(B_{2},N_{2}),(B_{1},N_{1})-(B_{2},N_{2})\right)>0,\;\;\forall(B_{1},N_{1})\neq(B_{2},N_{2})\;\mbox{in}\;\;\mathbf{H}.$

Also, $L$ is coercive if $\mbox{for some }\;\;(\psi_{0},\chi_{0})\in
K_{h}\times V_{h}$

(13)
$\displaystyle\frac{\left(L(B,N),(B,N)-(\psi_{0},\chi_{0})\right)}{\|(B,N)\|_{\mathbf{H}}}\longrightarrow+\infty\;\;\mbox{as}\;\;\|(B,N)\|_{\mathbf{H}}\rightarrow\infty.$

To show (12), we rewrite

(14)
$\left(L(B_{1},N_{1})-L(B_{2},N_{2}),(B_{1},N_{1})-(B_{2},N_{2})\right)=\|B_{1}-B_{2}\|_{0}^{2}\\\
+\Delta
t\left(D_{B}\nabla(B_{1}-B_{2}),\nabla(B_{1}-B_{2})\right)+\|N_{1}-N_{2}\|_{0}^{2}\\\
+\Delta t\left(D_{N}\nabla(N_{1}-N_{2}),\nabla(N_{1}-N_{2})\right)-\Delta
t\left(F(B_{1},N_{1})-F(B_{2},N_{2}),B_{1}-B_{2}\right)\\\ -\Delta
t\left(G(B_{1},N_{1})-G(B_{2},N_{2}),N_{1}-N_{2}\right).$

Using Assumption 1 part B, Cauchy-Schwarz inequality and the inequality (7),
the fifth term can be bounded as

(15) $-\Delta t(F(B_{1},N_{1})-F(B_{2},N_{2}),B_{1}-B_{2})\geq-\Delta
t(|F(B_{1},N_{1})-F(B_{2},N_{1})|,|B_{1}-B_{2}|)\\\ -\Delta
t(|F(B_{2},N_{1})-F(B_{2},N_{2})|,|B_{1}-B_{2}|)\geq-\frac{3\Delta
t}{2}M\|B_{1}-B_{2}\|_{0}^{2}-\frac{\Delta t}{2}M\|N_{1}-N_{2}\|_{0}^{2}.$

Similarly, the last term is estimated by

(16) $-\Delta t(G(B_{1},N_{1})-G(B_{2},N_{2}),N_{1}-N_{2})\geq-\Delta
t(|G(B_{1},N_{1})-G(B_{2},N_{1})|,|N_{1}-N_{2}|)\\\ -\Delta
t(|G(B_{2},N_{1})-G(B_{2},N_{2})|,|N_{1}-N_{2}|)\geq-\frac{3\Delta
t}{2}M\|N_{1}-N_{2}\|_{0}^{2}-\frac{\Delta t}{2}M\|B_{1}-B_{2}\|_{0}^{2}.$

Combining (14)–(16), and using Assumption 1 part A and (8), we finally obtain

(17)
$\displaystyle\left(L(B_{1},N_{1})-L(B_{2},N_{2}),(B_{1}-B_{2},N_{1}-N_{2})\right)\geq
A\left(\|B_{1}-B_{2}\|_{0}^{2}+\|N_{1}-N_{2}\|_{0}^{2}\right),$

with $A=(C_{PF}^{2}-2\Delta tMC_{PF}^{2}+\Delta t\gamma)$. This estimate
implies (12) as long as $A>0$. In turn, this is guaranteed if either (i) or
(ii) holds. In fact (i) requires that the diffusivities are large enough. If
this is not the case, (ii) requires a small enough $\Delta t$.

To show (13), set $(\psi_{0},\chi_{0})=(0,0)$, and apply Assumption 1 parts
A–B, and use inequality (8). Then for all $(B,N)\in\mathbf{H}$, we have

$\left(L(B,N),(B,N)-(\psi_{0},\chi_{0})\right)=\|B\|_{0}^{2}-\Delta
t(F(B,N),B)+\Delta t(D_{B}\nabla B,\nabla B)\\\ +\|N\|_{0}^{2}-\Delta
t(G(B,N),N)+\Delta t(D_{N}\nabla N,\nabla N)\\\
\geq\|B\|_{0}^{2}+\|N\|_{0}^{2}-\Delta tM(\|B\|_{0}^{2}+\|N\|_{0}^{2})+\Delta
t\gamma(\|B\|_{1}^{2}+\|N\|_{1}^{2})\\\ =(1-\Delta
tM)(\|B\|_{0}^{2}+\|N\|_{0}^{2})+\Delta t\gamma(\|B\|_{1}^{2}+\|N\|_{1}^{2}).$

If $\Delta t<\frac{1}{M}$, then we are done. Otherwise, apply (8), and we have

$\left(L(B,N),(B,N)-(\psi_{0},\chi_{0})\right)\geq(C_{PF}^{2}+\Delta t(\gamma-
MC_{PF}^{2}))\|(B,N)\|_{\mathbf{H}}^{2}\geq A\|(B,N)\|_{\mathbf{H}}^{2},$

where $A$ as above. Thus, (13) is satisfied if $\gamma>2MC_{PF}^{2}$, i.e.
condition (i) holds or if $\Delta t<\frac{C_{PF}^{2}}{2MC_{PF}^{2}-\gamma}$,
which is condition (ii), and the proof is complete. ∎

###  3.3. Assumptions on free boundary

The subsequent proof of convergence of the scheme (9) makes an assumption
about the behavior of the free boundary in time. This assumption seems
entirely reasonable for biofilm growth since we expect $\Omega^{*}$ to grow in
time and its boundary to “move forward” rather than “oscillate”.

Following [14] we define

(18) $\displaystyle D_{n}=\bigcup_{t\in
J_{n}}\Omega^{-}(t)\cup\Omega^{-}(t_{n+1})\backslash\overline{\Omega^{-}(t)\cap\Omega^{-}(t_{n+1})};\;\;n=0,\ldots
N_{T}-1,$

and assume as in ([14], p. 601, condition 2.3) that there exists a constant
$\delta$ such that

(19) $\sum_{n=0}^{N_{T}-1}m(D_{n})\leq\delta,\;\;$

where for a measurable $D\subset\mathbb{R}^{2}$, $m(D)$ denotes its Lebgesgue
measure.

We provide two examples to illustrate this assumption.

######  Example 3.1.

Consider $\Omega=(0,1)^{2}$, and assume that the interface
$\Gamma(t)=\partial\Omega^{-}(t)\cap\partial\Omega^{*}(t)$ is a set of points
$(x^{*}(t),y)$ with $x^{*}(t)$ decreasing in time starting from some
$x^{*}(0)$, and with $y\in(0,1)$. Then
$D_{n}=(x^{*}(t_{n+1}),x^{*}(t_{n}))\times(0,1)$. As $N_{T}\uparrow$ we find
that $\sum_{n}m(D_{n})=m\left((x^{*}(t_{n+1}),x^{*}(0))\right)$ is increasing
with $n$, but is bounded by $m(\Omega^{-}(0))$, thus (19) holds.

More generally, (19) asserts that $B(x,t)$ does not change too frequently from
reaching $B^{*}$ to lying strictly below $B^{*}$ or vice versa. In particular,
it avoids the following scenario.

######  Example 3.2.

Consider $\Omega$ as in Ex. 3.1, and the scenario for some partition
$\Upsilon$ of $J$ in which for every $n$, $D_{n}=(a,b)\times(0,1)$ with some
fixed $0<a<b<0$. This scenario corresponds to the free boundary $\Gamma(t)$
oscillating between $a$ and $b$. Then with $N_{T}\uparrow$ and for finer
partitions of $\Upsilon$ we find that $\sum_{n}m(D_{n})=N_{T}m((a,b))$ is
unbounded.

###  3.4. Error Estimate

In this section we prove an error estimate for the approximation of solutions
to (5). We follow the strategy of Johnson [14] which we adapt for the product
space $V_{h}\times V_{h}$ and our coupled system. The main challenge is to
handle the consistency error arising due to the nonlinearity of $F,G$ and the
coupled nature of the system.

We first state the main result. Next we state some auxiliary technical
results, and proceed with the proof. Throughout, $C$ denotes a generic
positive constant not depending on $\Delta t$ and $h$. We also define

$\displaystyle
e_{B}^{n}=B^{n}-B_{h}^{n},e_{N}^{n}=N^{n}-N_{h}^{n},n=0,\ldots,N_{T}.$

We first state an auxiliary result which follows from (19).

######  Lemma 3.2 (Result from [14], page 605; see also [1] for elaborated
details).

Let $\psi\in L^{\infty}(W^{2,p}),\;\;1\leq p<\infty$. and assume (19) holds.
Then $\forall$ $\epsilon>0$, and for $m=0,\ldots,N_{T}$

$\sum_{n=0}^{m-1}\Delta
t\int_{D_{n}}|\psi(x,t_{n+1})|\;dx\leq\frac{\epsilon}{2}\max_{n}\|\psi\|_{0}^{2}+\frac{\epsilon}{2}\sum_{n=0}^{m-1}\|\psi\|_{1}^{2}\Delta
t+C(\log\Delta t^{-1})^{1/2}\Delta t^{3/2}.$

######  Theorem 3.3.

Assume that Assumption 1 and condition (19) hold, and that $\Delta
t<\frac{1}{4M}$. Then there exists a constant $C$ independent of $\Delta t$
and $h$ such that

(20)
$\max_{n}\left(\|e_{B}^{n}\|_{0}+\|e_{N}^{n}\|_{0}\right)+\sqrt{\sum_{n}\left(\|e_{B}^{n}\|_{1}^{2}+\|e_{N}^{n}\|_{1}^{2}\right)\Delta
t}\leq C[(\log(\Delta t)^{-1})^{1/4}\Delta t^{3/4}+h].$

####  3.4.1. Auxiliary results adapted from [14]

In the proof we will use auxiliary technical results following ([14], page
602). In particular, we use the approximation operator $I_{h}$ constructed
therein which applies to functions not necessarily defined pointwise. One
smooths them out first, then interpolates. We only need formal properties of
the operator $I_{h}$, ([14], page 602) which we restate here.

######  Definition 3.1.

( [14], Lemmas 1 and 2). For each $h>0$, let
$I_{h}\mathrel{\mathop{\ordinarycolon}}H_{0}^{1}(\Omega)\rightarrow V_{h}$
denote a linear operator with the following properties:

  1. (i)

$\|\psi-I_{h}\psi\|_{j}\leq Ch^{k-j}\|\psi\|_{k},\;j=0,1,\;k=1,2,$

  2. (ii)

$I_{h}\psi\in K_{h}$ if $\psi\in K.$

This operator and approximation properties are needed, e.g., to handle terms
involving $\frac{\partial B}{\partial t}$. We mention that a similar goal of
approximating non-smooth functions is applied in the context of a-posteriori
error estimates; see, e.g. ([12] page 71).

Now we define $\eta(t)=B(\cdot,t)-I_{h}B(\cdot,t)$,
$\xi(t)=N(\cdot,t)-I_{h}N(\cdot,t)$ for $t\in J$. The operator $I_{h}$
commutes with the time differentiation, thus the following approximation
properties can be stated, as adapted from ([14], page 603).

######  Lemma 3.4.

([14], page 603)

  1. (i)

$\max_{n}\|\eta^{n}\|_{0}\leq Ch\|B\|_{L^{\infty}(J;H^{1}(\Omega))},$ and
$\max_{n}\|\xi^{n}\|_{0}\leq Ch\|N\|_{L^{\infty}(J;H^{1}(\Omega))}.$

  2. (ii)

$\left\|\frac{\partial\eta}{\partial
t}\right\|_{L^{2}(J_{n};L^{2}(\Omega))}\leq Ch\left\|\frac{\partial
B}{\partial t}\right\|_{L^{2}(J_{n};H^{1}(\Omega))},$ and
$\left\|\frac{\partial\xi}{\partial
t}\right\|_{L^{2}(J_{n};L^{2}(\Omega))}\leq Ch\left\|\frac{\partial
N}{\partial t}\right\|_{L^{2}(J_{n};H^{1}(\Omega))}.$

  3. (iii)

$\|\partial\eta^{n}\|_{0}\leq C(\Delta t)^{-1/2}h\left\|\frac{\partial
B}{\partial t}\right\|_{L^{2}(J_{n};H^{1}(\Omega))},$ and
$\|\partial\xi^{n}\|_{0}\leq C(\Delta t)^{-1/2}h\left\|\frac{\partial
N}{\partial t}\right\|_{L^{2}(J_{n};H^{1}(\Omega))}.$

Now we proceed to the proof of Theorem 3.3.

###### Proof.

We begin as usual deriving the error equations. Using Assumption 1 part A, we
have for $n=0,\ldots,N_{T}-1$ for the solutions of (9) that

(21) $(\partial e_{B}^{n},e_{B}^{n+1})+(\partial
e_{N}^{n},e_{N}^{n+1})+\mu_{1}(\nabla e_{B}^{n+1},\nabla
e_{B}^{n+1})+\nu_{1}(\nabla e_{N}^{n+1},\nabla e_{N}^{n+1})\\\ \leq(\partial
e_{B}^{n},e_{B}^{n+1})+(\partial e_{N}^{n},e_{N}^{n+1})+(D_{B}(x)\nabla
e_{B}^{n+1},\nabla e_{B}^{n+1})+(D_{N}(x)\nabla e_{N}^{n+1},\nabla
e_{N}^{n+1})\\\ =(\partial e_{B}^{n},\eta^{n+1})+(\partial
e_{N}^{n},\xi^{n+1})+\mu_{2}(\nabla
e_{B}^{n+1},\nabla\eta^{n+1})+\nu_{2}(\nabla e_{N}^{n+1},\nabla\xi^{n+1})\\\
+(\partial B^{n},I_{h}B^{n+1}-B_{h}^{n+1})-(\partial
B_{h}^{n},I_{h}B^{n+1}-B_{h}^{n+1})\\\ +(\partial
N^{n},I_{h}N^{n+1}-N_{h}^{n+1})-(\partial
N_{h}^{n},I_{h}N^{n+1}-N_{h}^{n+1})\\\ +(D_{B}(x)\nabla
B^{n+1},\nabla(I_{h}B^{n+1}-B_{h}^{n+1}))-(D_{B}(x)\nabla
B_{h}^{n+1},\nabla(I_{h}B^{n+1}-B_{h}^{n+1}))\\\ +(D_{N}(x)\nabla
N^{n+1},\nabla(I_{h}N^{n+1}-N_{h}^{n+1}))-(D_{N}(x)\nabla
N_{h}^{n+1},\nabla(I_{h}N^{n+1}-N_{h}^{n+1})).$

Taking $t=t_{n+1}$, $\psi=B_{h}^{n+1}$, $\chi=N_{h}^{n+1}-N^{n+1}$ in (5), and
$\psi=I_{h}B^{n+1}$, $\chi=I_{h}N^{n+1}-N_{h}^{n+1}$ in (9), and adding the
obtained inequalities and equations to (21), we have

(22) $(\partial e_{B}^{n},e_{B}^{n+1})+(\partial
e_{N}^{n},e_{N}^{n+1})+\mu_{1}\|e_{B}^{n+1}\|_{1}^{2}+\nu_{1}\|e_{N}^{n+1}\|_{1}^{2}\leq\sum_{j=1}^{10}P_{j}^{n},$

where

(23) $\displaystyle P_{1}^{n}$ $\displaystyle=$ $\displaystyle(\partial
e_{B}^{n},\eta^{n+1}),$ (24) $\displaystyle P_{2}^{n}$ $\displaystyle=$
$\displaystyle(\partial e_{N}^{n},\xi^{n+1}),$ (25) $\displaystyle P_{3}^{n}$
$\displaystyle=$ $\displaystyle\mu_{2}(\nabla e_{B}^{n+1},\nabla\eta^{n+1}),$
(26) $\displaystyle P_{4}^{n}$ $\displaystyle=$ $\displaystyle\nu_{2}(\nabla
e_{N}^{n+1},\nabla\xi^{n+1}),$ (27) $\displaystyle P_{5}^{n}$ $\displaystyle=$
$\displaystyle-(D_{B}(x)\nabla B^{n+1},\nabla\eta^{n+1})-(\partial
B^{n},\eta^{n+1})+(f^{n+1},\eta^{n+1}),$ (28) $\displaystyle P_{6}^{n}$
$\displaystyle=$ $\displaystyle-(D_{N}(x)\nabla
N^{n+1},\nabla\xi^{n+1})-(\partial N^{n},\xi^{n+1})+(g^{n+1},\xi^{n+1}),$ (29)
$\displaystyle P_{7}^{n}$ $\displaystyle=$ $\displaystyle(\partial
B^{n}-\frac{\partial B}{\partial t}(t_{n+1}),e_{B}^{n+1}),$ (30)
$\displaystyle P_{8}^{n}$ $\displaystyle=$ $\displaystyle(\partial
N^{n}-\frac{\partial N}{\partial t}(t_{n+1}),e_{N}^{n+1}),$ (31)
$\displaystyle P_{9}^{n}$ $\displaystyle=$
$\displaystyle(F(B^{n+1},N^{n+1}),e_{B}^{n+1})-(F(B_{h}^{n+1},N_{h}^{n+1}),I_{h}B^{n+1}-B_{h}^{n+1}),$
(32) $\displaystyle P_{10}^{n}$ $\displaystyle=$
$\displaystyle(G(B^{n+1},N^{n+1}),e_{N}^{n+1})-(G(B_{h}^{n+1},N_{h}^{n+1}),I_{h}N^{n+1}-N_{h}^{n+1}).$

Multiplying (22) by $\Delta t$ and summing over $n=0,\ldots,m-1$;
$m=1,\ldots,N_{T}$, we obtain

(33)
$\sum_{n=0}^{m-1}(e_{B}^{n+1}-e_{B}^{n},e_{B}^{n+1})+\sum_{n=0}^{m-1}(e_{N}^{n+1}-e_{N}^{n},e_{N}^{n+1})\\\
+\mu_{1}\sum_{n=0}^{m-1}\|e_{B}^{n+1}\|_{1}^{2}\Delta
t+\nu_{1}\sum_{n=0}^{m-1}\|e_{N}^{n+1}\|_{1}^{2}\Delta
t\leq\sum_{j=1}^{10}S_{j};$

Here we have defined $S_{j}=\sum_{n=0}^{m-1}|P_{j}^{n}|\Delta t$, for
$j=1,\ldots,10$.

We also have

(34) $\displaystyle 2\sum_{n=0}^{m-1}(e_{B}^{n+1}-e_{B}^{n},e_{B}^{n+1})$
$\displaystyle=$
$\displaystyle\sum_{n=0}^{m-1}\|e_{B}^{n+1}-e_{B}^{n}\|_{0}^{2}+\|e_{B}^{m}\|_{0}^{2}-\|e_{B}^{0}\|_{0}^{2},$

and

(35) $\displaystyle 2\sum_{n=0}^{m-1}(e_{N}^{n+1}-e_{N}^{n},e_{N}^{n+1})$
$\displaystyle=$
$\displaystyle\sum_{n=0}^{m-1}\|e_{N}^{n+1}-e_{N}^{n}\|_{0}^{2}+\|e_{N}^{m}\|_{0}^{2}.$

Multiplying (33) by 2 and using (34)–(35) give

(36)
$\sum_{n=0}^{m-1}\|e_{B}^{n+1}-e_{B}^{n}\|_{0}^{2}+\|e_{B}^{m}\|_{0}^{2}+\sum_{n=0}^{m-1}\|e_{N}^{n+1}-e_{N}^{n}\|_{0}^{2}+\|e_{N}^{m}\|_{0}^{2}\\\
+2\mu_{1}\sum_{n=0}^{m-1}\|e_{B}^{n+1}\|_{1}^{2}\Delta
t+2\nu_{1}\sum_{n=0}^{m-1}\|e_{N}^{n+1}\|_{1}^{2}\Delta
t\leq\|e_{B}^{0}\|_{0}^{2}+2\sum_{j=1}^{10}S_{j}.$

We shall estimate each one of the $S_{j}$’s. Many of these estimates are
direct analogues of the estimates in [14]. Other estimates handle the
consistency and coupling terms.

Estimation of $S_{1}$, and $S_{2}$, analogously as in [14]. Applying summation
by parts, Lemma 3.4, and inequalities (8), (7), we obtain

$\displaystyle S_{1}$ $\displaystyle\leq$ $\displaystyle\frac{\epsilon
C_{PF}^{2}}{2}\sum_{n=0}^{m-1}\|e_{B}^{n+1}\|_{1}^{2}\Delta
t+Ch^{2}\left\|\frac{\partial B}{\partial
t}\right\|_{L^{2}(J;H^{1}(\Omega))}^{2}+\frac{\epsilon}{2}\|e_{B}^{m}\|_{0}^{2}+\epsilon\|e_{B}^{0}\|_{0}^{2}+Ch^{2}\|B\|_{L^{\infty}(J;H^{1}(\Omega))}^{2}.$

By choosing an appropriate $\epsilon$, we have

(37) $S_{1}\leq\frac{\mu_{1}}{14}\sum_{n=0}^{m-1}\|e_{B}^{n+1}\|_{1}^{2}\Delta
t+\frac{1}{8}\|e_{B}^{m}\|_{0}^{2}+\|e_{B}^{0}\|_{0}^{2}+Ch^{2}.$

Similarly,

(38) $S_{2}\leq\frac{\nu_{1}}{10}\sum_{n=0}^{m-1}\|e_{N}^{n+1}\|_{1}^{2}\Delta
t+\frac{1}{4}\|e_{N}^{m}\|_{0}^{2}+Ch^{2}.$

Estimation of $S_{3}$, and $S_{4}$, analogously as in [14].

$S_{3}\leq\frac{\mu_{2}\epsilon}{2}\sum_{n=0}^{m-1}\|e_{B}^{n+1}\|_{1}^{2}\Delta
t+Ch^{2}\|B\|_{L^{\infty}(J;H^{2}(\Omega))}^{2},$

where we use Lemma 3.4. In particular, we have

(39) $S_{3}\leq\frac{\mu_{1}}{14}\sum_{n=0}^{m-1}\|e_{B}^{n+1}\|_{1}^{2}\Delta
t+Ch^{2}.$

Similarly,

(40) $S_{4}\leq\frac{\nu_{1}}{10}\sum_{n=0}^{m-1}\|e_{N}^{n+1}\|_{1}^{2}\Delta
t+Ch^{2}.$

Estimation of $S_{5}$ and $S_{6}$. These estimations are different from those
in [14], since they apply to nonconstant diffusivities $D_{B}$ and $D_{N}$. We
use Green’s theorem to rewrite the first terms of $P_{5}^{n}$ and $P_{6}^{n}$
as $(\nabla D_{B}(x)\cdot\nabla B^{n+1},\eta^{n+1})+(D_{B}(x)\Delta
B^{n+1},\eta^{n+1})$ and $(\nabla D_{N}(x)\cdot\nabla
N^{n+1},\xi^{n+1})+(D_{N}(x)\Delta N^{n+1},\xi^{n+1})$ respectively. Then we
apply Assumption 1 parts A and C, and we use Lemma 3.4 to obtain

(41) $\displaystyle
S_{5}\leq\sum_{n=0}^{m-1}\left(\|B^{n+1}\|_{1}+\|B^{n+1}\|_{2}+\|\partial
B^{n}\|_{0}+\|f^{n+1}\|_{0}\right)\|\eta^{n+1}\|_{0}\Delta t\leq Ch^{2}.$

Similarly,

(42) $\displaystyle S_{6}$ $\displaystyle\leq$
$\displaystyle\sum_{n=0}^{m-1}\left(\|N^{n+1}\|_{1}+\|N^{n+1}\|_{2}+\|\partial
N^{n}\|_{0}+\|g^{n+1}\|_{0}\right)\|\xi^{n+1}\|_{0}\Delta t\leq Ch^{2}.$

Estimation of $S_{7}$. In estimating $S_{7}$, we follow Johnson’s technique
[14]. However, we need to deal with coupling term. We can write $P_{7}^{n}$ as

$P_{7}^{n}=(\partial B^{n}-\frac{\partial B}{\partial
t}(t_{n+1}),e_{B}^{n+1})\\\ =\frac{1}{\Delta
t}\int_{\Omega}(B^{n+1}-B^{n}-\Delta t\frac{\partial B}{\partial
t}(t_{n+1}))e_{B}^{n+1}\;dx\\\ =\frac{1}{\Delta
t}\int_{\Omega}\left(\int_{t_{n}}^{t_{n+1}}\left(\frac{\partial B}{\partial
t}(s)-\frac{\partial B}{\partial
t}(t_{n+1})\right)\;ds\right)e_{B}^{n+1}\;dx\\\ =\frac{1}{\Delta
t}\int_{t_{n}}^{t_{n+1}}\left(\int_{\Omega^{-}}\left(\nabla\cdot(D_{B}\nabla
B(s))-\nabla\cdot(D_{B}\nabla B(t_{n+1}))\right)e_{B}^{n+1}\;dx\right)\;ds\\\
+\frac{1}{\Delta
t}\int_{t_{n}}^{t_{n+1}}\left(\int_{\Omega}\left(\widetilde{F}(B(s),N(s))-\widetilde{F}(B^{n+1},N^{n+1})\right)e_{B}^{n+1}\;dx\right)\;ds=q_{n}^{1}+q_{n}^{2},$

where

$\widetilde{F}(B(t),N(t))=\left\\{\begin{array}[]{ll}F(B(t),N(t))+f(x,t)&\mbox{if
$x\in\Omega^{-}(t)$},\\\ \min(F(B,N)+f(x,t),0)&\mbox{ if
$x\in\Omega^{\ast}(t)$,}\end{array}\right.$

with obvious notation for $q_{n}^{1}$ and $q_{n}^{2}$. Since
$\nabla\cdot(D_{B}\nabla B)=0$ a.e on $\Omega^{\ast}$, by using Green’s
theorem we obtain

$q_{n}^{1}=\frac{1}{\Delta
t}\int_{t_{n}}^{t_{n+1}}\left(\int_{\Omega}\nabla\cdot(D_{B}\nabla\left(B(s)-B(t_{n+1})\right)e_{B}^{n+1}\;dx\right)\;ds\\\
=\frac{1}{\Delta
t}\int_{t_{n}}^{t_{n+1}}\int_{s}^{t_{n+1}}\left(\int_{\Omega}D_{B}\nabla\frac{\partial
B}{\partial t}(t)\cdot\nabla e_{B}^{n+1}\;dx\right)\;dt\;ds.$

From the Cauchy-Schwarz inequality, we obtain

$|q_{n}^{1}|\leq\frac{\mu_{2}}{\Delta
t}\int_{t_{n}}^{t_{n+1}}\int_{s}^{t_{n+1}}\left\|\frac{\partial B}{\partial
t}(t)\right\|_{1}\|e_{B}^{n+1}\|_{1}\;dt\;ds\\\ \leq\mu_{2}(\Delta
t)^{1/2}\|e_{B}^{n+1}\|_{1}\left\|\frac{\partial B}{\partial
t}\right\|_{L^{2}(J_{n};H^{1}(\Omega))}\leq\frac{\epsilon}{2}\|e_{B}^{n+1}\|_{1}^{2}+\Delta
t\frac{\mu_{2}^{2}}{2\epsilon}\left\|\frac{\partial B}{\partial
t}\right\|_{L^{2}(J_{n};H^{1}(\Omega))}^{2}.$

Thus, $\sum_{n=0}^{m-1}|q_{n}^{1}|\Delta
t\leq\frac{\epsilon}{2}\sum_{n=0}^{m-1}\|e_{B}^{n+1}\|_{1}^{2}\Delta
t+C(\Delta t)^{2}\|\frac{\partial B}{\partial
t}\|_{L^{2}(J;H^{1}(\Omega))}^{2}.$ In particular,
$\sum_{n=0}^{m-1}|q_{n}^{1}|\Delta
t\leq\frac{\mu_{1}}{14}\sum_{n=0}^{m-1}\|e_{B}^{n+1}\|_{1}^{2}\Delta
t+C(\Delta t)^{2}.$ Now, to estimate $q_{n}^{2}$, we first notice that if
$x\in\Omega\backslash D_{n}$, then either
$x\in\Omega^{-}(t)\cap\Omega^{-}(t_{n+1})$ or
$x\in\Omega^{\ast}(t)\cap\Omega^{\ast}(t_{n+1})$ for all $t\in J_{n}$, so we
have

$|\widetilde{F}(B(t),N(t))-\widetilde{F}(B^{n+1},N^{n+1})|\leq|F(B(t),N(t))-F(B^{n+1},N^{n+1})|\\\
+|f(x,t)-f(x,t^{n+1})|\;\;\mbox{for all}\;x\in\Omega\backslash D_{n}.$

Thus,

$|q_{n}^{2}|\leq\frac{1}{\Delta t}\int_{t_{n}}^{t_{n+1}}\int_{\Omega\backslash
D_{n}}|F(B(s),N(s))-F(B^{n+1},N^{n+1})||e_{B}^{n+1}|\;dx\;\;ds\\\
+\frac{1}{\Delta t}\int_{t_{n}}^{t_{n+1}}\int_{\Omega\backslash
D_{n}}|f(x,s)-f(x,s)||e_{B}^{n+1}|\;dx\;\;ds\\\ +\frac{1}{\Delta
t}\int_{t_{n}}^{t_{n+1}}\int_{D_{n}}(\|F\|_{L^{\infty}(\mathbb{R}^{2})}+\|f\|_{L^{\infty}(J_{n};L^{\infty}(D_{n}))})|e_{B}^{n+1}|\;dx\;\;ds=k_{n}^{1}+k_{n}^{2}+k_{n}^{3}.$

Next we estimate the terms $k_{n}^{1}+k_{n}^{2}+k_{n}^{3}$ separately.

The first term $k_{n}^{1}$ contains the coupling and nonlinearity in $F$ which
are not present in [14], To handle that, we use Assumption 1 part B, and we
have

(43) $k_{n}^{1}\leq M\frac{1}{\Delta
t}\int_{t_{n}}^{t_{n+1}}\int_{\Omega\backslash
D_{n}}|B(s)-B(t_{n+1})||e_{B}^{n+1}|\;dx\;\;ds\\\ +M\frac{1}{\Delta
t}\int_{t_{n}}^{t_{n+1}}\int_{\Omega\backslash
D_{n}}|N(s)-N(t_{n+1})||e_{B}^{n+1}|\;dx\;\;ds\\\ \leq M\frac{1}{\Delta
t}\int_{t_{n}}^{t_{n+1}}\int_{\Omega\backslash
D_{n}}\int_{s}^{t_{n+1}}\left|\frac{\partial B}{\partial
t}(t)\right||e_{B}^{n+1}|\;dt\;dx\;ds+M\frac{1}{\Delta
t}\int_{t_{n}}^{t_{n+1}}\int_{\Omega\backslash
D_{n}}\int_{s}^{t_{n+1}}\left|\frac{\partial N}{\partial
t}(t)\right||e_{B}^{n+1}|\;dt\;dx\;ds\\\ \leq M\frac{1}{\Delta
t}\int_{t_{n}}^{t_{n+1}}\int_{s}^{t_{n+1}}\left\|\frac{\partial B}{\partial
t}(t)\right\|_{0}\|e_{B}^{n+1}\|_{0}\;dt\;ds+M\frac{1}{\Delta
t}\int_{t_{n}}^{t_{n+1}}\int_{s}^{t_{n+1}}\left\|\frac{\partial N}{\partial
t}(t)\right\|_{0}\|e_{B}^{n+1}\|_{0}\;dt\;ds\\\ \leq M(\Delta
t)^{1/2}\|e_{B}^{n+1}\|_{0}\left(\left\|\frac{\partial B}{\partial
t}\right\|_{L^{2}(J_{n};L^{2}(\Omega))}+\left\|\frac{\partial N}{\partial
t}\right\|_{L^{2}(J_{n};L^{2}(\Omega))}\right)\\\
\leq\frac{\epsilon}{2}C_{PF}^{2}\|e_{B}^{n+1}\|_{1}^{2}+\frac{M^{2}}{\epsilon}\Delta
t\left(\left\|\frac{\partial B}{\partial
t}\right\|_{L^{2}(J_{n};L^{2}(\Omega))}^{2}+\left\|\frac{\partial N}{\partial
t}\right\|_{L^{2}(J_{n};L^{2}(\Omega))}^{2}\right).$

Next

(44) $k_{n}^{2}=\frac{1}{\Delta t}\int_{t_{n}}^{t_{n+1}}\int_{\Omega\backslash
D_{n}}|f(x,s)-f(x,s)||e_{B}^{n+1}|\;dx\;\;ds\\\ \leq\frac{1}{\Delta
t}\int_{t_{n}}^{t_{n+1}}\int_{\Omega\backslash
D_{n}}\int_{s}^{t_{n+1}}\left|\frac{\partial f}{\partial
t}(x,t)\right||e_{B}^{n+1}|\;dt\;dx\;\;ds\leq\frac{1}{\Delta
t}\|e_{B}^{n+1}\|_{0}\int_{t_{n}}^{t_{n+1}}\int_{s}^{t_{n+1}}\left\|\frac{\partial
f}{\partial t}(x,t)\right\|_{0}\;dt\;\;ds\\\ \leq(\Delta
t)^{1/2}\|e_{B}^{n+1}\|_{0}\left\|\frac{\partial f}{\partial
t}\right\|_{L^{2}(J_{n};L^{2}(\Omega))}\leq\frac{\epsilon}{2}C_{PF}^{2}\|e_{B}^{n+1}\|_{1}^{2}+\frac{1}{2\epsilon}\Delta
t\left\|\frac{\partial f}{\partial
t}\right\|_{L^{2}(J_{n};L^{2}(\Omega))}^{2}.$

By (43) and (44), we have

(45) $\sum_{n=0}^{m-1}|k_{n}^{1}|\Delta t+\sum_{n=0}^{m-1}|k_{n}^{2}|\Delta
t\leq\frac{\mu_{1}}{14}\sum_{n=0}^{m-1}\|e_{B}^{n+1}\|_{1}^{2}\Delta
t+C(\Delta t)^{2}.$

Finally

$\displaystyle k_{n}^{3}$ $\displaystyle=$
$\displaystyle(\|F\|_{L^{\infty}(\mathbb{R}^{2})}+\|f\|_{L^{\infty}(J_{n};L^{\infty}(D_{n}))})\int_{D_{n}}|e_{B}^{n+1}|\;dx=(\|F\|_{L^{\infty}(\mathbb{R}^{2})}+\|f\|_{L^{\infty}(J_{n};L^{\infty}(D_{n}))})r_{n},$

By Lemma 3.2, $r_{n}$ can be estimated as

$\sum_{n=0}^{m-1}r_{n}\Delta
t\leq\frac{\epsilon}{2}\max_{n}\|e_{B}^{n+1}\|_{0}^{2}+\frac{\epsilon}{2}\sum_{n=0}^{m-1}\|e_{B}^{n+1}\|_{1}^{2}\Delta
t+C(\log\Delta t^{-1})^{1/2}\Delta t^{3/2}.$

Since at the end we want to kick the term $\max_{n}\|e_{B}^{n+1}\|_{0}^{2}$ to
the left hand side of inequality (36), we estimate as follows.

Let $\|e_{B}^{l+1}\|_{0}=\max_{n}\|e_{B}^{n+1}\|_{0}$ for some
$n,l\in\\{0,\ldots,m-1\\}$. Then

$\|e_{B}^{l+1}\|_{0}^{2}=\sum_{j=l+1}^{m-1}\|e_{B}^{j+1}-e_{B}^{j}\|_{0}^{2}+\|e_{B}^{m}\|_{0}^{2}\leq\sum_{n=0}^{m-1}\|e_{B}^{n+1}-e_{B}^{n}\|_{0}^{2}+\|e_{B}^{m}\|_{0}^{2}.$

Hence, we get

$\displaystyle\sum_{n=0}^{m-1}r_{n}\Delta t$ $\displaystyle\leq$
$\displaystyle\frac{\epsilon}{2}\left(\sum_{n=0}^{m-1}\|e_{B}^{n+1}-e_{B}^{n}\|_{0}^{2}+\|e_{B}^{m}\|_{0}^{2}\right)+\frac{\epsilon}{2}\sum_{n=0}^{m-1}\|e_{B}^{n+1}\|_{1}^{2}\Delta
t+C(\log\Delta t^{-1})^{1/2}\Delta t^{3/2},$

which yields

(46) $\sum_{n=0}^{m-1}k_{n}^{3}\Delta
t\leq\frac{1}{4}\sum_{n=0}^{m-1}\|e_{B}^{n+1}-e_{B}^{n}\|_{0}^{2}+\frac{1}{8}\|e_{B}^{m}\|_{0}^{2}\\\
+\frac{\mu_{1}}{14}\sum_{n=0}^{m-1}\|e_{B}^{n+1}\|_{1}^{2}\Delta
t+C(\log\Delta t^{-1})^{1/2}\Delta t^{3/2}.$

Estimation of $S_{8}$. Using Cauchy-Schwarz inequality, we have

$S_{8}\leq\frac{\epsilon}{2}\sum_{n=0}^{m-1}\|e_{N}^{n+1}\|_{1}^{2}\Delta
t+\frac{1}{2\epsilon}\sum_{n=0}^{m-1}\left\|\partial N^{n}-\frac{\partial
N}{\partial t}(t_{n+1})\right\|_{0}^{2}\Delta t.$

By Assumption 1 part F on $N$, we have $\frac{\partial^{2}N}{\partial
t^{2}}\in L^{2}(Q)$, thus the second term of the above inequality can be
estimated as

$\left\|\partial N^{n}-\frac{\partial N}{\partial
t}(t_{n+1})\right\|_{0}^{2}\Delta
t\leq\int_{\Omega}\int_{t_{n}}^{t_{n+1}}\left(\frac{\partial N}{\partial
t}(s)-\frac{\partial N}{\partial t}(t_{n+1})\right)^{2}\;ds\;dx\\\ \leq\Delta
t\int_{\Omega}\int_{t_{n}}^{t_{n+1}}\int_{s}^{t_{n+1}}\left(\frac{\partial^{2}N}{\partial
t^{2}}(t)\right)^{2}\;dt\;ds\;dx=(\Delta
t)^{2}\left\|\frac{\partial^{2}N}{\partial
t^{2}}\right\|_{L^{2}(J_{n};L^{2}(\Omega))},$

which implies $\sum_{n=0}^{m-1}\left\|\partial N^{n}-\frac{\partial
N}{\partial t}(t_{n+1})\right\|_{0}^{2}\Delta t\leq(\Delta
t)^{2}\left\|\frac{\partial^{2}N}{\partial
t^{2}}\right\|_{L^{2}(J;L^{2}(\Omega))}.$ Therefore we have
$S_{8}\leq\frac{\epsilon}{2}\sum_{n=0}^{m-1}\|e_{N}^{n+1}\|_{1}^{2}\Delta
t+C(\Delta t)^{2}.$ In particular,

(47) $S_{8}\leq\frac{\nu_{1}}{10}\sum_{n=0}^{m-1}\|e_{N}^{n+1}\|_{1}^{2}\Delta
t+C(\Delta t)^{2}.$

Estimation of $S_{9}$, and $S_{10}$. These are the consistency terms. By
Assumption 1 part B, $P_{9}^{n}$ can be estimated as

$P_{9}^{n}=(F(B^{n+1},N^{n+1})-F(B^{n+1},N_{h}^{n+1}),e_{B}^{n+1})+(F(B^{n+1},N_{h}^{n+1})-F(B_{h}^{n+1},N_{h}^{n+1}),e_{B}^{n+1})\\\
-(F(B_{h}^{n+1},N_{h}^{n+1}),\eta^{n+1})\\\ \leq
M\|e_{N}^{n+1}\|_{0}\|e_{B}^{n+1}\|_{0}+M\|e_{B}^{n+1}\|_{0}^{2}+\|F\|_{L^{\infty}(\mathbb{R}^{2})}\|\eta^{n+1}\|_{0}\\\
\leq\frac{M^{2}\epsilon}{2}\|e_{N}^{n+1}\|_{0}^{2}+\frac{1}{2\epsilon}\|e_{B}^{n+1}\|_{0}^{2}+M\|e_{B}^{n+1}\|_{0}^{2}+\|F\|_{L^{\infty}(\mathbb{R}^{2})}m(\Omega)^{1/2}\|\eta^{n+1}\|_{0}.$

Hence, we obtain

(48) $S_{9}\leq\frac{\nu_{1}}{10}\sum_{n=0}^{m-1}\|e_{N}^{n+1}\|_{1}^{2}\Delta
t+\frac{\mu_{1}}{14}\sum_{n=0}^{m-1}\|e_{B}^{n+1}\|_{1}^{2}\Delta
t+M\sum_{n=0}^{m-1}\|e_{B}^{n+1}\|_{0}^{2}\Delta t+Ch^{2}.$

Similarly,

(49)
$S_{10}\leq\frac{\nu_{1}}{10}\sum_{n=0}^{m-1}\|e_{N}^{n+1}\|_{1}^{2}\Delta
t+\frac{\mu_{1}}{14}\sum_{n=0}^{m-1}\|e_{B}^{n+1}\|_{1}^{2}\Delta
t+M\sum_{n=0}^{m-1}\|e_{N}^{n+1}\|_{0}^{2}\Delta t+Ch^{2}.$

Now we collect all the above estimates from (36)–(49) and we get

$\sum_{n=0}^{m-1}\|e_{B}^{n+1}-e_{B}^{n}\|_{0}^{2}+\sum_{n=0}^{m-1}\|e_{N}^{n+1}-e_{N}^{n}\|_{0}^{2}+\|e_{B}^{m}\|_{0}^{2}\\\
+\|e_{N}^{m}\|_{0}^{2}+\mu_{1}\sum_{n=0}^{m-1}\|e_{B}^{n+1}\|_{1}^{2}\Delta
t+\nu_{1}\sum_{n=0}^{m-1}\|e_{N}^{n+1}\|_{1}^{2}\Delta t\\\ \leq
2\|e_{B}^{0}\|_{0}^{2}+C\left(h^{2}+(\log\Delta t^{-1})^{1/2}\Delta
t^{3/2}\right)+\sum_{n=0}^{m-1}2M\left(\|e_{B}^{n+1}\|_{0}^{2}+\|e_{N}^{n+1}\|_{0}^{2}\right)\Delta
t.$

Thus we have

$\|e_{B}^{m}\|_{0}^{2}+\|e_{N}^{m}\|_{0}^{2}+\mu_{1}\sum_{j=1}^{m}\|e_{B}^{j}\|_{1}^{2}\Delta
t+\nu_{1}\sum_{j=1}^{m}\|e_{N}^{j}\|_{1}^{2}\Delta t\\\ \leq
2\|e_{B}^{0}\|_{0}^{2}+C\left(h^{2}+(\log\Delta t^{-1})^{1/2}\Delta
t^{3/2}\right)+\sum_{j=1}^{m}2M\left(\|e_{B}^{j}\|_{0}^{2}+\|e_{N}^{j}\|_{0}^{2}\right)\Delta
t.$

By the assumption on the smallness of $\Delta t$, we have $2M\Delta
t\left(\|e_{B}^{m}\|_{0}^{2}+\|e_{N}^{m}\|_{0}^{2}\right)<\frac{1}{2}\left(\|e_{B}^{m}\|_{0}^{2}+\|e_{N}^{m}\|_{0}^{2}\right)$,
and hence we have

(50)
$\|e_{B}^{m}\|_{0}^{2}+\|e_{N}^{m}\|_{0}^{2}+2\mu_{1}\sum_{j=1}^{m}\|e_{B}^{j}\|_{1}^{2}\Delta
t+2\nu_{1}\sum_{j=1}^{m}\|e_{N}^{j}\|_{1}^{2}\Delta t\\\ \leq
C\left(h^{2}+(\log\Delta t^{-1})^{1/2}\Delta
t^{3/2}\right)+\sum_{j=1}^{m-1}4M\left(\|e_{B}^{j}\|_{0}^{2}+\|e_{N}^{j}\|_{0}^{2}\right)\Delta
t.$

We apply the discrete Gronwall’s Lemma ([4], inequality (3.27), page 1114)

(51)
$(r^{0})^{2}+(s^{0})^{2}\leq(p^{0})^{2},\;\;\text{and}\;\;(r^{m})^{2}+(s^{m})^{2}\leq\sum_{j=0}^{m-1}(y^{j})^{2}(r^{j})^{2}+\sum_{j=0}^{m}(p^{j})^{2},\;\;m\geq
1\\\
\Longrightarrow(r^{m})^{2}+(s^{m})^{2}\leq\exp\left(\sum_{j=0}^{m-1}(y^{j})^{2}\right)\sum_{j=0}^{m}(p^{j})^{2},\;\;m\geq
1.$

with $(r^{m})^{2}=\|e_{B}^{m}\|_{0}^{2}+\|e_{N}^{m}\|_{0}^{2}$,
$(s^{m})^{2}=2\mu_{1}\sum_{j=1}^{m}\|e_{B}^{j}\|_{1}^{2}\Delta
t+2\nu_{1}\sum_{j=1}^{m}\|e_{N}^{j}\|_{1}^{2}\Delta t$,
$(p^{j})^{2}=(t_{j}-t_{j-1})C\left(h^{2}+(\log\Delta t^{-1})^{1/2}\Delta
t^{3/2}\right)$, $(y^{j})^{2}=4M\Delta t$ to (50) to get

$\|e_{B}^{m}\|_{0}^{2}+\|e_{N}^{m}\|_{0}^{2}+2\mu_{1}\sum_{j=1}^{m}\|e_{B}^{j}\|_{1}^{2}\Delta
t+2\nu_{1}\sum_{j=1}^{m}\|e_{N}^{j}\|_{1}^{2}\Delta t\\\ \leq
C\left(h^{2}+(\log\Delta t^{-1})^{1/2}\Delta
t^{3/2}\right)\exp(\sum_{j=0}^{m-1}4M\Delta t)\leq C\left(h^{2}+(\log\Delta
t^{-1})^{1/2}\Delta t^{3/2}\right)\exp(4MT).$

Therefore, there exists a constant $C$ independent of $h$ and $\Delta t$ such
that

$\|e_{B}^{m}\|_{0}^{2}+\|e_{N}^{m}\|_{0}^{2}+\sum_{j=1}^{m}\|e_{B}^{j}\|_{1}^{2}\Delta
t+\sum_{j=1}^{m}\|e_{N}^{j}\|_{1}^{2}\Delta t\\\ \leq C\left(h^{2}+(\log\Delta
t^{-1})^{1/2}\Delta t^{3/2}\right),$

and the proof is complete. ∎

######  Remark 3.

The order of convergence is actually close to first order in $h$ if
$h=O(\Delta t)$.

##  4\. Numerical Solver

The problem (9) is a system of nonlinear equations to be solved for
$(B_{h}^{n},N_{h}^{n})$ at each time step $n$. Below we formulate this
algebraic system and discuss the nonlinear solver which has the structure of
Newton-Raphson iterations complemented with the constraint implicit in the
variational inequality (9). We work in the framework of nonlinear
complementarity constraints (NCC), and with a Lagrange multiplier
$\Lambda_{h}^{n}$. Given some $B_{*}<B^{*}$ from
$\mathbb{R}\cup\\{\infty,-\infty\\}$ we define the Evans function

(52) $\displaystyle
P_{[B_{\ast},B^{\ast}]}(\psi)=\max\\{B_{\ast},\min(\psi,B^{\ast})\\}.$

This function is piecewise differentiable. In the model considered in this
paper $B_{*}=-\infty$, but the solver can be applied easily to the general
“double-obstacle” case.

To show how the Newton solver is applied to the solution of PVI, we rewrite
(9) in an equivalent algebraic form. Let
$V_{h}=\mathrm{span}\\{\phi_{i}\\}_{i=1}^{q}$ where $q$ is the number of
degrees of freedom. Define the vectors of degrees of freedom
$\mathbf{B}^{n},\mathbf{N}^{n}$ of $B_{h}^{n},N_{h}^{n}$. Define the matrices
$\mathbf{M},\;\mathbf{R},\;\mathbf{A}_{B},\;\mathbf{A}_{N}$, each with entries
with subscript $(ij)$

$\displaystyle M_{ij}$ $\displaystyle=$
$\displaystyle\int_{\Omega}\phi_{i}\phi_{j},\;\;R_{ij}=\int_{\Omega}P(N^{n}_{h})\phi_{i}\phi_{j},$
$\displaystyle(A_{B})_{ij}$ $\displaystyle=$
$\displaystyle\int_{\Omega}D_{B}({B^{n}_{h}})\nabla\phi_{i}\cdot\nabla\phi_{j},\;\;(A_{N})_{ij}=\int_{\Omega}D_{N}({B^{n}_{h}})\nabla\phi_{i}\cdot\nabla\phi_{j},$

where we applied time-lagging. For $F$ and $G$ we use (2) and define vectors
$\mathbf{F_{B}}$ and $\mathbf{F_{N}}$ with the entries
$(F_{B})_{i}=\int_{\Omega}f(x,t_{n+1})\phi_{i}$,
$(F_{N})_{i}=\int_{\Omega}g(x,t_{n+1})\phi_{i}$. We also seek the vector
$\Lambda^{n+1}$ of pointwise values on the grid.

In summary (9) can be written in the form

(53a) $\displaystyle(\mathbf{M}+\;\Delta t\;\mathbf{A_{B}}-\kappa_{B}\Delta
t\;\mathbf{R}\;)\;\mathbf{B}^{n+1}-\Delta t\;\mathbf{M}\Lambda^{n+1}-\Delta
t\mathbf{F_{B}}-\mathbf{M}\mathbf{B}^{n}$ $\displaystyle=0,$ (53b)
$\displaystyle\mathbf{B}-P_{[B_{\ast},B^{\ast}]}(\mathbf{B}^{n+1}-\Lambda^{n+1})$
$\displaystyle=0;$ (53c) $\displaystyle(\mathbf{M}+\;\Delta
t\;\mathbf{A_{N}})\mathbf{N}^{n+1}+\kappa_{N}\Delta
t\;\mathbf{R}\mathbf{B}^{n+1}-\Delta t\mathbf{F_{N}}-\mathbf{M}\mathbf{N}^{n}$
$\displaystyle=0.$

######  Remark 4.

The equations (53a) and (53c) are the discrete counterparts of the
diffusion–reaction equations at each time step. The middle equation (53b)
enforces pointwise that $B_{\ast}\leq B_{h}^{n+1}\leq B^{\ast}$. The Lagrange
multiplier is in fact a selection $-\Lambda_{h}^{n+1}\in\partial
I_{[B_{\ast},B^{\ast}]}(B_{h}^{n+1})$.

In summary, (53) can be written together with
$\mathcal{F}\mathrel{\mathop{\ordinarycolon}}\mathcal{R}^{3q}\rightarrow\mathcal{R}^{3q}$

$\mathcal{F}(\mathbf{B}^{n+1},\Lambda^{n+1},\mathbf{N}^{n+1})=0,$

where $\mathcal{F}$ is semismooth because of the mixed complementarity
condition (53b) expressed with the semismooth Evans function (52); see [23].
We solve this problem by iteration for
$(\overline{\mathbf{B}^{n+1}},\overline{\Lambda^{n+1}},\overline{\mathbf{N}^{n+1}})\approx(\mathbf{B}^{n+1},\Lambda^{n+1},\mathbf{N}^{n+1})$
only up to certain tolerance $tol$. The algorithm starts with
$(\mathbf{B}^{0},0,\mathbf{N}^{0})$ from initial data.

######  Algorithm 4.1.

(Semismooth Newton Method for (53), at every time step $n\geq 0$)

Given
$(\overline{\mathbf{B}^{n}},\overline{\Lambda^{n}},\overline{\mathbf{N}^{n}}$),
solve for
$(\overline{\mathbf{B}^{n+1}},\overline{\Lambda^{n+1}},\overline{\mathbf{N}^{n+1}})$
as shown below.

Denote the iterative guesses by
$({\mathbf{B}^{(k)}},\Lambda^{(k)},{\mathbf{N}^{(k)}})$.

  * Step (0).

Choose an initial guess ${\mathbf{B}^{(0)}}=\overline{\mathbf{B}^{n}}$,
$\Lambda^{(0)}=\overline{\Lambda^{n}}$, and
${\mathbf{N}^{(0)}}=\overline{\mathbf{N}^{n}}$, and set $k=0$.

  * Step (1).

Evaluate current residual
$\mathcal{F}^{(k)}=\mathcal{F}({\mathbf{B}^{(k)}},\Lambda^{(k)},{\mathbf{N}^{(k)}})$.
If $\|\mathcal{F}^{(k)}\|_{\infty}<tol,$ let $K=k$, and go to Step 5.

  * Step (2).

Evaluate current Jacobian. Select some
$\mathcal{J}^{(k)}\in\partial\mathcal{F}(\mathbf{B}^{(k)},\Lambda^{(k)},\mathbf{N}^{(k)})$.

  * Step (3).

Solve for the correction $\mathcal{J}^{(k)}s_{k}=-\mathcal{F}^{(k)}$; where
$s^{k}=(s_{B}^{k},s_{\Lambda}^{k},s_{N}^{k})$.

  * Step (4).

Correct the current guess: set

$\displaystyle(\mathbf{B}^{(k+1)},\Lambda^{(k+1)},\mathbf{N}^{(k+1)})=(\mathbf{B}^{(k)},\Lambda^{(k)},\mathbf{N}^{(k)})+(s_{B}^{(k)},s_{\Lambda}^{(k)},s_{N}^{(k)}).$

Increment $k$ by one, and go to Step 1.

  * Step (5).

Set $\overline{\mathbf{B}^{n+1}}=\mathbf{B}^{(K)}$,
$\overline{\Lambda^{n+1}}=\Lambda^{(K)}$ and
$\overline{\mathbf{N}^{n+1}}=\mathbf{N}^{(K)}$. STOP.

The Algorithm 4.1 is very efficient and performs in a way superior to other
solvers such as, e.g. relaxation methods given in [13]. See [1] for our
comparison of these solvers for a model PVI. We use $tol=10^{-6}$. Since the
average number of Newton iterations is between $2$ and $3$, we skip the
detailed report. In general, the solver can be proven to converge locally
$q$-superlinearly [23]. In our implementation we use MATLAB, and its built-in
linear solver for sparse systems.

##  5\. Numerical Experiments

In this section we present numerical experiments designed to show convergence
at the rates predicted by Theorem 3.3 for $d$=$2$. We also show convergence in
the cases not covered by the theory; in particular, we consider $d$=$1$,
$d$=$3$, as well as Neumann boundary conditions.

Further, we consider some simulations with nonlinear diffusivities as well as
in irregular geometries similar to those encountered at the porescale. Our
experiments are motivated by the experimental set-up in [18, 7]. We illustrate
the pointwise behavior of the biofilm and nutrient and of the total amount

(54) $\displaystyle\bar{B}(t)=\int_{\Omega}B(x,t)\;dx$

of biomass depending on modeling assumptions. For some simulations we present
both the approximations to $B(x,t)$ as well as of $N(x,t)$. In some other
simulations, the illustrations of $N(x,t)$ are predictable and are skipped. In
captions, we abbreviate “boundary conditions” to “bc”.

We examine the approximation error in two norms analyzed in Theorem 3.3. The
first ERR10=$\max_{n}(\|e_{B}^{n}\|_{0}+\|e_{N}^{n}\|_{0})$, and the second
error quantity
ERR20=$\sqrt{\sum_{n}\left(\|e_{B}^{n}\|_{1}^{2}+\|e_{N}^{n}\|_{1}^{2}\right)\Delta
t}$ for a scalar problem is shown in [14] to be of the same order as ERR10.

######  Remark 5.

In practice, we are unable to verify the convergence exactly using ERR10 and
ERR20, because the true solutions $B(x,t)$ and $N(x,t)$ to our coupled system
are not known. Manufacturing solutions under constraints is nontrivial, and
doing so for the coupled system is even more complicated except in scalar
examples in [1] or physically unrealistic examples. Thus we choose fine grid
solutions $B_{fine},N_{fine}$ as surrogates for $B,N$, with some $h_{fine}$
significantly smaller than $h$ considered in the convergence study.
Unfortunately, this approach requires also an appropriately small $\Delta
t_{fine}$, resulting in a very large number of time steps.
Therefore, our ability to actually check the convergence over all time steps
$n=1,2,\ldots$ as indicated in ERR10 and ERR20 with these large numbers of
time steps is limited. Instead, we limit ourselves to the sampling of the
spatial errors in time only over a limited set of a selected few $K$ time
steps $\Upsilon=\\{T_{1},T_{2},\ldots T_{K}\\}$ which correspond to some
selected indices $\\{N_{1},N_{2},\ldots N_{K}\\}$, different for each $\Delta
t$. In what follows we report

$\displaystyle ERR1^{\Upsilon}$ $\displaystyle=$
$\displaystyle\max_{n\in\\{N_{1},N_{2},\ldots
N_{K}\\}}(\|e_{B}^{n}\|_{0}+\|e_{N}^{n}\|_{0}),$ $\displaystyle
ERR2^{\Upsilon}$ $\displaystyle=$
$\displaystyle\sqrt{\sum_{n\in\\{N_{1},N_{2},\ldots
N_{K}\\}}\left(\|e_{B}^{n}\|_{1}^{2}+\|e_{N}^{n}\|_{1}^{2}\right)\Delta t},$

and in each instance we indicate which $T_{1},T_{2},\ldots T_{K}$ are in
$\Upsilon$.

######  Remark 6.

It is well known that using fine grid solution may somewhat overpredict the
convergence rate. In addition, in our examples the sampling of the error in
time is quite sparse, therefore we expect to see convergence rate higher than
that predicted by the theorem.

###  5.1. 1D experiments with Dirichlet Boundary Conditions

We start with a simple model problem in 1d with homogeneous Dirichlet boundary
conditions. The data in this model problem satisfies exactly the conditions in
Assumption 1, but strictly speaking the case $d$=$1$ is not covered by Theorem
3.3.

######  Example 5.1 (1D Simulation).

Let $\Omega=(0,1)$, and let the diffusion coefficients be constant
$D_{B}$=$0.5$, $D_{N}$=$0.1$. We set the initial biofilm
$B_{init}(x)=0.01|sin(\pi x)|$, the initial nutrient
$N_{init}(x)=0.02\chi_{(0.25,0.75)}$. We use Monod functions
$F(B,N)=\frac{2500N}{N+0.7}B$ and $G(B,N)=-\frac{100N}{N+0.7}B$. We also set
$B^{\ast}=0.02$, and use boundary condition $B(0,t)=B(1,t)=0=N(0,t)=N(1,t)$.

Figure 1. Evolution in Ex. 5.1 with $h=0.02$, and $\Delta t=0.005$.

#### 5.1.1.

Fig. 1 shows the typical growth of biofilm and the decay of nutrient over time
up until about $T\approx 0.02$. Since the nutrient is initially concentrated
in the middle of the domain, the majority of growth of $B$ and decay of $N$
occurs there. Eventually the nutrient diffuses away, and the biomass starts
also growing elsewhere.

At $t=0.02$ the biomass reaches the maximum density $B^{*}$, and the biofilm
“phase” in $\Omega^{*}$ forms. The Lagrange multiplier $\Lambda$ is active
when needed in $\Omega^{*}(t)$ to enforce the constraint. The Lagrange
multiplier takes always negative values, since its role (when moved to the
right hand-side) is to “push the solution down” so that (3) is satisfied. The
biofilm starts growing through the interface moving as a free boundary; this
scenario is similar to that described in Ex. 3.1. This shows that (19) is a
reasonable assumption.

The nutrient is not consumed in $\Omega^{*}$, but it slowly diffuses away
towards the external boundaries at $x=0$ and $x=1$ through which it escapes,
and towards the region $\Omega^{-}$ where it is consumed by the growing
biomass.

#### 5.1.2.

To test convergence of the numerical model, we use a fine grid solution with
$h=0.001$ and $\Delta t=0.0001$ as a surrogate for an analytical solution.
Since it would be difficult to set up $\Delta t$ to make the logarithmic terms
$[(\log(\Delta t)^{-1})^{1/4}\Delta t^{3/4}]$ conform to $h$, we choose
$\Delta t=O(h)$. We present the errors with $\Upsilon=\\{0.05,0.1\\}$. Table 1
shows $ERR1^{\Upsilon}$ and $ERR2^{\Upsilon}$.

The results demonstrate an essentially first order of convergence in ERR1Υ in
$h$, a bit higher than that predicted by the Theorem 3.3 for $d$=$2$, likely
thanks to an additional regularity of the solution exhibited by a modest size
of the region $D_{n}$ discussed in (18), and due to our error sampling
strategy discussed in Remark 6. On the other hand, we see that the convergence
order in ERR2 is about $O(h^{3/2})$. These results can be compared to those
theoretically predicted in [14] as well as our numerical results for that
theoretical case shown in the Appendix.

$h$ | $\Delta t$ | $ERR1$ | $ERR2$ |  $ERR1$ order |  $ERR2$ order
---|---|---|---|---|---
0.01 | 0.01 | 0.00026 | 0.00028 |  |
0.005 | 0.005 | 0.00014 | 0.00010 | 0.95433 | 1.4365
0.0025 | 0.0025 | 6.6251e-05 | 3.6292e-05 | 1.0347 | 1.4961
Table 1. Convergence for Ex. 5.1; Dirichlet bc

###  5.2. Experiments in $d=1$, with homogeneous Neumann bc and various
choices of diffusivity

Here we focus on the choice of diffusivity $D_{B}$, since this choice is the
primary difference between the various models in the literature. In
particular, we set up experiments with Neumann conditions to model the growth
in an isolated system, the only difference among the experiments being the
choice of $D_{B}(B)$.

######  Example 5.2 (1D Simulation, growth with Neumann conditions).

We fix the domain $\Omega=(0,1)$, $N_{init}=\chi_{\Omega}$,
$B_{init}=0.2\chi_{(0.4,0.6)}$, the diffusivity of nutrient $D_{N}=0.5$, the
maximum density $B^{*}=0.03$, the growth and utilization functions
$F(B,N)=\frac{10N}{N+0.007}B$, $G(B,N)=\frac{5N}{N+0.007}B$, respectively. We
vary the diffusivity of biofilm $D_{B}$. We consider four cases (i)
$D_{B}=0.1$, (ii) $D_{B}=0.001$, (iii)
$D_{B}=(D_{\max}-D_{\min})(B/B^{*})+D_{\min}$, (iv)
$D_{B}=(D_{\max}-D_{\min})(B^{8}/B^{*})+D_{\min}$; $D_{\max}=0.1$,
$D_{\min}=0.001$.

Figure 2. Evolution in Ex. 5.2 with $D_{B}=0.1$, $h=0.01$, $\Delta t=0.002$

Figure 3. Evolution in Ex. 5.2 with $D_{B}=0.001$, $h=0.01$, $\Delta t=0.0002$

Figure 4. Evolution in Ex. 5.2 with
$D_{B}=(D_{\max}-D_{\min})(B/B^{*})+D_{\min}$; $D_{\max}=0.1$,
$D_{\min}=0.001$, $h=0.01$, $\Delta t=0.0002$

Figure 5. Evolution of $B(x,t)$, $N(x,t)$ in Ex. 5.2 with
$D_{B}=(D_{\max}-D_{\min})(B^{8}/B^{*})+D_{\min}$; $D_{\max}=0.1$,
$D_{\min}=0.001$, $h=0.01$, $\Delta t=0.0002$

See the evolution of biofilm and nutrient in Fig. 2–5. The difference in
evolution patterns is substantial; this remains true also when $\bar{B}(t)$
alone is considered; this is explained below.

In fact $\bar{B}(t)$ is the only quantity that can be verified experimentally
without imaging. Consider momentarily the unconstrained and constrained models
governed by (1) and (5), respectively. With unlimited nutrient supply, the
growth of $\bar{B}(t)$ solving (1) in a closed system without the constraint
(3) should be exponential. (See the simulation in the Appendix.) However, the
experiments in [7] demonstrate that the microbial growth within the biofilm
phase is not exponential except in the beginning when
$\Omega^{*}\approx\emptyset$. Our simulation of (5) confirms that once the
constraint is active and the set $\Omega^{*}$ is growing, the growth ceases to
be exponential. These effects are stronger when we use small or nonlinear
diffusivities.

We explore these variants in simulations. Fig. 6 shows the evolution of the
total amount of biofilm $\bar{B}(t)$ depending on the variant (i) through
(iv); these particular models are entirely heuristic. We see that the growth
of $\bar{B}(t)$ before reaching the maximum density is exponential in the case
of constant $D_{B}$ up to the time the free boundary starts forming. The
growth is faster with larger $D_{B}$, but the free boundary moves faster with
smaller $D_{B}$. In the case of nonlinear diffusivity (iii) as in [18], as
well as in the more singular variant (iv) similar to that used in phase field
models, e.g., [29], or hybrid models in [9], we see again that the growth of
$\bar{B}(t)$ is exponential up to when $\Omega^{*}$ forms, but a more singular
diffusivity gives a more pronounced tapering effect. Here by “more singular”
$D_{B}(B)$ we mean a function with a large gradient close to $B^{*}$. We do
not use a truly singular $D_{B}(B)$ which blows up at $B\uparrow B^{*}$ such
as in [9].

Figure 6. Total $\bar{B}(t)$ on log-scale in Ex. 5.2 with (i) $D_{B}=0.1$ (top
left), (ii) $D_{B}=0.001$ (top right), (iii)
$D_{B}=(D_{\max}-D_{\min})(B/B^{*})+D_{\min}$ (bottom left), (iv)
$D_{B}=(D_{\max}-D_{\min})(B^{8}/B^{*})+D_{\min}$ (bottom right);
$D_{\max}=0.1$, $D_{\min}=0.001$. The tapering of exponential growth is more
pronounced for more singular diffusivity, with a high gradient close to
$B^{*}$.

###  5.3. Simulations in $d=2$

In this section we verify the theoretical result of Theorem 3.3. We also show
that the behavior of the coupled biofilm-nutrient dynamics depends
significantly on the geometry of the domain $\Omega$, the boundary conditions,
and the model for diffusivities $D_{B}(B)$. The examples are designed to show
the growth through interface, starting from an initial biomass concentrated in
a disk $D(r)$ centered at the origin with radius $r$.

######  Example 5.3 (Simulation in $d=2$, with Dirichlet boundary conditions).

Consider the square domain $\Omega=(-1,1)^{2}$ with $D_{B}=0.01$, and
$D_{N}=0.5$. We set $B_{{init}}=0.2\chi_{D(0.5)}$, and ${B^{\ast}}=0.3$. We
also set $N_{init}=\chi_{D(0.75)}$. The Monod functions are
$F=\frac{5N}{N+0.7}B$, $G=-\frac{0.5N}{N+0.7}B$.

Figure 7. Evolution of $B(x,t)$ in Ex. 5.3; Dirichlet bc.

Figure 8. Consumption of $N(x,t)$ in Ex 5.3; Dirichlet dc.

#### 5.3.1.

Figures 7 and 8 show the evolution of $B$ and $N$ over time; here $h=0.1$. We
use 494 nodes, and 899 elements.

The growth of $B(x,t)$ for $0<t<0.4$ is vigorous and concentrated near its
initial position, and then for $0.4<t<1$ its spread through the interface.
Nutrient is consumed most substantially where $B$ grows. Around $t=1$ both $B$
and $N$ start to decay, and the evolution is dominated by the escape of $B$
and $N$ through the boundary due to the homogeneous Dirichlet conditions.

#### 5.3.2.

For convergence test, we compute the fine grid solution $(B_{fine},N_{fine})$
as a proxy for $(B,N)$, triangulating the domain with
$\mathcal{T}_{fine}=22887$ triangles, with 11660 nodes and 34546 edges, where
the maximum length of each side of the triangles is $h_{fine}=0.02$. We use
$\Delta t_{fine}=5\times 10^{-5}$, and for the purposes of convergence testing
we store the numerical solution $(B_{fine},N_{fine})$ at
$\Upsilon=\\{0.1,0.2,0.3\\}$.

We compute the solution at coarse grids and compare it to $B_{fine}$ and
$N_{fine}$, calculating the associated values of the approximation error. The
errors are given in Table 2. It seems that $ERR1^{\Upsilon}$ is essentially of
the first order, whereas $ERR2^{\Upsilon}$ is of $O(h^{3/2})$, similarly as in
the case of $d=1$. Again we see that the free boundary is smoothly expanding,
and the size of $\sum_{n}m(D_{n})$ is only mildly increasing.

$h$ | $\Delta t$ |  $\\#$nodes |  $\\#$ elements | ERR1 | ERR2 | ERR1 ord. | ERR2 ord.
---|---|---|---|---|---|---|---
0.15 | 0.006 | 219 | 378 | 0.0499347 | 0.0474 |  |
0.1 | 0.004 | 494 | 899 | 0.0282514 | 0.0274 | 1.4047 | 1.3517
0.05 | 0.002 | 1906 | 3638 | 0.0113125 | 0.0087 | 1.3204 | 1.6551

$h$ | $\Delta t$ |  $\\#$nodes |  $\\#$ elements | ERR1 | ERR2 | ERR1 ord. | ERR2 ord.
---|---|---|---|---|---|---|---
0.15 | 0.006 | 219 | 378 | 0.0411 | 0.0550 |  |
0.1 | 0.004 | 494 | 899 | 0.0221 | 0.0371 | 1.5302 | 0.9710
0.05 | 0.002 | 1906 | 3638 | 0.0108 | 0.0141 | 1.0330 | 1.3957

Table 2. Convergence test for Ex. 5.3; Dirichlet bc (top) and for Ex. 5.4;
Neumann bc. (bottom). We use $\Upsilon=\\{0.1,0.2,0.3\\}$.

###  5.4. Experiments in $d=2$ with homogeneous Neumann boundary conditions

######  Example 5.4 (2D Simulation).

In this example, we consider data as in Ex. 5.3 except the boundary and
initial conditions. We choose homogeneous Neumann conditions to model an
isolated system, a nonsymmetric initial condition
$B_{init}=0.3\chi_{(-0.75,0)\times(-0.5,0.5)}$, and provide abundant nutrient
$N_{init}\equiv 1$.

We note that the case of Neumann boundary conditions is not covered by the
theory. However, the errors shown in Table 2 demonstrate that the order of
convergence is similar to that we obtained for Dirichlet boundary conditions.

Furthermore, we track the cumulative growth, i.e., $\bar{B}(t)$; see Fig. 10,
in which we compare the growth to that in more complex geometry and in $d$=$3$
addressed below.

###  5.5. Simulation in complex porescale geometry and in $d$=$3$.

######  Example 5.5 (Simulation in porescale geometry).

We use the data in Table 3, and we consider a domain $\Omega$ with geometry
motivated by [18]. We use nonlinear diffusivity $D_{B}=D_{B}(B)$,
$F(B,N)=\kappa_{B}\frac{N}{N+\eta}B$, $G(B,N)=-\kappa_{N}\frac{N}{N+\eta}B$,
$\Omega_{b}$ is the initial biofilm domain as shown in the illustration.

The evolution of biofilm is shown in Fig. 9 and the total amount $\bar{B}(t)$
in Fig. 10. We see the effects of reaching $B^{*}$ as well as of the limited
propagation through interface due to the complex geometry.

Maximum density of biomass | $B^{\ast}=0.12kg/m^{3}$
---|---
Growth constant | $\kappa_{B}=1.8/s$
Utilization constant | $\kappa_{N}=1.8\cdot 10/s$
Monod constant | $\eta=0.16kg/m^{3}$
Nutrient diffusivity | $D_{N}(x)=20m/s$
Biomass diffusivity | $D_{B}(B)=(D^{\ast}-D_{\ast})(B/B^{\ast})+D_{\ast}$
| $D^{\ast}=0.01,\;D_{\ast}=10^{-4}D^{\ast}$
Initial nutrient | $N_{0}(x)=\chi_{\Omega}$
Initial biomass | $B_{0}(x)=0.03\chi_{\Omega_{b}}$
Table 3. Data in Ex. 5.5 and 5.6 from [18] scaled by $10^{5}$

Figure 9. Growth of $B(x,t)$ in Ex. 5.5 simulated with $h=0.1$ and $899$
elements; Neumann bc.

Figure 10. Total $\bar{B}(t)$ for simulations with Neumann bc for Ex. 5.4
(simple $d$=$2$ geometry, porescale geometry in Ex. 5.5, and in $d$=$3$ in Ex.
5.6.

######  Example 5.6 (3D Simulation).

We use the data in Table 3. In this example, $\Omega$ is a ball with two holes
which represent solid grains. The initial biofilm is chosen to adhere to the
solid surfaces as is shown in Figure 11.

As time goes, biofilm keeps growing in its initial domain until $T\approx 1$
when it reaches its maximum $B^{*}=0.12$. After that, biofilm starts spreading
through the interface. The evolution of biofilm is illustrated in Fig. 11.

Figure 11. Evolution of biofilm in Ex. 5.6, with $h=0.2$ and 1397 elements.

We conclude with a comment on the behavior of $\bar{B}(t)$ shown in Fig. 10.
The qualitative behavior for the Ex. 5.5 and 5.6 seems different from than
that in the generic bulk 2d domain of Ex. 5.4. The growth of $\bar{B}(t)$ in
Fig. 10 appears initially fast in all cases, and a plateau of exponential
growth is achieved about the time the growth through interface begins. We do
not see this effect as strongly in Ex. 5.4 where constant diffusivity is used,
but we see this effect pronounced in the other two examples. Furthermore, the
growth through interface is further restricted in complex porescale geometry.

##  6\. Conclusions and summary

In this paper we proved and verified error estimates for a nonlinear coupled
system of parabolic variational inequalities modeling biofilm–nutrient
dynamics at microscale.

We derived error estimates of rate $O(h+(log\Delta t^{-1})^{1/4}\Delta
t^{3/4})$ in $l^{\infty}(L^{2})$ and $l^{2}(H^{1})$\- norms. The rate of
convergence was validated experimentally in 1D and 2D. Although the
theoretical analysis dealt with Dirichlet boundary conditions only, the
numerical results showed the same rate of convergence also with Neumann
boundary conditions. Moreover, simulations in 2D and 3D with irregular
geometries analogous to those obtained from imaging at porescale and with data
motivated by realistic simulations in [18] showed typical behavior of biofilm
and nutrient in porous media.

Further studies and analyses of the importance of modeling choice of
$D_{B}(B),D_{N}(N)$ are underway. We are also considering rigorous analysis of
the approximation with advection and coupling to the flow, and with other than
Galerkin FE. Last but not least, we hope to be able to address further
properties of the solutions such as non-negativity, e.g., through some form of
Discrete Maximum Principle.

### Acknowledgements

First, we would like to thank the anonymous reviewers for their suggestions
which helped to improve our paper.

Next, some of the work reported in this paper is part of the PhD thesis [1] by
Dr Azhar Alhammali at Oregon State University. She would like to gratefully
acknowledge the support of Saudi Arabian Cultural Mission to the US. In
addition, this research was partially supported by the NSF grant DMS-1522734
“Phase transitions in porous media across multiple scales”, and DMS-1912938
“Modeling with Constraints and Phase Transitions in Porous Media”, as well as
by M. Peszynska’s NSF IRD plan 2019-20.

Last but not least we would like to thank the collaborators on [18] and in
particular Dr Anna Trykozko for providing initial inspiration for the modeling
efforts on this project.

## Appendix

Here we provide two examples which provide context for and can be compared to
our convergence and simulations reported in Sec.5.

In particular, in Ex. 6.1 we show convergence of FE method for a scalar linear
PVI; this should be compared with that for the coupled constrained system
consider in Sec. 5.1.

######  Example 6.1.

In this example we approximate solutions to a scalar PVI

(55a) $\displaystyle B_{t}-0.5{B_{xx}}+\partial_{[-0.04,0.06]}B\ni$
$\displaystyle\pi^{2}\sin(x)B+3H(0.5-x)-3H(x-0.5)\;\mbox{on}\;(0,1),\;t>0,$
(55b) $\displaystyle B(0,t)=0=$ $\displaystyle B(1,t),\;t>0,$ (55c)
$\displaystyle B(x,0)=$ $\displaystyle 0.04\sin(\pi x)$ .

The approximation error in Fig. 12 for this scalar problem shows convergence
order similar to that in our 1d examples for the coupled system in Sec. 5.1.

In turn, in Ex. 6.2 we show convergence rate $O(\Delta t+h^{2})$ of an
unconstrained coupled system analyzed in [1].

Figure 12. Left: rate of convergence of FE solution for scalar PVI of Ex. 6.1,
$ERR1=\max_{n}\|B^{n}-B_{h}^{n}\|_{0}$.
$ERR2=\sqrt{\sum_{n}\|B^{n}-B_{h}^{n}\|_{1}^{2}\Delta t}$. Right: rate for
convergence of Ex. 6.2; here ERR1 and ERR2 are in a product space.

######  Example 6.2 (Unconstrained Coupled System in 1D).

Consider the same data in Example 5.1, but $B$ in this case is unconstrained,
i.e. $B^{*}=\infty$; see Fig. 13. The convergence error seems to be the usual
$O(h^{2}+\Delta t)$ in $ERR1=l^{\infty}(L^{2})$ and $ERR2=l^{2}(H_{0}^{1})$ as
in Fig. 12. This should be compared with the system under constraints in Sec
5.1 with rate close to $O(\Delta t+h)$.

Figure 13. Solutions of Ex. 6.2 without constraints

## References

  * [1] Azhar Alhammali.  Numerical Analysis of a System of Parabolic Variational Inequalities with Application to Biofilm Growth.  PhD thesis, Oregon State University, 2019.
  * [2] Claudio Baiocchi.  Discretization of evolution variational inequalities.  In Partial differential equations and the calculus of variations, Vol. I, volume 1 of Progr. Nonlinear Differential Equations Appl., pages 59–92. Birkhäuser Boston, Boston, MA, 1989.
  * [3] Viorel Barbu.  Nonlinear differential equations of monotone types in Banach spaces.  Springer Monographs in Mathematics. Springer, New York, 2010.
  * [4] John W. Barrett and Klaus Deckelnick.  Existence, uniqueness and approximation of a doubly-degenerate nonlinear parabolic system modelling bacterial evolution.  Math. Models Methods Appl. Sci., 17(7):1095–1127, 2007.
  * [5] Haïm Brézis.  Problèmes unilatéraux.  J. Math. Pures Appl. (9), 51:1–168, 1972.
  * [6] Franco Brezzi, William W. Hager, and P.-A. Raviart.  Error estimates for the finite element solution of variational inequalities.  Numer. Math., 28(4):431–443, 1977.
  * [7] FS Colwell, C Verba, AR Thurber, Y Alleau, D Koley, M Peszynska, and ME Torres.  Feasibility of biogeochemical sealing of wellbore cements: Lab and simulation tests.  In AGU Fall Meeting Abstracts, 2014.
  * [8] Timothy B Costa, Kenneth Kennedy, and Malgorzata Peszynska.  Hybrid three-scale model for evolving pore-scale geometries.  Computational Geosciences, 22:925–950, 2018.
  * [9] HJ Eberl, C Picioreanu, JJ Heijnen, and MCM van Loosdrecht.  A three-dimensional numerical study on the correlation of spatial structure, hydrodynamic conditions, and mass transfer and conversion in biofilms.  CHEMICAL ENGINEERING SCIENCE, 55(24):6209–6222, 2000.
  * [10] Anozie Ebigbo, A Phillips, Robin Gerlach, Rainer Helmig, Alfred B Cunningham, Holger Class, and Lee H Spangler.  Darcy-scale modeling of microbially induced carbonate mineral precipitation in sand columns.  Water Resources Research, 48(7), 2012.
  * [11] Charles M. Elliott.  On the finite element approximation of an elliptic variational inequality arising from an implicit time discretization of the Stefan problem.  IMA J. Numer. Anal., 1(1):115–125, 1981.
  * [12] Alexandre Ern and Jean-Luc Guermond.  Theory and practice of finite elements, volume 159 of Applied Mathematical Sciences.  Springer-Verlag, New York, 2004.
  * [13] Roland Glowinski.  Numerical methods for nonlinear variational problems.  Springer Series in Computational Physics. Springer-Verlag, New York, 1984\.
  * [14] Claes Johnson.  A convergence estimate for an approximation of a parabolic variational inequality.  SIAM J. Numer. Anal., 13(4):599–606, 1976.
  * [15] C. T. Kelley and Jim Rulla.  Solution of the time discretized Stefan problem by Newton’s method.  Nonlinear Anal., 14(10):851–872, 1990.
  * [16] F. A. MacLeod, H. M. Lappin-Scott, and J. W. Costerton.  Plugging of a model rock system by using starved bacteria.  APPLIED AND ENVIRONMENTAL MICROBIOLOGY, 54(6):1365–1372, 1988.
  * [17] R. H. Nochetto, M. Paolini, and C. Verdi.  Self-adaptive mesh modification for parabolic FBPs: theory and computation.  In Free boundary value problems (Oberwolfach, 1989), volume 95 of Internat. Ser. Numer. Math., pages 181–206. Birkhäuser, Basel, 1990.
  * [18] Malgorzata Peszynska, Anna Trykozko, Gabriel Iltis, Steffen Schlueter, and Dorthe Wildenschild.  Biofilm growth in porous media: Experiments, computational modeling at the porescale, and upscaling.  Advances in Water Resources, 95:288–301, 2016.
  * [19] R Schulz, N Ray, F Frank, HS Mahato, and P Knabner.  Strong solvability up to clogging of an effective diffusion–precipitation model in an evolving porous medium.  European Journal of Applied Mathematics, 28(2):179–207, 2017.
  * [20] R. E. Showalter.  Monotone operators in Banach space and nonlinear partial differential equations, volume 49 of Mathematical Surveys and Monographs.  American Mathematical Society, Providence, RI, 1997.
  * [21] Youneng Tang, Albert J. Valocchi, Charles J. Werth, and Haihu Liu.  An improved pore-scale biofilm model and comparison with a microfluidic flow cell experiment.  WATER RESOURCES RESEARCH, 49(12):8370–8382, 2013.
  * [22] Vidar Thomée.  Galerkin finite element methods for parabolic problems, volume 25 of Springer Series in Computational Mathematics.  Springer-Verlag, Berlin, 1997.
  * [23] Michael Ulbrich.  Semismooth Newton methods for variational inequalities and constrained optimization problems in function spaces, volume 11 of MOS-SIAM Series on Optimization.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA; Mathematical Optimization Society, Philadelphia, PA, 2011.
  * [24] Tycho L van Noorden, IS Pop, Anozie Ebigbo, and Rainer Helmig.  An upscaled model for biofilm growth in a thin strip.  Water resources research, 46(6), 2010.
  * [25] C Verba, AR Thurber, Y Alleau, D Koley, F Colwell, and ME Torres.  Mineral changes in cement-sandstone matrices induced by biocementation.  International Journal of Greenhouse Gas Control, 49:312–322, 2016\.
  * [26] C. Verdi.  On the numerical approach to a two-phase Stefan problem with nonlinear flux, volume 372 of Istituto di Analisi Numerica del Consiglio Nazionale delle Ricerche [Institute of Numerical Analysis of the National Research Council].  Istituto di Analisi Numerica del Consiglio Nazionale delle Ricerche, Pavia, 1984.
  * [27] C. Vuik.  An $L^{2}$-error estimate for an approximation of the solution of a parabolic variational inequality.  Numer. Math., 57(5):453–471, 1990.
  * [28] Mary Fanett Wheeler.  A priori $L_{2}$ error estimates for Galerkin approximations to parabolic partial differential equations.  SIAM J. Numer. Anal., 10:723–759, 1973.
  * [29] T. Zhang and I. Klapper.  Mathematical model of biofilm induced calcite precipitation.  WATER SCIENCE AND TECHNOLOGY, 61(11):2957–2964, 2010.