# Robust time-discretisation and linearisation schemes for singular and
degenerate evolution systems modelling biofilm growth

R.K.H. Smeets email<EMAIL_ADDRESS>University of Amsterdam, Korteweg-de
Vries Institute for Mathematics, The Netherlands K. Mitra Eindhoven University
of Technology, Department of Mathematics and Computer Science, The Netherlands
I.S. Pop Hasselt University, Faculty of Science, Belgium S. Sonner Radboud
University, IMAPP - Mathematics, The Netherlands

###### Abstract

We propose and analyse numerical schemes for a system of quasilinear,
degenerate evolution equations modelling biofilm growth as well as other
processes such as flow through porous media and the spreading of wildfires.
The first equation in the system is parabolic and exhibits degenerate and
singular diffusion, while the second is either uniformly parabolic or an
ordinary differential equation. First, we introduce a semi-implicit time
discretisation that has the benefit of decoupling the equations. We prove the
positivity, boundedness, and convergence of the time-discrete solutions to the
time-continuous solution. Then, we introduce an iterative linearisation scheme
to solve the resulting nonlinear time-discrete problems. Under weak
assumptions on the time-step size, we prove that the scheme converges
irrespective of the space discretisation and mesh. Moreover, if the problem is
non-degenerate, the convergence becomes faster as the time-step size
decreases. Finally, employing the finite element method for the spatial
discretisation, we study the behaviour of the scheme, and compare its
performance to other commonly used schemes. These tests confirm that the
proposed scheme is robust and fast.

Keywords: degenerate diffusion; time discretisation; linearisation;
unconditional convergence; stability; biofilm models; porous medium equation

MSC: 65M12, 65M22, 35K51, 35K65

##  1 Introduction

###  1.1 Motivation

Let $T>0$ be a maximal time and $\Omega\subset\mathds{R}^{d}$
($d\in\mathbb{N}$) be a bounded Lipschitz domain. With $Q=\Omega\times(0,T]$
denoting the parabolic space-time cylinder, we consider the following class of
degenerate quasilinear parabolic systems

$\displaystyle\partial_{t}u=\Delta\Phi(u)+f(v)u,$ (1.1a)
$\displaystyle\partial_{t}v=\mu\nabla\cdot\left(D(u)\nabla v\right)+g(u,v)$
(1.1b)

in $Q$. The first equation (1.1a) describes evolution of a population density
$u$ and exhibits degenerate and possibly also singular diffusion leading to
the formation of free boundaries propagating at a finite speed. More
specifically, the monotone function $\Phi$ vanishes as $u$ tends to zero and
can in addition have a singularity as $u$ tends to its maximum value. The
diffusion coefficient $D$ in the second equation which describes evolution of
a substrate concentration $v$ is assumed non-degenerate, i.e. it is bounded
above and below by positive constants. However, the mobility coefficient $\mu$
appearing in (1.1b) may be either 0 or 1, leading to either a coupled system
of a parabolic equation and an ordinary differential equation (if $\mu=0$), or
two parabolic equations (if $\mu=1$). The growth and spreading of the
population $u$ might depend on multiple substrates. Nevertheless, for
simplicity, we consider only one substrate $v$ in (1.1b), as an extension to
multiple substrates is straightforward. The system is completed by the initial
conditions $u(\cdot,0)=u_{0}$ and $v(\cdot,0)=v_{0}$ for given functions
$u_{0},v_{0}$, and by homogeneous Dirichlet boundary conditions for $u$ and,
if $\mu\neq 0$, also for $v$. Analysing system (1.1) analytically and
numerically is challenging due to the degenerate and singular diffusion in the
first equation which leads to free boundaries and steep gradients, and the
nonlinear coupling with the second equation.

The main motivation for this work comes from the biofilm growth models in [13,
15], where the solution component $u$ in (1.1a) describes the (normalized)
biomass density whose evolution is dictated by the diffusion operator

$\displaystyle\Delta\Phi(u)=\nabla\cdot\left(\frac{u^{\alpha}}{(1-u)^{\beta}}\nabla
u\right),\qquad\text{ for some }\qquad\alpha,\beta\geq 1.$ (1.2)

With this, (1.1a) is coupled to a reaction-diffusion partial differential
equation (PDE) or an ordinary differential equation (ODE) modelling the
evolution of the growth-limiting nutrient concentration $v$. In the resulting
model, the biofilm occupies the region where $\\{u>0\\}$. Observe that the
biomass diffusion in (1.2) shows a degeneracy of porous-medium type as $u$
approaches $0$, ensuring a finite speed of propagation of the interface
between the biofilm and the surrounding region, as well as a singularity as
$u$ approaches $1$. The latter implies that the solution $u$ remains bounded
by a constant strictly less than 1 despite the growth term $f$ in the
equation. The second equation, (1.1b), describes the evolution of the nutrient
concentration. The case $\mu=0$ corresponds to immobile substrates (e.g., in
the case of cellulolytic biofilms, [13]), while $\mu=1$ corresponds to
diffusive substrates (e.g., whenever biofilms grow in an aqueous medium,
[15]). The biofilm growth model was also extended to take multiple substrates
into account, both mobile and immobile, like in [17, 27]. The results of our
paper generalise directly to these cases. More complex multi-species biofilm
models including cross-diffusion have been studied in [11, 37].

Systems of the form (1.1) are not limited to models for biofilm growth. For
instance, coupled systems of parabolic and ordinary differential equations
appear in the modelling of two-phase or unsaturated flow through porous media
when effects like hysteresis or dynamic capillarity are taken into account
[27, 30] (degenerate, but non-singular diffusion). They also appear in
wildfire models [38] (nonlinear but non-degenerate diffusion), and reaction,
diffusion, and adsorption/desorption models in a porous medium [21], to name a
few other applications.

The aim of this paper is to develop robust, efficient, and structure-
preserving time discretisation and linearisation methods for System (1.1)
relying on minimal regularity assumptions, and converging even for degenerate
and singular diffusion. In what follows, the time discretisation is
introduced, as well as the linear iterative schemes. Furthermore, the main
results concerning the stability and convergence of these schemes are stated,
their proofs being given in the subsequent sections.

###  1.2 Time discretisation

To define the time discretisation, we take $N\in\mathbb{N}$ and let
$\tau=\frac{T}{N}$ be the time-step size, which is chosen uniform for the ease
of presentation. With $n\in\\{0,\dots,N\\}$, let $t_{n}=n\tau$ be the time-
step and denote by $u_{n}$ the approximation of $u$ at $t=t_{n}$, and
similarly for $v_{n}$. Then, we use an Euler semi-implicit approach for the
time discretisation. All terms in (1.1) are discretised implicitly except for
the reaction functions $f$ and $g$, which are discretised semi-implicitly.

######  Problem 1.1 (Semi-implicit time discretisation).

Given the approximate solutions $u_{n-1}$ and $v_{n-1}$ at time $t_{n-1}$,
find the approximate solution pair $(u_{n},v_{n})$ at the next time step
$t_{n}$ by solving the following system

$\displaystyle\frac{1}{\tau}(u_{n}-u_{n-1})=\Delta\Phi(u_{n})+f(v_{n-1})u_{n},$
(1.3a)
$\displaystyle\frac{1}{\tau}(v_{n}-v_{n-1})=\mu\nabla\cdot\left(D(u_{n})\nabla
v_{n}\right)+g(u_{n},v_{n-1}).$ (1.3b)

This approach has several advantages over explicit and implicit
discretisations. Explicit methods lead to a loss of regularity of the time-
discrete solutions, which results in instability due to the already low
regularity triggered by degenerate diffusion. On the other hand, implicit
schemes have the advantage that only a weak restriction on the time-step size,
independent of the space-discretisation, is needed to guarantee stability.
Fully implicit time-discretisation schemes for the mentioned biofilm models
were analysed in [4, 12, 19, 20]. However, for fully implicit schemes, the two
time-discrete equations originating from (1.1) are coupled which makes them
challenging to solve especially using an iterative method. In the semi-
implicit approach (1.3) the two equations are decoupled. This allows us to
solve them sequentially instead of iteratively, i.e. we first solve for
$u_{n}$ and then update $v_{n}$ using the known $u_{n}$. In fact, given
$u_{n}$, (1.3b) is a linear problem for $v_{n}$ which can be solved rather
easily. Moreover, with the proposed semi-implicit discretisation we retain the
same accuracy and stability one would expect from fully implicit
discretisations. Generalising the results in [12] for the scalar equation
(1.1a), we show that under weak assumptions on $\tau$, the discretisation
(1.3) is well-posed, the time-discrete solutions preserve positivity, remain
bounded, and converge to the time-continuous solutions as $\tau\to 0$. The
exact results are stated in Theorems 3.1 and 3.2. Below, we summarise them
omitting technical details.

######  Theorem (Well-posedness, boundedness, and convergence of the time-
discrete solutions).

For all sufficiently small time steps $\tau$, there exists a unique weak
solution of (1.3). Moreover, the time-discrete solutions $u_{n},v_{n}$ are
positive and $u_{n}$ is bounded almost everywhere in $\Omega$. In particular,
if $\Phi$ is singular, $u_{n}$ is bounded by a constant strictly less than the
singularity. Lastly, the time-discrete solutions converge to the time-
continuous solutions as $\tau\to 0$.

Concerning space discretisations, we also restrict ourselves to mentioning
works addressing specifically System (1.1). In this respect, a finite
difference method was used in [14]. The finite volume method was considered in
[20], and the convergence of the space-time dicretisation scheme was proven
using an entropy formulation. Convergence results of mixed finite elements for
a variation of the PDE-PDE model of biofilms were shown in [1], and for a PDE-
ODE model in the context of porous media flow in [8]. For the numerical
results presented here, we use finite elements. However, the numerical
analysis is done in a time-discrete, but continuous-in-space setting.
Therefore, the results are independent of the chosen spatial discretisation.

###  1.3 Linearisation

The time-discrete system (1.3) is nonlinear, degenerate and singular. For
approximating its numerical solution, stable iterative linearisation schemes
are needed. Most of the works addressing such type of problems are focusing on
a direct approach. More precisely, for solving (1.3a), $w=\Phi(u)$ is
considered as the primary unknown yielding $u=\Phi^{-1}(w)$ which is then
linearised by expanding in terms of the last iterate. With reference to
(1.3a), this approach is convenient because there are no spatial derivatives
applied to the nonlinear terms. However, this requires that $\Phi$ is an
invertible function; if this is not the case, then a regularisation step is
required. Alternatively, one can use $u$ as the primary variable, and avoid
inverting the function $\Phi$. In this case, the nonlinearity appears under
the Laplace operator, which makes the construction of robust linearisation
schemes and, in particular, proving the convergence a complex task.

Following [7] here we consider instead a split formulation involving two
primary unknowns, $u$ and $w$, which are related through the algebraic
relationship $w=\Phi(u)$. We reformulate the time-discrete version of (1.1a)
as a system of a linear elliptic equation and an algebraic one. For this
reformulated system, we construct linear iterative schemes having a stable and
robust behaviour. They all fit in the general framework given below.

######  Problem 1.2 (Linearisation scheme).

For $i\in\mathbb{N}$, let $u^{i}_{n}$ be the $i^{\rm th}$ iterate of the
$n^{\rm th}$ time-step, and let $u^{0}_{n}:=u_{n-1}$ be given. To compute
$u^{i}_{n}$ from $u^{i-1}_{n}$, first solve for
$(\tilde{u}^{i}_{n},w^{i}_{n})$ satisfying

$\displaystyle\frac{1}{\tau}(\tilde{u}^{i}_{n}-u_{n-1})=\Delta
w^{i}_{n}+f(v_{n-1})\tilde{u}^{i}_{n},$ (1.4a) $\displaystyle
L^{i}_{n}(\tilde{u}^{i}_{n}-u^{i-1}_{n})=w^{i}_{n}-\Phi(u^{i-1}_{n}).$ (1.4b)
The factors $L^{i}_{n}$ will be specified below. The way they are chosen is
defining the different schemes used here. Finally, we take the positive part
of $\tilde{u}^{i}_{n}$, $\displaystyle
u^{i}_{n}=\left[\tilde{u}^{i}_{n}\right]_{+}.$ (1.4c)

If $w_{n}^{i}$ and $\tilde{u}_{n}^{i}$ (and, consequently, also $u^{i}_{n}$)
converge to $w_{n}$ and $u_{n}$ respectively, then the limits satisfy
$w_{n}=\Phi(u_{n})$, and $u_{n}$ is a (weak) solution of (1.3a). Observe that
the formulation used in 1.2, where (1.3a) is split into a linear elliptic
equation and a nonlinear algebraic one, is well suited for degenerate
problems. In particular, no regularisation is needed in the slow diffusion
case, i.e. when $u^{i}_{n}$ takes values for which $\Phi^{\prime}$ vanishes.

As mentioned, the choice of the factors $L^{i}_{n}$ in (1.4b) leads to
different linearisation schemes. With $L^{i}_{n}=\Phi^{\prime}(u^{i-1}_{n})$,
one obtains the Newton Scheme (NS) in the context of the splitting formulation
(1.4). The convergence is guaranteed rigorously in the fully discrete case,
but this depends strongly on the spatial discretisation and mesh size [6, 32,
36]. For time-dependent problems, since the initial guess is often the
solution at the previous time step, this implies that the time-step size
should be sufficiently small, which may cancel the advantages brought by the
implicit discretisation. In the same category, we mention the modified Picard
scheme [9] and the Jäger–Kačur scheme [24, 25], for which the linear
convergence can be proven rigorously under similar conditions as for the NS.

Ideally, one works with a scheme that is unconditionally convergent, i.e. the
time-step $\tau$ does not depend on the spatial discretisation and the mesh-
size $h$, which was one of the main reasons for choosing an implicit time
discretisation. This can be achieved by using the L-scheme (LS) [26, 33],
which is nothing but choosing $L^{i}_{n}=L$ (a suficiently large constant) in
(1.4b). This scheme has a guaranteed but linear convergence, irrespective of
the initial guess and the spatial discretisation. These results were extended
in [35] to doubly-degenerate problems, where a Hölder continuous, not
necessarily strictly increasing nonlinearity appears under the time
derivative. The drawback of this scheme is its significantly slower
convergence when compared to the NS (whenever the latter converges) [26, 41].
Improvements can be made by choosing $L$ adaptively in each step, or by
performing first a number of LS steps, and then switching to the NS entirely
when the iterations are close enough to the solution [26, 44].

Such observations have lead to the Modified L-scheme, or M-scheme (MS) for
short, as introduced in [28], and on which we mainly focus here. We define
$L^{i}_{n}\coloneqq\max\\{\Phi^{\prime}(u^{i-1}_{n})\\\
+M\tau^{\gamma},\,2M\tau^{\gamma}\\}$, for some $\gamma\in(0,1]$ and
sufficiently large $M>0$. Hence, the MS can be viewed as a combination of the
NS and the LS where the LS is a first order global method ensuring stability,
while the NS is a first order local method speeding up convergence. In
particular, the MS has a first order local term to speed up convergence but
regularizes it with a second order global term which captures the evolution of
$u_{n}$ to ensure stability. In an earlier work [28], it has been shown for
the direct formulation that the MS is indeed unconditionally stable while
achieving much better convergence rates than the LS.

In this work, we apply the LS and MS to the splitting formulation in 1.2,
allowing us to handle systems with porous medium type degeneracies, as well as
singular diffusion. Similarly to [28], we prove that the MS converges
unconditionally and that, in the non-degenerate case, it has a contraction
rate that scales with $\tau$. The two main results for the LS and MS in
Section 4 are summarised in the following theorems.

######  Theorem (Convergence of the L-scheme).

For a sufficiently small time step size $\tau$ independent of the mesh size
$h$, the L-scheme converges to a function $u_{n}$ that is the weak solution of
our time-discretised eq. 1.3a. In the non-degenerate case, the L-scheme
converges with a contraction rate $\alpha<1$.

######  Theorem (Convergence of the M-scheme).

For a sufficiently small time step size $\tau$ independent of the mesh size
$h$ and under certain boundedness conditions on the iterates, the M-scheme
converges to a function $u_{n}$ that is the weak solution of our time-
discretised solution eq. 1.3a. In the non-degenerate case, the M-scheme
converges with a contraction rate $\alpha<1$, that scales with $\tau^{\gamma}$
for some $\gamma\in(0,1]$.

The outline of our paper is as follows: In Section 2 we provide the required
functional setting and background and state the structural assumptions. In
Section 3 we prove the results for the time discretisation and in Section 4 we
show the convergence results for both the LS and the MS. In Section 5 we
perform numerical simulations and compare the performances of the NS and MS
for a porous medium equation and both cases of the biofilm model (PDE-PDE
system and PDE-ODE system). Finally, in Section 6 we summarize our results and
discuss potential future research.

##  2 Preliminaries

###  2.1 Functional setting and background

Let $\Omega\subset\mathds{R}^{d}$ be a bounded Lipschitz domain. The
corresponding $L^{2}$ inner product and norm are denoted by
$\left(\cdot,\cdot\right)$ and $\norm{\cdot}$, and norms with respect to other
Banach spaces $V$ by $\norm{\cdot}_{V}$. By
$\left(W^{1,p}(\Omega),\|\cdot\|_{W^{1,p}(\Omega)}\right)$, $1\leq p<\infty$,
we denote the Sobolev spaces and use the short-hand notation
$H^{1}(\Omega)\coloneqq W^{1,2}(\Omega)$. The space $H^{1}_{0}(\Omega)$ is the
closure of $C^{\infty}_{c}(\Omega)$ in $H^{1}(\Omega)$, which is equipped with
the equivalent norm $\norm{u}_{H^{1}_{0}(\Omega)}\coloneqq\norm{\nabla u}$ due
to the Poincaré inequality

$\norm{u}\leq C_{\Omega}\norm{\nabla u}\quad\text{for all }u\in
H^{1}_{0}(\Omega),$ (2.1)

where $C_{\Omega}>0.$ The dual space of $H^{1}_{0}\left(\Omega\right)$ is
denoted by $H^{-1}\coloneqq\left(H^{1}_{0}\left(\Omega\right)\right)^{\ast}$
with the norm

$\norm{u}_{H^{-1}\Omega}\coloneqq\sup_{\phi\in
H^{1}_{0}\left(\Omega\right)}\frac{\langle u,\phi\rangle}{\norm{\nabla\phi}},$
(2.2)

where $\langle\cdot,\cdot\rangle$ denotes the duality paring. Since we
consider homogeneous Dirichlet boundary conditions, we will mainly use
$H^{1}_{0}(\Omega)$. Lastly, we consider the Bochner spaces $L^{p}(0,T;V)$,
with $V$ a Banach space, equipped with the norm

$\norm{u}_{L^{p}(0,T;V)}\coloneqq\left(\int_{0}^{T}\norm{u(t)}_{V}^{p}dt\right)^{1/p}<\infty.$
(2.3)

We will frequently use Young’s inequality

$uv\leq\frac{1}{2\rho}u^{2}+\frac{\rho}{2}v^{2},\quad\text{for }\rho>0\text{
and }u,v\in\mathds{R},$ (2.4)

the Cauchy-Schwarz inequality

$\left|\int_{\Omega}uv\right|\leq\norm{u}\norm{v},\qquad u,v\in
L^{2}(\Omega),$ (2.5)

and the discrete Gronwall Lemma: Let
$\\{u_{n}\\}_{n\in\mathbb{N}},\\{a_{n}\\}_{n\in\mathbb{N}},\\{b_{n}\\}_{n\in\mathbb{N}}$
be non-negative sequences such that $u_{n}\leq
a_{n}+\sum_{k=1}^{n-1}b_{k}u_{k}$, then

$u_{n}\leq
a_{n}+\sum_{k=1}^{n-1}a_{k}b_{k}\exp\left(\sum_{k<j<n}b_{j}\right).$ (2.6)

Lastly, for convex functions $\Psi\in C\left(\mathds{R}^{+}\right)$ with
$\Psi(0)=0$, we have the following inequalities

$\displaystyle\text{Jensen's
inequality:}\quad\Psi\left(\frac{1}{|\Omega|}\int_{\Omega}|f|\right)\leq\frac{1}{|\Omega|}\int_{\Omega}\Psi(|f|)\quad\text{for
}f\in L^{2}(\Omega),$ (2.7a) $\displaystyle\text{Super-
additivity:}\quad\Psi(a)+\Psi(b)\leq\Psi(a+b)\quad\text{for all }a,b\geq 0.$
(2.7b)

To denote the positive and negative part we write
$\left[u\right]_{+}\coloneqq\max\\{0,u\\}$ and
$\left[u\right]_{-}\coloneqq\min\\{0,u\\}$ and use the notation $a\lesssim b$
if $a\leq Cb$ for some constant $C>0$. Moreover,
$\mathds{R}^{+}=\\{x\in\mathds{R}\;\colon\;x\geq 0\\}$ and
$\mathds{R}^{+}_{\ast}=\\{x\in\mathds{R}\;\colon\;x>0\\}$. In some proofs we
will also use the notation

$I(a,b)=\\{x\;\colon\;\min\\{a,b\\}\leq x\leq\max\\{a,b\\}\\},\qquad
a,b\in\mathds{R}.$ (2.8)

If $u,v$ are two given functions then $I(u,v)$ should be considered pointwise
a.e.

Finally, throughout this work, $C>0$ will denote a generic constant that might
change in each occurrence and from line to line but will always be independent
of $\tau$.

###  2.2 Structural assumptions

Dependng on the function $\Phi$ in (1.2), the problems may be singular-
degenerate, or have a structure resembling the porous medium equations. To
cover both cases, we introduce a maximum density/concentration value $b$,
where $b=1$ in the former case (the biofilm model) and $b=\infty$ for porous
medium type equations.

With $b$ as above, We make the following structural assumptions for Equation
1.1

  1. (P1)

$\Phi\colon[0,b)\to\mathds{R}^{+}$ is an increasing function with locally
Lipschitz continuous derivatives, satisfying

$\displaystyle\quad\Phi(0)=0,\quad\lim_{u\nearrow b}\Phi(u)=\infty,\text{ and
}\Phi^{\prime}(u)>0\text{ for }u\in(0,b),$

with $\inf_{[0,b]}\Phi^{\prime}=:\phi_{\mathrm{m}}\in[0,\infty)$, and
$\sup_{[0,b]}\Phi^{\prime}=:\phi_{\mathrm{M}}\in(0,\infty]$. We furthermore
require that $\Phi^{\prime}$ is strictly increasing in $[0,\varepsilon_{0})$
for some $\varepsilon_{0}\in(0,b)$.

  2. (P2)

$f\colon\mathds{R}^{+}\to\mathds{R}$ is Lipschitz continuous and bounded, with
$\|f\|_{\infty}=f_{\mathrm{M}}$ for some constant $f_{\mathrm{M}}\geq 0$.
$g\colon[0,b)\times\mathds{R}^{+}\to\mathds{R}$ is Lipschitz continuous with
Lipschitz constant $g_{\mathrm{M}}>0$. Moreover, we assume that
$g(\cdot,0)\geq 0$.

  1. (P3)

$D\colon[0,b)\to\mathds{R}$ is a continuous function. There exist constants
$D_{\mathrm{m}},D_{\mathrm{M}}$ s.t. $0<D_{\mathrm{m}}\leq D(u)\leq
D_{\mathrm{M}}<\infty$ for all $u\in[0,b)$.

The functions $\Phi$, $f$, $g$, $D$ are extended for arguments
$u\in\mathds{R}^{-}$ by their values at $u=0$. If $b=1$, the functions that
are bounded at $u=b$ are also extended for $u\in[b,\infty)$ by their values at
$u=b$.

######  Remark 2.0.1 (Validity of the assumptions (P1) \- (P3)).

The biofilm models [13, 15] (see Equation (5.8) in Section 5) and the porous
medium equation satisfy the assumptions (P1) \- (P3), with $b=1$ and
$b=\infty$ respectively. Note that we only consider non-negative solutions as
$u$ and $v$ denote densities and/or concentrations.

For the initial data we assume the following.

  1. (P4)

The initial conditions $u_{0},v_{0}\colon\Omega\to[0,\infty)$ are s.t.
$v_{0}\in L^{2}(\Omega)$, $u_{0}\in L^{\infty}(\Omega)$ and
$\|u_{0}\|_{\infty}<b$.

Finally, for the ease of presentation we consider homogeneous Dirichlet
boundary conditions for $u$, and also for $v$ if $\mu=1$. The results here can
be the extended to mixed Dirichlet-Neumann boundary conditions, following the
ideas in [23, 29].

###  2.3 Weak formulation of the time continuous problem

We consider weak solutions of Equation 1.1 with homogeneous Dirichlet boundary
conditions and initial data satisfying (P4).

######  Definition 2.1 (Weak formulation).

A weak solution of (1.1) is a pair $(u,v)\in\\\
C([0,T];L^{2}(\Omega))^{2}\allowbreak\cap\,H^{1}(0,T;H^{-1}(\Omega))^{2}$ s.t.
$(\Phi(u),\mu v)\in L^{2}(0,T;H^{1}_{0}\left(\Omega\right))^{2}$,
$(u,v)(0)=(u_{0},v_{0})$, and

$\displaystyle\int_{0}^{T}\left\langle\partial_{t}u,\phi\right\rangle+\int_{0}^{T}(\nabla\Phi(u),\nabla\phi)$
$\displaystyle=\int_{0}^{T}(f(v)u,\phi),$
$\displaystyle\int_{0}^{T}\left\langle\partial_{t}v,\eta\right\rangle+\int_{0}^{T}\mu(D(u)\nabla
v,\nabla\eta)$ $\displaystyle=\int_{0}^{T}(g(u,v),\eta)$

hold for all $\phi,\eta\in L^{2}(0,T;H^{1}_{0}\left(\Omega\right))$.

For the well-posedness of Problem (1.1) with $b=\infty$ and $\mu=1$, we refer
to [3, 31]. If $\mu=0$, the existence and uniqueness of solutions for the
coupled system follow by $L^{1}$-contraction similarly as in [29]. Well-
posedness results for the system with $b=1$ and either Dirichlet or mixed
Dirichlet-Neumann boundary conditions were obtained in [23, 29]. In
particular, uniqueness can be shown if $D$ in (P3) is independent of $u$ [23]
or if $\mu=0$ [29]. Under these assumptions, local well-posedness was also
shown for homogeneous Neumann boundary conditions in [29]. Furthermore, it was
shown that solutions are non-negative, and if (P4) holds, then the solution
$u$ is bounded by a constant strictly less than $1$, i.e. the singularity in
the diffusion coefficient is not attained. The local Hölder continuity of
solutions of such systems was studied in [22]. In this direction, we also
mention [16] where the specific PDE-PDE biofilm model [15] (corresponding to
$\mu=1$ with diffusion coefficient (1.2)) was analyzed, [27] where the
existence of solutions for a similar degenerate PDE-ODE system was studied,
and [5] where a doubly degenerate PDE-PDE system was analysed.

Lastly, we remark that, for simplicity, we assume homogeneous Dirichlet
boundary conditions. Extending the results to mixed Dirichlet-Neumann boundary
conditions is possible following the arguments in [23, 29].

##  3 Time discretisation

In this section we analyse the following weak form of the time discretised
system (1.3). We first state it in a weak form.

######  Problem (Weak formulation of the time-discretised system).

Let $n\in\\{1,\dots,N\\}$ and $u_{n-1},\,v_{n-1}\in L^{2}(\Omega)$ given. Find
$(u_{n},v_{n})\subset L^{2}(\Omega)^{2}$ such that $\Phi(u_{n}),\,\mu v_{n}\in
H^{1}_{0}\left(\Omega\right)$, and for all $\phi,\eta\in
H^{1}_{0}\left(\Omega\right)$ it holds

$\displaystyle\left(\frac{1}{\tau}(u_{n}-u_{n-1}),\phi\right)+(\nabla\Phi(u_{n}),\nabla\phi)$
$\displaystyle=(f(v_{n-1})u_{n},\phi),$ (3.1a)
$\displaystyle\left(\frac{1}{\tau}(v_{n}-v_{n-1}),\eta\right)+\mu(D(u_{n})\nabla
v_{n},\nabla\eta)$ $\displaystyle=(g(u_{n},v_{n-1}),\eta).$ (3.1b)

Throughout this paper, we write $w_{n}=\Phi(u_{n})$ and use the shorthand
notation

$\displaystyle h_{n-1}\coloneqq 1-\tau f(v_{n-1}).$ (3.2)

Note that $h_{n-1}$ is positive if $\tau<1/f_{\mathrm{M}}$.

######  Remark 3.0.1 (The decoupling of the equations).

Observe that the solution $u_{n}\in L^{2}(\Omega)$ of (3.1a) does not depend
on the solution $v_{n}\in L^{2}(\Omega)$ of (3.1b). Hence, the system (3.1)
can be solved sequentially, i.e. we first solve the nonlinear problem (3.1a)
and subsequently the linear problem (3.1b).

In this section we prove the following results, already briefly mentioned in
Section 1.2.

######  Theorem 3.1 (Well-posedness and boundedness of the time-discrete
solutions).

For $\tau<1/f_{\mathrm{M}}$, there exists a unique weak solution
$(u_{n},v_{n})$ of (3.1). Moreover, there exist $\tau_{\rm
disc}\coloneqq\min\\{1/f_{\mathrm{M}},1/g_{\mathrm{M}}\\}>0$ and
$\breve{u}\in[0,b)$ independent of $n$, such that

$\displaystyle 0\leq u_{n}\leq\breve{u},\quad\text{ and }\quad 0\leq
v_{n}\quad\text{ a.e. in }\Omega\quad\text{ for all }1\leq n\leq N\text{ and
}\tau<\tau_{\rm disc}.$ (3.3)

In fact, $\breve{u}$ is given by

$\displaystyle\breve{u}=\begin{cases}\|u_{0}\|_{L^{\infty}}\exp\left(\frac{Tf_{\mathrm{M}}}{1-\tau
f_{\mathrm{M}}}\right)&\text{ if }b=\infty,\\\\[4.5pt]
\Phi^{-1}\left(\|\Phi(u_{0})\|_{L^{\infty}}+\frac{{\rm
diam}(\Omega)^{2}}{2d}f_{\mathrm{M}}\right)<1&\text{ if }b=1,\end{cases}$
(3.4)

######  Remark 3.1.1 (Computable upper bound for $u_{n}$).

Observe that (3.4)–(3.3) provides a uniform upper bound for $u_{n}$ that is a
priori computable. This will be used in Section 4 to show the convergence of
the iterative linearisation scheme.

######  Theorem 3.2 (Convergence of the time-discrete solutions).

Let $(u,v)\in C([0,T];L^{2}(\Omega))^{2}$ be the unique weak solution of
Equation 1.1. For a time-step size $\tau=\frac{T}{N_{\tau}}>0$,
$N_{\tau}\in\mathbb{N}$, let
$\\{(u_{n},v_{n})\\}_{n\in\mathbb{N}}\subset(L^{2}(\Omega))^{2}$ be the time-
discrete solution of (3.1) with $\\{w_{n}\\}_{n\in\mathbb{N}}\subset
H^{1}_{0}\left(\Omega\right)$. Then, in addition to (P4), if $u_{0}\in
H^{1}_{0}(\Omega)$, then along any sequence of $\tau$ converging to 0 we have

$\displaystyle\sum_{n=1}^{N_{\tau}}\int_{(n-1)\tau}^{n\tau}\left[\|u_{n}-u(t)\|^{2}+\|w_{n}-\Phi(u(t))\|^{2}+\|v_{n}-v(t)\|^{2}\right]\mathrm{d}t\to
0.$ (3.5)

Moreover, if $u_{0}\not\in H^{1}_{0}(\Omega)$, then consider an approximation
$u^{\varepsilon}_{0}\in H^{1}_{0}(\Omega)$ of the initial data such that
$\|u^{\varepsilon}_{0}-u_{0}\|\leq\varepsilon$, for fixed $\varepsilon>0$, and
let
$\\{(u_{n}^{\varepsilon},v_{n}^{\varepsilon})\\}_{n\in\mathbb{N}}\subset(L^{2}(\Omega))^{2}$
be the corresponding time-discrete solutions with
$\\{w^{\varepsilon}_{n}\\}_{n\in\mathbb{N}}\subset
H^{1}_{0}\left(\Omega\right)$. Then, along any sequence of
$(\varepsilon,\tau)$ converging to $(0,0)$ one has

$\displaystyle\sum_{n=1}^{N_{\tau}}\int_{(n-1)\tau}^{n\tau}\left[\|u^{\varepsilon}_{n}-u(t)\|^{2}+\|w_{n}^{\varepsilon}-\Phi(u(t))\|^{2}+\|v^{\varepsilon}_{n}-v(t)\|^{2}\right]\mathrm{d}t\to
0.$ (3.6)

The proofs of Theorems 3.1 and 3.2 are based on several lemmas.

###  3.1 Proof of Theorem 3.1: well-posedness, positivity, and boundedness

####  3.1.1 Existence and uniqueness

We first prove the existence and uniqueness of solutions of the system of
equations (3.1).

######  Lemma 3.3 (Well-posedness for (3.1)).

For $\tau<1/f_{\mathrm{M}}$, there exists a unique weak solution of the time
discretised system (3.1).

###### Proof.

(Step 1) Existence of $u_{n}$: As the equations are decoupled, we can first
prove the existence of the solution $u_{n}$ of (3.1a). To this end, we use
arguments in [34]. We consider the function
$\Psi\coloneqq\Phi^{-1}:\mathds{R}^{+}\to[0,b)$ which satisfies

$\displaystyle\Psi(0)=0,\quad\Psi^{\prime}=\frac{1}{\Phi^{\prime}\circ\Psi}\geq
0$ (3.7)

by (P1). Then, the energy $J\colon H^{1}_{0}\left(\Omega\right)\to\mathds{R}$,
defined by

$\displaystyle
J(w)\coloneqq\int_{\Omega}\left[h_{n-1}\int_{0}^{w}\Psi+\frac{\tau}{2}|\nabla
w|^{2}-u_{n-1}w\right]$ (3.8)

is convex and coercive for $\tau<1/f_{\mathrm{M}}$. Hence, a minimizer
$w_{n}\in H^{1}_{0}\left(\Omega\right)$ of $J$ exists, and $u_{n}=\Psi(w_{n})$
solves the corresponding Euler-Lagrange equation (3.1a). Then using (P1), for
an arbitrary $\varepsilon\in(0,b)$, we have
$0\leq\Psi\leq\Psi(\varepsilon)<\infty$ in $[0,\varepsilon]$ and $\Psi$ is
Lipschitz in $(\varepsilon,b)$. Hence, by (3.7) it follows that, for a.e.
$x\in\Omega$ we have $0\leq u_{n}^{2}(x)=\Psi^{2}(w_{n})(x)$, so $u_{n}$ is
measurable. Integrating over $\Omega$, one gets

$\displaystyle 0\leq\int_{\Omega}u_{n}^{2}$
$\displaystyle=\int_{\Omega}\Psi^{2}(w_{n})=\int_{\\{0\leq
w_{n}\leq\varepsilon\\}}\Psi^{2}(w_{n})+\int_{\\{w_{n}>\varepsilon\\}}\Psi^{2}(w_{n})$
$\displaystyle\leq\Psi(\varepsilon)^{2}|\Omega|+C\left(1+\norm{w_{n}}^{2}\right)\leq
C+C\norm{\nabla w_{n}}_{L^{2}(\Omega)}^{2},$

for some constant $C>0$, where, in the last estimate, we used the Poincaré
inequality (2.1). Since $w_{n}\in H^{1}_{0}\left(\Omega\right)$, this shows
that the integral is finite, so $u_{n}\in L^{2}(\Omega)$.

(Step 2) Uniqueness of $u_{n}$: Assume that for a given $v_{n-1}\in
L^{2}(\Omega)$, there are two solutions $u_{n},\,\tilde{u}_{n}\in
L^{2}(\Omega)$ of (3.1a) with $w_{n}=\Phi(u_{n})$ and
$\tilde{w}_{n}=\Phi(\tilde{u}_{n})$ in $H^{1}_{0}(\Omega)$. For their
difference we obtain

$\displaystyle\frac{1}{\tau}(h_{n-1}(u_{n}-\tilde{u}_{n}),\varphi)+(\nabla(\Phi(u_{n})-\Phi(\tilde{u}_{n})),\nabla\varphi)=0\qquad\forall\varphi\in
H^{1}_{0}(\Omega).$ (3.9)

Note that $\phi=\left[\Phi(u_{n})-\Phi(\tilde{u}_{n})\right]_{+}\in
H^{1}_{0}\left(\Omega\right)$, see e.g. [10]. Choosing $\phi$ as a test
function in (3.9) leads to

$\frac{1}{\tau}\left(h_{n-1}\left(u_{n}-\tilde{u}_{n}\right),\left[\Phi(u_{n})-\Phi(\tilde{u}_{n})\right]_{+}\right)+\norm{\nabla\left[\Phi(u_{n})-\Phi(\tilde{u}_{n})\right]_{+}}_{L^{2}(\Omega)}^{2}=0.$
(3.10)

As $\Phi$ is an increasing function, we note that
$\left(u_{n}-\tilde{u}_{n}\right)\left[\Phi(u_{n})-\Phi(\tilde{u}_{n})\right]_{+}\geq
0$. Hence, both terms in (3.10) are non-negative and therefore, have to be
equal to 0. We conclude that
$\norm{\nabla\left[\Phi(u_{n})-\Phi(\tilde{u}_{n})\right]_{+}}_{L^{2}(\Omega)}^{2}=0$,
which results in
$\norm{\left[\Phi(u_{n})-\Phi(\tilde{u}_{n})\right]_{+}}_{L^{2}(\Omega)}^{2}=0$
by the Poincaré inequality (2.1). This implies that
$\Phi(\tilde{u}_{n})\geq\Phi(u_{n})$ a.e. in $\Omega$, but as $\Phi$ is an
increasing function, we also find that $\tilde{u}_{n}\geq u_{n}$ a.e. in
$\Omega$. Due to the symmetry in the arguments, it follows in the same way
that $u_{n}\geq\tilde{u}_{n}$ a.e. which implies that $u_{n}=\tilde{u}_{n}$
a.e. in $\Omega$. The uniqueness of $u_{n}$ can also be shown via the
$L^{1}$-contraction principle [46].

(Step 3) Existence-uniqueness of $v_{n}$: We now prove the existence and
uniqueness for the solution $v_{n}$ of (3.1b). In the PDE-ODE case, i.e.
$\mu=0$, we have an explicit expression for $v_{n}$,

$v_{n}=v_{n-1}+\tau g(u_{n},v_{n-1}).$ (3.11)

Hence, the existence and uniqueness of $v_{n}$ follows from the existence and
uniqueness of $u_{n}$.

In the PDE-PDE case, i.e. $\mu=1$, existence and uniqueness follows from the
Lax-Milgram theorem [18]. Indeed, the weak form can be rewritten as

$a(v_{n},\eta)=l(\eta)\quad\forall\eta\in H^{1}_{0}\left(\Omega\right),$
(3.12)

where the bilinear form $a\colon H^{1}_{0}\left(\Omega\right)\times
H^{1}_{0}\left(\Omega\right)\to\mathds{R}$ is given by
$a(v_{n},\eta)=(v_{n},\eta)+\tau\mu(D(u_{n})\nabla v_{n},\nabla\eta)$. It is
bounded and coercive since $0<D_{\mathrm{m}}\leq D\leq D_{\mathrm{M}}<\infty$
by (P3). Moreover, $l(\eta)=(\tau g(v_{n-1},u_{n})+v_{n-1},\eta)$ is a bounded
linear functional on $H^{1}_{0}\left(\Omega\right)$. Consequently, there
exists a unique solution $v_{n}\in H^{1}_{0}\left(\Omega\right)$, which
concludes the proof. ∎

####  3.1.2 Positivity and boundedness in $L^{\infty}(\Omega)$

For the time-continuous biofilm models it was shown in [23, 29] that the
solutions $u$ and $v$ are non-negative, and that $u<1$. We aim to prove that
these properties also hold for the time-discrete solutions. First, we derive
bounds for $u_{n}$ in the general case, i.e. including porous medium type
diffusion.

######  Lemma 3.4 (Positivity and boundedness of $u_{n}$).

Let $u_{n-1}\in L^{\infty}(\Omega)$ be positive a.e. in $\Omega$. Then for
$\tau<1/f_{M}axi$, the solution $u_{n}\in L^{2}(\Omega)$ of (3.1a) is positive
and bounded a.e. in $\Omega$. More precisely, we have

$0\leq u_{n}\leq\sup\left\\{\frac{u_{n-1}}{1-\tau
f(v_{n-1})}\right\\}\quad\text{ a.e. in }\Omega$ (3.13)

for all $1\leq n\leq N.$

###### Proof.

To prove that $u_{n}$ is bounded from above, we use the test function
$\left[\Phi(u_{n})-\Phi(a)\right]_{+}$, for some $a\in\mathds{R}^{+}$ in
(3.1a). Note that $[\Phi(u_{n})-\Phi(a)]_{\pm}\in H^{1}(\Omega)$ and
$\left[\Phi(u_{n})\mp\Phi(a)\right]_{\pm}$ = 0 if $u_{n}=0$ as $\Phi(0)=0$,
and thus $\left[\Phi(u_{n})\mp\Phi(a)\right]_{\pm}\in H^{1}_{0}(\Omega)$. We
find

$\displaystyle\int_{\Omega}h_{n-1}\left(u_{n}-a\right)\left[\Phi(u_{n})-\Phi(a)\right]_{+}+\int_{\Omega}h_{n-1}\left(a-\frac{1}{h_{n-1}}u_{n-1}\right)\left[\Phi(u_{n})-\Phi(a)\right]_{+}$
$\displaystyle=-\tau\int_{\Omega}\nabla\Phi(u_{n})\cdot\nabla\left[\Phi(u_{n})-\Phi(a)\right]_{+}=-\tau\int_{\Omega}\nabla\left[\Phi(u_{n})-\Phi(a)\right]_{+}^{2}\leq
0.$ (3.14)

Let $\tau<1/f_{\mathrm{M}}$ and $a\coloneqq\sup\frac{1}{h_{n-1}}u_{n-1}$,
which implies that $a\geq 0$ as $h_{n-1},u_{n-1}\geq 0$. Then, the second term
on the left hand side is positive. The first term is also positive as
$\left(u_{n}-a\right)\left[\Phi(u_{n})-\Phi(a)\right]_{+}\geq 0$ since $\Phi$
is increasing. We conclude that the inequality in (3.14) must be an equality.
This is only possible if $\left[\Phi(u_{n})-\Phi(a)\right]_{+}=0$, and thus
$\Phi(u_{n})\leq\Phi(a)$. As $\Phi$ is an increasing function, this implies
that $u_{n}\leq a$ and hence,

$u_{n}\leq
a=\sup\left\\{\frac{1}{h_{n-1}}u_{n-1}\right\\}=\sup\left\\{\frac{u_{n-1}}{1-\tau
f(v_{n-1})}\right\\}.$ (3.15)

We use the same arguments to prove that $u_{n}\geq 0$, but this time with
$a=0$ i.e. $\Phi(a)=\Phi(0)=0$. Using the test function
$\phi=\left[\Phi(u_{n})\right]_{-}$ we conclude that
$\left[\Phi(u_{n})\right]_{-}=0$, and thus $u_{n}\geq 0$. ∎

An explicit upper bound for $u_{n}$ can be given in terms of the initial
conditions and $f_{\mathrm{M}}$, which we provide in the following result.

######  Lemma 3.5 (Explicit upper bound $u_{n}$).

Let $u_{0}\in L^{\infty}(\Omega)$ satisfy assumption (P4). Then, for the
sequence $\\{(u_{n},v_{n})\\}_{n=1}^{N}\subset L^{2}(\Omega)^{2}$ solving
(3.1), one has

$\|u_{n}\|_{L^{\infty}(\Omega)}\leq\|u_{0}\|_{L^{\infty}(\Omega)}\exp\left(\frac{n\tau
f_{\mathrm{M}}}{1-\tau f_{\mathrm{M}}}\right).$ (3.16)

######  Remark 3.5.1 (Upper bound as $\tau\to 0$).

Assuming that $u_{n}\to u(t)$ in $L^{\infty}(\Omega)$ when $\tau=t/n\to 0$,
the upper bound (3.16) implies that

$\|u(t)\|_{L^{\infty}(\Omega)}\leq\|u_{0}\|_{L^{\infty}(\Omega)}\exp\left(tf_{M}\right).$
(3.17)

###### Proof.

By (3.15) in the proof of Lemma 3.4, we have using
$\frac{1}{(1-x)}\leq\exp(\frac{x}{1-x})$ for $|x|<1$ that,

$\begin{split}\|u_{n}\|_{L^{\infty}(\Omega)}&\leq\sup\left\\{\frac{u_{n-1}}{1-\tau
f(v_{n-1})}\right\\}\leq\frac{\|u_{n-1}\|_{L^{\infty}(\Omega)}}{1-\tau
f_{\mathrm{M}}}\leq\frac{\|u_{0}\|_{L^{\infty}(\Omega)}}{\left(1-\tau
f_{\mathrm{M}}\right)^{n}}\\\
&\leq\|u_{0}\|_{L^{\infty}(\Omega)}\exp\left(\frac{n\tau
f_{\mathrm{M}}}{1-\tau f_{M}}\right).\end{split}$

∎

In the biofilm case, i.e. $b=1$, we can improve the upper bound. As shown e.g.
in [23, 29], the solution $u$ of the time continuous system is strictly less
than 1, so we aim to prove this also for the approximate solutions $u_{n}$.

######  Lemma 3.6 (Upperbound $u_{n}$ if $b=1$).

Consider the biofilm case, i.e. $b=1$, and let $u_{0}\in L^{\infty}(\Omega)$
satisfy assumption (P4) and $\tau<1/f_{\mathrm{M}}$. Then for the sequence
$\\{(u_{n},v_{n})\\}_{n=1}^{N}\\\ \subset L^{2}(\Omega)^{2}$ solving (3.1),
one has

$0\leq u_{n}\leq 1-\delta\qquad\text{a.e. in}\ \Omega,$ (3.18)

for all $1\leq n\leq N$, and some $\delta>0$.

###### Proof.

By Lemma 3.4 we have $0\leq u_{n}\leq C$, for some constant $C>0$. Let
$\tilde{\omega}\in H^{1}_{0}(\Omega)+\|\Phi(u_{0})\|_{L^{\infty}(\Omega)}$ be
the solution of

$\begin{split}\left(\nabla\tilde{\omega},\nabla\phi\right)=\left(Cf_{\mathrm{M}},\phi\right)\quad\text{for
all }\phi\in H^{1}_{0}(\Omega).\end{split}$ (3.19)

As $Cf_{\mathrm{M}}\in\mathds{R}^{+}_{\ast}$, by properties of the Poisson
equation, we know that $\tilde{\omega}\in L^{\infty}(\Omega)$. Further, since
$\tilde{\omega}$ is superharmonic, the maximum principle implies that
$\Phi(u_{0})\leq\tilde{\omega}$. We will prove that
$\Phi(u_{n})\leq\tilde{\omega}$ for all $1\leq n\leq N$ by induction. Assuming
it holds for $n-1$, we subtract (3.19) from (3.1a) and multiply both sides by
$\tau$. We then choose the test function
$\phi=\left[\Phi(u_{n})-\tilde{\omega}\right]_{+}\in
H^{1}_{0}\left(\Omega\right)$ to find

$\begin{split}&\left(u_{n}-\Phi^{-1}(\tilde{\omega})+\Phi^{-1}(\tilde{\omega})-u_{n-1},\left[\Phi(u_{n})-\tilde{\omega}\right]_{+}\right)+\tau\left(\nabla\left(\Phi(u_{n})-\tilde{\omega}\right),\nabla\left[\Phi(u_{n})-\tilde{\omega}\right]_{+}\right)\\\
&+\tau\left(f_{M}C-f(v_{n-1})u_{n},\left[\Phi(u_{n})-\tilde{\omega}\right]_{+}\right)=0.\end{split}$
(3.20)

By the induction hypothesis, we have $\Phi^{-1}(\tilde{\omega})-u_{n-1}\geq
0$, while
${\left(u_{n}-\Phi^{-1}(\tilde{\omega})\right)\left[\Phi(u_{n})-\tilde{\omega}\right]_{+}}$
${\geq 0}$ as $\Phi$ is an increasing function. The other terms are also
positive as

$\tau\left(\nabla\left(\Phi(u_{n})-\tilde{\omega}\right),\nabla\left[\Phi(u_{n})-\tilde{\omega}\right]_{+}\right)=\tau\norm{\nabla\left[\Phi(u_{n})-\tilde{\omega}\right]_{+}}^{2}\text{
and }f_{\mathrm{M}}C-f(v_{n-1})u_{n}\geq 0$

by definition of $f_{\mathrm{M}}$ and $C$. Hence, the Poincare inequality
implies that

$\norm{\left[\Phi(u_{n})-\tilde{\omega}\right]_{+}}_{L^{2}(\Omega)}^{2}\leq
C_{\Omega}\norm{\nabla\left[\Phi(u_{n})-\tilde{\omega}\right]_{+}}_{L^{2}(\Omega)}^{2}=0,$
(3.21)

and thus $\Phi(u_{n})\leq\tilde{\omega}$. To conclude the proof, we recall
that $\tilde{\omega}$ is bounded and hence,

$\begin{split}0\leq
u_{n}\leq\Phi^{-1}(\tilde{\omega})=1-\delta,\quad\delta>0.\end{split}$ (3.22)

∎

######  Remark 3.6.1 (Effective Lipschitz continuity of $\Phi$).

By Lemma 3.6, in the biofilm case, i.e. $b=1$, we can effectively restrict the
domain of $\Phi$ to $[0,1-\delta]\subset[0,1)$. Within this interval,
$\Phi^{\prime}$ is Lipschitz continuous as stated in assumption (P1).

We have shown that $u_{n}\leq 1-\delta$ for some $\delta>0$, but we aim to
derive an explicit bound. Such a bound will be useful in Section 4 when we
propose the linearisation scheme and is provided in the following lemma.

######  Lemma 3.7 (Explicit upper bound $u_{n}$ if $b=1$).

Consider the biofilm case, i.e. $b=1$, and let $u_{0}\in L^{\infty}(\Omega)$
satisfy assumption (P4) and $\tau<1/f_{\mathrm{M}}$. Then for the sequence
$\\{(u_{n},v_{n})\\}_{n=1}^{N}\subset L^{2}(\Omega)^{2}$ solving (3.1), one
has

$0\leq u_{n}\leq\Phi^{-1}(\tilde{C})<1\qquad\text{a.e. in}\ \Omega,$ (3.23)

for all $1\leq n\leq N$, where

$\tilde{C}=\norm{\Phi(u_{0})}_{L^{\infty}(\Omega)}+\frac{{\rm
diam}(\Omega)^{2}}{2d}f_{\mathrm{M}},$ (3.24)

and $d$ is the spatial dimension of $\Omega\subset\mathds{R}^{d}$.

###### Proof.

Let $\tilde{\omega}$ be the solution of Equation 3.19 as in the proof of Lemma
3.6. Then, $\tilde{\omega}$ is a classical solution to the elliptic problem

$\begin{split}\begin{cases}-\Delta\tilde{\omega}=f_{\mathrm{M}}&\text{on
$\Omega$},\\\ \tilde{\omega}=\norm{\Phi(u_{0})}_{L^{\infty}(\Omega)}&\text{on
$\partial\Omega$},\end{cases}\end{split}$ (3.25)

as we can take $C=1$ by Lemma 3.6.

As $f_{\mathrm{M}}>0$, $\tilde{\omega}$ is superharmonic, and we conclude that
$\tilde{\omega}\geq\norm{\Phi(u_{0})}_{L^{\infty}(\Omega)}\geq 0$ by the
maximum principle. We then define the function
$w=f_{\mathrm{M}}\norm{x-\bar{x}}^{2}/(2d)\geq 0$, where $d$ is the spatial
dimension and $\bar{x}$ is given by
$\bar{x}_{i}=|\Omega|^{-1}\int_{\Omega}x_{i}$. It is straightforward to see
that $\Delta w=f_{\mathrm{M}}$.

If we consider $z=\tilde{\omega}+w$, it satisfies $-\Delta
z=-\Delta\tilde{\omega}-\Delta w=0$. Hence, by the maximum principle it
follows that

$\begin{split}\norm{\tilde{\omega}}_{L^{\infty}(\Omega)}\leq\norm{z}_{L^{\infty}(\Omega)}\leq\norm{z}_{L^{\infty}(\partial\Omega)}&=\norm{\tilde{\omega}}_{L^{\infty}(\partial\Omega)}+\frac{\norm{x-\bar{x}}^{2}_{L^{\infty}(\partial\Omega)}}{2d}f_{\mathrm{M}}\\\
&=\norm{\Phi(u_{0})}_{L^{\infty}(\Omega)}+\frac{{\rm
diam}(\Omega)^{2}}{2d}f_{\mathrm{M}},\end{split}$ (3.26)

where we used that $w,\tilde{w},z\geq 0$, which implies that
$0\leq\tilde{w}\leq z$. This concludes the proof. ∎

We now show that $v_{n}$ is positive, similarly as in Lemma 3.4 for $u_{n}$,
for both cases $b=1$ and $b=\infty$.

######  Lemma 3.8 (Positivity of $v_{n}$).

Let the assumptions of Lemma 3.4 hold. Let $v_{n-1}\in L^{2}(\Omega)$ be
positive a.e. in $\Omega$. Then for $\tau<1/g_{\mathrm{M}}$ the solution
$v_{n}\in L^{2}(\Omega)$, $\mu v_{n}\in H^{1}_{0}\left(\Omega\right)$ of
(3.1b) is positive a.e. in $\Omega$.

###### Proof.

By Lemmas 3.4, 3.5 and 3.6, $u_{n}$ is positive and bounded in $[0,b)$
implying that $g(u,\cdot)$ is well-defined. For the positivity of $v$, first
observe that $G(u,t):=t+\tau g(u,t)$ is an increasing function in $t$ for
$\tau<1/g_{\mathrm{M}}$. For $\mu=1$, inserting $\eta=[v_{n}]_{-}$ in (3.1b)
gives

$\displaystyle(v_{n},[v_{n}]_{-})+\tau\mu(D(u_{n})\nabla
v_{n},\nabla[v_{n}]_{-})=\tau(g(u_{n},v_{n-1}),[v_{n}]_{-})+(v_{n-1},[v_{n}]_{-})$
$\displaystyle=$ $\displaystyle\
\tau(g(u_{n},v_{n-1})-g(u_{n},0),[v_{n}]_{-})+\tau(g(u_{n},0),[v_{n}]_{-})+(v_{n-1},[v_{n}]_{-})$
$\displaystyle=$ $\displaystyle\
\tau(G(u_{n},v_{n-1})-G(u_{n},0),[v_{n}]_{-})+\tau(g(u_{n},0),[v_{n}]_{-})\leq
0,$

since $G(u_{n},v_{n-1})\geq G(u_{n},0)$ due to $v_{n-1}$ being positive, and
$g(u_{n},0)\geq 0$ from (P2). Using a similar test function in the case
$\mu=0$, we conclude that

$\int_{\Omega}[v_{n}]_{-}^{2}+\tau\mu\int_{\Omega}D(u_{n})\nabla[v_{n}]_{-}^{2}\leq
0$ (3.27)

from which we conclude that $[v_{n}]_{-}=0$ a.e. in $\Omega$, or in other
words, $v_{n}\geq 0$. ∎

We now have all the necessary results to prove Theorem 3.1.

###### Proof of Theorem 3.1.

The existence and uniqueness of weak solutions of (3.1) is provided by Lemma
3.3. The positivity and boundedness of $u_{n}$ are proven in Lemma 3.4, while
the explicit bounds are given in Lemma 3.5 and Lemma 3.7. Finally, the
positivity of $v_{n}$ is the result of Lemma 3.8. ∎

###  3.2 Proof of Theorem 3.2: convergence of the time-discrete solutions

Here Rothe’s method is used to prove the convergence of the time-discrete
solutions $\\{(u_{n},v_{n})\\}_{n\in\mathbb{N}}\subset(L^{2}(\Omega))^{2}$ of
(3.1). For a time-step size $\tau>0$ with time-steps $t_{n}:=n\tau$ (recall
that $T=N_{\tau}\tau$ is fixed), and a sequence
$\\{z_{n}\\}_{n\in\mathbb{N}}\subset L^{2}(\Omega)$, we construct the piece-
wise constant and affine time-interpolations
$\hat{z}_{\tau},\,\bar{z}_{\tau}\in L^{2}(\Omega\times[0,T])$ as

$\displaystyle\hat{z}_{\tau}(t):=z_{n},\quad\bar{z}_{\tau}(t):=z_{n-1}+\frac{t-t_{n-1}}{\tau}(z_{n}-z_{n-1})\quad\text{
if }t\in(t_{n-1},t_{n}]\text{ for some }n\in\mathbb{N}.$ (3.28)

####  3.2.1 Uniform boundedness of the interpolates in Bochner spaces

######  Lemma 3.9 (Uniform boundedness of $\hat{w}_{\tau}$, $\bar{w}_{\tau}$,
$\bar{v}_{\tau}$ with respect to $\tau$).

For a time-step size $\tau>0$, let
$\\{(u_{n},w_{n},v_{n})\\}_{n\in\mathbb{N}}$ satisfy the assumptions in
Theorem 3.2, and let $\\{\hat{w}_{\tau}\\}$, $\\{\hat{v}_{\tau}\\}$,
$\\{\bar{w}_{\tau}\\}$, $\\{\bar{v}_{\tau}\\}$ be the piecewise constant,
respectively piecewise linear time-interpolations introduced in (3.28). Then,
there exists a constant $\bar{C}>0$ independent of $\tau>0$ such that for
both, $z=\bar{w}_{\tau},$ and $z=\hat{w}_{\tau}$, it holds

$\displaystyle\sup_{0\leq t\leq
T}\|z(t)\|^{2}_{L^{\infty}(\Omega)}+\int_{0}^{T}\|\nabla z\|^{2}\leq\bar{C},$
(3.29) $\displaystyle\sup_{0\leq t\leq T}\|\nabla
z(t)\|+\int_{0}^{T}\|\partial_{t}\bar{w}_{\tau}\|^{2}\leq\bar{C}[1+\|\nabla\Phi(u_{0})\|^{2}].$
(3.30)

Additionally, for $\mu\in\\{0,1\\}$ it holds

$\displaystyle\sup_{0\leq t\leq
T}\|\bar{v}_{\tau}(t)\|^{2}+\int_{0}^{T}\|\partial_{t}\bar{v}_{\tau}\|^{2}_{H^{-1}(\Omega)}+(1-\mu)\sup_{0\leq
t\leq
T}\|\partial_{t}\bar{v}_{\tau}(t)\|^{2}+\mu\int_{0}^{T}\|\nabla\bar{v}_{\tau}\|^{2}\leq\bar{C}.$
(3.31)

###### Proof.

Observe that Theorem 3.1, specially (3.4) directly yields

$\displaystyle\sup_{0\leq t\leq
T}\|\hat{w}_{\tau}(t)\|_{L^{\infty}(\Omega)}\leq\sup_{1\leq n\leq
N_{\tau}}\|w_{n}\|_{L^{\infty}(\Omega)}\overset{\eqref{eq:uup}}{\leq}\Phi(\breve{u})<\infty.$

Similarly, $\|\bar{w}_{\tau}(t)\|_{L^{\infty}(\Omega)}<\infty$ since
$\bar{w}_{\tau}(t)$ is a convex combination of $\\{w_{n}\\}_{n\in\mathbb{N}}$.
The other estimates follow closely the Rothe method, see e.g. [29] for an
identical context, or [27].

(Step 1) Bound (3.31): Inserting $\eta=v_{n}$ in (3.1b) one has

$\displaystyle\frac{1}{\tau}(v_{n}-v_{n-1},v_{n})+\mu(D(u_{n})\nabla
v_{n},\nabla v_{n})=(g(u_{n},v_{n-1}),v_{n}).$ (3.32)

To rewrite the first term we use the identity
$a(a-b)=\frac{1}{2}[a^{2}-b^{2}+(a-b)^{2}]$,

$\displaystyle\frac{1}{\tau}(v_{n}-v_{n-1},v_{n})=\frac{1}{2\tau}[\|v_{n}\|^{2}-\|v_{n-1}\|^{2}+\|v_{n}-v_{n-1}\|^{2}],$
(3.33a) and for the second term, (P3) implies that
$\displaystyle\mu(D(u_{n})\nabla v_{n},\nabla v_{n})\geq\mu
D_{\mathrm{m}}\|\nabla v_{n}\|^{2}.$ (3.33b) For the third term, notice that
$|g(u_{n},v_{n-1})|\leq|g(u_{n},v_{n-1})-g(u_{n},0)|+|g(u_{n},0)|\leq
C[1+|v_{n-1}|]$ for some constant $C>0$, which follows from the Lipschitz
continuity of $g$ in (P2) and (3.4). Then, one has
$\displaystyle(g(u_{n},v_{n-1}),v_{n})=(g(u_{n},v_{n-1}),v_{n}-v_{n-1})+(g(u_{n},v_{n-1}),v_{n-1})$
(3.33c) $\displaystyle\overset{\ref{ass: 2},\eqref{eq:
Young}}{\leq}C[1+\|v_{n-1}\|^{2}]+\|v_{n}-v_{n-1}\|^{2},$ (3.33d)

and summing up the estimates above from $n=1$ to $n=N_{\tau}$, we obtain

$\displaystyle\|v_{N_{\tau}}\|^{2}+(1-2\tau)\sum_{n=1}^{N_{\tau}}\|v_{n}-v_{n-1}\|^{2}+2\mu
D_{\mathrm{m}}\sum_{n=1}^{N_{\tau}}\|\nabla
v_{n}\|^{2}\tau\leq\|v_{0}\|^{2}+2C\sum_{n=1}^{N_{\tau}}[1+\|v_{n-1}\|^{2}]\tau.$

For $\tau<1/2$, applying the discrete Gronwall lemma (2.6) to the above
inequality reveals that $\|v_{n}\|$ is uniformly bounded with respect to
$\tau$ provided $1\leq n\leq N_{\tau}$. Substituting this back into the above
inequality, one obtains

$\displaystyle\|v_{N_{\tau}}\|^{2}+\sum_{n=1}^{N_{\tau}}\|v_{n}-v_{n-1}\|^{2}+\mu
D_{\mathrm{m}}\sum_{n=1}^{N_{\tau}}\|\nabla v_{n}\|^{2}\tau\leq C.$ (3.34)

Observe that $\hat{v}_{\tau}(t)=v_{n}$ for $t\in(t_{n-1},t_{n}]$, and
$\bar{v}_{\tau}$ is a convex combination of $\\{v_{n}\\}_{n\in\mathbb{N}}$.
Hence, the above inequality implies that $\|\hat{v}_{\tau}(t)\|$ and
$\|\bar{v}_{\tau}(t)\|$ are uniformly bounded with respect to $\tau$.
Likewise,
$\int_{0}^{T}\|\nabla\hat{v}_{\tau}\|^{2}=\sum_{n=1}^{N_{\tau}}\|\nabla
v_{n}\|^{2}\tau$ which is uniformly bounded due to (3.34) if $\mu=1$, and the
same also holds for $\int_{0}^{T}\|\nabla\bar{v}_{\tau}\|^{2}$. Observe that
for $t\in(t_{n-1},t_{n}]$,

$\displaystyle\|\partial_{t}\bar{v}_{\tau}(t)\|_{H^{-1}(\Omega)}:=\sup_{\begin{subarray}{c}\eta\in
H^{1}_{0}(\Omega)\\\
\|\nabla\eta\|=1\end{subarray}}\left\langle\frac{1}{\tau}(v_{n}-v_{n-1}),\eta\right\rangle$
$\displaystyle\qquad\overset{\eqref{eq: time-discrete weak
1}}{=}\sup_{\begin{subarray}{c}\eta\in H^{1}_{0}(\Omega)\\\
\|\nabla\eta\|=1\end{subarray}}\left[-\mu(D(u_{n})\nabla
v_{n},\nabla\eta)+(g(u_{n},v_{n-1}),\eta)\right]$
$\displaystyle\qquad\overset{\eqref{eq:Poincare},\ref{ass: 2}}{\leq}\mu
D_{\mathrm{M}}\|\nabla
v_{n}\|+C_{\Omega}\|g(u_{n},v_{n-1})\|\overset{\eqref{eq:v-disc-
bound},\eqref{eq:uup}}{\leq}\mu D_{\mathrm{M}}\|\nabla\hat{v}_{\tau}\|+C.$
(3.35a) This implies that
$\int_{0}^{T}\|\partial_{t}\bar{v}_{\tau}(t)\|^{2}_{H^{-1}(\Omega)}$ is
uniformly bounded with respect to $\tau$ since
$\int_{0}^{T}\|\nabla\hat{v}_{\tau}\|^{2}$ is. If in addition $\mu=0$, then
$\displaystyle\|\partial_{t}\bar{v}_{\tau}(t)\|=\|g(u_{n},v_{n-1})\|\overset{\eqref{eq:v-disc-
bound},\eqref{eq:uup}}{\leq}C.$ (3.35b)

Combining (3.34) and (3.35) we obtain (3.31).

(Step 2) Bounds (3.29)–(3.30): Proving (3.29), requires taking $\phi=w_{n}$ as
a test function in (3.1a). The arguments are identical to Step 2 in the proof
of Lemma 4.3 in [29] and hence, will be omitted for the sake of brevity. For
obtaining (3.30), we insert $\phi=w_{n}-w_{n-1}=\Phi(u_{n})-\Phi(u_{n-1})$ in
(3.1a) to get

$\displaystyle\left(\frac{1}{\tau}(u_{n}-u_{n-1}),\Phi(u_{n})-\Phi(u_{n-1})\right)+(\nabla
w_{n},\nabla(w_{n}-w_{n-1}))=(f(v_{n-1})u_{n},w_{n}-w_{n-1}).$ (3.36)

Noting that $\partial_{t}\bar{w}_{\tau}=(\Phi(u_{n})-\Phi(u_{n-1}))/\tau$ for
$t\in(t_{n-1},t_{n}]$ and
$L_{\Phi}:=\sup_{u\in[0,\breve{u}]}\\{\Phi^{\prime}(u)\\}<\infty$ from (3.4),
the first term in (3.36) is estimated as

$\displaystyle\left(\frac{1}{\tau}(u_{n}-u_{n-1}),\Phi(u_{n})-\Phi(u_{n-1})\right)$
$\displaystyle\overset{\eqref{eq:uup}}{\geq}\frac{\tau}{\sup\limits_{u\in[0,\breve{u}]}\Phi^{\prime}(u)}\left\|\frac{\Phi(u_{n})-\Phi(u_{n-1})}{\tau}\right\|^{2}$
$\displaystyle=\frac{\tau}{L_{\Phi}}\|\partial_{t}\bar{w}_{\tau}\|^{2}.$
(3.37a) Using the identity $a(a-b)=\frac{1}{2}[a^{2}-b^{2}+(a-b)^{2}]$, the
second-term is estimated as $\displaystyle(\nabla
w_{n},\nabla(w_{n}-w_{n-1}))=$ $\displaystyle\frac{1}{2}[\|\nabla
w_{n}\|^{2}-\|\nabla w_{n-1}\|^{2}+\|\nabla(w_{n}-w_{n-1})\|^{2}],$ (3.37b)
Similarly as in Step 1, using that $\|u_{n}\|_{L^{\infty}(\Omega)}<C$ by
(3.4), $\|v_{n}\|<C$ by (3.34), and that $f$ is a Lipschitz function by (P2),
we have that $\|f(v_{n-1})u_{n}\|<C$. Then, the final term is estimated as
$\displaystyle(f(v_{n-1})u_{n},w_{n}-w_{n-1}))$
$\displaystyle\overset{\eqref{eq:
Young}}{\leq}\frac{L_{\Phi}}{2}\|f(v_{n-1})u_{n}\|^{2}\tau+\frac{\tau}{2L_{\Phi}}\left\|\frac{w_{n}-w_{n-1}}{\tau}\right\|^{2}$
$\displaystyle\leq
C\tau+\frac{\tau}{2L_{\Phi}}\|\partial_{t}\bar{w}_{\tau}\|^{2}$ (3.37c)

Combining the above estimates and summing from $n=1$ to $n=N_{\tau}$ we get

$\displaystyle\frac{1}{L_{\Phi}}\sum_{n=1}^{N_{\tau}}\|\partial_{t}\bar{w}_{\tau}\|^{2}\tau+\|\nabla
w_{N_{\tau}}\|^{2}+\sum_{n=1}^{N_{\tau}}\|\nabla(w_{n}-w_{n-1})\|^{2}\leq
CT+\|\nabla\Phi(u_{0})\|^{2}.$ (3.38)

Since $\hat{w}_{\tau}(t)=w_{n}$ for $t\in(t_{n-1},t_{n}]$, and
$\bar{w}_{\tau}$ is a convex combination of $w_{n}$ and $w_{n-1}$, similarly
as in Step 1, we conclude that $\|\nabla\hat{w}_{t}(t)\|$ and
$\|\nabla\bar{w}_{\tau}(t)\|$ are bounded for all $t\in[0,T]$. Finally,
observing that
$\int_{0}^{T}\|\partial_{t}\bar{w}_{\tau}\|^{2}=\sum_{n=1}^{N_{\tau}}\|\partial_{t}\bar{w}_{\tau}\|^{2}\tau$,
we have (3.30). ∎

####  3.2.2 Convergence to the time-continuous solution if $u_{0}\in
H^{1}_{0}(\Omega)$

We first prove the following result which will be used frequently:

######  Lemma 3.10 (An important convergence result).

Let $\psi\in C^{1}(\mathds{R}^{+})$ be strictly increasing, convex in
$[0,\varepsilon_{0})$ for some $\varepsilon_{0}>0$, and assume that for
$\psi_{\mathrm{m}}:=\inf_{[\varepsilon_{0},\infty)}\psi^{\prime}$ one has
$\psi_{\mathrm{m}}>0$. For a measurable set $\omega\subset\mathds{R}^{d}$, let
$\\{\varphi_{n}\\}_{n\in\mathbb{N}}\subset L^{1}(\omega)$ be a sequence of
non-negative functions such that
$\|\psi(\varphi_{n})-\psi(\varphi)\|_{L^{1}(\omega)}\to 0$ for a fixed (non-
negative) $\varphi\in L^{1}(\omega)$. Then,
$\|\varphi_{n}-\varphi\|_{L^{1}(\omega)}\to 0$ as $n\to\infty$.

###### Proof.

Let $\bar{\Psi}\in C(\mathds{R}^{+})$ be defined as

$\bar{\Psi}(\varphi)=\begin{cases}\psi(\varphi)-\psi(0)&\text{ for
}\varphi\in[0,\varepsilon_{0}),\\\
\psi(\varepsilon_{0})-\psi(0)+\psi^{\prime}(\varepsilon_{0})(\varphi-\varepsilon_{0})&\text{
for }\varphi\geq\varepsilon_{0}.\end{cases}$

It is straightforward to verify that $\bar{\Psi}$ is convex, strictly
increasing, $\bar{\Psi}(0)=0$, and

$\displaystyle|\bar{\Psi}(\varphi_{1})-\bar{\Psi}(\varphi_{2})|\leq(\psi^{\prime}(\varepsilon_{0})/\psi_{\mathrm{m}})|\psi(\varphi_{1})-\psi(\varphi_{2})|\text{
for all }\varphi_{1/2}\geq 0.$ (3.39)

The inequality above follows from considering separately the cases
$\varphi_{1/2}\leq\varepsilon_{0}$ which gives
$|\bar{\Psi}(\varphi_{1})-\bar{\Psi}(\varphi_{2})|=|\psi(\varphi_{1})-\psi(\varphi_{2})|$;
$\varphi_{1/2}\geq\varepsilon_{0}$ which gives
$|\bar{\Psi}(\varphi_{1})-\bar{\Psi}(\varphi_{2})|=\psi^{\prime}(\varepsilon_{0})|\varphi_{1}-\varphi_{2}|\leq(\psi^{\prime}(\varepsilon_{0})/\psi_{\mathrm{m}})|\psi(\varphi_{1})-\psi(\varphi_{2})|$;
and $\varphi_{1}$, $\varphi_{2}$ being on different sides of $\varepsilon_{0}$
which gives also (3.39) by combining the estimates for the other two cases.
Moreover, using the super-additivity property (2.7b) one has for
$\varphi_{n}>\varphi$ that
$\bar{\Psi}(\varphi_{n}-\varphi)\leq\bar{\Psi}(\varphi_{n})-\bar{\Psi}(\varphi)$,
and by symmetry, we conclude that
$\bar{\Psi}(|\varphi_{n}-\varphi|)\leq|\bar{\Psi}(\varphi_{n})-\bar{\Psi}(\varphi)|$.
Consequently,

$\begin{split}\bar{\Psi}\left(\frac{1}{|\omega|}\int_{\omega}|\varphi_{n}-\varphi|\right)&\overset{\eqref{eq:
jensen}}{\leq}\frac{1}{|\omega|}\int_{\omega}\bar{\Psi}\left(|\varphi_{n}-\varphi|\right)\leq\frac{1}{|\omega|}\int_{\omega}|\bar{\Psi}(\varphi_{n})-\bar{\Psi}(\varphi)|\\\
&\leq C\|\psi(\varphi_{n})-\psi(\varphi)\|_{L^{1}(\omega)}\longrightarrow
0.\end{split}$

Since $\bar{\Psi}$ is strictly increasing, it follows that
$\|\varphi_{n}-\varphi\|_{L^{1}(\omega)}\to 0$. ∎

The above result has previously been used in Lemma 3.3 of [29] to prove strong
convergence of solutions, see also [23]. Here, we use it in a similar way.

######  Proof of (3.5) in Theorem 3.2.

Observe that $\bar{w}_{\tau}\in H^{1}(Q)$ is uniformly bounded with respect to
$\tau$ if $u_{0}\in H^{1}_{0}(\Omega)$ by Lemma 3.9 since $\Phi(u_{0})\in
H^{1}_{0}(\Omega)$ in this case. Hence, by the compact embedding
$H^{1}(Q)\hookrightarrow L^{2}(Q)$, there exists $w\in H^{1}(Q)$ such that
along a sub-sequence of $\tau$ converging to $0$,

$\displaystyle\bar{w}_{\tau}\rightharpoonup w\ \text{ weakly in }H^{1}(Q),$
(3.40a) $\displaystyle\bar{w}_{\tau}\longrightarrow w\ \text{ strongly in
}L^{2}(Q).$ (3.40b) Define $u:=\Phi^{-1}(w)$ which is bounded in
$[0,\breve{u}]$ a.e. in $\Omega$ for all $t>0$ due to (3.4). We will prove
that $\displaystyle\hat{w}_{\tau}\longrightarrow w\ \text{ strongly in
}L^{2}(Q)$ (3.40c) $\displaystyle\hat{u}_{\tau}\longrightarrow u\ \text{
strongly in }L^{2}(Q).$ (3.40d)

The convergence (3.40c) follows from (3.28) and (3.40b) since

$\displaystyle\int_{0}^{T}\|\hat{w}_{\tau}-\hat{w}_{\tau}\|^{2}$
$\displaystyle\overset{\eqref{eq:time_interpolants}}{=}\sum_{n=0}^{N}\int_{t_{n-1}}^{t_{n}}\left(\frac{t-t_{n-1}}{\tau}\right)^{2}\|w_{n}-w_{n-1}\|^{2}=\frac{1}{3}\sum_{n=0}^{N}\|w_{n}-w_{n-1}\|^{2}\tau$
$\displaystyle\overset{\eqref{eq:Poincare}}{\leq}\frac{C_{\Omega}\tau}{3}\sum_{n=0}^{N}\|\nabla(w_{n}-w_{n-1})\|^{2}\overset{\eqref{eq:w-disc-
bound}}{\leq}C\tau\longrightarrow 0.$

To show (3.40d), noting that $\Phi(\hat{u}_{\tau})=\hat{w}_{\tau}$, we have

$\|\Phi(\hat{u}_{\tau})-\Phi(u)\|_{L^{2}(Q)}=\|\hat{w}_{\tau}-w\|^{2}_{L^{2}(Q)}\longrightarrow
0,$

which also implies that $\|\Phi(\hat{u}_{\tau})-\Phi(u)\|_{L^{1}(Q)}\to 0$.
Hence, using Lemma 3.10 with $\psi=\Phi$ gives that
$\|\hat{u}_{\tau}-u\|_{L^{1}(Q)}\to 0$ and since both $\hat{u}_{\tau},u\in
L^{\infty}(Q)$, we have (3.40d).

For the convergence of $v$, note that if $\mu=1$ then (3.31) implies that
$\bar{v}_{\tau}\in H^{1}(0,T;H^{-1}(\Omega))\cap
L^{2}(0,T;H^{1}_{0}(\Omega))=:{\cal W}$ is uniformly bounded with respect to
$\tau$. The space ${\cal W}$ is compactly embedded into $L^{2}(Q)$ and
continuously into $C([0,T];L^{2}(\Omega))$ (Aubin-Lions lemma). Hence, for
$\mu=1$, there exists $v\in{\cal W}\subset C([0,T];L^{2}(\Omega))$ such that

$\displaystyle\bar{v}_{\tau}\longrightarrow v\ \text{ strongly in }L^{2}(Q),$
(3.41a) $\displaystyle\hat{v}_{\tau}\longrightarrow v\ \text{ strongly in
}L^{2}(Q),$ (3.41b)

For $\mu=0$, let $v\in C([0,T];L^{2}(\Omega))$ be the solution of
$\partial_{t}v=g(u,v)$ with $v(0)=v_{0}$. Then,

$\displaystyle\frac{1}{2}\|(\bar{v}_{\tau}-v)(T)\|^{2}$
$\displaystyle=\int_{0}^{T}(\bar{v}_{\tau}-v,\partial_{t}(\bar{v}_{\tau}-v))\leq\|\bar{v}_{\tau}-v\|_{L^{2}(Q)}\|\partial_{t}(\bar{v}_{\tau}-v)\|_{L^{2}(Q)}$
$\displaystyle\overset{\eqref{eq:Uni_bnd_v}}{\leq}C\|\partial_{t}(\bar{v}_{\tau}-v)\|_{L^{2}(Q)}.$

Using that
$\partial_{t}\bar{v}_{\tau}(t)=(v_{n}-v_{n-1})/\tau=g(u_{n},v_{n-1})=g(\hat{u}_{\tau}(t),v_{n-1})$
for $t\in(t_{n-1},t_{n}]$, one further estimates

$\displaystyle\|\partial_{t}(\bar{v}_{\tau}-v)\|^{2}_{L^{2}(Q)}=\sum_{n=1}^{N_{\tau}}\int_{t_{n-1}}^{t_{n}}\|g(\hat{u}_{\tau},v_{n-1})-g(u,v)\|^{2}$
$\displaystyle\overset{\ref{ass:
2}}{\leq}C\sum_{n=1}^{N_{\tau}}\int_{t_{n-1}}^{t_{n}}\left[\|\hat{u}_{\tau}-u\|^{2}+\|v_{n-1}-v\|^{2}\right]$
$\displaystyle\leq
C\left(\|\hat{u}_{\tau}-u\|^{2}_{L^{2}(Q)}+\int_{0}^{T}\|\bar{v}_{\tau}-v\|^{2}+\sum_{n=1}^{N_{\tau}}\int_{t_{n-1}}^{t_{n}}\|v_{n-1}-\bar{v}_{\tau}\|^{2}\right).$

Note that $\|\hat{u}_{\tau}-u\|^{2}_{L^{2}(Q)}\to 0$ from (3.40d), and

$\begin{split}\sum_{n=1}^{N_{\tau}}\int_{t_{n-1}}^{t_{n}}\|v_{n-1}-\bar{v}_{\tau}\|^{2}&\overset{\eqref{eq:time_interpolants}}{=}\sum_{n=1}^{N_{\tau}}\int_{t_{n-1}}^{t_{n}}\left(\frac{t-t_{n-1}}{\tau}\right)^{2}\|v_{n}-v_{n-1}\|^{2}\\\
&=\frac{\tau}{3}\sum_{n=1}^{N_{\tau}}\|v_{n}-v_{n-1}\|^{2}\overset{\eqref{eq:v-disc-
bound}}{\leq}C\tau\longrightarrow 0,\end{split}$

Hence, applying Gronwall’s lemma we get that $\|(\bar{v}_{\tau}-v)(T)\|\to 0$
which proves the strong convergence result in (3.41a). The convergence of
$\hat{v}_{\tau}$ in (3.41b) follows from (3.34) similar to (3.40c).

It is straightforward to show that $(u,v)$ is indeed a weak solution of (1.1),
a detailed proof can be found in Theorem 3.1 in [27]. Since
$(\hat{u}_{\tau},\hat{w}_{\tau},\hat{w}_{\tau})$ is bounded uniformly
componentwise in $L^{2}(Q)$ for $\tau$ small, and every converging subsequence
of it converges to the unique limit $(u,v,w)$ weakly solving (1.1), along
every sequence of $\tau\to 0$ this limit is obtained. ∎

####  3.2.3 Convergence to the time-continuous solution if $u_{0}\not\in
H^{1}_{0}(\Omega)$

For less regular initial data we need to use the $L^{1}$-contraction
principle, see [31] for the general idea, and [23] for a proof for this
specific case.

######  Lemma 3.11 ($L^{1}$-contraction principle).

Let $(u_{1},v_{1})$ and $(u_{2},v_{2})$ be the weak solutions of (1.1)
corresponding to the initial data $u_{1}(0)=u_{1,0}$ and $u_{2}(0)=u_{2,0}$
and let $u_{1,0},u_{2,0}$ satisfy (P4). Then, for any $t>0$

$\displaystyle\|(u_{1}-u_{2})(t)\|_{L^{1}(\Omega)}\leq\|u_{1,0}-u_{2,0}\|_{L^{1}(\Omega)}+\int_{0}^{t}\|f(v_{1})u_{1}-f(v_{2})u_{2}\|_{L^{1}(\Omega)}.$
(3.42)

######  Lemma 3.12 (Convergence of the continuous solutions as $\varepsilon\to
0$).

Let $(u,v)$ and $(u^{\varepsilon},v^{\varepsilon})$ be the weak solutions of
(1.1) corresponding to the initial conditions $u(0)=u_{0}$ and
$u^{\varepsilon}(0)=u_{0}^{\varepsilon}$, where $u_{0}^{\varepsilon}$ is as in
Theorem 3.2, and let $w=\Phi(u),w^{\varepsilon}=\Phi(u^{\varepsilon})$. Then,
for any $t>0$, along any sequence of $\varepsilon$ converging to 0 we have

$\displaystyle\|(u^{\varepsilon}-u)(t)\|_{L^{2}(\Omega)}+\|(w^{\varepsilon}-w)(t)\|_{L^{2}(\Omega)}+\|v^{\varepsilon}-v\|_{L^{2}(Q)}\longrightarrow
0.$ (3.43)

###### Proof.

Observe that the uniform bound in (3.31) holds also for $z=v^{\varepsilon}$
with the constant $\bar{C}$ independent of $\varepsilon$. Hence, similar to
(3.41a), along a subsequence of $\varepsilon\to 0$, one has
$\|v^{\varepsilon}-v\|_{L^{2}(Q)}\to 0$. Moreover, noting that $0\leq
u(t),\,u^{\varepsilon}(t)\leq\breve{u}<C$ a.e. in $\Omega$ due to (3.29), one
has by Lemma 3.11 that

$\displaystyle\|(u^{\varepsilon}-u)(t)\|_{L^{1}(\Omega)}$
$\displaystyle\leq\|u^{\varepsilon}_{0}-u_{0}\|_{L^{1}(\Omega)}+\int_{0}^{t}\|f(v^{\varepsilon})u^{\varepsilon}-f(v)u\|_{L^{1}(\Omega)}$
$\displaystyle\leq|\Omega|^{\frac{1}{2}}\|u^{\varepsilon}_{0}-u_{0}\|+\int_{0}^{t}\|(f(v^{\varepsilon})-f(v))u^{\varepsilon}\|_{L^{1}(\Omega)}+\int_{0}^{t}\|f(v)(u^{\varepsilon}-u)\|_{L^{1}(\Omega)}$
$\displaystyle\overset{\ref{ass:
2}}{\leq}\varepsilon|\Omega|^{\frac{1}{2}}\overset{|u^{\varepsilon}|<C}{+}C\int_{0}^{t}\|v^{\varepsilon}-v\|_{L^{1}(\Omega)}+f_{\mathrm{M}}\int_{0}^{t}\|u^{\varepsilon}-u\|_{L^{1}(\Omega)}.$
(3.44)

Applying Gronwall’s lemma (2.6) along with
$\|v^{\varepsilon}-v\|_{L^{1}(Q)}\to 0$ we get that
$\|(u^{\varepsilon}-u)(t)\|_{L^{1}(\Omega)}\to 0$ for all $t>0$, which further
implies that $\|(u^{\varepsilon}-u)(t)\|_{L^{2}(\Omega)}\to 0$ since $0\leq
u(t),\,u^{\varepsilon}(t)\leq\breve{u}<C$. It also implies that
$\|\Phi(u^{\varepsilon})-\Phi(u)\|_{L^{2}(\Omega)}\to 0$ since $\Phi$ is
Lipschitz in $[0,\breve{u}]$. This proves the result. ∎

######  Proof of (3.6) in Theorem 3.2.

We choose $\varepsilon>0$ small enough such that along the subsequence in
Lemma 3.12 we have

$\displaystyle\int_{0}^{T}\left[\|u^{\varepsilon}-u\|^{2}+\|w^{\varepsilon}-w\|^{2}+\|v^{\varepsilon}-v\|^{2}\right]\leq\tfrac{1}{2}\delta.$
(3.45a) for some arbitrary $\delta>0$. For this fixed $\varepsilon>0$, noting
that $u^{\varepsilon}_{0}\in H^{1}_{0}(\Omega)$, one can choose a time-step
$\tau>0$ small enough such that by (3.5) one has
$\displaystyle\sum_{n=0}^{N_{\tau}}\int_{t_{n-1}}^{t_{n}}\left[\|u_{n}^{\varepsilon}-u^{\varepsilon}(t)\|^{2}+\|w^{\varepsilon}_{n}-w^{\varepsilon}(t)\|^{2}+\|v^{\varepsilon}_{n}-v^{\varepsilon}(t)\|^{2}\right]\mathrm{d}t\leq\tfrac{1}{2}\delta.$
(3.45b) Combining these estimates, one finds the desired subsequence
$(\varepsilon,\tau)\to(0,0)$ such that (3.6) holds.

∎

##  4 Linearisation

We have shown that the time-discretised system (3.1) is well-posed and that
its solutions possess the qualitative behaviour we expect from the time-
continuous system. In this section, we propose linearisation schemes and prove
their well-posedness and convergence. Recall that $\Phi$ is possibly not
Lipschitz continuous if $b=1$. However, $u_{n}$ takes values in
$[0,\breve{u}]$ where $\breve{u}<b$ is a uniform a priori computable upper
bound (see Remark 3.1.1), and $\Phi$ is Lipschitz in $[0,\breve{u}]$. Hence,
we can regularize $\Phi$ as follows.

######  Definition 4.1 (Regularization of $\Phi$).

If $b=1$ and with $\breve{u}>0$ given in (3.4), let the function
$\breve{\Phi}\colon\mathds{R}^{+}\to\mathds{R}^{+}$ be defined as

$\begin{split}\breve{\Phi}(u)=\begin{cases}\Phi(u),&\text{if
}u\leq\breve{u},\\\
\Phi^{\prime}(\breve{u})(u-\breve{u})+\Phi(\breve{u}),&\text{if
}u\geq\breve{u}.\end{cases}\end{split}$ (4.1)

If $b=\infty$, we set $\breve{\Phi}=\Phi$.

Recalling that $\Phi(u)=w$ posseses space regularity, we propose an iterative
linearisation scheme to solve (3.1a) which splits the equation into two
coupled equations. The iterations are obtained by solving the following.

######  Problem (The splitting linearisation).

Let $n\in\mathbb{N}$ and $i\in\mathbb{N}_{0}$ be fixed, and assume
$u_{n-1},v_{n-1}\in L^{2}(\Omega)$ and $(u^{i-1}_{n},w^{i-1}_{n})\in
L^{2}(\Omega)\times H^{1}_{0}\left(\Omega\right)$ be given, satisfying
$u_{n-1},u^{i-1}_{n}\leq\breve{u}$. Find the pair $(u^{i}_{n},w^{i}_{n})\in
L^{2}(\Omega)\times H^{1}_{0}\left(\Omega\right)$ such that, for all $\phi\in
H^{1}_{0}\left(\Omega\right)$ and $\xi\in L^{2}(\Omega)$ it holds that

$\displaystyle\left(\frac{1}{\tau}\left(\tilde{u}^{i}_{n}-u_{n-1}\right),\phi\right)+\left(\nabla
w^{i}_{n},\nabla\phi\right)=\left(f(v_{n-1})\tilde{u}^{i}_{n},\phi\right),$
(4.2a)
$\displaystyle\left(L^{i}_{n}(\tilde{u}^{i}_{n}-u^{i-1}_{n}),\xi\right)=\left(w^{i}_{n}-\breve{\Phi}\left(u^{i-1}_{n}\right),\xi\right),$
(4.2b) $\displaystyle\ u^{i}_{n}=\left[\tilde{u}^{i}_{n}\right]_{+}\text{ a.e.
in }\Omega,$ (4.2c)

for some specific choice of a bounded function
$L^{i}_{n}:\Omega\to\mathds{R}^{+}$, which depends only on iterates up to
$u^{i-1}_{n}$ but not on $u^{i}_{n}$. The iteration starts with the initial
guess $u^{0}_{n}=u_{n-1}$.

Such a splitting method was first proposed in [7, Section 4.2] for the
$L$-scheme assuming that $\Phi$ is Lipschitz. Here, we generalize the results.

######  Remark 4.1.1 (Positivity of $u^{i}_{n}$).

To shorten the proofs, throughout this section, we will simplify (4.2) (where
we first determine $\tilde{u}^{i}_{n}$ and then set
$u^{i}_{n}=\left[\tilde{u}^{i}_{n}\right]_{+}$) by referring to
$\tilde{u}^{i}_{n}$ interchangeably as $u^{i}_{n}$. All inequalities and
results in this section remain valid since

$\norm{u^{i}_{n}-u_{n}}_{L^{p}(\Omega)}=\norm{\left[\tilde{u}^{i}_{n}\right]_{+}-u_{n}}_{L^{p}(\Omega)}\leq\norm{\tilde{u}^{i}_{n}-u_{n}}_{L^{p}(\Omega)},$
(4.3)

for all $p\geq 1$. Indeed, as $u_{n}\geq 0$ which gives
$\left[\tilde{u}^{i}_{n}\right]_{+}-u_{n}=\tilde{u}^{i}_{n}-u_{n}$ if
$\tilde{u}^{i}_{n}\geq 0$ and
$\tilde{u}^{i}_{n}-u_{n}<\left[\tilde{u}^{i}_{n}\right]_{+}-u_{n}\leq 0$ if
$\tilde{u}^{i}_{n}<0$. The reason we introduce the formulation (4.2) is that
it guarantees that $u^{i}_{n}\geq 0$, which is important for the numerical
implementation.

Before proving results for particular linearisation schemes, we show that the
regularization $\breve{\Phi}$ does not alter the solution $u_{n}$. This is
obvious for $b=\infty$, while for $b=1$ it follows from the proposition below.

######  Proposition 4.2 (Consistency of the regularized $\breve{\Phi}$).

Let $b=1$, and $\breve{\Phi}$ the regularized approximation of $\Phi$ given in
Definition 4.1. Then, the solution of (4.2) coincides with the solution to
(3.1).

###### Proof.

Suppose $u_{n}$ and $\tilde{u}_{n}$ are the weak solution of

$\displaystyle\frac{1}{\tau}(\tilde{u}_{n}-u_{n-1})$
$\displaystyle=\Delta\breve{\Phi}(\tilde{u}_{n})+f(v_{n-1})\tilde{u}_{n},$
(4.4) $\displaystyle\frac{1}{\tau}(u_{n}-u_{n-1})$
$\displaystyle=\Delta\Phi(u_{n})+f(v_{n-1})u_{n}.$ (4.5)

Since $u_{n}\leq\breve{u}$ by Theorem 3.1, one has
$\breve{\Phi}(u_{n})=\Phi(u_{n})$, i.e. $u_{n}$ is a solution of Equation
(4.4). However, this solution is unique due to Theorem 3.1 which implies that
$\tilde{u}_{n}=u_{n}$. ∎

To show that the linearisation scheme is well-defined, we prove that if it
converges, the limit is indeed a solution of the time-discretised equation
(3.1a).

######  Proposition 4.3 (Consistency of the linearisation scheme).

Let $u^{i}_{n}$ be uniformly bounded with respect to $i\in\mathbb{N}$ in
$L^{2}(\Omega)$, $u^{i}_{n}\to\tilde{u}_{n}$ strongly in $L^{1}(\Omega)$, and
$w^{i}_{n}\to\tilde{w}_{n}$ strongly in $H^{1}_{0}(\Omega)$. Then
$\tilde{u}_{n}$ is the weak solution to the time-discretised equation (3.1a)
and $\tilde{w}_{n}=\Phi(\tilde{u}_{n})$ a.e. in $\Omega$.

Theorems 4.4 and 4.5 will show that the hypotheses of Proposition 4.3 are
indeed satisfied for the L- and M-schemes. Hence, the iterates
$(u^{i}_{n},w^{i}_{n})$ converge to the time-discrete solutions.

###### Proof.

First we observe that $(u_{n}^{i},\phi)\to(\tilde{u}_{n},\phi)$ for all
$\phi\in L^{2}(\Omega)$, since $u^{i}_{n}$ is bounded in $L^{2}(\Omega)$ and
$u^{i}_{n}\to\tilde{u}_{n}$ strongly in $L^{1}(\Omega)$. Hence, taking the
limit in (4.2a) implies that

$\left(h_{n-1}\tilde{u}_{n}-u_{n-1},\phi\right)+\tau\left(\nabla\tilde{w}_{n},\nabla\phi\right)=0\qquad\forall\phi\in
H^{1}_{0}\left(\Omega\right),$ (4.6)

where $h_{n-1}$ is defined in (3.2). Similarly, taking the limit $i\to\infty$
in (4.2b) we get

$\left(\tilde{w}_{n}-\breve{\Phi}(\tilde{u_{n}}),\xi\right)=0\qquad\forall\xi\in
L^{2}(\Omega).$ (4.7)

Here, we used that $\breve{\Phi}$ is Lipschitz continuous implying that
$(\breve{\Phi}(u^{i}_{n}),\xi)\to(\breve{\Phi}(\tilde{u}_{n}),\xi)$, and that
$L^{i}_{n}$ is bounded which yields
$\left(L^{i}_{n}\left(u^{i}_{n}-u^{i-1}_{n}\right),\xi\right)\to 0$ for all
$\xi\in L^{2}(\Omega)$. We conclude that
$\tilde{w}_{n}=\breve{\Phi}(\tilde{u}_{n})$ a.e., which allows us to
substitute it back into Equation (4.6) and hence,

$\left(h_{n-1}\tilde{u}_{n}-u_{n-1},\phi\right)+\tau\left(\nabla\breve{\Phi}\left(\tilde{u}_{n}\right),\nabla\phi\right)=0\qquad\forall\phi\in
H^{1}_{0}\left(\Omega\right).$ (4.8)

This coincides with the time-discretised equation for $u_{n}$, and thus
$\tilde{u}_{n}=u_{n}$, as solutions are unique. ∎

We can identify the linearisation schemes mentioned in Section 1 as special
cases of (4.2):

$\displaystyle\textbf{Newton
scheme}:L^{i}_{n}:=\breve{\Phi}^{\prime}(u^{i-1}_{n}),$ (4.9a)
$\displaystyle\textbf{L-scheme}:\qquad\quad L^{i}_{n}:=L,\quad\text{for a
constant }L>0,$ (4.9b)
$\displaystyle\textbf{M-scheme}:\qquad\;\;L^{i}_{n}:=\max\\{\breve{\Phi}^{\prime}(u^{i-1}_{n})+M\tau^{\gamma},\,2M\tau^{\gamma}\\},\text{
for constants }M>0,\;\gamma\in(0,1].$ (4.9c)

In the sequel, we consider the L- and M-schemes as our main focus will be on
degenerate problems. We denote the errors of the iterates at the
$n^{\text{th}}$ time step by

$\displaystyle e^{i}_{u}=u^{i}_{n}-u_{n},\quad e^{i}_{w}=w^{i}_{n}-w_{n},$
(4.10)

where $u_{n},w_{n},u_{n}^{i},w_{n}^{i}$ are the solutions of (1.3) and (1.4)
respectively.

The following theorems provide the main convergence results for the L- and
M-scheme, their proofs are given in Subsections 4.1 and 4.2. The results for
both schemes are similar, but the proofs for the L-scheme are more
straightforward.

######  Theorem 4.4 (Convergence of the L-scheme).

For $\tau<1/f_{\mathrm{M}}$ there exist unique solutions
$\\{(u^{i}_{n},w^{i}_{n})\\}_{i\in\mathbb{N}}\subset L^{2}(\Omega)\times
H^{1}_{0}\left(\Omega\right)$ of (4.2) with (4.9b), i.e. $L^{i}_{n}:=L$.
Furthermore, if $L>\sup\breve{\Phi}^{\prime}$, then $u^{i}_{n}\to u_{n}$ in
$L^{1}(\Omega)$ and $w^{i}_{n}\to w_{n}$ in $H^{1}_{0}\left(\Omega\right)$.
For non-degenerate problems, i.e. if $\inf\Phi^{\prime}=\phi_{\mathrm{m}}>0$,
the error-norm is a strict contraction

$\|(e^{i}_{u},e^{i}_{w})\|_{L}\leq\alpha\|(e^{i-1}_{u},e^{i-1}_{w})\|_{L}$

with rate $\alpha=\sqrt{\frac{L}{L+\phi_{\mathrm{m}}}},$ where

$\norm{\left(e^{i}_{u},e^{i}_{w}\right)}_{L}^{2}\coloneqq\int_{\Omega}h_{n-1}|e^{i}_{u}|^{2}+\frac{2\tau}{L+\phi_{\mathrm{m}}}\norm{\nabla
e^{i}_{w}}^{2}_{L^{2}(\Omega)}.$ (4.11)

For the M-scheme, we need to impose an additional assumption.

  1. (A1)

For a given $n\in\mathbb{N}$, there exists $\Lambda\geq 0$ and
$\gamma\in(0,1]$ such that
$\|u_{n}-u_{n-1}\|_{L^{\infty}(\Omega)}\leq\Lambda\tau^{\gamma}$.

######  Remark 4.4.1 (Assumption (A1)).

Assumption (A1) was used in [28] with $\gamma=1$ in the context of nonlinear
diffusion problems, and this property was proven for a particular case in
Proposition 3.1, but not for porous medium type diffusion. Note that (A1) with
$\gamma=1$ is the time-discrete counterpart of the regularity assumption
$\partial_{t}u\in L^{\infty}(\Omega)$. But for degenerate problems this is
typically not satisfied. However, solutions of porous medium type equations
are Hölder continuous, and for degenerate and singular systems of the form
(1.1), the Hölder continuity of solutions was shown in [22]. Hence, Assumption
(A1) is expected to hold as a time-discrete counterpart of the Hölder
continuity with exponent $\gamma\in(0,1]$.

######  Theorem 4.5 (Convergence of the M-scheme).

For $\tau<1/f_{\mathrm{M}}$ there exist unique solutions
$\\{(u^{i}_{n},w^{i}_{n})\\}_{i\in\mathbb{N}}\subset L^{2}(\Omega)\times
H^{1}_{0}\left(\Omega\right)$ of (4.2) (4.2) with (4.9c). Furthermore, assume
that (A1) holds, take $M>M_{0}\coloneqq\|\breve{\Phi}^{\prime}\|_{\rm
Lip}\Lambda$, and let $\\{u_{n}^{i}\\}_{i\in\mathbb{N}}$ satisfy

$\displaystyle\|u^{i}_{n}-u_{n}\|_{L^{\infty}(\Omega)}\leq\Lambda\tau^{\gamma}\text{
for all }i\in\mathbb{N}.$ (4.12)

Then, $u^{i}_{n}\to u_{n}$ in $L^{1}(\Omega)$ and $w^{i}_{n}\to w_{n}$ in
$H^{1}_{0}\left(\Omega\right)$.
For non-degenerate problems, i.e. if $\inf\Phi^{\prime}=\phi_{\mathrm{m}}>0$,
and if $\tau<\left(\phi_{\mathrm{m}}/M\right)^{\frac{1}{\gamma}}$, then the
error-norm is a strict contraction,

$\|(e^{i}_{u},e^{i}_{w})\|_{M}\leq\alpha\|(e^{i-1}_{u},e^{i-1}_{w})\|_{M},$

with rate $\alpha=\frac{2M\tau^{\gamma}}{\phi_{\mathrm{m}}+M\tau^{\gamma}},$
where

$\norm{\left(e^{i}_{u},e^{i}_{w}\right)}_{M}^{2}\coloneqq\int_{\Omega}h_{n-1}|e^{i}_{u}|^{2}+\frac{2\tau}{\phi_{\mathrm{m}}+M\tau^{\gamma}}\norm{\nabla
e^{i}_{w}}^{2}_{L^{2}(\Omega)}.$ (4.13)

For the convergence of $(u^{i}_{n},w^{i}_{n})$ in the $L^{\infty}$-norm, see
Proposition 4.9.

######  Remark 4.5.1 (Boundedness condition (4.12) and contraction in
$L^{\infty}$).

In the non-degenerate case, i.e. if $\phi_{\mathrm{m}}>0$, the boundedness
condition (4.12) follows from (A1) for time-step sizes
$\tau\leq(\phi_{\mathrm{m}}/3M)^{\frac{1}{\gamma}}$, as stated in Proposition
4.9. In fact, Proposition 4.9 even provides linear convergence of $u^{i}_{n}$
to $u_{n}$ in $L^{\infty}(\Omega)$ with a contraction rate that scales with
$\tau^{\gamma}$. For the case of singular diffusion, a proof of (4.12) was
given in [28, Lemma 3.1]. We expect that the result also holds in our case.
However. since it is not the main focus of this work, we state it as an
assumption.

######  Remark 4.5.2 (Comparison L- and M-scheme).

Note that the extra assumptions (A1) and (4.12) are not required for the
L-scheme, and hence, the L-scheme is expected to be more robust than the
M-scheme. However, this comes at the cost of being considerably slower than
the Newton scheme. On the other hand, the assumptions required for M-scheme
are expected to hold for problems such as (1.1). Since the contraction rate
for the M-scheme scales with $\tau$, for practical purposes the M-scheme
results in a more competitive iterative solver than the L-scheme.

We first prove the existence and uniqueness results stated in Theorems 4.4 and
4.5. Recall that $L^{i}_{n}=L$ is constant for the L-scheme and
$\sup\breve{\Phi}^{\prime}$ is bounded due to assumption (P1) and the
construction of $\breve{\Phi}$ in Definition 4.1. The proof of the following
lemma applies to both schemes.

######  Lemma 4.6 (Existence-uniqueness).

For $\tau<1/f_{M}$, the system of equations (4.2) with

$L^{i}_{n}\coloneqq L>\sup\breve{\Phi}^{\prime}\qquad\text{or}\qquad
L^{i}_{n}\coloneqq\max\\{\breve{\Phi}^{\prime}(u^{i-1}_{n})+M\tau^{\gamma},\,2M\tau^{\gamma}\\}$

has a unique solution.

###### Proof.

We eliminate $\tilde{u}^{i}_{n}$ in (4.2a) through (4.2b) and find

$\left(\frac{h_{n-1}}{L^{i}_{n}}w_{n}^{i},\phi\right)+\tau(\nabla
w_{n}^{i},\nabla\phi)=(g_{n}^{i},\phi)\quad\text{ for all $\phi\in
H_{0}^{1}(\Omega)$},$ (4.14)

where
$g_{n}^{i}=\frac{h_{n-1}}{L^{i}_{n}}\breve{\Phi}(u^{i-1}_{n})-(h_{n-1}u^{i-1}_{n}-u_{n-1})$.
Consider the bilinear form
$B(w,\phi)=((h_{n-1}/L^{i}_{n})\,w,\phi)+\tau(\nabla w,\nabla\phi)$ and the
linear functional $l(\phi)=(g_{n}^{i},\phi)$. We observe that $L^{i}_{n}$ is
constant, or bounded from above and below by positive constants in case of the
M-scheme, see (4.9), and $0<h_{n-1}<1$ due to $\tau<1/f_{\mathrm{M}}$. Hence,
using the Cauchy-Schwarz and Poincaré inequality implies that $B$ is coercive
and bounded and $l$ is a bounded linear functional on $H_{0}^{1}(\Omega)$. The
Lax-Milgram theorem now provides the existence and uniqueness of a solution
$w^{i}_{n}\in H^{1}_{0}\left(\Omega\right)$. The existence and uniqueness of
$\tilde{u}_{n}^{i}\in L^{2}(\Omega)$ then follows from (4.2b), while
$u^{i}_{n}$ can be found through (4.2c). ∎

###  4.1 L-scheme

First, we show that the solutions of the L-scheme converge to the time-
discrete solutions $u_{n}$ and $w_{n}$.

######  Lemma 4.7 (Convergence of the L-scheme).

Under the assumptions of Theorem 4.4, the stated convergence results hold.

###### Proof.

We use ideas from the proof of Lemma 2.6 in [7]. Subtracting (3.1a) from
Equation (4.2a) we find

$\left(h_{n-1}e^{i}_{u},\phi\right)+\tau\left(\nabla
e^{i}_{w},\nabla\phi\right)=0,$ (4.15)

where $e^{i}_{u}=u^{i}_{n}-u_{n}$ and $e^{i}_{w}=w^{i}_{n}-w_{n}$, see (4.10).
By adding and subtracting $Lu_{n}$, and adding and subtracting
$w_{n}=\breve{\Phi}(u_{n})$ on the right hand side of (4.2b), we can rewrite
it as

$\left(L\left(e^{i}_{u}-e^{i-1}_{u}\right),\xi\right)=\left(e^{i}_{w}-\delta\breve{\Phi}^{i-1},\xi\right)$
(4.16)

where
$\delta\breve{\Phi}^{i-1}\coloneqq\breve{\Phi}\left(u^{i-1}_{n}\right)-\breve{\Phi}(u_{n})$.
Choosing $\phi=e^{i}_{w}\in H^{1}_{0}\left(\Omega\right)$ and
$\xi=h_{n-1}e^{i}_{u}\in L^{2}(\Omega)$, we combine the two equations and
obtain

$\left(h_{n-1}L\left(e^{i}_{u}-e^{i-1}_{u}\right),e^{i}_{u}\right)+\left(h_{n-1}\delta\breve{\Phi}^{i-1},e^{i}_{u}\right)+\tau\norm{\nabla
e^{i}_{w}}_{L^{2}(\Omega)}^{2}=0.$ (4.17)

This is rewritten as

$\left(h_{n-1}\left(L-\left(\frac{\delta\breve{\Phi}}{\delta
u}\right)^{i-1}\right)\left(e^{i}_{u}-e^{i-1}_{u}\right),e^{i}_{u}\right)+\left(h_{n-1}\left(\frac{\delta\breve{\Phi}}{\delta
u}\right)^{i-1}e^{i}_{u},e^{i}_{u}\right)+\tau\norm{\nabla
e^{i}_{w}}_{L^{2}(\Omega)}^{2}=0,$ (4.18)

where

$\left(\frac{\delta\breve{\Phi}}{\delta
u}\right)^{i-1}\coloneqq\frac{\breve{\Phi}\left(u^{i-1}_{n}\right)-\breve{\Phi}(u_{n})}{u^{i-1}_{n}-u_{n}}=\frac{\delta\breve{\Phi}^{i-1}}{e^{i-1}_{u}}.$
(4.19)

Using the identity
$(a-b)a=\frac{1}{2}\left(a^{2}-b^{2}+\left(a-b\right)^{2}\right)$ with
$a=e^{i}_{u}$ and $b=e^{i-1}_{u}$ we rewrite the first term in (4.18) and
obtain

$\begin{split}&\frac{1}{2}\int_{\Omega}h_{n-1}\left(L+\left(\frac{\delta\breve{\Phi}}{\delta
u}\right)^{i-1}\right)|e^{i}_{u}|^{2}+\frac{1}{2}\int_{\Omega}h_{n-1}\left(L-\left(\frac{\delta\breve{\Phi}}{\delta
u}\right)^{i-1}\right)|e^{i}_{u}-e^{i-1}_{u}|^{2}\\\ &+\tau\norm{\nabla
e^{i}_{w}}_{L^{2}(\Omega)}^{2}=\frac{1}{2}\int_{\Omega}h_{n-1}\left(L-\left(\frac{\delta\breve{\Phi}}{\delta
u}\right)^{i-1}\right)|e^{i-1}_{u}|^{2}.\end{split}$ (4.20)

Note that $L>\sup\breve{\Phi}^{\prime}$ and assumption (P1) imply that

$0\leq\phi_{\mathrm{m}}\leq\left(\frac{\delta\breve{\Phi}}{\delta
u}\right)^{i-1}<L.$ (4.21)

Combining this with equation (4.20) we find that

$\frac{L+\phi_{\mathrm{m}}}{2}\int_{\Omega}h_{n-1}|e^{i}_{u}|^{2}+\frac{\varepsilon}{2}\int_{\Omega}h_{n-1}|e^{i}_{u}-e^{i-1}_{u}|^{2}+\tau\norm{\nabla
e^{i}_{w}}^{2}_{L^{2}(\Omega)}\leq\frac{L}{2}\int_{\Omega}h_{n-1}|e^{i-1}_{u}|^{2},$
(4.22)

where

$\varepsilon\coloneqq
L-\sup\breve{\Phi}^{\prime}\leq\left(L-\left(\frac{\delta\breve{\Phi}}{\delta
u}\right)^{i-1}\right).$ (4.23)

Note that the norm $\sqrt{\int_{\Omega}h_{n-1}|e^{i}_{u}|^{2}}$ is equivalent
to $\norm{e^{i}_{u}}_{L^{2}(\Omega)}$, since $0<h_{n-1}<1$ for
$\tau<1/\norm{f}_{L^{\infty}}$.

In the non-degenerate case, i.e. if $\phi_{\mathrm{m}}>0$, we obtain a
contraction as the second term in (4.22) is positive,
$\norm{\left(e^{i}_{u},e^{i}_{w}\right)}_{L}\leq\sqrt{\frac{L}{L+\phi_{\mathrm{m}}}}\norm{\left(e^{i-1}_{u},e^{i-1}_{w}\right)}_{L},$
with the norm defined in (4.11). Consequently, $u^{i}_{n}\to u_{n}$ in
$L^{2}(\Omega)$ and $w^{i}_{n}\to w_{n}$ in $H^{1}_{0}\left(\Omega\right)$ by
Banach’s fixed-point theorem.

In the degenerate case, we can sum up both sides of Equation (4.22) to find

$\begin{split}0&\leq\frac{\varepsilon}{2}\sum_{i=1}^{N}\left(\int_{\Omega}h_{n-1}|e^{i}_{u}-e^{i-1}_{u}|^{2}\right)+\tau\sum_{i=1}^{N}\norm{\nabla
e^{i}_{w}}_{L^{2}(\Omega)}^{2}\\\
&\leq\frac{L}{2}\int_{\Omega}h_{n-1}|e^{0}_{u}|^{2}-\frac{L}{2}\int_{\Omega}h_{n-1}|e^{N}_{u}|^{2}<\infty.\end{split}$
(4.24)

Hence, taking the limit $N\to\infty$ we conclude that of both sums must go to
0, yielding

$\displaystyle\norm{\nabla e^{i}_{w}}_{L^{2}(\Omega)}\to 0,$ (4.25a)
$\displaystyle\norm{e^{i}_{u}-e^{i-1}_{u}}_{L^{2}(\Omega)}\to 0.$ (4.25b)

From (4.25a) it follows that $w^{i}_{n}\to w_{n}$ in
$H^{1}_{0}\left(\Omega\right)$. Moreover, we can rewrite (4.25b) and use the
strong form of (4.2b) to find

$\norm{\frac{1}{L}\left(w^{i}_{n}-\breve{\Phi}(u^{i-1}_{n})\right)}_{L^{2}(\Omega)}=\norm{u^{i}_{n}-u^{i-1}_{n}}_{L^{2}(\Omega)}=\norm{e^{i}_{u}-e^{i-1}_{u}}_{L^{2}(\Omega)}\to
0.$ (4.26)

Hence, $\breve{\Phi}(u^{i}_{n})\to w_{n}=\breve{\Phi}(u_{n})$ in
$L^{2}(\Omega)$, as $L>0$. Finally, to prove that $u^{i}_{n}\to u_{n}$ in
$L^{1}(\Omega)$ we use Lemma 3.10 replacing the function $\psi$ with
$\breve{\Phi}$.

The result for the positive part of $u^{i}_{n}$ follows from Remark 4.1.1,
which completes the proof. ∎

###### Proof of Theorem 4.4.

The existence and uniqueness of solutions
$\\{(u^{i}_{n},w^{i}_{n})\\}_{i\in\mathbb{N}}\subset L^{2}(\Omega)\times
H^{1}_{0}\left(\Omega\right)$ follows from Lemma 4.6 and the convergence from
Lemma 4.7. ∎

We have shown that the L-scheme converges and that the error is a strict
contraction in the non-degenerate case. Unfortunately, the contraction rate is
very close to 1 if $\phi_{\mathrm{m}}>0$ is small compared to $L$, as reported
in [26, 28, 44]. A closer inspection indicates that setting
$L^{i}_{n}>\sup\breve{\Phi}^{\prime}$ everywhere in the domain is superfluous
and that this is the main reason for the slow convergence rate. To overcome
this drawback we aim to modify the L-scheme such that it is stable but
converges fast. This leads us to the M-scheme, first introduced in [28].

###  4.2 M-scheme

Note that in contrast to the L-scheme, now $L^{i}_{n}$ is a function of the
previous iterate $u^{i-1}_{n}$. We first derive two useful estimates that are
needed to prove the main convergence results.

######  Lemma 4.8 (Some useful inequalities).

Let (P1) and (4.12) hold and $M\geq M_{0}=\|\breve{\Phi}^{\prime}\|_{\rm
Lip}\Lambda$. Then, the following inequalities hold:

$\displaystyle L^{i}_{n}\geq 2M\tau^{\gamma},$ (4.27a) $\displaystyle
0\leq(M-M_{0})\tau^{\gamma}\leq
L^{i}_{n}-\left(\frac{\delta\breve{\Phi}}{\delta u}\right)^{i-1}\leq
2M\tau^{\gamma},$ (4.27b)

where $\left(\frac{\delta\breve{\Phi}}{\delta u}\right)^{i-1}$ was defined in
(4.19).

###### Proof.

Note that (4.27a) is an immediate consequence of the definition of
$L^{i}_{n}$. To prove (4.27b), we first note that

$\left(\frac{\delta\breve{\Phi}}{\delta
u}\right)^{i-1}=\breve{\Phi}^{\prime}(\zeta),$ (4.28)

for some $\zeta\in I(u^{i-1}_{n},u_{n})$ by the mean-value theorem. Moreover,
for any $\zeta\in I(u^{i-1}_{n},u_{n})$, we have

$|\breve{\Phi}^{\prime}(u^{i-1}_{n})-\breve{\Phi}^{\prime}(\zeta)|\leq\|\breve{\Phi}^{\prime}\|_{\rm
Lip}\|u^{i-1}_{n}-\zeta\|_{L^{\infty}(\Omega)}\leq\|\breve{\Phi}^{\prime}\|_{\rm
Lip}\Lambda\tau=M_{0}\tau^{\gamma},$ (4.29)

where we used (4.12) in the last inequality. If
$L^{i}_{n}=\breve{\Phi}^{\prime}(u^{i-1}_{n})+M\tau^{\gamma}$, then
$L^{i}_{n}-\breve{\Phi}^{\prime}(\zeta)\geq(M-M_{0})\tau^{\gamma}\geq 0$. On
the other hand, if $L^{i}_{n}=2M\tau^{\gamma}$, then
$\breve{\Phi}^{\prime}(u^{i-1}_{n})\leq M\tau^{\gamma}$ by the definition of
$L^{i}_{n}$. Together with (4.29) we conclude that
$\breve{\Phi}^{\prime}(\zeta)\leq\breve{\Phi}^{\prime}(u^{i-1}_{n})+M_{0}\tau^{\gamma}\leq(M+M_{0})\tau^{\gamma}$,
and thus again
$L^{i}_{n}-\breve{\Phi}^{\prime}(\zeta)\geq(M-M_{0})\tau^{\gamma}\geq 0$.

To derive the upper bound we argue analogously. If
$L^{i}_{n}=\breve{\Phi}^{\prime}(u^{i-1}_{n})+M\tau^{\gamma}$, we find

$L^{i}_{n}-\breve{\Phi}^{\prime}(\zeta)\leq\|\breve{\Phi}^{\prime}\|_{\rm
Lip}(u^{i-1}_{n}-\zeta)+M\tau^{\gamma}\leq\|\breve{\Phi}^{\prime}\|_{\rm
Lip}\Lambda\tau^{\gamma}+M\tau^{\gamma}\leq 2M\tau^{\gamma}.$ (4.30)

If $L^{i}_{n}=2M\tau^{\gamma}$, then we have
$L^{i}_{n}-\breve{\Phi}^{\prime}(\zeta)\leq L^{i}_{n}-\phi_{\mathrm{m}}\leq
2M\tau^{\gamma}$. Hence, combining all estimates we find

$0\leq(M-M_{0})\tau^{\gamma}\leq L^{i}_{n}-\breve{\Phi}^{\prime}(\zeta)\leq
2M\tau^{\gamma}.$ (4.31)

As the estimates hold for any $\zeta\in I(u^{i-1}_{n},u_{n})$, the statement
follows from (4.28).

∎

Next, we prove the first convergence result in $L^{\infty}$ for non-degenerate
problems.

######  Proposition 4.9 ($L^{\infty}$ convergence of $u^{i}_{n}$).

Assume that $\inf\Phi^{\prime}=\phi_{\mathrm{m}}>0$ and (A1) holds. Then, for
$M>M_{0}\coloneqq\|\breve{\Phi}^{\prime}\|_{\rm Lip}\Lambda$, one has

$\|w^{i}_{n}-w_{n}\|_{L^{\infty}(\Omega)}\leq
2M\tau^{\gamma}\|u^{i-1}_{n}-u_{n}\|_{L^{\infty}(\Omega)}.$ (4.32)

Moreover, if $\tau<(\phi_{\mathrm{m}}/M)^{\frac{1}{\gamma}}$, then

$\|u^{i}_{n}-u_{n}\|_{L^{\infty}(\Omega)}\leq\frac{4M\tau^{\gamma}}{\phi_{\mathrm{m}}+M\tau^{\gamma}}\|u^{i-1}_{n}-u_{n}\|_{L^{\infty}(\Omega)}.$
(4.33)

Therefore, if $\tau<(\phi_{\mathrm{m}}/(3M))^{\frac{1}{\gamma}}$, then
$u^{i}_{n}$ converges linearly in $L^{\infty}(\Omega)$ to $u_{n}$ and the
uniform boundedness of the iterates
$\|u^{i}_{n}-u_{n}\|_{L^{\infty}(\Omega)}\leq\Lambda\tau^{\gamma}$ in (4.12)
holds.

###### Proof.

We prove the statement by induction in $i\in\mathbb{N}$. For $i=1$ it is
satisfied by assumption (A1) and Remark 4.1.1. Remark 4.1.1 will also be used
in the following estimates, i.e. taking the positive part of $u_{n}$ does not
alter the inequalities. For the induction step we assume
$\|u^{i-1}_{n}-u_{n}\|_{L^{\infty}(\Omega)}\leq\Lambda\tau^{\gamma}$, which
allows us to use Lemma 4.8. We split the proof into two parts. First, we show
that (4.32) holds and subsequently, we deduce from it (4.33).

(Step 1:) Note that (4.32) is equivalent to showing that $[e^{i}_{w}-a]_{+}=0$
and $[e^{i}_{w}+a]_{-}=0$ for a specific $a>0$. We subtract (3.1a) from
Equation (4.2a) and rewrite (4.2b) as in the proof of Lemma 4.7, which yields

$\displaystyle(h_{n-1}e^{i}_{u},\phi)+\tau(\nabla e^{i}_{w},\nabla\phi)$
$\displaystyle=0,$ (4.34a)
$\displaystyle\left(L^{i}_{n}\left(e^{i}_{u}-e^{i-1}_{u}\right),\xi\right)$
$\displaystyle=\left(e^{i}_{w}-\delta\breve{\Phi}^{i-1},\xi\right).$ (4.34b)

Choosing $\phi=[e^{i}_{w}-a]_{+}\in H^{1}_{0}\left(\Omega\right)$ in (4.34a)
yields

$(h_{n-1}e^{i}_{u},[e^{i}_{w}-a]_{+})+\tau(\nabla
e^{i}_{w},\nabla[e^{i}_{w}-a]_{+})=0,$ (4.35)

and since the second term is positive, we find that

$(h_{n-1}e^{i}_{u},[e^{i}_{w}-a]_{+})\leq 0.$ (4.36)

To eliminate $e^{i}_{u}$ we observe that equation (4.34b) implies that

$e^{i}_{u}=\frac{L^{i}_{n}-\left(\frac{\delta\breve{\Phi}}{\delta
u}\right)^{i-1}}{L^{i}_{n}}e^{i-1}_{u}+\frac{e^{i}_{w}}{L^{i}_{n}},$ (4.37)

almost everywhere. Combining (4.36) and (4.37) yields

$\int_{\Omega}\frac{h_{n-1}}{L^{i}_{n}}[e^{i}_{w}-a][e^{i}_{w}-a]_{+}+\int_{\Omega}\frac{h_{n-1}}{L^{i}_{n}}\left(a+\left(L^{i}_{n}-\left(\frac{\delta\breve{\Phi}}{\delta
u}\right)^{i-1}\right)e^{i-1}_{u}\right)[e^{i}_{w}-a]_{+}\leq 0.$ (4.38)

The first term is positive and the second term can be made positive by
choosing

$a=2M\tau^{\gamma}\norm{u^{i-1}_{n}-u_{n}}_{L^{\infty}(\Omega)},$ (4.39)

as $0\leq L^{i}_{n}-\left(\frac{\delta\breve{\Phi}}{\delta u}\right)^{i-1}\leq
2M\tau^{\gamma}$ by Lemma 4.8. Hence, we find that $[e^{i}_{w}-a]_{+}=0$. The
proof for $\phi=[e^{i}_{w}+a]_{-}\in H^{1}_{0}\left(\Omega\right)$ is
analogous which proves (4.32).

(Step 2:) To show (4.33) we again note that
$L^{i}_{n}-\left(\frac{\delta\breve{\Phi}}{\delta u}\right)^{i-1}\leq
2M\tau^{\gamma}$ by Lemma 4.8 and that in the non-degenerate case, we have

$\frac{1}{L^{i}_{n}}\leq\min\left\\{\frac{1}{2M\tau^{\gamma}},\frac{1}{\phi_{\mathrm{m}}+M\tau^{\gamma}}\right\\}.$
(4.40)

Hence, using (4.32) in Equation (4.37) implies that

$\begin{split}\|e^{i}_{u}\|_{L^{\infty}(\Omega)}&\leq\min\left\\{\frac{1}{2M\tau^{\gamma}},\frac{1}{\phi_{\mathrm{m}}+M\tau^{\gamma}}\right\\}(2M\tau^{\gamma}\|e^{i-1}_{u}\|_{L^{\infty}(\Omega)}+2M\tau^{\gamma}\|e^{i-1}_{u}\|_{L^{\infty}(\Omega)})\\\
&=\min\left\\{2,\frac{4M\tau^{\gamma}}{\phi_{\mathrm{m}}+M\tau^{\gamma}}\right\\}\|e^{i-1}_{u}\|_{L^{\infty}(\Omega)}.\end{split}$
(4.41)

Consequently, if $\tau<(\frac{\phi_{\mathrm{m}}}{M})^{\frac{1}{\gamma}}$, we
get (4.33), and

Finally, note that the linear convergence and uniform $L^{\infty}$-bound of
the iterates (4.12) follows if
$\tau<(\frac{\phi_{\mathrm{m}}}{3M})^{\frac{1}{\gamma}}$. Indeed, with the
contraction rate
$\bar{\alpha}={4M\tau^{\gamma}}/({\phi_{\mathrm{m}}+M\tau^{\gamma}})<1$ one
has

$\|e^{i}_{u}\|_{L^{\infty}(\Omega)}<\bar{\alpha}\|e^{i-1}_{u}\|_{L^{\infty}(\Omega)}\leq\dots\leq\bar{\alpha}^{i}\|e^{0}_{u}\|_{L^{\infty}(\Omega)}=\bar{\alpha}^{i}\|u_{n-1}-u_{n}\|_{L^{\infty}(\Omega)}\leq\Lambda\tau^{\gamma},$

the last inequality resulting from (A1). ∎

Finally, we prove the convergence result for the M-scheme similar to Lemma 4.7
for the L-scheme. The proof is analogous, but we obtain a better contraction
rate in the non-degenerate case as the time step $\tau$ is made smaller.

######  Lemma 4.10 (Convergence of the M-scheme).

Under the assumptions of Theorem 4.5 with
$M>M_{0}:=\|\breve{\Phi}^{\prime}\|_{\rm Lip}\Lambda$, the stated convergence
results hold.

###### Proof.

As in the proof of Proposition 4.9, we subtract (3.1a) from Equation (4.2a)
and rewrite (4.2b) which yields

$\displaystyle(h_{n-1}e^{i}_{u},\phi)+\tau(\nabla e^{i}_{w},\nabla\phi)=0,$
(4.42a)
$\displaystyle(L^{i}_{n}e^{i}_{u},\xi)=(e^{i}_{w},\xi)+((L^{i}_{n}-\breve{\Phi}^{\prime}(\zeta))e^{i-1}_{u},\xi),$
(4.42b)

where $\zeta\in I(u^{i-1}_{n},u_{n})$. Choosing $\phi=e^{i}_{w}\in
H_{0}^{1}(\Omega)$ and $\xi=h_{n-1}e_{u}^{i}\in L^{2}(\Omega)$ we combine the
equations and obtain

$(h_{n-1}L_{n}^{i}e^{i}_{u},e^{i}_{u})+\tau(\nabla e_{w}^{i},\nabla
e^{i}_{w})=(h_{n-1}(L_{n}^{i}-\breve{\Phi}^{\prime}(\zeta))e^{i-1}_{u},e^{i}_{u}).$
(4.43)

We estimate the right hand side using Young’s inequality (2.4) and $0\leq
L^{i}_{n}-\Phi^{\prime}(\zeta)\leq 2M\tau^{\gamma}$, as proven in Lemma 4.8,
to find that for any $\rho>0$,

$\begin{split}(h_{n-1}(L_{n}^{i}-\breve{\Phi}^{\prime}(\zeta))e^{i-1}_{u},e^{i}_{u})&\leq
2M\tau^{\gamma}\int_{\Omega}\sqrt{h_{n-1}}e^{i-1}_{u}\sqrt{h_{n-1}}e^{i}_{u},\\\
&\leq\frac{M\tau^{\gamma}}{\rho}\int_{\Omega}h_{n-1}|e^{i-1}_{u}|^{2}+\rho
M\tau^{\gamma}\int_{\Omega}h_{n-1}|e^{i}_{u}|^{2}.\end{split}$ (4.44)

In the non-degenerate case, i.e. $\phi_{m}>0$, we estimate the left hand of
(4.43) similarly using that $L^{i}_{n}\geq\phi_{\mathrm{m}}+M\tau^{\gamma}$
and obtain

$(\phi_{\mathrm{m}}+M\tau^{\gamma})\int_{\Omega}h_{n-1}|e^{i}_{u}|^{2}+\tau\|\nabla
e^{i}_{w}\|_{L^{2}(\Omega)}^{2}\leq(h_{n-1}L_{n}^{i}e^{i}_{u},e^{i}_{u})+\tau\|\nabla
e^{i}_{w}\|_{L^{2}(\Omega)}^{2}.$ (4.45)

Combining (4.44) and (4.45) it follows that

$(\phi_{\mathrm{m}}+(1-\rho)M\tau^{\gamma})\int_{\Omega}h_{n-1}(e^{i}_{u})^{2}+\tau\|\nabla
e^{i}_{w}\|_{L^{2}(\Omega)}^{2}\leq\frac{M\tau^{\gamma}}{\rho}\int_{\Omega}h_{n-1}(e^{i-1}_{u})^{2}.$
(4.46)

which implies that

$\norm{(e^{i}_{u},e^{i}_{w})}_{M,\rho}\leq\sqrt{\frac{M\tau^{\gamma}}{\rho(\phi_{\mathrm{m}}+(1-\rho)M\tau^{\gamma})}}\norm{(e^{i-1}_{u},e^{i-1}_{w})}_{M,\rho},$
(4.47)

where

$\norm{(e^{i}_{u},e^{i}_{w})}_{M,\rho}^{2}\coloneqq\int_{\Omega}h_{n-1}|e^{i}_{u}|^{2}+\frac{\tau}{\phi_{\mathrm{m}}+(1-\rho)M\tau^{\gamma}}\|\nabla
e^{i}_{w}\|_{L^{2}(\Omega)}^{2}.$ (4.48)

Equation (4.48) defines a norm if
$0<\rho<1+\frac{\phi_{\mathrm{m}}}{M\tau^{\gamma}}$, and choosing
$0<\rho=\rho^{\ast}=\frac{1}{2}(1+\frac{\phi_{\mathrm{m}}}{M\tau^{\gamma}})<1+\frac{\phi_{\mathrm{m}}}{M\tau^{\gamma}}$,
minimizes the contraction rate. Hence,

$\norm{(e^{i}_{u},e^{i}_{w})}_{M,\rho^{\ast}}\eqqcolon\norm{(e^{i}_{u},e^{i}_{w})}_{M}\leq\frac{2M\tau^{\gamma}}{\phi_{\mathrm{m}}+M\tau^{\gamma}}\norm{(e^{i-1}_{u},e^{i-1}_{w})}_{M},$
(4.49)

which is a contraction if $\tau<({\phi_{\mathrm{m}}}/{M})^{\frac{1}{\gamma}}$,
and the contraction rate scales with $\tau^{\gamma}$. As in the proof of Lemma
4.7 we conclude that $u^{i}_{n}\to u_{n}$ in $L^{2}(\Omega)$ and $w^{i}_{n}\to
w_{n}$ in $H^{1}(\Omega)$ by Banach’s fixed-point theorem.

The degenerate case is dealt with in the same manner as in the proof of Lemma
4.7. Analogous to (4.22), we find using $L^{i}_{n}\geq 2M\tau^{\gamma}$ and
$\varepsilon:=\inf\left(L^{i}_{n}-\left(\frac{\delta\breve{\Phi}}{\delta
u}\right)^{i-1}\right)\overset{\eqref{ineq: 2}}{\geq}(M-M_{0})\tau>0$ that

$M\tau^{\gamma}\int_{\Omega}h_{n-1}|e^{i}_{u}|^{2}+\frac{\varepsilon}{2}\int_{\Omega}h_{n-1}|e^{i}_{u}-e^{i-1}_{u}|^{2}+\tau\norm{\nabla
e^{i}_{w}}^{2}_{L^{2}(\Omega)}\leq
M\tau^{\gamma}\int_{\Omega}h_{n-1}|e^{i-1}_{u}|^{2}.$ (4.50)

Summing both sides of (4.50) yields

$\begin{split}0&\leq\frac{\varepsilon}{2}\sum_{i=1}^{N}\left(\int_{\Omega}h_{n-1}|e^{i}_{u}-e^{i-1}_{u}|^{2}\right)+\tau\sum_{i=1}^{N}\norm{\nabla
e^{i}_{w}}_{L^{2}(\Omega)}^{2}\\\ &\leq
M\tau^{\gamma}\left(\int_{\Omega}h_{n-1}|e^{0}_{u}|^{2}-\int_{\Omega}h_{n-1}|e^{N}_{u}|^{2}\right)<\infty.\end{split}$
(4.51)

This implies that $w^{i}_{n}\to w_{n}$ in $H^{1}_{0}\left(\Omega\right)$,
$\breve{\Phi}(u^{i}_{n})\to w_{n}=\breve{\Phi}(u_{n})$ in $L^{2}(\Omega)$ and
$u^{i}_{n}\to u_{n}$ in $L^{1}(\Omega)$. The result for the positive part of
$u^{i}_{n}$ follows from Remark 4.1.1 which completes the proof. ∎

######  Remark 4.10.1.

Following the proof for the L-scheme in Lemma 4.7 we would obtain the
contraction rate

$\norm{(e^{i}_{u},e^{i}_{w})}_{M}\leq\sqrt{\frac{2M\tau^{\gamma}}{\phi_{\mathrm{m}}+M\tau^{\gamma}}}\norm{(e^{i-1}_{u},e^{i-1}_{w})}_{M},$
(4.52)

which is a larger than the rate in Lemma 4.10 if
$\tau<(\phi_{\mathrm{m}}/M)^{\frac{1}{\gamma}}$. However, this is the range of
$\tau$ where the M-scheme provides a contraction, and therefore Lemma 4.10,
specifically (4.49), provides a sharper result. Furthermore, the contraction
rate stated in Lemma 4.10 is half the contraction rate obtained for the
$L^{\infty}$-norm in Proposition 4.9.

###### Proof of Theorem 4.5.

The existence and uniqueness of solutions
$\\{(u^{i}_{n},w^{i}_{n})\\}_{i\in\mathbb{N}}\subset L^{2}(\Omega)\times
H^{1}_{0}\left(\Omega\right)$ for (4.2) follows from Lemma 4.6 and the
convergence results from Lemma 4.10. ∎

##  5 Numerical Results

We use the finite element method (FEM) to compute the solutions as it directly
links to the weak form of the Equations (3.1) and (4.2). The FEniCSx package
in Python is used to solve the finite-dimensional problems [2, 39, 40], and
all the code is made available on GitHub111Link to the GitHub repository:
https://github.com/Rsmeets99/M_scheme_biofilm_PDE. Let $\mathcal{T}$ denote
the triangulation of $\Omega$ and let $\mathcal{P}_{p}(\mathcal{T})$ be the
space of element-wise polynomials of degree up to $p\in\mathds{N}$. We define
the FEM solutions to be $\tilde{u}_{n,h}^{i}\in
U_{h}\coloneqq\mathcal{P}_{0}(\mathcal{T})$ and $w_{n,h},w^{i}_{n,h}\in
V_{h}\coloneqq\mathcal{P}_{1}(\mathcal{T})\cap H^{1}_{0}(\Omega)$. For
$\mu=1$, we take the spatial approximation $v_{n,h}\in V_{h}$ of $v_{n}$ since
in this case $v_{n}$ is differentiable, and $v_{n,h}\in U_{h}$ otherwise. This
leads us to the following problem:

######  Problem 5.1 (Finite element system).

Given $u_{n-1,h},u^{i-1}_{n,h},v_{n-1,h}$, find
$\left(\tilde{u}^{i}_{n,h},w^{i}_{n,h}\right)\in Z$, such that

$\begin{split}&\left(h_{n-1}\tilde{u}^{i}_{n,h}\right),\phi_{h}+\tau\left(\nabla
w^{i}_{n,h},\nabla\phi_{h}\right)=\left(u_{n-1,h},\phi_{h}\right)\\\
&\left(L^{i}_{n}\tilde{u}^{i}_{n,h}-w^{i}_{n,h},\xi_{h}\right)=\left(L^{i}_{n}u^{i-1}_{n,h}-\Phi\left(u^{i-1}_{n,h}\right),\xi_{h}\right)\end{split}$
(5.1)

for all $\left(\xi_{h},\phi_{h}\right)\in Z$, where $Z$ denotes the (mixed)
finite element space $Z=U_{h}\times V_{h}$. Afterwards, set
$u^{i}_{n,h}=\left[\tilde{u}^{i}_{n,h}\right]_{+}$.

We iteratively solve $u^{i}_{n,h}$ and $w^{i}_{n,h}$ until the following
stopping criteria is met:

$\norm{\left(e^{i}_{u,h},e^{i}_{w,h}\right)}=\int_{\Omega}L^{i}_{n}|e^{i}_{u,h}|^{2}\text{d}x+\tau\norm{\nabla
e^{i}_{w,h}}^{2}_{L^{2}(\Omega)}<\text{tol},$ (5.2)

where $e^{i}_{u,h}\coloneqq u^{i}_{n,h}-u^{i-1}_{n,h}$, $e^{i}_{w,h}\coloneqq
w^{i}_{n,h}-w^{i-1}_{n,h}$ and $\text{tol}\in\mathds{R}^{+}_{\ast}$ is some
tolerance. Once the tolerance is reached, we set $u_{n,h}=u^{i}_{n,h}$ and
calculate $v_{n,h}\in V_{h}$ by solving

$(v_{n,h},\eta_{h})+\tau\mu(D(u_{n,h})\nabla
v_{n,h},\nabla\eta_{h})=\tau(g(u_{n,h},v_{n-1,h}),\eta_{h})+(v_{n-1,h},\eta_{h}),\qquad\forall\eta_{h}\in
V_{h},$ (5.3)

in the PDE-PDE case ($\mu=1$), or by solving $v_{n,h}\in U_{h}$

$(v_{n,h},\eta_{h})=\tau(g(u_{n,h},v_{n-1,h}),\eta_{h})+(v_{n-1,h},\eta_{h}),\qquad\forall\eta_{h}\in
U_{h},$ (5.4)

in the PDE-ODE case ($\mu=0$).

Depending on the specific function $g$ we could update $v_{n,h}$ in the PDE-
ODE case explicitly through

$v_{n,h}=v_{n-1,h}+\tau g(u_{n,h},v_{n-1,h}).$ (5.5)

However, in general we cannot guarantee that $v_{n,h}\in U_{h}$, while
Equation (5.4) provides a projection onto the correct space. The full
algorithm is summarised in Algorithm 1. Instead of solving the full system
(4.2), it is possible to eliminate $u$ from the equations and solve only for
$w$. While faster to solve due to the reduced dimension of the resulting
linear system, it does require correct projection operators and this leads to
a modification of $M$ making it dependent on the mesh size $h$. More details
are given in [42, Section 5.1.2].

$t=t_{\text{start}}$ ;

for  _$t <T$_ do

error = 1 ;

while  _$\text{error} >\text{tol}$_ do

Solve system (5.1) from 5.1 for $\tilde{u}^{i}_{n,h},w^{i}_{n,h}$ ;

Set $u^{i}_{n,h}=\left[\tilde{u}^{i}_{n,h}\right]_{+}$ ;

Compute new error;

end while

set $u_{n,h}=u^{i}_{n,h}$ ;

set $w_{n,h}=w^{i}_{n,h}$ ;

if  _$\mu=1$_ then

Compute $v_{n,h}$ in the PDE case using $u_{n,h}$ through solving equation
(5.3)

else

Compute $v_{n,h}$ in the ODE case using $u_{n,h}$ through solving equation
(5.4)

end if

$t=t+\tau$ ;

end for

Algorithm 1 M-scheme algorithm

We test our scheme for 3 different problems: a porous medium equation, the
biofilm PDE-PDE model and the biofilm PDE-ODE model. The goal is to get
numerical convergence results, as well as to compare the performance of the
M-scheme to the Newton.

######  Remark 5.1.1 (Newton scheme).

The ‘true’ Newton scheme with $M=0$ may not converge in the degenerate case
without regularization. To overcome this problem we use the M-scheme with a
very small $M$ (e.g. $M=10^{-7}\ll\text{tol}$) which is still large enough so
that the scheme converges in most cases. This is a form of a regularized
Newton scheme (4.9).

######  Remark 5.1.2 (L-scheme).

For the porous medium equation, the L-scheme is at least an order of magnitude
slower than the M-scheme, while for the biofilm models the L required for
convergence is so large, that it becomes multiple orders of magnitude slower.
Therefore, we will not show the results for the L-scheme in our comparison.

######  Remark 5.1.3 (Adaptive M-scheme).

When using the M scheme, in practice, it is beneficial to choose $M$
adaptively in each step using a posteriori estimators. In the biofilm case,
with an adaptive scheme, the required $M$ increases if $u_{n}$ approaches 1
which ensures convergence, while $M$ is small and therefore the scheme is fast
when $u_{n}$ is bounded away from 1. Similar work for the L-scheme has been
done in [44]. This is however beyond the scope of our current work.

###  5.1 Porous medium equation

As a first test case, we consider the 1D porous medium equation with a
reaction term

$\partial_{t}u=\Delta\left(u^{m}\right)+\beta u,$ (5.6)

where $m>1$ and $\beta\in\mathds{R}$ in a bounded interval
$\Omega\subset\mathds{R}$ with homogeneous Dirichlet boundary conditions. In
our notation, this corresponds to $b=\infty$, $\Phi(u)=u^{m}$ and
$f(v)=\beta$. Note here that $h_{n-1}>0$ if $\beta<1/\tau$. It serves as a
good benchmark problem as (5.6) has the exact solution $u(x,t)=e^{\beta
t}z(x,s)$, where $s=\frac{1}{\beta(m-1)}e^{\beta(m-1)t}$ and $z$ is the
Barenblatt solution [46] given by

$z(x,t)=t^{-\frac{d}{d(m-1)+2}}\left[C-\frac{m-1}{2m(d(m-1)+2)}\left|xt^{-\frac{1}{d(m-1)+2}}\right|^{2}\right]_{+}^{\frac{1}{m-1}}.$
(5.7)

The exact solution $u$ is Hölder continuous in time with exponent
$\gamma=1/(m-1)$.

We first verify the consistency of the time discretisation stated in Theorem
3.2 by computing the error in the left-hand side of (3.5) for different values
of $\tau$. As initial condition, we take the exact solution $u$ evaluated at
$t=0.5$. The results are shown in Figure 1 exhibiting an order of convergence
between 0.5 and 1. Note that the results are independent of the choice of $M$
and $\gamma$ as long as the linearisation scheme converges.

Figure 1: Error in (3.5) against time step size $\tau$ for $h=10^{-4}$, $m=4$,
$\gamma=1/3$, time $0.5\leq t\leq 1$, $\text{tol}=10^{-7}$, $d=1$, $\beta=1$.

A convergence study of the iterative schemes for different time-steps
$\tau\in\\{10^{-1},10^{-1.5},$ $10^{-2},10^{-2.5}\\}$ was also performed with
the mesh size $h$ ranging between $0.1$ and $0.005$. Then the average number
of iterations needed to get to a final time $T=1.1$ was determined for
different values of $M\in\\{10^{-1},10^{-2},10^{-3},10^{-7}\\}$. The results
are displayed in Figure 2. We first note that for each $\tau$ there is an
optimal value $M$ and the M-scheme out-performs the Newton-scheme in this
case, which is most apparent for smaller mesh sizes and larger time-steps. It
is expected that the convergence of the Newton-scheme is conditioned by
restrictions on the time step size depending on the mesh size [36]. For
instance, for $\tau=10^{-1}$ the M-scheme performs significantly better than
the Newton-scheme for small mesh sizes, while the schemes are equivalent for
$\tau=10^{-2.5}$. Secondly, we note that the optimal value $M$ stays optimal
for all mesh sizes. Hence, we can find the optimal $M$ for a coarse mesh, and
use it then for computations on finer meshes [45]. Lastly, the number of
iterations required decreases with decreasing $\tau$. The reason is two-fold:
first, we expect the convergence rate to increase as $\tau$ gets smaller by
Theorem 4.5, and secondly, the difference between the solutions of two
consecutive time steps $u_{n}$ and $u_{n-1}$ decreases when $\tau$ does, and
therefore the iterations start with a better initial guess
$u^{0}_{n}=u_{n-1}$.

(a) $\tau=10^{-1}$

(b) $\tau=10^{-1.5}$

(c) $\tau=10^{-2}$

(d) $\tau=10^{-2.5}$

Figure 2: Average iterations required for solving (5.6) in 1D for varying mesh
size $h$ and time steps $\tau$, with $m=4$, $\gamma=1/3$, for time $0.5\leq
t\leq 1.1$, $\text{tol}=10^{-5}$.

Having found an optimal $M$ for $\gamma=1/(m-1)$ (which is $M=10^{-3}$), we
next test scaling of the contraction rates predicted by Theorem 4.5. Observe
that not all assumptions are satisfied as the problem is degenerate.
Nevertheless, we find a scaling of the contraction rate with some power of
$\tau$, as shown in Figure 3. The contraction rate is calculated as the
geometric mean of
$\norm{(e^{i}_{u},e^{i}_{w})}/\norm{(e^{i-1}_{u},e^{i-1}_{w})}$ over the first
3 iterations. Note that the convergence rate $\alpha$ appears to scale
linearly with $\tau^{0.42}$ instead of $\tau^{\gamma}=\tau^{0.33}$. This
‘super-convergence’ can be explained due to the fact that the challenging part
of the numerical solution is the free boundary, while the solution is much
more regular in the rest of the domain. For the test case in 1D, the free
boundary consists of just two points. Hence, it does not play a deciding role
in the convergence behaviour.

Figure 3: Convergence rate $\alpha$ against time step size $\tau$ for
$h=10^{-4}$, $m=4$, $\gamma=1/3$, $t=0.5$, $d=1$.

###  5.2 Biofilm equations

In this section, we investigate the robustness of the M-schem for the more
challenging biofilm models (1.1) which are coupled systems involving a
singular-degenerate diffusion equation. We consider the case $\mu=1$
corresponding to a PDE-PDE coupling [15], and $\mu=0$ corresponding to a PDE-
ODE coupling [13]. The corresponding functions in (1.1) are as follows:

$\displaystyle\Phi^{\prime}(u)=d_{1}\frac{u^{\alpha}}{\left(1-u\right)^{\beta}},$
$\displaystyle f(v)=k_{3}\frac{v}{v+k_{2}}-k_{4},$ (5.8a) $\displaystyle
D(u)=d_{2},$ $\displaystyle g(u,v)=-k_{1}\frac{uv}{v+k_{2}},$ (5.8b)

for some given constants $k_{1},k_{2},k_{3},k_{4},d_{1},d_{2}>0$ and
$\alpha,\beta>1$. For our comparison of the M-scheme and (regularized) Newton-
scheme, we will use the same parameters as in [13], which are $k_{1}=0.4$,
$k_{2}=0.01$, $k_{3}=1$, $k_{4}=0.42$, $d_{1}=10^{-6}$, $\alpha=4$, $\beta=4$.
For fixed $\alpha=\beta=4$, we can calculate $\Phi(u)$ explicitly,

$\begin{split}\Phi(u)=10^{-6}\int_{0}^{u}\frac{s^{4}}{\left(1-s\right)^{4}}ds=10^{-6}\left(\frac{18u^{2}-30u+13}{3\left(1-u\right)^{3}}+u+4\ln\left(1-u\right)-\frac{13}{3}\right).\end{split}$

Note that we cannot simply use $\Phi$ but have to use its regularized form
$\breve{\Phi}$ as given in Definition 4.1. To define $\breve{\Phi}$ we use
Theorem 3.1 to calculate the upper bound $\breve{u}$. As the initial condition
we take

$u_{0}(x)=\frac{h}{r}\left(\sqrt{\max(0,r^{2}-(x-x_{1})^{2})}+\sqrt{\max(0,r^{2}-(x-x_{2})^{2})}\right),$
(5.9)

with a maximum height $h=0.9$, radius $r=0.2$, $x_{1}=-0.3$, $x_{2}=0.3$. For
the domain $\Omega=(-1,1)$, this yields $\breve{u}=0.992$. For $\gamma$ we
take $\gamma=1/\alpha=1/4$ since the regularity of solutions is expected to be
similar as for the porous medium equation with the diffusion coefficient
$\Phi^{\prime}\approx u^{\alpha}$ when $u$ is small (close to the free
boundary). We assume homogeneous Neumann boundary conditions for $u$ and, if
$\mu=1$, mixed boundary conditions for $v$. Namely, at the boundary $x=-1$ we
specify the Dirichlet condition $v=1$ and at $x=1$ homogeneous Neumann
boundary conditions. While homogeneous Neumann conditions for $u$ are not
covered by our theory, we still expect the results to hold as in the
simulations, the biofilm region marked by the support of $u$, never reaches
the boundary. The well-posedness results for time-continuous models in [23,
29] also apply to inhomogeneous Dirichlet and mixed boundary conditions, and
to homogeneous Neumann conditions under certain time restrictions. We impose
these boundary conditions as they were chosen for the numerical results in
[13, 15], from where we also took our parameter values.

The results of the M-scheme and (regularized) Newton-scheme are given in
Figure 4 for the PDE-ODE case. For the chosen parameters, the behaviour of the
iterative schemes for the PDE-PDE case is almost identical. As in Figure 2 we
see that smaller time-steps $\tau$ require fewer iterations for the reasons
mentioned before. We only note that the amount of required iterations for the
biofilm system is considerably higher. When solutions of the biofilm model
approach 1, the diffusion coefficient $\breve{\Phi}^{\prime}$ becomes very
large, and therefore $L^{i}_{n}$ as well. This slows down the convergence of
the iterative scheme.

A noticeable difference between the two figures is that if the mesh size is
too small for the corresponding time-step, the Newton scheme becomes unstable
in the biofilm case and often does not converge. As the time-step size
decreases, the Newton scheme starts converging for smaller mesh sizes.
However, we see a gap in the performance between the M-scheme and the Newton-
scheme for these larger mesh sizes. The reason is that it is impossible to
choose an optimal $M$ for the entire time range. A smaller $M$ would have a
similar performance as the Newton-scheme for large mesh sizes, but as $u_{n}$
approaches values closer to 1, the diffusion coefficient blows up and
convergence is no longer guaranteed for smaller mesh sizes. The choice of $M$
is therefore dictated by how close $u_{n}$ gets to 1. A larger $M$ improves
stability at the cost of convergence speed.

(a) $\tau=10^{-1}$

(b) $\tau=10^{-1.5}$

(c) $\tau=10^{-2}$

(d) $\tau=10^{-2.5}$

Figure 4: Average iterations required for solving (5.8) in 1D for varying mesh
size $h$ and time steps $\tau$, with $m=4$, $\gamma=1/4$, for time $0\leq
t\leq 1.2$, $\mu=0$ and $\text{tol}=10^{-5}$.

Similarly to Figure 3, we calculate the contraction rate as the geometric mean
over the first three iterations for different values of $\tau$ in Figure 5. We
find the contraction rate to be approximately $\tau^{0.25}$, which aligns with
our predicted value $\gamma$.

Figure 5: Convergence rate $\alpha$ against time step size $\tau$ for
$h=10^{-4}$, $m=4$, $\gamma=1/4$.

A test example for both the PDE-PDE and PDE-ODE cases in two dimensions is
given in Figure 6 and Figure 7 respectively. For the initial conditions, we
have chosen two hemispheres as in the 1D case and used similar parameter
values as disclosed in the caption of both figures. Computationally, these two
problems are challenging. We find values of $u_{n}$ very close to $\breve{u}$
($\breve{u}=0.989$ and $\breve{u}=0.988$ for the PDE-PDE and PDE-ODE cases
respectively), which leads to a blow-up of the diffusion coefficient. On top
of that, the two blobs possess sharp interfaces that merge at some point
creating additional singularities, see Figure 7. For these reasons, the mesh
size is kept relatively small to accurately resolve this merging. Despite
these challenges, the numerical methods perform robustly, and we recover the
expected qualitative behaviour of the solutions of both models. In the PDE-PDE
simulation, we see that since the nutrients $v_{n}$ diffuse and are constantly
added through the Dirichlet boundary conditions on the top boundary, the
biofilm expands towards the top, while slowly dying off at the bottom. For the
PDE-ODE simulation, the biofilm expands in the radial direction as it consumes
the nutrients while dying off in places where nutrients have been depleted.
This leads to crater-like structures and inverse colony formation as seen in
experiments, e.g. see [13]. Each 2D simulation required a long run-time, and
due to the limitation of computational resources, a thorough comparison of
iterative schemes could not be conducted in 2D.

Figure 6: Simulation of the PDE-PDE model using the M-scheme ($M=10^{-2})$,
with $k_{1}=5,k_{2}=0.01,k_{3}=1,k_{4}=0.42,d_{1}=5\cdot
10^{-6},d_{2}=0.2,\tau=0.01,h=0.02$ and $\gamma=0.5$. The first and third row
picture $u_{n}$ while the second and fourth row $v_{n}$. For $v_{n}$ we have
homogeneous Neumann boundary conditions at the sides and bottom, and the
Dirichlet boundary condition $v=1$ at the top. Figure 7: Simulation of the
PDE-ODE model using the M-scheme ($M=10^{-2})$, with
$k_{1}=0.8,k_{2}=0.01,k_{3}=1,k_{4}=0.42,d_{1}=8\cdot
10^{-6},\tau=0.01,h=0.02$ and $\gamma=0.5$. The first and third row picture
$u_{n}$ while the second and fourth row $v_{n}$.

##  6 Conclusion

We introduced a semi-implicit time-discretisation scheme for solving a class
of degenerate quasi-linear parabolic problems of porous medium type with
diffusion coefficients that can also be singular. Such systems model biofilm
growth and other nonlinear diffusion processes with sharp interfaces that
propagate at a finite speed. The well-posedness of the time-discrete solutions
was shown, as well as explicit upper bounds were proved for both the biofilm
model ($b=1$) and porous medium equations ($b=\infty$). We then introduced the
L/M-scheme as an iterative linearisation method for solving the quasi-linear
elliptic PDEs that resulted from the time discretisation. It was shown that
these schemes converge irrespective of the spatial discretisation. In the non-
degenerate case, these schemes will even show a contraction with a contraction
rate that scales with some power of $\tau$ for the M-scheme provided $\tau$ is
small. Finally, the schemes were implemented numerically using a finite
element method and it was shown that for larger time steps $\tau$ and finer
mesh sizes $h$, the M-scheme outperforms the Newton scheme.

The schemes can be generalised to systems that allow for additional substrates
and admit terms for an advective flow field in these equations, see [29]. In a
future work, we are considering including a nonlinearity in the time
derivative as well which makes the problem doubly degenerate. Such problems
are commonly found in multiphase flow through porous medium. Furthermore, one
can consider multi-species biofilm models with or without cross-diffusion that
comprise multiple degenerate equations that are strongly coupled through the
diffusion operator, e.g. see [19, 43].

## Acknowledgements

K. Mitra acknowledges the support of Research Foundation - Flanders (FWO), the
Junior Postdoctoral Fellowship grant 1209322N. K. Mitra and S. Sonner thank
the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO) for
financial support (grant OCENW.KLEIN.358). The work of I.S. Pop was supported
by FWO through the Odysseus programme (Project G0G1316N) and the German
Research Foundation (DFG) through the SFB 1313, Project Number 327154368.

## References

  * [1] A. Alhammali, M. Peszynska, and C. Shin, Numerical analysis of a mixed finite element approximation of a coupled system modeling biofilm growth in porous media with simulations, International Journal of Numerical Analysis and Modeling, 21 (2024), pp. 20–64.
  * [2] M. S. Alnaes, A. Logg, K. B. Ølgaard, M. E. Rognes, and G. N. Wells, Unified form language: A domain-specific language for weak formulations of partial differential equations, ACM Transactions on Mathematical Software, 40 (2014).
  * [3] H. W. Alt and S. Luckhaus, Quasilinear elliptic-parabolic differential equations, Math. z, 183 (1983), pp. 311–341.
  * [4] E. Balsa-Canto, A. López-Núñez, and C. Vázquez, Numerical methods for a nonlinear reaction–diffusion system modelling a batch culture of biofilm, Applied Mathematical Modelling, 41 (2017), pp. 164–179.
  * [5] J. W. Barrett and K. Deckelnick, Existence, uniqueness and approximation of a doubly-degenerate nonlinear parabolic system modelling bacterial evolution, Mathematical Models and Methods in Applied Sciences, 17 (2007), pp. 1095–1127.
  * [6] K. Brenner and C. Cances, Improving newton’s method performance by parametrization: the case of the richards equation, SIAM Journal on Numerical Analysis, 55 (2017), pp. 1760–1785.
  * [7] C. Cancès, J. Droniou, C. Guichard, G. Manzini, M. B. Olivares, and I. S. Pop, Error estimates for the gradient discretisation method on degenerate parabolic equations of porous medium type, in Polyhedral methods in geosciences, Springer, 2021, pp. 37–72.
  * [8] X. Cao and K. Mitra, Error estimates for a mixed finite element discretization of a two-phase porous media flow model with dynamic capillarity, Journal of Computational and Applied Mathematics, 353 (2019), pp. 164–178.
  * [9] M. A. Celia, E. T. Bouloutas, and R. L. Zarba, A general mass-conservative numerical solution for the unsaturated flow equation, Water Resources Research, 26 (1990), pp. 1483–1496.
  * [10] D. L. Cohn, Measure Theory: Second Edition, Birkhäuser Advanced Texts Basler Lehrbücher, Springer New York, 2013.
  * [11] E. S. Daus, P. Milišić, and N. Zamponi, Analysis of a degenerate and singular volume-filling cross-diffusion system modeling biofilm growth, SIAM Journal on Mathematical Analysis, 51 (2019), pp. 3569–3605.
  * [12] A. Duvnjak and H. J. Eberl, Time-discretisation of a degenerate reaction-diffusion equation arising in biofilm modeling, Electronic Transactions on Numerical Analysis, 23 (2006), pp. 15–38.
  * [13] H. J. Eberl, E. M. Jalbert, A. Dumitrache, and G. M. Wolfaardt, A spatially explicit model of inverse colony formation of cellulolytic biofilms, Biochemical Engineering Journal, 122 (2017), pp. 141–151.
  * [14] H. J. Eberl and D. Laurent, A finite difference scheme for a degenerated diffusion equation arising in microbial ecology, Electronic Journal of Differential Equations, CS15 (2007).
  * [15] H. J. Eberl, D. F. Parker, and M. C. M. van Loosdrecht, A new deterministic spatio-temporal continuum model for biofilm development, Journal of Theoretical Medicine, 3 (2001), pp. 161–175.
  * [16] M. A. Efendiev, S. Zelik, and H. J. Eberl, Existence and longtime behavior of a biofilm model, Communications on Pure & Applied Analysis, 8 (2009), pp. 509–531.
  * [17] B. O. Emerenini, S. Sonner, and H. J. Eberl, Mathematical analysis of a quorum sensing induced biofilm dispersal model and numerical simulation of hollowing effects, Mathematical Biosciences & Engineering, 14 (2016), pp. 625–653.
  * [18] L. C. Evans, Partial Differential Equations: Second Edition, American Mathematical Society, Providence, R.I., 2010.
  * [19] M. Ghasemi, S. Sonner, and H. J. Eberl, Time adaptive numerical solution of a highly non-linear degenerate cross-diffusion system arising in multi-species biofilm modelling, European Journal of Applied Mathematics, 29 (2018), p. 1035–1061.
  * [20] C. Helmer, A. Jüngel, and A. Zurek, Analysis of a finite-volume scheme for a single-species biofilm model, Applied Numerical Mathematics, 185 (2023), pp. 386–405.
  * [21] R. Helmig, Multiphase Flow and Transport Processes in the Subsurface, Environmental Science and Engineering, Springer Berlin, Heidelberg, 1997.
  * [22] V. Hissink Muller, Interior Hölder continuity for singular-degenerate porous medium type equations with an application to a biofilm model, Journal of Evolution Equations, 22 (2022), pp. 1–92.
  * [23] V. Hissink Muller and S. Sonner, Well-posedness of singular-degenerate porous medium type equations and application to biofilm models, Journal of Mathematical Analysis and Applications, 509 (2022), p. 125894.
  * [24] W. J. Jäger and J. Kacur, Solution of porous medium type systems by linear approximation schemes., Numerische Mathematik, 60 (1991/92), pp. 407–428.
  * [25] W. J. Jäger and J. Kacur, Solution of doubly nonlinear and degenerate parabolic problems by relaxation schemes, ESAIM: M2AN, 29 (1995), pp. 605–627.
  * [26] F. List and F. A. Radu, A study on iterative methods for solving richards‘ equation, Computational Geosciences, 20 (2016).
  * [27] K. Mitra, Existence and properties of solutions of the extended play-type hysteresis model, Journal of Differential Equations, 288 (2021), pp. 118–140.
  * [28] K. Mitra and I. S. Pop, A modified l-scheme to solve nonlinear diffusion problems, Computers & Mathematics with Applications, 77 (2019), pp. 1722–1738.  7th International Conference on Advanced Computational Methods in Engineering (ACOMEN 2017).
  * [29] K. Mitra and S. Sonner, Well-posedness and qualitative properties of quasilinear degenerate evolution systems, arXiv preprint, 2304.00175 (2023).
  * [30] K. Mitra and C. J. van Duijn, Wetting fronts in unsaturated porous media: The combined case of hysteresis and dynamic capillary pressure, Nonlinear Analysis: Real World Applications, 50 (2019), pp. 316–341.
  * [31] F. Otto, L1-contraction and uniqueness for quasilinear elliptic–parabolic equations, Journal of Differential Equations, 131 (1996), pp. 20–38.
  * [32] E. J. Park, Mixed finite elements for nonlinear second-order elliptic problems, SIAM Journal on Numerical Analysis, 32 (1995), pp. 865–885.
  * [33] I. S. Pop, F. A. Radu, and P. Knabner, Mixed finite elements for the richards’ equation: linearization procedure, Journal of Computational and Applied Mathematics, 168 (2004), pp. 365–373.  Selected Papers from the Second International Conference on Advanced Computational Methods in Engineering (ACOMEN 2002).
  * [34] I. S. Pop and B. Schweizer, Regularization schemes for degenerate richards equations and outflow conditions, Mathematical Models and Methods in Applied Sciences, 21 (2011), pp. 1685–1712.
  * [35] F. A. Radu, K. Kumar, J. M. Nordbotten, and I. S. Pop, A robust, mass conservative scheme for two-phase flow in porous media including hölder continuous nonlinearities, IMA Journal of Numerical Analysis, 38 (2017), pp. 884–920.
  * [36] F. A. Radu, I. S. Pop, and P. Knabner, Newton—type methods for the mixed finite element discretization of some degenerate parabolic equations, in Numerical Mathematics and Advanced Applications, A. B. de Castro, D. Gómez, P. Quintela, and P. Salgado, eds., Berlin, Heidelberg, 2006, Springer Berlin Heidelberg, pp. 1192–1200.
  * [37] K. Rahman, R. Sudarsan, and H. J. Eberl, A mixed-culture biofilm model with cross-diffusion, Bulletin of Mathematical Biology, 77 (2015), pp. 2086–2124.
  * [38] C. Reisch, A. Navas-Montilla, and I. Özgen-Xian, Analytical and numerical insights into wildfire dynamics: Exploring the advection-diffusion-reaction model, arXiv preprint arXiv:2307.16174, (2023).
  * [39] M. W. Scroggs, I. A. Baratta, C. N. Richardson, and G. N. Wells, Basix: a runtime finite element basis evaluation library, Journal of Open Source Software, 7 (2022), p. 3982.
  * [40] M. W. Scroggs, J. S. Dokken, C. N. Richardson, and G. N. Wells, Construction of arbitrary order finite element degree-of-freedom maps on polygonal and polyhedral cell meshes, ACM Transactions on Mathematical Software, 48 (2022).
  * [41] D. Seus, K. Mitra, I. S. Pop, F. A. Radu, and C. Rohde, A linear domain decomposition method for partially saturated flow in porous media, Computer Methods in Applied Mechanics and Engineering, 333 (2018), pp. 331–355.
  * [42] R. K. H. Smeets, Numerical schemes for singular porous medium type equations, Master’s thesis, Radboud University Nijmegen & Hasselt University, 2023.
  * [43] S. Sonner, M. A. Efendiev, and H. J. Eberl, On the well-posedness of mathematical models for multicomponent biofilms, Math. Methods Appl. Sci., 38 (2015), pp. 3753–3775.
  * [44] J. S. Stokke, K. Mitra, E. Storvik, J. Both, and F. A. Radu, An adaptive solution strategy for Richards’ equation, Computers & Mathematics with Applications, 152 (2023), pp. 155–167.
  * [45] E. Storvik, J. W. Both, K. Kumar, J. M. Nordbotten, and F. A. Radu, On the optimization of the fixed-stress splitting for biot’s equations, International Journal for Numerical Methods in Engineering, 120 (2019), pp. 179–194.
  * [46] J. L. Vazquez, The Porous Medium Equation: Mathematical Theory, Oxford Mathematical Monographs, Clarendon Press, 2007.