# Well-posedness and qualitative properties of quasilinear degenerate
evolution systems

K. Mitra111email<EMAIL_ADDRESS>Faculty of Science, Hasselt
University, Hasselt, Belgium  S. Sonner  Faculty of Science, Radboud
University, Nijmegen, The Netherlands

###### Abstract

We analyze nonlinear degenerate coupled PDE-PDE and PDE-ODE systems that
arise, for example, in the modelling of biofilm growth. One of the equations,
describing the evolution of a biomass density, exhibits degenerate and
singular diffusion. The other equations are either of advection-reaction-
diffusion type or ordinary differential equations. Under very general
assumptions the existence of weak solutions is proven by considering
regularized systems, deriving uniform bounds and using fixed point arguments.
Assuming additional structural assumptions we also prove the uniqueness of
solutions.

Global-in-time well-posedness is established for Dirichlet and mixed boundary
conditions, whereas, only local well-posedness can be shown for homogeneous
Neumann boundary conditions. Using a suitable barrier function and comparison
theorems we formulate sufficient conditions for finite-time blow-up or uniform
boundedness of solutions. Finally, we show that solutions of the degenerate
parabolic equation inherit additional global spatial regularity if the
diffusion coefficient has a power-law growth.

Keywords: degenerate diffusion $\bullet$ biofilm models $\bullet$ quasilinear
parabolic systems $\bullet$ PDE-ODE systems $\bullet$ well-posedness $\bullet$
regularity $\bullet$ finite time blow up

MSC: 35K65, 35K59, 35A01, 35A02, 35B44, 35B45, 35B50

###### Contents

  1. 1 Introduction
  2. 2 Problem formulation
    1. 2.1 Preliminaries
    2. 2.2 Assumptions on the data
    3. 2.3 Weak solutions
  3. 3 Well-posedness for Dirichlet and mixed boundary conditions
    1. 3.1 A regularized problem
    2. 3.2 Existence of solutions of the degenerate parabolic problem
    3. 3.3 A contraction argument for proving Theorem 3.2
    4. 3.4 A fixed point argument for proving Theorem 3.1
  4. 4 Homogeneous Neumann boundary conditions
    1. 4.1 Existence of weak solutions
      1. 4.1.1 Well-posedness of backward Euler time–discretizations
      2. 4.1.2 Interpolations in time
      3. 4.1.3 Proof of Lemma 3.1
    2. 4.2 Finite time blow-up
  5. 5 Spatial regularity of the biomass density
    1. 5.1 Some auxiliary functions
    2. 5.2 Boundedness of $M$ in $L^{2}(0,T;H^{r}(\Omega))$

##  1 Introduction

This paper investigates the well-posedness and qualitative properties of weak
solutions of a wide class of quasilinear parabolic systems where one of the
equations shows degenerate and singular diffusion. We also consider couplings
of such degenerate parabolic equations with ordinary differential equations
(ODEs). The motivation for our work is models describing the growth of
spatially heterogeneous biofilms in dependence of growth limiting substrates.
The models are either formulated as systems of partial differential equations
(PDEs) or as coupled PDE-ODE systems, e.g. see [7, 6]. Their characteristic
and challenging features are the degenerate and singular diffusion effects in
the equation for the biomass density and the nonlinear coupling of this
equation to additional ODEs and/or PDEs for the substrates.

Let $\Omega\subset{\mathbb{R}}^{d},\,d\in{\mathbb{N}},$ be a bounded Lipschitz
domain and $T>0$. We denote the parabolic cylinder by $Q:=\Omega\times(0,T]$.
Throughout this study, for a fixed $k\in{\mathbb{N}}$, $j\in\\{1,\dots,k\\}$
will denote an integer, and $\vec{w}=(w_{1},\dots,w_{k})$ a $k$-dimensional
vector. We consider the following problem in $Q$,

$\displaystyle\partial_{t}M=\nabla\cdot[D_{0}(M)\nabla M]+f_{0}(M,\vec{S}),$
(1.1a)
$\displaystyle\partial_{t}S_{j}=\nu_{j}\nabla\cdot[D_{j}(M,\vec{S})\nabla
S_{j}+\bm{v}_{j}S_{j}]+f_{j}(M,\vec{S}),$ (1.1b) for $j=1,\dots,k$, where
$M:Q\to{\mathbb{R}}$ denotes the biomass density and the vector-valued
function $\vec{S}:Q\to{\mathbb{R}}^{k}$ the substrate concentrations. The
biomass density $M$ is normalized with respect to the maximum biomass density
and hence, it takes values in $[0,1)$. The biomass diffusion coefficient
$D_{0}:[0,1)\to[0,\infty)$ is degenerate, it satisfies $D_{0}(0)=0$ and
$\lim_{m\nearrow 1}D_{0}(m)=\infty$. Although, we remark that large parts of
our analysis are also valid for non-degenerate functions $D_{0}$. The
diffusion coefficients of the substrates
$D_{j}:[0,1]\times{\mathbb{R}}^{k}\to[0,\infty)$ are non-degenerate, i.e. they
are bounded from above and below by positive constants. The constants
$\nu_{j}\geq 0$ will be referred to as the mobility coefficients of the
substrates. It is important to point out that the case of immobilized
substrates ($\nu_{j}=0$) is included in our setting which leads to a coupling
of Equation (1.1a) with ODEs in (1.1b). Moreover,
$\bm{v}_{j}:Q\to{\mathbb{R}}^{d}$ is a given flow-field. Finally, the reaction
terms $f_{0},\,f_{j}:{\mathbb{R}}^{k+1}\to{\mathbb{R}}$ describe the complex
interplay between the substrates and biomass.

In biofilm modelling applications, it is important to allow for mixed
Dirichlet-Neumann or homogeneous Neumann boundary conditions for $M$. To this
end, we divide the boundary $\partial\Omega$ into two disjoint parts
$\Gamma_{1}$ and $\Gamma_{2}$ that are both Lipschitz boundaries. We
complement (1.1a)–(1.1b) with the following initial and boundary conditions
for $M$ and $\vec{S}$,

$\displaystyle M(0)$ $\displaystyle=M_{0},\quad\vec{S}(0)=\vec{S}_{0},$ (1.1c)
$\displaystyle M|_{\Gamma_{1}}$ $\displaystyle=h_{0},\quad[\nabla
M\cdot\bm{\hat{n}}]|_{\Gamma_{2}}=0,\quad\nu_{j}S_{j}|_{\partial\Omega}=\nu_{j}h_{j},$
(1.1d)

where $\bm{\hat{n}}$ denotes the outward unit normal to $\partial\Omega$ and
$M_{0}:\Omega\to[0,1)$, $\vec{S}_{0}:\Omega\to{\mathbb{R}}^{k}$,
$h_{0}:\Gamma_{1}\to[0,1)$ and $h_{j}:\partial\Omega\to{\mathbb{R}}$ are
given. We remark that the case $\Gamma_{1}=\emptyset$ is allowed in our
setting which corresponds to homogeneous Neumann boundary conditions for $M$.
The case $\Gamma_{2}=\emptyset$ is also included which corresponds to
Dirichlet boundary conditions for $M$. Note that in (1.1c) we do not prescribe
boundary conditions for immobilized substrates $S_{j}$, i.e. if $\nu_{j}=0$
for some $j\in\\{1,\dots,k\\}.$ To simplify the presentation of our results we
assume Dirichlet boundary conditions for the substrates, but the analysis
remains valid if we impose mixed boundary conditions for the substrates, see
Remark 2.4.

In models for biofilm growth, the actual biofilm is described by the region
where $M$ is positive,

$\Omega^{+}(t)=\\{x\in\Omega:M(t,x)>0\\}.$

Due to the degeneracy of the biomass diffusion coefficient, $D_{0}(0)=0$,
there is a sharp interface between the biofilm and the surrounding region, and
the interface propagates at a finite speed. The additional singularity in the
diffusion coefficient, $\lim_{m\nearrow 1}D_{0}(m)=\infty$, ensures that the
biomass density does not exceed its maximum value, i.e. $M$ remains bounded by
a constant strictly less than $1$.

In Figure 1 typical situations modelled by (1.1) are sketched for biofilm
colonies depending on a single substrate. In the left figure the substrate is
dissolved in the spatial domain $\Omega$ and transported by diffusion and
convection. The biofilm colony grows into the aqueous phase. In the right
figure the substrate is immobilized and contained in the spatial domain
$\Omega$. The bacteria consume and degrade the substrate, a biofilm front
develops and propagates through the substratum.

A system of the form (1.1) with a single dissolved substrate $S=S_{1}$, i.e.
$k=1$ and $\nu_{1}>0$, was first proposed in [7] to model biofilm growth in an
aqueous medium. In this case,

$\displaystyle D_{0}(M)$ $\displaystyle=d_{2}\frac{M^{a}}{(1-M)^{b}},$
$\displaystyle D_{1}(M,S)$ $\displaystyle=d_{1},$ (1.2a) $\displaystyle
f_{0}(M,S)$ $\displaystyle=k_{3}\frac{SM}{k_{4}+S}-k_{2}M,$ $\displaystyle
f_{1}(M,S)$ $\displaystyle=-k_{1}\frac{SM}{k_{4}+S},$ (1.2b)

for some constants $k_{1},k_{2},k_{3},k_{4},d_{1},d_{2}>0$ and $a,b\geq 1$,
and $\mathbf{v_{1}}$ is a given flow field. An ODE-PDE system of the form
(1.1) with a single substrate was used in [6] to model cellulolytic biofilms
degrading an immobilized cellulose material. In this case, the functions
$D_{0},f_{0}$ and $f_{1}$ are as in (1.2) and $D_{1}\equiv 0$,
$\mathbf{v_{1}}\equiv 0$.

(a) Coupled PDE-PDE systems

(b) Coupled PDE-ODE systems

Figure 1:  Schematic figures illustrating biofilm growth in dependence of a
single nutrient $S$ in an aqueous medium (a) and in an immobilized medium (b).
The biofilm is represented by the region where $M(\bm{x},t)>0$, which is
separated by the surrounding region by a sharp interface. Nutrients are
consumed by bacteria resulting in the production of biomass. The parts of the
boundary where homogeneous Dirichlet and Neumann conditions are specified are
also marked in the diagrams, $\Gamma_{1}$ (Dirichlet) in blue and $\Gamma_{2}$
(Neumann) in red.
(a) PDE-PDE systems [7]: The biofilm colonies grow in a liquid containing
substrates. The substrates diffuse and are transported by a flow field
$\bm{v}$. The diffusion coefficient of the substrate might depend on $M$, i.e.
it differs inside and outside the biofilm.
(b) PDE-ODE systems [6]: The bacteria degrade and consume an immobilized
medium which is the case, e.g. for cellulolytic biofilms. The biofilm colony
propagates consuming the immobile cellulose, leaving at its wake a region of
low substrate concentrations.

The existence of weak solutions of scalar nonlinear degenerate parabolic
equations such as (1.1a) was shown in the seminal papers [1, 2], however, for
bounded diffusion coefficients $D_{0}$. Uniqueness of solutions was proven in
[21] using $L^{1}$-contraction. The existence of weak solutions for the
biofilm model [7] with the diffusion coefficients and reaction functions in
(1.2) and $\mathbf{v_{1}}\equiv 0$ was proven in [9] under the assumptions of
homogeneous Dirichlet boundary conditions for $M$, i.e. $\Gamma_{2}=\emptyset$
and $h_{0}\equiv 0$. The existence of the global attractor for the generated
semigroup in $L^{1}(\Omega)$ was also shown. The well-posedness theory was
generalized in [15] where more general functions $D_{0}$, $f_{0}$ and $f_{1}$
and mixed Dirichlet-Neumann boundary conditions were considered. The Hölder
continuity of solutions was studied in [14].

Several extensions and variations of the single species biofilm growth model
[7] have been proposed and analyzed. Most works are simulation studies and
only few analytical results have been obtained. The well-posedness of multi-
substrate biofilm models with $k>1$, $\nu_{j}>0$ in (1.2), appearing in
antibiotic disinfection and quorum sensing applications, was established in
[25, 10]. A PDE–ODE system with an immobile substrate, i.e. $k=1$ and
$\nu_{1}=0$, was proposed and numerically studied in [6]. The simulations
reproduced many experimentally observed features of cellulolytic biofilms. The
existence and stability of travelling wave solutions for this model were shown
in [19], but the well-posedness of the model remained an open problem. Many
examples of semilinear coupled PDE-ODE models appearing in biology are
discussed in [20, Chapter 13]. For a PDE–ODE model for hysteretic flow through
porous media with a diffusion coefficient $D_{0}$ depending on both $M$ and
$\vec{S}$, the existence of solutions was shown in [18].

We aim to develop a unifying solution theory for a large class of systems with
degenerate diffusion that is motivated by models for biofilm growth, but the
analysis is not limited to these applications. In fact, we expect that such
models can also be used, e.g. to describe cancer cell invasion or the spread
of wildfires. In our paper we extend previous well-posedness results in the
following directions:

(a) Well-posedness results for PDE-PDE systems: our results extend the theory
developed for systems with one substrate in [15] to systems with an arbitrary
number of substrates $k\in{\mathbb{N}}$. Moreover, the existence of weak
solutions is proven for a broad class of diffusion coefficients $D_{0}$ and
$D_{j}$, reaction terms $f_{j}$, and allows for flow-fields $\bm{v}_{j}$ which
has not been considered in earlier works.

(b) Well-posedness of PDE–ODE systems ($\nu_{j}=0$): The well-posedness of
PDE-ODE systems of the form (1.1) with a degenerate and/or singular diffusion
coefficient $D_{0}$ has been an open problem. The theory we develop applies to
the cellulolytic biofilm model [6] and implies its local well-posedness.

(c) Mixed as well as homogeneous Neumann conditions for $M$: Global well-
posedness is shown for mixed Dirichlet-Neumann boundary conditions and a local
well-posedness result is established assuming homogeneous Neumann boundary
conditions for $M$. Moreover, apart from well-posedness results we also
analyze qualitative properties such as boundedness or blow-up of solutions.

(d) Global spatial regularity of M: We further show that under certain porous
medium type growth conditions on $D_{0}$ close to zero, the biomass
concentration $M$ inherits some global spatial regularity.

The outline of our paper is as follows: In Section 2 we introduce notation,
state our assumptions on the data and introduce the concept of weak solutions.
In Section 3 we prove global well-posedness for systems with Dirichlet or
mixed Dirichlet-Neumann boundary conditions for $M$. In Section 4 we establish
local well-posedness for systems with homogeneous Neumann conditions for $M$.
We also derive criteria ensuring finite-time blow-up of the model and discuss
some important examples. In Section 5 we show that even in the degenerate
case, the biomass density $M$ possesses some global spatial regularity.

##  2 Problem formulation

In this section, we introduce notation and a suitable functional framework. We
state the properties of the coefficient functions and the boundary and initial
data for system (1.1) that will be assumed throughout the paper. Moreover, we
introduce weak solutions of the problem.

###  2.1 Preliminaries

##### Functional setting:

Let $\Omega\subset{\mathbb{R}}^{d}$ be a bounded Lipschitz domain. The
boundary $\partial\Omega$ is divided into two regular open subsets
$\Gamma_{1}$ and $\Gamma_{2}$ that are both Lipschitz boundaries and such that
$\partial\Omega=\overline{\Gamma}_{1}\cup\overline{\Gamma}_{2}$ and
$\Gamma_{1}\cap\Gamma_{1}=\emptyset$, e.g. see [22]. We denote by
$(\cdot,\cdot)$ and $\|\cdot\|$ the $L^{2}(\Omega)$ inner product and norm.
The norm of any other Banach space $V$ will be denoted by $\|\cdot\|_{V}$. For
$1\leq p\leq\infty$, let $W^{1,p}(\Omega)$ denote the Sobolev space of
functions $u\in L^{p}(\Omega)$ such that the weak derivative $\nabla u$ exists
and $\nabla u\in(L^{p}(\Omega))^{d}$. For $r\in(0,1)$ and $p\in[0,\infty)$,
the Sobolev-Slobodeckij space $W^{r,p}(\Omega)$ is the set of functions $u\in
L^{p}(\Omega)$ such that

$\displaystyle\|u\|_{W^{r,p}(\Omega)}:=\|u\|_{L^{p}(\Omega)}+\int_{\Omega}\int_{\Omega}\dfrac{|u(\bm{x})-u(\bm{y})|^{p}}{|\bm{x}-\bm{y}|^{d+rp}}d\bm{x}d\bm{y}<\infty.$
(2.1)

We define $H^{r}(\Omega):=W^{r,2}(\Omega)$ for $r\in(0,1]$. Let
$H^{1}_{0}(\Omega)$ denote the closure of $C^{\infty}_{c}(\Omega)$ in
$H^{1}(\Omega)$, which is equipped with the norm
$\|u\|_{H^{1}_{0}(\Omega)}:=\|\nabla u\|$. Similarly, we define

$\displaystyle{\cal H}^{1}:=\\{u\in H^{1}(\Omega):\;\mathrm{tr}(u)=0\text{ in
}\Gamma_{1}\\}\quad\text{with the norm}\quad\|u\|_{{\cal
H}^{1}}:=\|u\|_{H^{1}(\Omega)}.$ (2.2a) The dual spaces of $H^{1}_{0}(\Omega)$
and ${\cal H}^{1}$ are defined as: $\displaystyle
H^{-1}:=(H^{1}_{0}(\Omega))^{*}\quad\text{ and }\quad{\cal H}^{-1}=({\cal
H}^{1})^{*}.$ (2.2b) Observe that, $\displaystyle\text{ if
}\Gamma_{1}=\emptyset\text{ then }{\cal H}^{1}=H^{1}(\Omega),\quad\text{ if
}\Gamma_{2}=\emptyset\text{ then }{\cal H}^{1}=H^{1}_{0}(\Omega).$ (2.2c)

Let $\langle\cdot,\cdot\rangle$ denote the duality pairing of ${\cal H}^{1}$
and ${\cal H}^{-1}$. The duality pairing of any other Sobolev space $V$ will
be denoted by $\langle\cdot,\cdot\rangle_{V,V^{*}}$.

Finally, we introduce the following Bochner spaces that are important for our
analysis:

$\displaystyle{\bf{\cal W}}:=L^{\infty}(0,T;L^{\infty}(\Omega))\cap
H^{1}(0,T;{\cal H}^{-1})\cap C([0,T];L^{2}(\Omega)),$ (2.3a)
$\displaystyle{\bf{\cal X}}:=L^{2}(0,T;{\cal H}^{1})\cap H^{1}(0,T;{\cal
H}^{-1}),$ (2.3b) $\displaystyle{\bf{\cal
Y}}:=L^{\infty}(0,T;H^{1}(\Omega))\cap H^{1}(0,T;L^{2}(\Omega)),$ (2.3c)
$\displaystyle{\bf{\cal Z}}:=C([0,T];(L^{2}(\Omega))^{k}).$ (2.3d)

Note that we have the continuous embedding ${\bf{\cal X}}\hookrightarrow
C([0,T];L^{2}(\Omega))$.

##### Inequalities:

Note that the Poincar$\acute{\text{e}}$ inequality, i.e. $\|u\|\leq
C_{\Omega}\|\nabla u\|$ for $u\in H^{1}_{0}(\Omega)$, where $C_{\Omega}>0$
denotes the Poincaré constant, also holds for functions $u\in{\cal H}^{1}$ if
$\Gamma_{1}\not=\emptyset$.

We recall Young’s inequality stating that for any $\sigma>0$ one has

$ab\leq\frac{1}{2\sigma}a^{2}+\frac{\sigma}{2}b^{2}\qquad\forall
a,b\in{\mathbb{R}}.$ (2.4)

We will also frequently use Gronwall’s Lemma stating that if $u,\,a,\,b\in
C({\mathbb{R}})$ are non-negative, then $u(t)\leq
a(t)+\int_{0}^{t}u(\varrho)\,b(\varrho)\,\mathrm{d}\varrho$ implies that

$u(t)\leq
a(t)+\int_{0}^{t}a(\varrho)\,b(\varrho)e^{\int_{\varrho}^{t}b(\tau)d\tau}\,\mathrm{d}\varrho$
(2.5a) for all $t>0$; and the discrete counterpart of the 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.5b)

Finally, for a convex $\eta\in C({\mathbb{R}}^{+})$ with $\eta(0)=0$ we will
use Jensen’s inequality and the super-additivity property:

Jensen’s inequality:
$\displaystyle\eta\left(\tfrac{1}{|\Omega|}\int_{\Omega}|f|\right)\leq\tfrac{1}{|\Omega|}\int_{\Omega}\eta(|f|)\qquad\text{
for }f\in L^{1}(\Omega);$ (2.6a) Super-additivity:
$\displaystyle\eta(a)+\eta(b)\leq\eta(a+b)\qquad\text{ for all }a,\,b\geq 0.$
(2.6b)

##### Further notation:

We denote by $[\cdot]_{+}$ and $[\cdot]_{-}$ the positive and negative part of
functions, i.e. $[\cdot]_{+}:=\max\\{\cdot,0\\}$ and
$[\cdot]_{-}:=\min\\{\cdot,0\\}$, respectively. By $C>0$ we refer to an
undisclosed constant in the estimates that may vary in each occurrence and
from line to line. Finally, the notation

$\displaystyle a\lesssim b\quad\text{ implies that }a\leq Cb\quad\text{for
some constant }C>0$ (2.7)

which does not depend on a parameter $\varepsilon>0$ (to be specified later).

###  2.2 Assumptions on the data

We specify the hypotheses on the data associated with (1.1).

  1. (P1)

The diffusion coefficient $D_{0}:[0,1)\mapsto[0,\infty)$ is a continuous
function that is strictly increasing in $[0,\epsilon_{0})$ for some
$\epsilon_{0}\in(0,1]$, and satisfies

$D_{0}(0)=0,\;\;\lim\limits_{m\nearrow 1}D_{0}(m)=\infty\;\text{ and
}\;\;D_{0}(m)>0\text{ for all }m\in(0,1).$

The primitive of $D_{0}$, expressed by the Kirchhoff transform function
$\Phi:[0,1)\to[0,\infty)$,
$\Phi(m)=\int_{0}^{m}D_{0}(\varrho)\,\mathrm{d}\varrho,$ satisfies

$\displaystyle\lim\limits_{m\nearrow 1}\Phi(m)=\infty.$ (2.8)

  2. (P2)

The diffusion coefficients
$D_{j}:[0,1]\times{\mathbb{R}}^{k}\to[D_{\min},D_{\max}]$, $j=1,\dots,k,$ with
constants $0<D_{\min}<D_{\max}<\infty$, are Lipschitz continuous with respect
to both variables.

######  Remark 2.1 (Biofilm models).

In models for biofilm growth, see e.g. [6, 9], the diffusion coefficient
$D_{0}$ is given by the function in (1.2), and the diffusion coefficient
$D_{1}$ for the (single) substrate is assumed to be constant. These functions
satisfy all assumptions in (P1)–(P2).

######  Remark 2.2 (Generalizations of the assumptions (P1)–(P2)).

Our analysis can be extended to systems where the diffusion coefficient
$D_{0}$ is piecewise constant, non-degenerate, and/or has a porous media type
degeneracy, e.g., $D_{0}(m)=m^{a}$, for some constant $a>1$. To keep the
analysis uniform and self-contained, we only analyze the case (P1) which is
more involved and arises in models for biofilm growth.

We could also allow for degenerate diffusion coefficients $D_{j}$ if they only
depend on the substrate $S_{j}$. Some additional assumptions are required to
cover this case which are discussed in Corollary 3.2.1.

For the flow-field and reaction terms we make the following assumptions:

  1. (P3)

The flow-field satisfies $\bm{v}_{j}\in(L^{\infty}(Q))^{d}$, $j=1,\dots,k$.

  2. (P4)

The functions $f_{0},\,f_{j}\in C([0,1]\times{\mathbb{R}}^{k})$ are uniformly
Lipschitz continuous. They can be extended to uniformly Lipschitz continuous
functions on ${\mathbb{R}}^{k+1}$ which (to simplify notation) we will also
denote by $f_{0},\,f_{j}$. The constant $C_{L}\geq 0$ is the maximum of the
Lipschitz constants of $f_{0},f_{1},\dots,f_{k}$. Moreover,
$f_{0}(0,\vec{s})\geq 0$ for all $\vec{s}\in{\mathbb{R}}^{k}$.

  1. (P4enumi)

There exists a non-negative and locally Lipschitz continuous function
$f_{max}\in C({\mathbb{R}})$ such that $f_{0}(\cdot,\vec{s})\leq
f_{\max}(\cdot)$ for all $\vec{s}\in{\mathbb{R}}^{k}$.

######  Remark 2.3 (Assumptions (P4)–(P4enumi)).

Assumption (P4) admits reaction functions $f_{0}(\cdot,\vec{s})$,
$f_{j}(\cdot,\vec{s})$ (for $\vec{s}\in{\mathbb{R}}^{k}$) that have
superlinear growth with respect to their first argument as long as they are
Lipschitz continuous within the interval $[0,1]$. This is because the
physically relevant solutions satisfy $M\in[0,1]$. However, before proving the
upper bound for $M$ (in Lemma 3.2), we need the functions $f_{0}$ and $f_{j}$
to be defined in ${\mathbb{R}}^{k+1}$ (in Lemma 3.1) which is why we introduce
the extensions.

Assumption (P4enumi) is needed to derive the $L^{\infty}$ bound for the
solution $M$. This is important in our setting since physically relevant
solutions take values in $[0,1)$, otherwise, the models are not valid. Such
$L^{\infty}$ bounds may not be required in other applications, for example for
porous medium-type equations. We therefore explicitly state in all theorems
and lemmas where this assumption is required and where it can be omitted.

Under additional assumptions, the analysis can be generalized also to systems
with reaction functions that depend on $\bm{x}$ and $t$, e.g. see [15]. To
simplify the presentation of our results, we omit this dependency here.

Finally, we remark that the condition $f_{0}(\cdot,\vec{s})\leq
f_{\max}(\cdot)$ can be relaxed to $f_{0}(\cdot,\vec{s})\leq
g(|\vec{s}|)\,f_{\max}(\cdot)$ for some $g\in C({\mathbb{R}}^{+})$. Then, for
the proofs to go through, we need uniform $L^{\infty}$ bounds on the solution
$(M,\vec{S})$ which can be established for a certain class of functions
$f_{j}$. For an example, we refer to Corollary 3.2.1.

For the boundary and initial data we assume the following properties:

  1. (P5)

The initial data $M_{0}\in L^{\infty}(\Omega)$ satisfies

$\underline{M}:=\mathrm{ess}\inf\limits_{\bm{x}\in\Omega}\\{M_{0}\\}\geq
0\quad\text{and}\quad\overline{M}:=\mathrm{ess}\sup\limits_{\bm{x}\in\Omega}\\{M_{0}\\}<1.$

The initial data $\vec{S}_{0}\in(L^{\infty}(\Omega))^{k}$ satisfies

$\underline{S}:=\min\limits_{1\leq j\leq
k}\mathrm{ess}\inf\limits_{\bm{x}\in\Omega}\\{S_{0,j}\\}>-\infty\quad\text{and}\quad\overline{S}:=\max\limits_{1\leq
j\leq k}\mathrm{ess}\sup\limits_{\bm{x}\in\Omega}\\{S_{0,j}\\}<\infty.$

  2. (P6)

The Dirichlet boundary data $h_{0}:\Gamma_{1}\to[\underline{M},\overline{M}]$
is such that there exists $h^{e}_{0}\in H^{1}(\Omega)$ satisfying
$h^{e}_{0}|_{\Gamma_{1}}=h_{0}$ in a trace sense. If $\Gamma_{1}=\emptyset$
then set $h^{e}_{0}\equiv 0$. For the Dirichlet data
$h_{j}:\partial\Omega\to[\underline{S},\overline{S}]$ there also exist
functions $h^{e}_{j}\in H^{1}(\Omega)$ such that
$h^{e}_{j}|_{\partial\Omega}=h_{j}$ in a trace sense.

######  Remark 2.4 (Assumption (P6)).

Observe that, under the assumptions in (P6), it is always possible to choose
the extensions to $\Omega$ such that
$h^{e}_{0}\in[\underline{M},\overline{M}]$ and
$h^{e}_{j}\in[\underline{S},\overline{S}]$ a.e. in $\Omega$. For example,
consider $\bar{h}^{e}=\min\\{h^{e}_{0},\overline{M}\\}\in H^{1}(\Omega)$. Then
$\bar{h}^{e}\leq\overline{M}$ a.e. and $\bar{h}^{e}=h_{0}$ on
$\partial\Omega$. This choice will implicitly be used in the proofs that
follow. Similar arguments apply to the boundary conditions $h_{j}.$

To keep notations simple we only consider Dirichlet boundary conditions for
the substrates $\vec{S}$. Mixed Dirichlet-/Neumann boundary conditions with
different divisions of the boundary depending on $j$ can also be assumed
without major modifications in the subsequent arguments, see e.g. [15].

###  2.3 Weak solutions

We introduce the following notion of weak solutions.

######  Definition 1 (Weak solution).

The pair $(M,\vec{S})$ with $M\in{\bf{\cal W}}$ (see (2.3)), $\Phi(M)\in
L^{2}(0,T;H^{1}(\Omega))$, $S_{j}\in H^{1}(0,T;H^{-1}(\Omega))\cap
C([0,T];L^{2}(\Omega))$, and $\nu_{j}S_{j}\in L^{2}(0,T;H^{1}(\Omega))$,
$j=1,\dots,k,$ is a weak solution of (1.1) provided that $M$ is bounded in
$[0,1)$ a.e. in $Q$, $M(0)=M_{0}$ and $\vec{S}(0)=\vec{S}_{0}$ a.e. in
$\Omega$, $\Phi(M)=\Phi(h_{0})$ on $\Gamma_{1}$ and
$\nu_{j}S_{j}=\nu_{j}h_{j}$ on $\partial\Omega$ in the trace sense, and for
all $\varphi\in L^{2}(0,T;{\cal H}^{1})$, $\vec{\zeta}\in
L^{2}(0,T;(H^{1}_{0}(\Omega))^{k})$, we have

$\displaystyle\int_{0}^{T}\langle\varphi,\partial_{t}M\rangle+\int_{0}^{T}(\nabla\Phi(M),\nabla\varphi)=\int_{0}^{T}(f_{0}(M,\vec{S}),\varphi),$
(2.9a)
$\displaystyle\int_{0}^{T}\langle\zeta_{j},\partial_{t}S_{j}\rangle_{H^{1}_{0},H^{-1}}+\nu_{j}\int_{0}^{T}(D_{j}(M,\vec{S})\nabla
S_{j}+\bm{v}_{j}\,S_{j},\nabla\zeta_{j})=\int_{0}^{T}(f_{j}(M,\vec{S}),\zeta_{j}).$
(2.9b)

######  Remark 2.5.

In Equation 2.9 we take $\nu_{j}\,S_{j}\in L^{2}(0,T;H^{1}(\Omega))$ instead
of taking $S_{j}\in L^{2}(0,T;H^{1}(\Omega))$. This is required since
$\nu_{j}=0$ is allowed in our setting, and in this case, $S_{j}$ might not
possess any spatial regularity. Similarly, we see that $\Phi(M)$, and not $M$,
possesses spatial regularity. Therefore, the traces are also only defined for
functions with sufficient spatial regularity.

##  3 Well-posedness for Dirichlet and mixed boundary conditions

In this section, we prove the well-posedness of weak solutions for the case
when $\Gamma_{1}$ has non-zero measure, i.e. $M$ either satisfies Dirichlet
boundary condition or mixed Dirichlet-Neumann boundary condition. The main
results of this section are stated in the following theorems.

######  Theorem 3.1 (Existence and boundedness).

Let (P1)–(P6) and (P4enumi) be satisfied, and $\Gamma_{1}$ have non-zero
measure. Then, there exists a weak solution $(M,\vec{S})$ of(1.1) in the sense
of Equation 2.9. Furthermore, a constant $\delta\in(0,1)$ exists such that
$0\leq M\leq 1-\delta$ a.e. in $Q$.

######  Theorem 3.2 (Uniqueness).

Let the assumptions of Theorem 3.1 hold. In addition, for each
$j\in\\{1,\dots,k\\}$, assume that either $\nu_{j}=0$ or the diffusion
coefficient $D_{j}$ depends only on $S_{j}$. Then, a unique weak solution
$(M,\vec{S})$ of(1.1) exists in the sense of Equation 2.9.

The proof of Theorem 3.2 is based on a contraction argument which, along with
the existence of solutions, guarantees uniqueness. However, a different
argument based on Schauder’s fixed point theorem is required for proving the
existence of solutions in the more general setting of Theorem 3.1.

For the proof of these theorems, we initially focus on the first equation in
(2.9) for a given $\vec{S}$. In Section 3.1 we consider a non-degenerate
approximation of (1.1a) and then discuss existence (Lemma 3.1) and boundedness
(Lemma 3.2) of solutions. In Section 3.2, we pass the regularization parameter
to zero to show the existence of weak solutions of the original problem (Lemma
3.3). We then consider the coupled system. In Section 3.3, the $L^{1}$
contraction principle (Lemma 3.4) is applied to prove Theorem 3.2, and in
Section 3.4, a Schauder argument (Lemma 3.5) is used to prove Theorem 3.1.

We remark that Lemmas 3.1–3.5 hold for all boundary conditions including
homogeneous Neumann boundary condition. The proof of Lemma 3.1 is postponed to
Section 4, due to complications arising from homogeneous Neumann condition.
Lemmas 3.2–3.5 are proven in the general case.

###  3.1 A regularized problem

We introduce the following regularization of the Kirchhoff transform $\Phi$:
for $\varepsilon>0$, let $\Phi_{\varepsilon}\in C^{1}({\mathbb{R}})$ be a non-
degenerate approximation of $\Phi$ satisfying

$\displaystyle\varepsilon\leq{\Phi_{\varepsilon}}^{\prime}\leq\varepsilon^{-1}\quad\text{
and }\quad\lim\limits_{\varepsilon\to
0}\Phi_{\varepsilon}(m)=\Phi(m)\quad\text{ for all }m\in[0,1).$ (3.1)

A specific choice of $\Phi_{\epsilon}$ that will be used in the sequel and in
Section 5 is

$\displaystyle\Phi_{\varepsilon}(m):=\int_{0}^{m}\min\left\\{\max\\{\varepsilon,D(\varrho)\\},\varepsilon^{-1}\right\\}\,\mathrm{d}\varrho.$
(3.2)

Then, recalling the functional spaces defined in (2.3), the following lemma
holds.

######  Lemma 3.1 (Existence for a regularized problem).

Let (P1)–(P6) hold. Let $\vec{s}\in{\bf{\cal Z}}$ be given and
$\varepsilon\in(0,1)$ be sufficiently small. Then there exists a unique
$M_{s,\varepsilon}\in{\bf{\cal X}}+h^{e}_{0}$ which satisfies
$M_{s,\varepsilon}(0)=M_{0}$, and for all $\varphi\in L^{2}(0,T;{\cal
H}^{1})$,

$\displaystyle\int_{0}^{T}\langle\varphi,\partial_{t}M_{s,\varepsilon}\rangle+\int_{0}^{T}(\nabla\Phi_{\varepsilon}(M_{s,\varepsilon}),\nabla\varphi)=\int_{0}^{T}(f_{0}(M_{s,\varepsilon},\vec{s}),\varphi).$
(3.3)

Moreover, $M_{s,\varepsilon}\in C([0,T];L^{2}(\Omega))$, and for all
$t\in[0,T]$ we have

$\displaystyle\|M_{s,\varepsilon}(t)\|^{2}+\int_{0}^{T}\left[\|\nabla\Phi_{\varepsilon}(M_{s,\varepsilon})\|^{2}+\|\partial_{t}M_{s,\varepsilon}\|^{2}_{{\cal
H}^{-1}}\right]$ $\displaystyle\lesssim
1+\int_{0}^{T}\left(\|\vec{s}\|^{2}+\|\Phi_{\varepsilon}(M_{s,\varepsilon})\|^{2}\right)+\left(1+\|\Phi^{\prime}_{\varepsilon}(M_{s,\varepsilon})\|_{L^{\infty}(Q)}\right)\|h^{e}_{0}\|^{2}_{H^{1}(\Omega)}.$
(3.4)

Furthermore, if $M_{0}\in H^{1}(\Omega)$ in (P5), and
$M_{0}|_{\Gamma_{1}}=h_{0}$ in the trace sense then, in addition, it holds
that

$\displaystyle\|\nabla\Phi_{\varepsilon}(M_{s,\varepsilon}(t))\|^{2}+\int_{0}^{T}\int_{\Omega}\frac{|\partial_{t}\Phi(M_{s,\varepsilon})|^{2}}{\Phi^{\prime}_{\varepsilon}(M_{s,\varepsilon})}$
$\displaystyle\lesssim\|\nabla\Phi_{\varepsilon}(M_{0})\|^{2}+\|\Phi^{\prime}_{\varepsilon}(M_{s,\varepsilon})\|_{L^{\infty}(Q)}\left[1+\int_{0}^{T}\|\vec{s}\|^{2}\right].$
(3.5)

The proof of Lemma 3.1 is postponed to Section 4.1. For $\Gamma_{1}$ having
non-zero measure, the existence of $M_{s,\varepsilon}\in{\bf{\cal
X}}+h^{e}_{0}$ follows immediately from [1] since
${\Phi_{\varepsilon}}^{\prime}$ satisfies the uniform ellipticity condition.
However, the result in [1] does not cover the case of homogeneous Neumann
conditions, i.e. the case $\Gamma_{1}=\emptyset$.

The assumption (P4enumi) is not required in Lemma 3.1, but it is needed for
the next result.

######  Lemma 3.2 (Boundedness for the regularized problem).

In addition to the hypothesis of Lemma 3.1 we assume that (P4enumi) holds. Let
$M_{s,\varepsilon}\in{\bf{\cal X}}+h^{e}_{0}$ denote the solution in Lemma
3.1. Moreover, let $\hat{M}\in C^{1}({\mathbb{R}}^{+})$ denote the solution of
the integral equation

$\hat{M}(t)=\overline{M}+\int_{0}^{t}f_{\max}(\hat{M}(\varrho))\,\mathrm{d}\varrho,\qquad
t\in[0,T].$

  1. (a)

Then, $0\leq M_{s,\varepsilon}(t)\leq\hat{M}(t)$ a.e. in $\Omega$ for all
$t\in[0,T]$.

  2. (b)

If, in addition, $\Gamma_{1}$ has non-zero measure, then a constant
$\delta\in(0,1)$ exists such that

$0\leq M_{s,\varepsilon}(t)\leq 1-\delta\ \text{ a.e. in }\Omega\text{ for all
}t\in[0,T].$

Observe that the above lemma implies that $M_{s,\varepsilon}\in
L^{\infty}(0,T;L^{\infty}(\Omega))$ and the family $M_{s,\varepsilon}$ is
uniformly bounded with respect to $\varepsilon>0$ in
$L^{\infty}(0,T;L^{\infty}(\Omega))$.

###### Proof.

The existence of $\hat{M}$ follows from the Picard-Lindelöf Theorem since
$f_{\max}\in\mathrm{Lip}({\mathbb{R}})$. Moreover, it satisfies
$\partial_{t}\hat{M}=f_{\max}(\hat{M})$ and therefore,
$\hat{M}\geq\overline{M}$ since $f_{\max}$ was assumed to be non-negative in
(P4enumi).

(Step 1) $\mathbf{M_{s,\varepsilon}\geq 0}$: Inserting the test function
$\varphi=[M_{s,\varepsilon}]_{-}$ in (3.3) implies that

$\int_{0}^{T}\left[\partial_{t}\left(\tfrac{1}{2}\|[M_{s,\varepsilon}]_{-}\|^{2}\right)+\varepsilon\|\nabla[M_{s,\varepsilon}]_{-}\|^{2}\right]\overset{\ref{prop:fg1}}{\leq}C_{L}\int_{0}^{T}\|[M_{s,\varepsilon}]_{-}\|^{2}.$

Since $[M_{s,\varepsilon}(0)]_{-}=[M_{0}]_{-}=0$, we have
$\|[M_{s,\varepsilon}]_{-}(T)\|=0$ using Gronwall’s Lemma (2.5a).

(Step 2) $\mathbf{M_{s,\varepsilon}\leq\hat{M}}$: Inserting the test function
$\varphi=[M_{s,\varepsilon}-\hat{M}]_{+}\in L^{2}(0,T;{\cal H}^{1})$ in (3.3)
we obtain

$\displaystyle\int_{0}^{T}\langle[M_{s,\varepsilon}-\hat{M}]_{+},\partial_{t}M_{s,\varepsilon}\rangle=\int_{0}^{T}\langle[M_{s,\varepsilon}-\hat{M}]_{+},\partial_{t}[M_{s,\varepsilon}-\hat{M}]\rangle+\int_{0}^{T}\langle[M_{s,\varepsilon}-\hat{M}]_{+},\partial_{t}\hat{M}\rangle$
$\displaystyle\qquad\,=\int_{0}^{T}\partial_{t}\left(\frac{1}{2}\|[M_{s,\varepsilon}-\hat{M}]_{+}\|^{2}\right)+\int_{0}^{T}(\partial_{t}\hat{M},[M_{s,\varepsilon}-\hat{M}]_{+})$
$\displaystyle\qquad\overset{\ref{prop:IC}}{=}\frac{1}{2}\|[M_{s,\varepsilon}-\hat{M}]_{+}(T)\|^{2}+\int_{0}^{T}(\partial_{t}\hat{M},[M_{s,\varepsilon}-\hat{M}]_{+}),$
(3.6a)
$\displaystyle\int_{0}^{T}(\nabla\Phi_{\varepsilon}(M_{s,\varepsilon}),\nabla[M_{s,\varepsilon}-\hat{M}]_{+})=\int_{0}^{T}({\Phi_{\varepsilon}}^{\prime}(M_{s,\varepsilon})\nabla
M_{s,\varepsilon},\nabla[M_{s,\varepsilon}-\hat{M}]_{+})\geq 0,$ (3.6b) and
for the reaction term we obtain
$\displaystyle\int_{0}^{T}(f_{0}(M_{s,\varepsilon},\vec{s}),[M_{s,\varepsilon}-\hat{M}]_{+})=\int_{0}^{T}(f_{0}(M_{s,\varepsilon},\vec{s})-f_{0}(\hat{M},\vec{s})+f_{0}(\hat{M},\vec{s}),[M_{s,\varepsilon}-\hat{M}]_{+})$
$\displaystyle\overset{\ref{prop:fg1},\ref{prop:fg2}}{\leq}C_{L}\int_{0}^{T}\|[M_{s,\varepsilon}-\hat{M}]_{+}\|^{2}+\int_{0}^{T}f_{\max}(\hat{M})[M_{s,\varepsilon}-\hat{M}]_{+}.$
(3.6c)

Combining the estimates in (3.6) it follows that

$\displaystyle\frac{1}{2}\|[M_{s,\varepsilon}-\hat{M}]_{+}(T)\|^{2}+\int_{0}^{T}(\partial_{t}\hat{M}-f_{\max}(\hat{M}),[M_{s,\varepsilon}-\hat{M}]_{+})\leq
C_{L}\int_{0}^{T}\|[M_{s,\varepsilon}-\hat{M}]_{+}\|^{2}.$ (3.7)

The second term is zero by the definition of $\hat{M}$. Hence, using
Gronwall’s Lemma (2.5a) we have the result.

(Step 3) $\mathbf{M_{s,\varepsilon}\leq 1-\delta}$: This is a generalization
of Proposition 6 in [9] to the case of mixed or homogeneous Neumann boundary
conditions, see also the proof of Theorem 2.7 in [15]. For $f_{\max}(\cdot)$
introduced in (P4enumi), let $\hat{u}\in{\cal H}^{1}$ solve the elliptic
problem

$\displaystyle(\nabla\hat{u},\nabla\varphi)=(\hat{C},\varphi)\quad\text{for
all}\ \varphi\in{\cal H}^{1},\quad\text{where }\hat{C}:=\max_{0\leq t\leq
T}f_{\max}(\hat{M}(t)).$ (3.8)

The existence of a unique weak solution $\hat{u}$ directly follows from the
Lax-Milgram Lemma. If $d=1$ (one space-dimension), then we immediately have
$\hat{u}\in L^{\infty}(\Omega)$ from Morrey’s inequality [11, Chapter 5].
Hence, let $d\geq 2$. Set $q=2d/(d-2)$ for $d>2$, and $q>2$ for $d=2$. Then
for $m\geq 0$, inserting the test function $\varphi=[\hat{u}-m]_{+}\in{\cal
H}^{1}$ and denoting $A(m):=\\{\bm{x}\in\Omega:\hat{u}(\bm{x})>m\\}$ we have
the estimates,

$\displaystyle\|\nabla[\hat{u}-m]_{+}\|^{2}$
$\displaystyle\leq\hat{C}\,\|[\hat{u}-m]_{+}\|_{L^{1}(\Omega)},$
$\displaystyle\|[\hat{u}-m]_{+}\|_{L^{1}(\Omega)}$
$\displaystyle\leq|A(m)|^{1-\frac{1}{q}}\|[\hat{u}-m]_{+}\|_{L^{q}(\Omega)}\leq
C^{\prime}|A(m)|^{1-\frac{1}{q}}\|\nabla[\hat{u}-m]_{+}\|,$

where the last inequality follows from the Sobolev inequality [11, Chapter 5].
Hence, we have
$\|[\hat{u}-m]_{+}\|_{L^{1}(\Omega)}\leq{C^{\prime}}^{2}\hat{C}|A(m)|^{\gamma}$,
where $\gamma=2-\frac{2}{q}>1$. Thus, following the steps of [13, Lemma 7.3]
we conclude that $\hat{u}\in L^{\infty}(\Omega)$. Hence, using the comparison
principle [21], one has
$\Phi_{\varepsilon}(M_{s,\varepsilon}(t))\leq\hat{u}+\bar{M}<\infty$ a.e. in
$\Omega$ for all $\varepsilon>0$ and $t>0$, which concludes the proof. ∎

######  Remark 3.1 (Generalization to $\vec{s}\in(L^{2}(Q))^{k}$).

Although Lemmas 3.1 and 3.2, and the following Lemma 3.3, assume
$\vec{s}\in{\bf{\cal Z}}$ to simplify the presentation, the results remain
valid for all $\vec{s}\in(L^{2}(Q))^{k}$ as evident from the a-priori
estimates (3.4)–(3.1). This observation will become important in Lemma 3.5
which provides the setting for the proof of Theorem 3.1.

###  3.2 Existence of solutions of the degenerate parabolic problem

######  Lemma 3.3 (Existence for the degenerate problem).

Let (P1)–(P6) and (P4enumi) hold. Let $\vec{s}\in{\bf{\cal Z}}$ be given and
let $0<T^{*}\leq\infty$ denote the time, independent of $\vec{s}$ and
$\varepsilon>0$, such that the solutions $M_{s,\varepsilon}\in{\bf{\cal
X}}+h^{e}_{0}$ of (3.3) remain bounded in $[0,1)$ for all $t<T^{*}$. Let
$T<T^{*}$. Then there exists a unique $M_{s}\in{\bf{\cal W}}$ with
$\Phi(M_{s})\in L^{2}(0,T;H^{1}(\Omega))$ satisfying $M_{s}(0)=M_{0}$,
$\Phi(M_{s})=\Phi(h_{0})$ on $\Gamma_{1}$ in the trace sense, and

$\displaystyle\int_{0}^{T}\langle\varphi,\partial_{t}M_{s}\rangle+\int_{0}^{T}(\nabla\Phi(M_{s}),\nabla\varphi)=\int_{0}^{T}(f_{0}(M_{s},\vec{s}),\varphi),$
(3.9)

for all $\varphi\in L^{2}(0,T;{\cal H}^{1})$. Moreover, $0\leq M_{s}<1-\delta$
a.e. in $Q$ for some constant $\delta\in(0,1)$.

###### Proof.

In this proof, we will first assume that

$\displaystyle M_{0}\in H^{1}(\Omega).$ (3.10)

This constraint will later be dropped. Lemma 3.2 and the assumption $T<T^{*}$,
imply the existence of $\bar{\delta}\in(0,1)$ such that

$\displaystyle M_{s,\varepsilon}\in[0,1-\bar{\delta}],\ \ \text{and
consequently, }\ \
\phi_{\varepsilon}:=\Phi_{\varepsilon}(M_{s,\varepsilon})\in[0,\Phi(1-\bar{\delta})]\
\text{ a.e. in }Q,$ (3.11)

for small $\varepsilon>0$. The shorthand $\phi_{\varepsilon}$ will be used to
denote $\Phi_{\varepsilon}(M_{s,\varepsilon})$ for the rest of the proof.

(Step 1) Convergence of $\mathbf{M_{s,\varepsilon}}$ and
$\mathbf{\phi_{\varepsilon}}$, assuming (3.10): Taking (3.11) into account
which implies that ${\Phi_{\varepsilon}}^{\prime}(M_{s,\varepsilon})$ is
bounded above independent of $\varepsilon$, we conclude from (3.1) in Lemma
3.1 that $\phi_{\varepsilon}$ is uniformly bounded in ${\bf{\cal Y}}$ (see
(2.3)). Observe that ${\bf{\cal Y}}\subset H^{1}(Q)$. Using the compact
embedding $H^{1}(Q)\hookrightarrow\hookrightarrow L^{2}(Q)$, we conclude that
there exists $\phi\in H^{1}(Q)$ and a subsequence $\phi_{\varepsilon}$ such
that for $\varepsilon\to 0$,

$\displaystyle\phi_{\varepsilon}\rightharpoonup\phi\text{ weakly in
}H^{1}(Q),$ (3.12a) $\displaystyle\phi_{\varepsilon}\to\phi\text{ strongly in
}L^{2}(Q).$ (3.12b)

Setting

$\displaystyle M_{s}:=\Phi^{-1}(\phi)\in L^{\infty}(0,T;L^{\infty}(\Omega))$
(3.13)

we claim that for $\varepsilon\to 0$,

$\displaystyle\partial_{t}M_{s,\varepsilon}\rightharpoonup\partial_{t}M_{s}\text{
weakly in }L^{2}(0,T;{\cal H}^{-1}(\Omega)),$ (3.14a) $\displaystyle
M_{s,\varepsilon}\to M_{s}\text{ strongly in }L^{2}(Q).$ (3.14b)

To see this, we consider a convex strictly increasing function $\eta\in
C^{1}([0,1))$ such that

$\displaystyle\eta=\Phi\text{ in }[0,\epsilon_{0}),\quad\text{ and
}\quad\eta\circ\Phi^{-1}\in\mathrm{Lip}({\mathbb{R}}^{+}),$ (3.15)

where $\epsilon_{0}\in(0,1)$ was fixed in (P1). Recall that $D_{0}$ is
strictly increasing in $[0,\epsilon_{0})$, and therefore, $\Phi$ is convex in
$[0,\epsilon_{0})$. For $\varrho>\epsilon_{0}$,
$\Phi^{\prime}(\varrho)=D_{0}(\varrho)$ is bounded away from 0, implying that
$(\Phi^{-1})^{\prime}(\varrho)$ is bounded for $\varrho>\Phi(\epsilon_{0})$.
Hence, it is always possible to find such a function $\eta$.

Using (3.11) and that $\eta$ is strictly increasing and convex, we obtain

$\displaystyle\eta\left(\frac{1}{|Q|}\int_{Q}|M_{s,\varepsilon}-M_{s}|\right)\overset{\eqref{eq:ConvexFuncIneq1}}{\leq}\frac{1}{|Q|}\int_{Q}\eta(|M_{s,\varepsilon}-M_{s}|)\overset{\eqref{eq:ConvexFuncIneq2}}{\leq}\frac{1}{|Q|}\int_{Q}|\eta(M_{s,\varepsilon})-\eta(M_{s})|$
$\displaystyle\quad\overset{\eqref{eq:cond_eta}}{\lesssim}$
$\displaystyle\|\Phi(M_{s,\varepsilon})-\Phi(M_{s})\|_{L^{1}(Q)}\leq\|\Phi(M_{s,\varepsilon})-\phi_{\varepsilon}\|_{L^{1}(Q)}+\|\phi_{\varepsilon}-\Phi(M_{s})\|_{L^{1}(Q)}$
$\displaystyle\quad\leq\;\,$
$\displaystyle\|\Phi(M_{s,\varepsilon})-\Phi_{\varepsilon}(M_{s,\varepsilon})\|_{L^{1}(Q)}+\|\phi_{\varepsilon}-\phi\|_{L^{1}(Q)}\overset{\eqref{eq:propPhieps},\eqref{eq:convergencePhiL2}}{\longrightarrow}0,\quad\text{
for }\varepsilon\searrow 0.$

More specifically, the term
$\|\Phi(M_{s,\varepsilon})-\Phi_{\varepsilon}(M_{s,\varepsilon})\|_{L^{1}(Q)}$
vanishes due to (3.2) and (3.11) since for small $\varepsilon$ we observe that

$\displaystyle|\Phi(M_{s,\varepsilon})-\Phi_{\varepsilon}(M_{s,\varepsilon})|\overset{\eqref{eq:PhiepsDef},\eqref{eq:MPhiLess1}}{=}\left|\int_{0}^{M_{s,\varepsilon}}(D(\varrho)-\max\\{\varepsilon,D(\varrho)\\})\,\mathrm{d}\varrho\right|\leq
C\varepsilon.$

Then, the continuity of $\eta$ implies that $M_{s,\varepsilon}\to M_{s}$ in
$L^{1}(Q)$. Using (3.11) and (3.13), the convergence also holds in $L^{2}(Q)$.

Since $M_{s,\varepsilon}-h^{e}_{0}\in{\bf{\cal X}}$ (see (2.3)), (3.11) and
$(M_{s,\varepsilon}-h^{e}_{0})\in C([0,T];L^{2}(\Omega))$ from Lemma 3.3 also
imply that $M_{s,\varepsilon}\in{\bf{\cal W}}$. It follows from estimate (3.4)
in Lemma 3.1 and (3.11) that $M_{s,\varepsilon}$ is bounded in ${\bf{\cal
W}}$, and the bound is independent of $\varepsilon$. Hence,
$M_{s,\varepsilon}$ has a weak limit in $H^{1}(0,T;{\cal H}^{-1})\cap
L^{2}(Q)\supset{\bf{\cal W}}$. The strong convergence $M_{s,\varepsilon}\to
M_{s}$ in $L^{2}(Q)$ then implies by the uniqueness of weak limits that (3.14)
holds.

Next we show that $M_{s}\in{\bf{\cal W}}$, i.e. it remains to show that
$M_{s}\in C([0,T];L^{2}(\Omega))$. From the uniform boundedness of
$M_{s,\varepsilon}$ in ${\bf{\cal W}}$ we conclude that
$M_{s,\varepsilon}(t)\rightharpoonup M_{s}(t)$ weakly in ${\cal H}^{-1}$ for
almost all $t\in[0,T]$ since
$(M_{s,\varepsilon}(t)-M_{s}(t),\tilde{\varphi})=\int_{0}^{t}\langle\partial_{t}(M_{s,\varepsilon}-M_{s}),\tilde{\varphi}\rangle\to
0$ for all $\tilde{\varphi}\in{\cal H}^{1}$, see (3.14a). Hence, for all
$\varphi\in L^{2}(\Omega)$ one has

$\displaystyle|(M_{s,\varepsilon}(t)-M_{s}(t),\varphi)|$
$\displaystyle=|(M_{s,\varepsilon}(t)-M_{s}(t),\varphi-\tilde{\varphi})+(M_{s,\varepsilon}(t)-M_{s}(t),\tilde{\varphi})|$
$\displaystyle\leq(\|M_{s,\varepsilon}(t)\|+\|M_{s}(t)\|)\|\varphi-\tilde{\varphi}\|+|(M_{s,\varepsilon}(t)-M_{s}(t),\tilde{\varphi})|.$
(3.16)

Let $\mu>0$ be arbitrary. Since $M_{s,\varepsilon}(t)$ and $M_{s}(t)$ are
uniformly bounded in $L^{\infty}(\Omega)$ by Lemma 3.2 and (3.13), one can
choose $\tilde{\varphi}\in{\cal H}^{1}$ such that the first term on the right
hand side in (3.2) is less than $\mu/2$. Then, we choose $\varepsilon$ small
enough such that the second term is less than $\mu/2$ implying that

$\displaystyle M_{s,\varepsilon}(t)\rightharpoonup M_{s}(t)\text{ weakly in
}L^{2}(\Omega)\text{ for almost all }t\in[0,T].$ (3.17)

Finally to show the continuity of $M_{s}$ in time, we observe that for
$\tau>0$ one has

$\displaystyle\|M_{s}(t+\tau)-M_{s}(t)\|^{2}=(M_{s}(t+\tau)-M_{s,\varepsilon}(t+\tau),M_{s}(t+\tau)-M_{s}(t))$
$\displaystyle\quad+(M_{s,\varepsilon}(t+\tau)-M_{s,\varepsilon}(t),M_{s}(t+\tau)-M_{s}(t))+(M_{s,\varepsilon}(t)-M_{s}(t),M_{s}(t+\tau)-M_{s}(t)).$

The first and last term on the right side are arbitrarily small for all
sufficiently small $\varepsilon$ due to the weak convergence of
$M_{s,\varepsilon}(t)$ in $L^{2}(\Omega)$. The term in the middle vanishes as
$\tau$ tends to zero since $M_{s,\varepsilon}\in C([0,T];L^{2}(\Omega))$. This
proves that $M_{s}\in C([0,T];L^{2}(\Omega))$.

Using (3.12),(3.14) we can now pass to the limit in (3.3) and conclude that
$M_{s}$ is a solution of (3.9). This completes the proof for the case
$M_{0}\in H^{1}(\Omega)$.

(Step 2) Existence for $\mathbf{M_{0}\in L^{1}(\Omega)}$ satisfying (P5):  We
postulate that for a given $M_{0}\in L^{1}(\Omega)$ and $\mu>0$, there exists
$M_{0}^{\mu}\in H^{1}(\Omega)$ such that

$\displaystyle 0\leq M_{0}^{\mu}\leq\bar{M}<1\text{ a.e. in }Q,\text{ and
}\|M_{0}^{\mu}-M_{0}\|_{L^{1}(\Omega)}\leq\mu.$ (3.18)

The existence of $\tilde{M}_{0}^{\mu}\in C^{\infty}_{c}({\mathbb{R}}^{d})$
such that $\|\tilde{M}_{0}^{\mu}-M_{0}\|_{L^{1}(\Omega)}\leq\tfrac{1}{2}\mu$
follows from the fact that $C^{\infty}_{c}({\mathbb{R}}^{d})$ is dense in
$L^{1}({\mathbb{R}}^{d})$, see Theorem 4.3 in [4]. Define
$M_{0}^{\mu}=\min\\{\overline{M},\tilde{M}_{0}^{\mu}\\}\in H^{1}(\Omega)$,
where $\overline{M}$ was defined in (P5). Then $M_{0}^{\mu}$ satisfies (3.18)
since

$\displaystyle\|M_{0}^{\mu}-M_{0}\|_{L^{1}(\Omega)}$
$\displaystyle\leq\|M_{0}^{\mu}-\tilde{M}^{\mu}_{0}\|_{L^{1}(\Omega)}+\|\tilde{M}_{0}^{\mu}-M_{0}\|_{L^{1}(\Omega)}$
$\displaystyle=\|[\tilde{M}_{0}^{\mu}-\bar{M}]_{+}\|_{L^{1}(\Omega)}+\|\tilde{M}_{0}^{\mu}-M_{0}\|_{L^{1}(\Omega)}$
$\displaystyle\overset{\ref{prop:IC}}{\leq}\|\tilde{M}_{0}^{\mu}-M_{0}\|_{L^{1}(\Omega)}+\|\tilde{M}_{0}^{\mu}-M_{0}\|_{L^{1}(\Omega)}\leq\mu.$

Now, let $M_{s}^{\mu}\in{\bf{\cal W}}$ be the weak solution corresponding to
the initial data $M^{\mu}_{s}(0)=M^{\mu}_{0}\in H^{1}(\Omega)$, $\mu>0$, which
exists by Step 1. Consider a sequence
$\\{\mu_{n}\\}_{n\in{\mathbb{N}}}\subset{\mathbb{R}}^{+}$ converging to zero.
Then, the $L^{1}$-contraction result in [21] implies that there exists a
constant $C>0$ independent of $\mu$ such that

$\displaystyle\|(M_{s}^{\mu_{m}}-M_{s}^{\mu_{n}})(t)\|_{L^{1}(\Omega)}\leq
C\|M_{0}^{\mu_{m}}-M_{0}^{\mu_{n}}\|_{L^{1}(\Omega)}$ (3.19)

for all $m,n\in{\mathbb{N}}$, $t\in[0,T]$. Note that an $L^{1}$-contraction
result also holds for homogeneous Neumann boundary conditions, see [3]. Hence,
$\\{M_{s}^{\mu_{n}}\\}_{n\in{\mathbb{N}}}$ is a Cauchy sequence in
$L^{1}(\Omega)$, and since $M_{s}^{\mu_{n}}\in L^{\infty}(\Omega)$ is
uniformly bounded with respect to $\mu_{n}$ (see Lemma 3.2), it is also a
Cauchy sequence in $L^{2}(\Omega)$. Since $M_{s}^{\mu_{n}}\in
C([0,T];L^{2}(\Omega))$, we conclude that there exists $M_{s}\in
C([0,T];L^{2}(\Omega))$ such that

$\displaystyle\|M^{\mu_{n}}_{s}-M_{s}\|_{C([0,T];L^{2}(\Omega))}\to
0\quad\text{ as }n\to\infty.$

The uniform boundedness of $M^{\mu_{n}}_{s}$ in ${\bf{\cal W}}$ and
$\Phi(M^{\mu_{n}}_{s})$ in $L^{2}(0,T;H^{1}(\Omega))$ follow directly from
(3.4), Lemma 3.1. The strong convergence of $M^{\mu_{n}}_{s}\to M_{s}$ and its
uniform $L^{\infty}$-boundedness away from 1 (the singular point of $\Phi$),
implies that $\Phi(M^{\mu_{n}}_{s})$ also converges strongly to $\Phi(M_{s})$
in $L^{2}(\Omega)$ for all $t\in[0,T]$. Hence, similar to before, passing the
limit $n\to\infty$ it follows that $M_{s}$ solves (3.9). ∎

###  3.3 A contraction argument for proving Theorem 3.2

We first show the existence of a unique weak solution under the additional
assumptions stated in Theorem 3.2 compared to Theorem 3.1. The assumptions in
Theorem 3.2 demand that $D_{j}$ depends only on $S_{j}$, i.e.
$D_{j}:{\mathbb{R}}\to[D_{\min},D_{\max}]$. This is unless $\nu_{j}=0$ in
which case we can also define $D_{j}$ as such. Hence, similar to (2.8), we
introduce the function

$\displaystyle\Phi_{j}(S):=\int_{0}^{S}D_{j}(\varrho)\,\mathrm{d}\varrho,\qquad\text{for
}j\in\\{1,\dots,k\\}.$ (3.20)

Observe that due to (P2), $\Phi_{j}$ is Lipschitz continuous and strictly
increasing.

For a given $\vec{s}\in{\bf{\cal Z}}$ let $M_{s}\in{\bf{\cal W}}$ be the
corresponding solution in Lemma 3.3. Define the operator
$\mathfrak{A}:{\bf{\cal Z}}\to{\bf{\cal Z}}$ such that for all
$j\in\\{1,\dots,k\\}$, $\mathfrak{A}(\vec{s})_{j}$ satisfies
$\nu_{j}\,\mathfrak{A}(\vec{s})_{j}\in L^{2}(0,T;H^{1}(\Omega))$,
$\mathfrak{A}(\vec{s})_{j}\in H^{1}(0,T;H^{-1}(\Omega))$, and for all
$\zeta_{j}\in L^{2}(0,T;H^{1}_{0}(\Omega))$,

$\displaystyle\int_{0}^{T}\left[\langle\zeta_{j},\partial_{t}\mathfrak{A}(\vec{s})_{j}\rangle_{H^{1}_{0},H^{-1}}+\nu_{j}\left(\nabla\Phi_{j}(\mathfrak{A}(\vec{s})_{j})+\bm{v}_{j}\,\mathfrak{A}(\vec{s})_{j},\nabla\zeta_{j}\right)\right]=\int_{0}^{T}(f_{j}(M_{s},\vec{s}),\zeta_{j}),$
(3.21a) $\displaystyle\text{with }\mathfrak{A}(\vec{s})_{j}(0)=S_{0,j}\text{
and }\nu_{j}\mathfrak{A}(s)_{j}=\nu_{j}h_{j}\text{ on $\partial\Omega$ in the
trace sense}.$ (3.21b)

To prove Theorem 3.2 we need the following lemma.

######  Lemma 3.4 ($L^{1}$-contraction property of $\mathfrak{A}$).

Under the assumptions of Theorem 3.2, define
$\Phi_{j}:{\mathbb{R}}\to{\mathbb{R}}$ by (3.20). Assume that $T<T^{*}$ for
$T^{*}>0$ introduced in Lemma 3.3. Then the operator $\mathfrak{A}:{\bf{\cal
Z}}\to{\bf{\cal Z}}$, introduced in (3.21), is well-defined. Moreover, there
exists a strictly increasing function $\mathfrak{C}\in
C^{1}({\mathbb{R}}^{+})$ with $\mathfrak{C}(0)=0$ such that for all
$t\in[0,T]$ and $\vec{s}_{1},\,\vec{s}_{2}\in{\bf{\cal Z}}$,

$\int_{0}^{t}\|\mathfrak{A}(\vec{s}_{1})-\mathfrak{A}(\vec{s}_{2})\|_{(L^{1}(\Omega))^{k}}\leq\mathfrak{C}(t)\int_{0}^{t}\|\vec{s}_{1}-\vec{s}_{2}\|_{(L^{1}(\Omega))^{k}}.$

###### Proof.

Since $f_{j}(M_{s},\vec{s})\in C([0,T];L^{2}(\Omega))$, and $D_{j}$ is bounded
from above and below by a positive constant by (P2), the existence and
regularity results in [1] imply that $\mathfrak{A}(\vec{s})_{j}$ is well-
defined for $\nu_{j}>0$ (similar to Lemma 3.1). If $\nu_{j}=0$, then
$\mathfrak{A}(\vec{s})_{j}$ is simply the solution of an ODE with known right
hand side. From (3.9), using the $L^{1}$-contraction result in [21, 3] and the
Lipschitz continuity of $f_{0}$, it follows that for all $t\in[0,T]$,

$\displaystyle\|(M_{s_{1}}-M_{s_{2}})(t)\|_{L^{1}(\Omega)}$
$\displaystyle\leq\int_{0}^{t}\|f_{0}(M_{s_{1}},\vec{s}_{1})-f_{0}(M_{s_{2}},\vec{s}_{2})\|_{L^{1}(\Omega)}$
$\displaystyle\overset{\ref{prop:fg1}}{\leq}C_{L}\int_{0}^{t}\|\vec{s}_{1}-\vec{s}_{2}\|_{(L^{1}(\Omega))^{k}}+C_{L}\int_{0}^{t}\|M_{s_{1}}-M_{s_{2}}\|_{L^{1}(\Omega)}.$
(3.22)

Applying Gronwall’s Lemma (2.5a) we conclude that

$\displaystyle\|(M_{s_{1}}-M_{s_{2}})(t)\|_{L^{1}(\Omega)}\leq
C_{L}\exp(C_{L}\,t)\int_{0}^{t}\|\vec{s}_{1}-\vec{s}_{2}\|_{(L^{1}(\Omega))^{k}}.$
(3.23)

We now apply the $L^{1}$-contraction principle to (3.21) and use the Lischitz
continuity of $f_{j}$ and the previous estimate to get

$\displaystyle\|(\mathfrak{A}(\vec{s}_{1})_{j}-\mathfrak{A}(\vec{s}_{2})_{j})(t)\|_{L^{1}(\Omega)}$
$\displaystyle\leq\int_{0}^{t}\|f_{j}(M_{s_{1}},\vec{s}_{1})-f_{j}(M_{s_{2}},\vec{s}_{2})\|_{L^{1}(\Omega)}$
$\displaystyle\leq
C_{L}\int_{0}^{t}\|\vec{s}_{1}-\vec{s}_{2}\|_{(L^{1}(\Omega))^{k}}+C_{L}\int_{0}^{t}\|M_{s_{1}}-M_{s_{2}}\|_{L^{1}(\Omega)}$
$\displaystyle\leq
C_{L}(1+C_{L}\,t\,\exp(C_{L}\,t))\int_{0}^{t}\|\vec{s}_{1}-\vec{s}_{2}\|_{(L^{1}(Q))^{k}}.$
(3.24)

Note that this estimate also holds for the case $\nu_{j}=0$. Hence, setting
$\mathfrak{C}(t)=k\,C_{L}\,t\,(1+C_{L}\,t\,\exp(C_{L}\,t))$ the result
follows. ∎

###### Proof of Theorem 3.2.

Choosing $T>0$ small enough such that $\mathfrak{C}(T)<1$ the existence of a
unique weak solution $(M,\vec{S})$ of (1.1) follows from Lemma 3.4 and
Banach’s fixed point theorem. Since Lemma 3.2 implies that $T^{*}=\infty$
provided that $\Gamma_{1}$ has a non-zero measure, the argument can be
repeated and solutions can be patched together to cover the interval $[0,T]$
for an arbitrary $T>0$, thus concluding the proof. ∎

###  3.4 A fixed point argument for proving Theorem 3.1

In this section, we use Schauder’s fixed point theorem to prove the existence
of solutions for general diffusion coefficients $D_{j}$ satisfying (P2). Since
the case of ODE-PDE couplings, where $\nu_{j}=0$, is already covered by the
previous section, here we assume that $\nu_{j}>0$. Similarly as in the
previous section, we define the map
$\mathfrak{B}:(L^{2}(Q))^{k}\to(L^{2}(Q))^{k}$ such that for all
$j\in\\{1,\dots,k\\}$, $\mathfrak{B}(\vec{s})_{j}\in
L^{2}(0,T;H^{1}(\Omega))\cap H^{1}(0,T;H^{-1}(\Omega))$, and for all
$\zeta_{j}\in L^{2}(0,T;H^{1}_{0}(\Omega))$,

$\displaystyle\int_{0}^{T}[\langle\zeta_{j},\partial_{t}\mathfrak{B}(\vec{s})_{j}\rangle_{H^{1}_{0},H^{-1}}+\nu_{j}(D_{j}(M_{s},\vec{s})\nabla\mathfrak{B}(\vec{s})_{j}+\bm{v}_{j}\,\mathfrak{B}(\vec{s})_{j},\nabla\zeta_{j})]=\int_{0}^{T}(f_{j}(M_{s},\mathfrak{B}(\vec{s})),\zeta_{j}),$
(3.25a) $\displaystyle\text{with }\mathfrak{B}(\vec{s})_{j}(0)=S_{0,j},\text{
and }\mathfrak{B}(s)_{j}=h_{j}\text{ on $\partial\Omega$ in the trace sense}.$
(3.25b)

######  Lemma 3.5 (Schauder criteria for $\mathfrak{B}$).

Let $\nu_{j}>0$ for all $j\in\\{1,\dots,k\\}$. Then under the assumptions of
Lemma 3.3, the operator $\mathfrak{B}:(L^{2}(Q))^{k}\to(L^{2}(Q))^{k}$
introduced in (3.25) is well-defined, continuous, compact, and
$\|\mathfrak{B}(\vec{s})\|_{\bf{\cal Z}}$ is bounded for all
$\vec{s}\in(L^{2}(Q))^{k}$.

###### Proof.

(Step 1): Well-posedness, boundedness and compactness. Recalling Remark 3.1,
we have existence and boundedness of weak solutions $M_{s}\in{\bf{\cal W}}$
for any $\vec{s}\in(L^{2}(Q))^{k}$. We observe that $D_{j}$ satisfies the
ellipticity condition by (P2), $\bm{v}_{j}\in(L^{\infty}(\Omega))^{d}$,
$f_{j}(\cdot,s)$ is Lipschitz continuous, and for $h^{e}_{j}$ in (P6), we have

$|f_{j}(M_{s},h^{e}_{j})|\overset{\ref{prop:fg1}}{\leq}C(1+|h^{e}_{j}|+|M_{s}|)\in
L^{2}(\Omega)\text{ for a constant }C>0.$

Consequently, $\mathfrak{B}(\vec{s})_{j}\in L^{2}(0,T;H^{1}(\Omega))\cap
H^{1}(0,T;H^{-1}(\Omega))$ is well-defined. This is also evident from a
Schaefer’s fixed point argument [11, Chapter 9] using the a-priori estimate
obtained by inserting $\zeta_{j}=\mathfrak{B}(\vec{s})_{j}-h^{e}_{j}$ in
(3.25). For the first term this yields,

$\displaystyle\int_{0}^{T}\langle\mathfrak{B}(\vec{s})_{j}-h^{e}_{j},\partial_{t}\mathfrak{B}(\vec{s})_{j}\rangle_{H^{1}_{0},H^{-1}}=\frac{1}{2}\left[\|\mathfrak{B}(\vec{s}(T))_{j}-h^{e}_{j}\|^{2}-\|S_{0,j}-h^{e}_{j}\|^{2}\right],$

and for the diffusion and convection terms we obtain,

$\displaystyle\nu_{j}\int_{0}^{T}(D_{j}(M_{s},\vec{s})\nabla\mathfrak{B}(\vec{s})_{j},\nabla(\mathfrak{B}(\vec{s})_{j}-h^{e}_{j}))$
$\displaystyle\qquad\quad=\frac{\nu_{j}}{2}\int_{Q}D_{j}(M_{s},\vec{s})\left[|\nabla\mathfrak{B}(\vec{s})_{j}|^{2}-|\nabla
h^{e}_{j}|^{2}+|\nabla(\mathfrak{B}(\vec{s})_{j}-h^{e}_{j})|^{2}\right],$
$\displaystyle\left|\nu_{j}\int_{0}^{T}(\bm{v}_{j}\,\mathfrak{B}(\vec{s})_{j},\nabla(\mathfrak{B}(\vec{s})_{j}-h^{e}_{j}))\right|$
$\displaystyle\qquad\quad\overset{\eqref{Eq:YoungsIneq}}{\leq}\tfrac{\nu_{j}\|\bm{v}_{j}\|^{2}_{L^{\infty}(\Omega)}}{2D_{\min}}\int_{0}^{T}\|\mathfrak{B}(\vec{s})_{j}\|^{2}+\tfrac{\nu_{j}D_{\min}}{2}\int_{0}^{T}\|\nabla(\mathfrak{B}(\vec{s})_{j}-h^{e}_{j})\|^{2}.$

Finally, we can estimate the reaction term by

$\displaystyle\int_{0}^{T}(f_{j}(M_{s},\mathfrak{B}(\vec{s})),\mathfrak{B}(\vec{s})_{j}-h^{e}_{j})=\int_{0}^{T}(f_{j}(M_{s},\mathfrak{B}(\vec{s}))-f_{j}(M_{s},\vec{h}^{e})+f_{j}(M_{s},\vec{h}^{e}),\mathfrak{B}(\vec{s})_{j}-h^{e}_{j})$
$\displaystyle\qquad\overset{\eqref{Eq:YoungsIneq},\ref{prop:fg1}}{\leq}C\sum_{i=1}^{k}\int_{0}^{T}\|\mathfrak{B}(\vec{s})_{i}-h^{e}_{i}\|^{2}_{(L^{2}(\Omega))^{k}}+\int_{0}^{T}\|f_{j}(M_{s},\vec{h}^{e})\|^{2}.$

Combining the above estimates, using Young’s inequality and summing from $j=1$
to $k$, one has

$\displaystyle\sum_{j=1}^{k}\left[\|\mathfrak{B}(\vec{s}(T))_{j}-h^{e}_{j}\|^{2}+\nu_{j}D_{\min}\int_{0}^{T}\|\nabla\mathfrak{B}(\vec{s})_{j}\|^{2}\right]\lesssim
1+\int_{0}^{T}\sum_{j=1}^{k}\|\mathfrak{B}(\vec{s})_{j}-h^{e}_{j}\|^{2}.$

Thus, using Gronwall’s Lemma (2.5a), we conclude that $\mathfrak{B}(\vec{s})$
is uniformly bounded in $L^{2}(0,T;(H^{1}(\Omega))^{k})$ and
$L^{\infty}(0,T;(L^{2}(\Omega))^{k})$. Moreover, from (3.25), we obtain

$\displaystyle\|\partial_{t}\mathfrak{B}(\vec{s})_{j}\|_{L^{2}(0,T;H^{-1}(\Omega))}$
$\displaystyle=\sup_{\begin{subarray}{c}\|\zeta_{j}\|_{L^{2}(0,T;H^{1}_{0}(\Omega))}=1\end{subarray}}\int_{0}^{T}\left[\nu_{j}(D_{j}(M_{s},\vec{s})\nabla\mathfrak{B}(\vec{s})_{j}+\bm{v}_{j}\,\mathfrak{B}(\vec{s})_{j},\nabla\zeta_{j})-(f_{j}(M_{s},\mathfrak{B}\vec{s}),\zeta_{j})\right]$
$\displaystyle\leq\nu_{j}(D_{\max}\|\nabla\mathfrak{B}(\vec{s})_{j}\|_{L^{2}(Q)}+\|\bm{v}_{j}\|_{L^{\infty}(\Omega)}\|\mathfrak{B}(\vec{s})_{j}\|_{L^{2}(Q)})+C_{\Omega}\,\|f_{j}(M_{s},\mathfrak{B}(\vec{s}))\|_{L^{2}(Q)}.$

Hence $\mathfrak{B}(\vec{s})_{j}-h^{e}_{j}\in L^{2}(0,T;H^{1}_{0}(\Omega))$
and $\mathfrak{B}(\vec{s})_{j}\in H^{1}(0,T;H^{-1}(\Omega))$ is bounded
uniformly with respect to $\vec{s}\in(L^{2}(Q))^{k}$. Due to the compact
embedding $L^{2}(0,T;H^{1}_{0}(\Omega))\cap
H^{1}(0,T;H^{-1}(\Omega))\hookrightarrow\hookrightarrow L^{2}(Q)$ (see [23]),
the mapping $\mathfrak{B}$ is also compact. Moreover, due to the continuous
embedding $L^{2}(0,T;H^{1}_{0}(\Omega))\cap H^{1}(0,T;H^{-1}(\Omega))\subset
C([0,T];L^{2}(\Omega))$, $\|\mathfrak{B}(\vec{s})\|_{\bf{\cal Z}}$ is bounded
uniformly.

(Step 2): Continuity. Let the sequence $\\{\vec{s}^{i}\\}_{i\in{\mathbb{N}}}$
converge to $\vec{s}^{\,\star}$ in $(L^{2}(Q))^{k}$. Let $M_{s}^{i}$ and
$M^{\star}_{s}$ denote the solutions of (3.9) for $\vec{s}=\vec{s}^{i}$ and
$\vec{s}=\vec{s}^{\,\star}$ respectively. Then, by using the
$L^{1}$-contraction result (3.23) we have that
$\|(M_{s}^{i}-M_{s}^{\star})(t)\|_{L^{1}(\Omega)}\to 0$ for all $t\in(0,T]$.
Since, $M_{s}^{i}$ are bounded in $L^{\infty}(\Omega)$ (Lemma 3.3), one
further has that

$\displaystyle\|\vec{s}^{i}-\vec{s}^{\,\star}\|_{(L^{2}(Q))^{k}}+\|M_{s}^{i}-M_{s}^{\star}\|_{C([0,T];L^{2}(\Omega))}\rightarrow
0\ \text{ as }i\rightarrow\infty.$ (3.26)

Observe that, since $\mathfrak{B}(\vec{s}^{\,\star})_{j}\in
L^{2}(0,T;H^{1}(\Omega))$, for any given $\varepsilon>0$, there exists
$s_{j}^{\varepsilon,\star}\in C^{\infty}({\mathbb{R}}^{d})$ such that

$\displaystyle\|\mathfrak{B}(\vec{s}^{\,\star})_{j}-s_{j}^{\varepsilon,\star}\|_{L^{2}(0,T;H^{1}(\Omega))}\leq\varepsilon/4D_{\max}.$
(3.27)

We consider the difference of $\mathfrak{B}(\vec{s}^{i})_{j}$ and
$\mathfrak{B}(\vec{s}^{\,\star})_{j}$ by subtracting two versions of (3.25).
First, we split up the diffusion term,

$\displaystyle\int_{0}^{T}(D_{j}(M_{s}^{i},\vec{s}^{i})\nabla\mathfrak{B}(\vec{s}^{i})_{j}-D_{j}(M_{s}^{\star},\vec{s}^{\,\star})\nabla\mathfrak{B}(\vec{s}^{\,\star})_{j},\nabla\zeta_{j})$
$\displaystyle=\int_{0}^{T}(D_{j}(M_{s}^{i},\vec{s}^{i})\nabla(\mathfrak{B}(\vec{s}^{i})_{j}-\mathfrak{B}(\vec{s}^{\,\star})_{j})+(D_{j}(M_{s}^{i},\vec{s}^{i})-D_{j}(M_{s}^{\star},\vec{s}^{\,\star}))\nabla\mathfrak{B}(\vec{s}^{\,\star})_{j},\nabla\zeta_{j}),$
(3.28a) and use Hölder’s inequality to estimate the second term as follows
$\displaystyle\|(D_{j}(M_{s}^{i},\vec{s}^{i})-D_{j}(M_{s}^{\star},\vec{s}^{\,\star}))\nabla\mathfrak{B}(\vec{s}^{\,\star})_{j}\|_{L^{2}(Q)}$
$\displaystyle\leq\|(D_{j}(M_{s}^{i},\vec{s}^{i})-D_{j}(M_{s}^{\star},\vec{s}^{\,\star}))\nabla
s^{\varepsilon,\star}_{j}\|_{L^{2}(Q)}+\|(D_{j}(M_{s}^{i},\vec{s}^{i})-D_{j}(M_{s}^{\star},\vec{s}^{\,\star}))\nabla(\mathfrak{B}(\vec{s}^{\,\star})_{j}-s^{\varepsilon,\star}_{j})\|_{L^{2}(Q)}$
$\displaystyle\leq\|s^{\varepsilon,\star}_{j}\|_{C^{1}(Q)}\|D_{j}(M_{s}^{i},\vec{s}^{i})-D_{j}(M_{s}^{\star},\vec{s}^{\,\star})\|_{L^{2}(Q)}+2D_{\max}\|\nabla(\mathfrak{B}(\vec{s}^{\,\star})_{j}-s^{\varepsilon,\star}_{j}))\|_{L^{2}(Q)}\overset{\eqref{eq:MSto0inL2},\eqref{eq:Sstar_eps_def}}{\leq}\varepsilon,$
(3.28b) for all $i\geq i_{\varepsilon,1}$, where
$i_{\varepsilon,1}\in{\mathbb{N}}$ is large enough. Here we used the Lipschitz
continuity of $D_{j}$, see (P2). Furthermore (3.26) implies for $i\geq
i_{\varepsilon,2}$, where $i_{\varepsilon,2}\in{\mathbb{N}}$ is large enough,
that
$\displaystyle\|f_{j}(M_{s}^{i},\mathfrak{B}(\vec{s}^{i}))-f_{j}(M_{s}^{\star},\mathfrak{B}(\vec{s}^{\,\star}))\|_{L^{2}(\Omega)}\overset{\ref{prop:fg1}}{\lesssim}\varepsilon+\|\mathfrak{B}(\vec{s}^{i})-\mathfrak{B}(\vec{s}^{\,\star})\|_{(L^{2}(\Omega))^{k}}.$
(3.28c)

Hence, defining $\eth
s_{j}^{i}:=\mathfrak{B}(\vec{s}^{i})_{j}-\mathfrak{B}(\vec{s}^{\,\star})_{j}$,
it follows from (3.25) that for
$i\geq\max\\{i_{\varepsilon,1},i_{\varepsilon,2}\\}$,

$\displaystyle\left|\int_{0}^{T}[\langle\zeta_{j},\partial_{t}\,\eth
s^{i}_{j}\rangle_{H^{1}_{0},H^{-1}}+\nu_{j}(D_{j}(M_{s}^{i},\vec{s}^{i})\nabla\eth
s^{i}_{j}+\bm{v}_{j}\,\eth s^{i}_{j},\nabla\zeta_{j})]\right|$
$\displaystyle\overset{\eqref{eq:DTermsplit}}{\lesssim}\int_{0}^{T}\left[\sum_{j=1}^{k}\|\eth
s^{i}_{j}\|\|\zeta_{j}\|+\varepsilon\|\zeta_{j}\|+\varepsilon\|\nabla\zeta_{j}\|\right]$
$\displaystyle\ \,\lesssim\ \int_{0}^{T}\left[\sum_{j=1}^{k}\|\eth
s^{i}_{j}\|\|\zeta_{j}\|+\varepsilon+\varepsilon\left(\|\zeta_{j}\|^{2}+\|\nabla\zeta_{j}\|^{2}\right)\right].$

Finally, we insert the test function $\zeta_{j}=\eth s^{i}_{j}\in
L^{2}(0,T;H^{1}_{0}(\Omega))$ and sum up the resulting estimates from $j=1$ to
$k$. Then for $\varepsilon>0$ small enough, one obtains using Gronwall’s lemma
(2.5a)

$\displaystyle\sum_{j=1}^{k}\left[\|\eth
s^{i}_{j}\|^{2}+\int_{0}^{T}\|\nabla\eth
s^{i}_{j}\|^{2}\right]\lesssim\varepsilon.$

Hence, the right hand side can be made arbitrary small by choosing
$\varepsilon>0$ small enough which simply requires
$i\geq\max\\{i_{\varepsilon,1},i_{\varepsilon,2}\\}$. Passing to the limit
$\varepsilon\to 0$ we conclude that $\|\eth
s^{i}_{j}\|_{L^{2}(0,T;H^{1}(\Omega))}+\|\eth s^{i}_{j}\|_{\bf{\cal Z}}\to 0$
for all $j=1,\dots,k$. This shows that the operator $\mathfrak{B}$ is
continuous, thus, concluding the proof. ∎

###### Proof of Theorem 3.1.

If $\nu_{j}>0$ for all $j\in\\{1,\dots,k\\}$, then using Lemma 3.5 and
Schauder’s fixed point theorem, see [11, Chapter 9], we conclude that a fixed
point $\vec{s}=\vec{S}\in{\bf{\cal Z}}\subset(L^{2}(Q))^{k}$ exists of the
mapping $\mathfrak{B}$, i.e. $\mathfrak{B}(\vec{S})=\vec{S}$. It is easy to
verify that this fixed point $(M_{S},\vec{S})$ is a weak solution of (1.1).

If $\nu_{\bar{j}}=0$ for all $\bar{j}\in\mathcal{I}\subset\\{1,\dots,k\\}$,
then the theorem is proved by first applying the contraction mapping in Lemma
3.4 for (1.1a) and (1.1b) with $\bar{j}\in\mathcal{I}$, followed by applying
the fixed point argument for $j\in\\{1,\dots,k\\}\setminus\mathcal{I}$. The
details are left for the avid reader. ∎

The approach developed in this section can be extended to systems with
degenerate diffusion coefficients $D_{j}$ under some additional assumptions.
Below we discuss an example of such a case.

######  Corollary 3.2.1 (Existence of weak solutions for degenerate $D_{j}$).

Assume that (P1), (P3)–(P6) and (P4enumi) hold. For some
$\ell\in\\{1,\dots,k\\}$, instead of (P2), assume that the diffusion
coefficient $D_{\ell}$ satisfies $D_{\ell}(m,\vec{s})=D_{\ell}(s_{\ell})$
where

$\displaystyle D_{\ell}:[0,\infty)\to{\mathbb{R}}\text{ is continuous and
strictly increasing, with }D_{\ell}(0)=0.$

Moreover, let $\mathrm{ess}\inf\\{S_{0,\ell}\\}\geq
0,\mathrm{ess}\inf\\{h_{\ell}\\}\geq 0$, $\bm{v}_{\ell}=\bm{0}$ in $Q$,
$f_{\ell}(\cdot,\vec{s})\leq f^{\ell}_{\max}(s_{j})$ for some function
$f^{\ell}_{\max}\in\mathrm{Lip}({\mathbb{R}}^{+})$ and
$f_{\ell}(m,\vec{s})\geq 0$ if $s_{\ell}=0$.

Then a weak solution $(M,\vec{S})$ of (1.1) exists in the sense of Equation
2.9, but with $\nu_{\ell}\,S_{\ell}\in L^{2}(0,T;H^{1}(\Omega))$ replaced by
$\nu_{\ell}\int_{0}^{S_{\ell}}D_{\ell}\in L^{2}(0,T;H^{1}(\Omega))$. The
solution is unique if for all $j\in\\{1,\dots,k\\}$ either $\nu_{j}=0$ or
$D_{j}$ depends only on $s_{j}$. Furthermore, $S_{\ell}$ is non-negative and
bounded almost everywhere in $Q$.

If $D_{\ell}=D_{\ell}(s_{\ell})$ is degenerate for some
$\ell\in\\{1,\dots,k\\}$, without loss of generality we assume $\nu_{\ell}>0$.
We define $\Phi_{\ell}$ as in (3.20) and fix all $s_{j}\in
C([0,T];L^{2}(\Omega))$, $j\not=\ell$, for $s_{\ell}\in
C([0,T];L^{2}(\Omega))$. Let $M_{s}\in{\bf{\cal W}}$ be the solution of (3.9)
from Lemma 3.3. Then, $\tilde{s}_{\ell}\in H^{1}(0,T;H^{-1}(\Omega))\cap
C([0,T];L^{2}(\Omega))$ with $\Phi_{\ell}(\tilde{s}_{\ell})\in
L^{2}(0,T;H^{1}(\Omega))$ is defined as the solution of the following problem,
for all $\zeta_{\ell}\in L^{2}(0,T;H^{1}_{0}(\Omega))$,

$\displaystyle\int_{0}^{T}[\langle\zeta_{\ell},\partial_{t}\tilde{s}_{\ell}\rangle_{H^{1}_{0},H^{-1}}+\nu_{\ell}(\nabla\Phi_{\ell}(\tilde{s}_{\ell}),\nabla\zeta_{\ell})]=\int_{0}^{T}(f_{\ell}(M_{s},(s_{1},\dots,\tilde{s}_{\ell},\dots,s_{k})),\zeta_{\ell}),$
$\displaystyle\text{with }\tilde{s}_{\ell}(0)=S_{0,\ell}\ \text{ and }\
\Phi_{\ell}(\tilde{s}_{\ell})=\Phi_{\ell}(h_{\ell})\text{ on $\partial\Omega$
in the trace sense}.$

The existence and uniqueness of $\tilde{s}_{\ell}$ then follows from Lemmas
3.1 and 3.3. Defining
$\hat{S}_{\ell}(t):=\bar{S}+\int_{0}^{t}f^{\ell}_{\max}$, we have, similar to
Lemma 3.2 that $0\leq\tilde{s}_{\ell}(t)\leq\hat{S}_{\ell}(t)$ a.e. in
$\Omega$ for all $t>0$. Following Lemma 3.4, we further conclude that
$\tilde{s}_{\ell}$ satisfies an $L^{1}$-contraction result with respect to
$s_{\ell}$ since all other $s_{j}$, $j\not=\ell$, are fixed. Finally, the
arguments in the proof of Theorem 3.2 concludes the proof.

##  4 Homogeneous Neumann boundary conditions

In this section, we show the existence of solutions for homogeneous Neumann
boundary conditions and present the proof of Lemma 3.1. The global existence
of solutions cannot be guaranteed for homogeneous Neumann boundary conditions,
since the density $M$ might reach 1 in finite time. The local existence and
uniqueness of solutions is analyzed in Section 4.1 and the finite time blow-up
in Section 4.2.

###  4.1 Existence of weak solutions

######  Theorem 4.1 (Local well-posedness for homogeneous Neumann conditions).

Let $\Gamma_{1}=\emptyset$. We assume that (P1)–(P6) and (P4enumi) hold. Then,
there exists a positive time $T^{*}\geq\sup\\{t:\hat{M}(t)<1\\}>0$, where
$\hat{M}\in C^{1}({\mathbb{R}}^{+})$ is defined in Lemma 3.2, such that for
$T\in(0,T^{*})$, a weak solution $(M,\vec{S})$ of (1.1) exists in the sense of
Equation 2.9. Moreover, the solution is unique if either $\nu_{j}=0$, or
$D_{j}$ depends only on $S_{j}$ for all $j\in\\{1,\dots,k\\}$.

This result essentially follows from the proof of Theorems 3.2 and 3.1.
Indeed, note that Lemmas 3.2, 3.3, 3.4 and 3.5 were proven for the general
case, i.e., they also hold for homogeneous Neumann boundary conditions. Hence,
it remains to show Lemma 3.1 for the case that $\Gamma_{1}=\emptyset$. In this
subsection, we present the proof under more general assumptions that cover
mixed as well as homogeneous Neumann boundary conditions. The proof follows
the Rothe method [16] that is based on time-discrete approximations of the
solutions. To simplify the notation for different boundary conditions (see
(2.2c))

$\text{without loss of generality we {assume that} }h_{0}\equiv 0,\text{
implying }h_{0}^{e}\equiv 0.$

####  4.1.1 Well-posedness of backward Euler time–discretizations

We consider an equivalent formulation of (3.3) and discrtize it using the
backward Euler scheme. Following (3.1), we introduce

$\displaystyle\beta_{\varepsilon}:={\Phi_{\varepsilon}}^{-1}\quad\text{ such
that }\quad\varepsilon\leq{\beta_{\varepsilon}}^{\prime}\leq\varepsilon^{-1}.$
(4.1)

Then replacing $\Phi_{\varepsilon}(M_{s,\varepsilon})$ by $u$ and
$M_{s,\varepsilon}$ by $\beta_{\varepsilon}(u)$, we demand that $u\in
L^{2}(0,T;{\cal H}^{1})$ with $\beta_{\varepsilon}(u)\in H^{1}(0,T;{\cal
H}^{-1})$ and $\beta_{\varepsilon}(u(0))=M_{0}$ satisfies

$\displaystyle\int_{0}^{T}\langle\varphi,\partial_{t}\beta_{\varepsilon}(u)\rangle+\int_{0}^{T}(\nabla
u,\nabla\varphi)=\int_{0}^{T}(f_{0}(\beta_{\varepsilon}(u),\vec{s}),\varphi)\quad\text{forall}\
\varphi\in L^{2}(0,T;{\cal H}^{1})$ (4.2)

and a given $\vec{s}\in{\bf{\cal Z}}$. For $N\in{\mathbb{N}}$, we denote by
$\tau:=T/N$ the time-step size and set $t_{n}:=n\tau$ for
$n\in\\{0,1,\dots,N\\}$. Then we define the time–discrete sequence
$\\{u_{n}\\}_{n=1}^{N}\subset{\cal H}^{1}$ recursively as follows: setting
$u_{0}:=\Phi_{\varepsilon}(M_{0})$ (i.e., $\beta_{\varepsilon}(u_{0})=M_{0}$),
let $u_{n}\in{\cal H}^{1}$ be the solution of

$\displaystyle\tfrac{1}{\tau}(\beta_{\varepsilon}(u_{n})-\beta_{\varepsilon}(u_{n-1}),\zeta)+(\nabla
u_{n},\nabla\zeta)=(f_{0}(\beta_{\varepsilon}(u_{n}),\vec{s}(t_{n})),\zeta)\qquad\text{for
all }\zeta\in{\cal H}^{1}.$ (4.3)

The following lemma implies the well-posedness of the time–discrete
formulation.

######  Lemma 4.1 (Well-posedness a semilinear elliptic problem).

For a given $F\in L^{2}(\Omega)$, there exists a unique solution $w\in{\cal
H}^{1}$ of the elliptic problem

$\displaystyle(\beta_{\varepsilon}(w),\zeta)+(\nabla
w,\nabla\zeta)=(F,\zeta)\qquad\text{for all }\zeta\in{\cal H}^{1}.$ (4.4)

###### Proof.

The proof is based on monotonicity arguments. Let the operator
$\mathfrak{F}:{\cal H}^{1}\to{\cal H}^{-1}$ be defined by the inner product

$\displaystyle\langle\zeta,\mathfrak{F}(w)\rangle:=(\beta_{\varepsilon}(w),\zeta)+(\nabla
w,\nabla\zeta).$ (4.5)

Then $\mathfrak{F}$ is strongly monotone since

$\langle
w-v,\mathfrak{F}(w)-\mathfrak{F}(v)\rangle\geq\varepsilon\|w-v\|^{2}+\|\nabla(w-v)\|^{2}\geq\varepsilon\|w-v\|_{{\cal
H}^{1}}^{2}.$

Furthermore, $\mathfrak{F}$ is Lipschitz continuous since, using the Cauchy-
Schwarz inequality, we obtain

$\|\mathfrak{F}(w)-\mathfrak{F}(v)\|_{{\cal H}^{-1}}\leq\sup_{\zeta\in{\cal
H}^{1}}\left(\frac{\varepsilon^{-1}\|w-v\|\|\zeta\|+\|\nabla(w-v)\|\|\nabla\zeta\|}{\|\zeta\|_{{\cal
H}^{1}}}\right)\leq\varepsilon^{-1}\|w-v\|_{{\cal H}^{1}}.$

Hence, invoking the nonlinear Lax-Milgram Lemma [26, Theorem 2.G] completes
the proof. ∎

Observe that the operator $\widetilde{\mathfrak{F}}:{\cal H}^{1}\to{\cal
H}^{-1}$ defined by the inner product

$\displaystyle\langle\zeta,\widetilde{\mathfrak{F}}(w)\rangle:=(\beta_{\varepsilon}(w)-\tau\,f_{0}(\beta_{\varepsilon}(w),\vec{s}(t_{n})),\zeta)+(\nabla
w,\nabla\zeta).$ (4.6)

is strictly monotone with respect to $w$ if $\tau<C_{L}^{-1}$ by (P4) and
Lipschitz continuous. Hence, adjusting the arguments in the proof of Lemma 4.1
to (4.3) we obtain the existence and uniqueness of the time–discrete
solutions.

######  Lemma 4.2 (Well-posedness of the time–discrete solutions).

Let (P1)–(P6) hold. Then the sequence $\\{u_{n}\\}_{n=1}^{N}\subset{\cal
H}^{1}$ introduced in (4.3) is well-defined for $\tau<C_{L}^{-1}$.

####  4.1.2 Interpolations in time

For a fixed $N\in{\mathbb{N}}$ with $\tau=T/N$, we define the time
interpolates $\hat{u}_{\tau}\in L^{\infty}(0,T;{\cal H}^{1})$ and
$\bar{u}_{\tau}\in C([0,T];{\cal H}^{1})$ from the time–discrete solutions
$\\{u_{n}\\}_{n=1}^{N}\subset{\cal H}^{1}$ such that for
$t\in(t_{n-1},t_{n}]$, $n\in\\{1,\dots,N\\}$,

$\displaystyle\hat{u}_{\tau}:=u_{n},\quad\text{ and
}\quad\bar{u}_{\tau}:=\beta_{\varepsilon}^{-1}\left(\beta_{\varepsilon}(u_{n-1})+\tfrac{t-t_{n-1}}{\tau}(\beta_{\varepsilon}(u_{n})-\beta_{\varepsilon}(u_{n-1}))\right).$
(4.7)

Observe that $\bar{u}_{\tau}$ satisfies for all $n\in\\{1,\dots,N\\}$,

$\displaystyle\bar{u}_{\tau}(t_{n})=u_{n},\;\;\text{ and
}\;\;\partial_{t}\beta(\bar{u}_{\tau})=\frac{\beta_{\varepsilon}(u_{n})-\beta_{\varepsilon}(u_{n-1})}{\tau}\quad\text{
for }t\in(t_{n-1},t_{n}].$ (4.8)

######  Lemma 4.3 (Uniform boundedness of the time interpolates with respect
to $\tau$).

Let (P1)–(P6) hold. Then there exist constants $\tau^{*},C>0$, independent of
$\tau$, such that for $\tau<\tau^{*}$,

$\displaystyle\|\beta_{\varepsilon}(\hat{u}_{\tau})\|^{2}+\int_{0}^{T}\|\nabla\hat{u}_{\tau}\|^{2}$
$\displaystyle\leq
C+C\int_{0}^{T}\left(\|\vec{s}\|^{2}+\|\hat{u}_{\tau}\|^{2}\right),$ (4.9a)
$\displaystyle\|\beta_{\varepsilon}(\bar{u}_{\tau})\|^{2}+\int_{0}^{T}[\|\nabla\bar{u}_{\tau}\|^{2}+\|\partial_{t}\beta_{\varepsilon}(\bar{u}_{\tau})\|_{{\cal
H}^{-1}}^{2}]$ $\displaystyle\leq
C+C\int_{0}^{T}\left(\|\vec{s}\|^{2}+\|\hat{u}_{\tau}\|^{2}\right).$ (4.9b)

The above inequalities imply the uniform boundedness of
$\beta_{\varepsilon}(\hat{u}_{\tau}),\,\beta_{\varepsilon}(\bar{u}_{\tau})\in
L^{\infty}(0,T;L^{2}(\Omega))$, $\hat{u}_{\tau},\,\bar{u}_{\tau}\in
L^{2}(0,T;{\cal H}^{1})$ and $\beta_{\varepsilon}(\bar{u}_{\tau})\in
H^{1}(0,T;{\cal H}^{-1})$ with respect to $\tau<\tau^{*}$.

###### Proof.

(Step 1) Uniform boundedness of $\|\beta_{\varepsilon}(u_{\tau})\|$: We choose
the test function $\zeta=\beta_{\varepsilon}(u_{n})\in{\cal H}^{1}$ in (4.3),
yielding

$\displaystyle(\beta_{\varepsilon}(u_{n})-\beta_{\varepsilon}(u_{n-1}),\beta_{\varepsilon}(u_{n}))+\tau(\nabla
u_{n},\nabla\beta_{\varepsilon}(u_{n}))=\tau(f_{0}(\beta_{\varepsilon}(u_{n}),\vec{s}(t_{n})),\beta_{\varepsilon}(u_{n})).$
(4.10)

Observe from the identity $2a(a-b)=a^{2}-b^{2}+(a-b)^{2}$ that

$(\beta_{\varepsilon}(u_{n})-\beta_{\varepsilon}(u_{n-1}),\beta_{\varepsilon}(u_{n}))=\tfrac{1}{2}[\|\beta_{\varepsilon}(u_{n})\|^{2}-\|\beta_{\varepsilon}(u_{n-1})\|^{2}+\|\beta_{\varepsilon}(u_{n})-\beta_{\varepsilon}(u_{n-1})\|^{2}].$

Moreover, one has for some constant $C^{\prime}>0$ independent of
$\varepsilon$ and $\tau$ that

$\displaystyle(f_{0}(\beta_{\varepsilon}(u_{n}),\vec{s}(t_{n})),\beta_{\varepsilon}(u_{n}))$
$\displaystyle\overset{\ref{prop:fg1}}{\leq}C^{\prime}(1+\|\vec{s}(t_{n})\|^{2}+\|\beta_{\varepsilon}(u_{n})\|^{2}),$
$\displaystyle(\nabla u_{n},\nabla\beta_{\varepsilon}(u_{n}))$
$\displaystyle\geq\varepsilon\|\nabla u_{n}\|^{2}.$

Then, combining these inequalities we obtain

$\|\beta_{\varepsilon}(u_{n})\|^{2}+\|\beta_{\varepsilon}(u_{n})-\beta_{\varepsilon}(u_{n-1})\|^{2}+\varepsilon\tau\|\nabla
u_{n}\|^{2}\leq\|\beta_{\varepsilon}(u_{n-1})\|^{2}+\tau
C^{\prime}\left(1+\|s_{n}\|^{2}+\|\beta_{\varepsilon}(u_{n})\|^{2}\right).$

Applying the discrete Gronwall Lemma (2.5b) for small enough $\tau>0$, we have
for a constant $C>0$ independent of $N$ or $\varepsilon$ that

$\displaystyle\|\beta_{\varepsilon}(u_{N})\|^{2}+\sum_{n=0}^{N}[\varepsilon\|\nabla
u_{n}\|^{2}\tau+\|\beta_{\varepsilon}(u_{n})-\beta_{\varepsilon}(u_{n-1})\|^{2}]\leq\|\beta_{\varepsilon}(u_{0})\|^{2}+C+C\sum_{n=0}^{N}\|\vec{s}(t_{n})\|^{2}\tau.$
(4.11)

For $\tau>0$ small enough, one can estimate

$\displaystyle\sum_{n=0}^{N}\|\vec{s}(t_{n})\|^{2}\tau\leq\left(1+\int_{0}^{T}\|\vec{s}\|^{2}\right).$
(4.12)

Combining (4.11)-(4.12) we conclude that $\beta_{\varepsilon}(u_{N})$ and, in
extension of the method, all $\beta_{\varepsilon}(u_{n})$ are uniformly
bounded in $L^{2}(\Omega)$ with respect to $N$ and $\varepsilon$. Then, the
definition (4.7) implies that $\beta_{\varepsilon}(\hat{u}_{\tau})$
($=\beta_{\varepsilon}(u_{n})$ for $t\in(t_{n-1},t_{n}]$) and
$\beta_{\varepsilon}(\bar{u}_{\tau})$
($\leq\max\\{\beta_{\varepsilon}(u_{n}),\beta_{\varepsilon}(u_{n-1})\\}$ for
$t\in(t_{n-1},t_{n}]$) are uniformly bounded.

(Step 2) Uniform boundedness of $\|\nabla u_{\tau}\|_{L^{2}(0,T;{\cal
H}^{1})}$: Let us now test (4.3) with $\zeta=u_{n}\in{\cal H}^{1}$. This
yields

$\displaystyle(\beta_{\varepsilon}(u_{n})-\beta_{\varepsilon}(u_{n-1}),u_{n})+\tau\|\nabla
u_{n}\|^{2}=\tau(f_{0}(\beta_{\varepsilon}(u_{n}),\vec{s}(t_{n})),u_{n}).$
(4.13)

Now, from the convexity of the function $\int_{0}^{m}\Phi_{\varepsilon}$ (see
(3.1)), one has

$\displaystyle\int_{\beta_{\varepsilon}(u_{n-1})}^{\beta_{\varepsilon}(u_{n})}\Phi_{\varepsilon}\leq\Phi_{\varepsilon}(\beta_{\varepsilon}(u_{n}))(\beta_{\varepsilon}(u_{n})-\beta_{\varepsilon}(u_{n-1}))=u_{n}(\beta_{\varepsilon}(u_{n})-\beta_{\varepsilon}(u_{n-1})).$

For the last term, we observe that

$\displaystyle(f_{0},u_{n})\overset{\eqref{Eq:YoungsIneq}}{\leq}\frac{1}{2}[\|f_{0}(\beta_{\varepsilon}(u_{n}),\vec{s}(t_{n}))\|^{2}+\|u_{n}\|^{2}]\overset{\ref{prop:fg1}}{\leq}C[1+\|\beta_{\varepsilon}(u_{n})\|^{2}+\|\vec{s}(t_{n})\|^{2}+\|u_{n}\|^{2}].$
(4.14)

Hence, summing the inequalities from $n=0$ to $n=N$, using the uniform
boundedness of $\|\beta_{\varepsilon}(u_{n})\|$ from (4.11), and

$\sum_{n=1}^{N}\int_{\beta_{\varepsilon}(u_{n-1})}^{\beta_{\varepsilon}(u_{n})}\Phi_{\varepsilon}=\int^{\beta_{\varepsilon}(u_{N})}_{0}\Phi_{\varepsilon}-\int^{\beta_{\varepsilon}(u_{0})}_{0}\Phi_{\varepsilon},$

the estimate becomes

$\displaystyle\int_{\Omega}\int^{\beta_{\varepsilon}(u_{N})}_{0}\Phi_{\varepsilon}+\sum_{n=0}^{N}\|\nabla
u_{n}\|^{2}\tau\leq\int_{\Omega}\int^{M_{0}}_{0}\Phi_{\varepsilon}+C\,T+C\sum_{n=0}^{N}\left(\|\vec{s}(t_{n})\|^{2}+\|u_{n}\|^{2}\right)\,\tau.$
(4.15)

Hence, noting that
$\int_{0}^{T}\|\hat{u}_{\tau}\|^{2}=\sum_{n=0}^{N}\tau\|u_{n}\|^{2}$ we have
that $\int_{0}^{T}\|\nabla\hat{u}_{\tau}\|^{2}=\sum_{n=1}^{N}\|\nabla
u_{n}\|^{2}\,\tau$ is bounded as stated in (4.9a), and correspondingly,
following its definition, the other interpolate $\bar{u}_{\tau}$ is also
bounded in $L^{2}(0,T;{\cal H}^{1})$ as in (4.9b).

(Step 3) Uniform boundedness of
$\|\partial_{t}\beta_{\varepsilon}(\bar{u}_{\tau})\|_{L^{2}(0,T;{\cal
H}^{-1})}$: We have from (4.3) that

$\displaystyle\|\partial_{t}\beta_{\varepsilon}(\bar{u}_{\tau})\|_{{\cal
H}^{-1}}\overset{\eqref{eq:prop_ubar}}{=}\|\tfrac{1}{\tau}(\beta_{\varepsilon}(u_{n})-\beta_{\varepsilon}(u_{n-1}))\|_{{\cal
H}^{-1}}=\sup_{\zeta\in{\cal
H}^{1}}\frac{\frac{1}{\tau}\langle\beta_{\varepsilon}(u_{n})-\beta_{\varepsilon}(u_{n-1}),\zeta\rangle}{\|\zeta\|_{{\cal
H}^{1}}}$ $\displaystyle=\sup_{\zeta\in{\cal H}^{1}}\frac{-(\nabla
u_{n},\nabla\zeta)+(f_{0}(\beta_{\varepsilon}(u_{n}),\vec{s}(t_{n})),\zeta)}{\|\zeta\|_{{\cal
H}^{1}}}\leq\|\nabla
u_{n}\|+C(1+\|\vec{s}(t_{n})\|+\|\beta_{\varepsilon}(u_{n})\|).$

The bound in (4.9b) for $\partial_{t}\beta_{\varepsilon}(\bar{u}_{\tau})$ now
follows from Steps 1 and 2.

Observe that, since $\beta_{\varepsilon}$ satisfies (4.1), the
$\int_{0}^{T}\|\hat{u}_{\tau}\|^{2}$ terms on the right hand side of
(4.9a)–(4.9b) can be bounded above using the Gronwall Lemma in (4.9a), which
yields the uniform boundedness of the quantities stated in (4.3) in their
respective spaces with respect to $\tau$. However, note that the bounds may
still depend on $\varepsilon>0$. ∎

######  Lemma 4.4 (Higher regularity of the time interpolates for
$M_{0}\in{\cal H}^{1}$).

Let the assumptions of Lemma 4.3 hold. If, in addition $M_{0}\in{\cal H}^{1}$,
then for a constant $C>0$ independent of $\tau$, one has

$\displaystyle\|\nabla\hat{u}_{\tau}(t)\|^{2}+\tfrac{\varepsilon}{2}\int_{0}^{T}\|\partial_{t}\beta_{\varepsilon}(\bar{u}_{\tau})\|^{2}\leq\|\nabla\Phi_{\varepsilon}(M_{0})\|^{2}+C+C\int_{0}^{T}\|\vec{s}\|^{2}.$
(4.16)

###### Proof.

We insert the test function $\zeta=u_{n}-u_{n-1}$ in (4.3). This gives term-
wise

$\displaystyle(\beta_{\varepsilon}(u_{n})-\beta_{\varepsilon}(u_{n-1}),u_{n}-u_{n-1})$
$\displaystyle\overset{\eqref{def:BetaEps}}{\geq}\varepsilon\tau^{2}\int_{\Omega}\left|\frac{\beta_{\varepsilon}(u_{n})-\beta_{\varepsilon}(u_{n-1})}{\tau}\right|^{2}\overset{\eqref{eq:prop_ubar}}{=}\tau^{2}\varepsilon\|\partial_{t}\beta_{\varepsilon}(\bar{u}_{\tau})\|^{2},$
$\displaystyle\tau(\nabla u_{n},\nabla(u_{n}-u_{n-1}))$
$\displaystyle=\tfrac{\tau}{2}[\|\nabla u_{n}\|^{2}-\|\nabla
u_{n-1}\|^{2}+\|\nabla(u_{n}-u_{n-1})\|^{2}],$
$\displaystyle\tau(f_{0},u_{n}-u_{n-1})$
$\displaystyle\leq\tfrac{\tau^{2}}{2\varepsilon^{3}}\|f_{0}\|^{2}+\tfrac{\varepsilon^{3}}{2}\|u_{n}-u_{n-1}\|^{2}\overset{\eqref{def:BetaEps},\eqref{eq:prop_ubar}}{\leq}\tfrac{\tau^{2}}{2\varepsilon^{3}}\|f_{0}\|^{2}+\tfrac{\varepsilon\tau^{2}}{2}\|\partial_{t}\beta_{\varepsilon}(\bar{u}_{\tau})\|^{2}.$

Similarly to (4.14), we obtain

$\displaystyle\|f_{0}\|^{2}\leq
C(1+\|\vec{s}(t_{n})\|^{2}+\|\beta_{\varepsilon}(u_{n})^{2}\|)\leq
C\left[1+\|\vec{s}(t_{n})\|^{2}+\sum_{n=1}^{N}\,\|\vec{s}(t_{n})\|^{2}\tau\right],$

where we used Lemma 4.3. Finally, summing the resulting inequalities from
$n=1$ to $n=N$ and cancelling out $\tau$ one has

$\displaystyle\|\nabla
u_{N}\|^{2}+\tfrac{\varepsilon}{2}\sum_{n=1}^{N}\|\partial_{t}\beta_{\varepsilon}(\bar{u}_{\tau})\|^{2}\tau\leq\|\nabla\Phi_{\varepsilon}(M_{0})\|^{2}+C\sum_{n=1}^{N}(1+\|\vec{s}(t_{n})\|^{2})\tau,$
(4.17)

which proves the lemma. ∎

######  Remark 4.1 (Covering homogeneous Neumann condition).

The above lemmas cover both homogeneous mixed boundary conditions and
homogeneous Neumann conditions. In the latter case, ${\cal
H}^{1}=H^{1}(\Omega)$. To cover the case of inhomogeneous mixed boundary
conditions, we have to test with $\zeta=\beta_{\varepsilon}(u_{n})-h^{e}_{0}$
in Step 1 of Lemma 4.3 and with $\zeta=u_{n}-\Phi_{\varepsilon}(h^{e}_{0})$ in
Step 2. The details are straightforward, and hence, omitted.

####  4.1.3 Proof of Lemma 3.1

(Step 1) Existence: Note that $\beta_{\varepsilon}$ is Lipschitz and strictly
increasing by (4.1). Using this fact and applying Gronwall’s Lemma to (4.9a)
implies that $\hat{u}_{\tau},\,\beta_{\varepsilon}(\hat{u}_{\tau})$ are
uniformly bounded in $L^{\infty}(0,T;L^{2}(\Omega))$ with respect to $\tau$.
Consequently $\bar{u}_{\tau},\,\beta_{\varepsilon}(\bar{u}_{\tau})$ are
uniformly bounded as well by (4.9b). Thus, using (4.9b) we obtain the uniform
boundedness of $\beta_{\varepsilon}(\bar{u}_{\tau})\in{\bf{\cal X}}$. Due to
the compact embedding of $L^{2}(Q)$ in ${\bf{\cal X}}$ [23], there exists
$u\in{\bf{\cal X}}$ such that along a subsequence $\tau\to 0$,

$\displaystyle\beta_{\varepsilon}(\bar{u}_{\tau})\rightharpoonup\beta_{\varepsilon}(u)\text{
weakly in }{\bf{\cal X}}=L^{2}(0,T;{\cal H}^{1})\cap H^{1}(0,T;{\cal
H}^{-1}),$ (4.18a)
$\displaystyle\beta_{\varepsilon}(\bar{u}_{\tau})\to\beta_{\varepsilon}(u)\text{
strongly in }L^{2}(Q).$ (4.18b) Using (4.11), one has
$\displaystyle\int_{0}^{T}\|\beta_{\varepsilon}(\hat{u}_{\tau})-\beta_{\varepsilon}(\bar{u}_{\tau})\|^{2}=\sum_{n=1}^{N}\int_{t_{n-1}}^{t_{n}}\left\|\tfrac{t_{n}-t}{{\tau}}(\beta_{\varepsilon}(u_{n})-\beta_{\varepsilon}(u_{n-1}))\right\|^{2}\mathrm{d}t$
$\displaystyle=\sum_{n=1}^{N}\left\|(\beta_{\varepsilon}(u_{n})-\beta_{\varepsilon}(u_{n-1}))\right\|^{2}\int_{t_{n-1}}^{t_{n}}\left(\tfrac{t_{n}-t}{{\tau}}\right)^{2}\mathrm{d}t=\tfrac{{\tau}}{3}\sum_{n=1}^{N}\|\beta_{\varepsilon}(u_{n})-\beta_{\varepsilon}(u_{n-1})\|^{2}\to
0.$

This, along with the uniform bound with respect to $\tau$ of
$\hat{u}_{\tau},\,\bar{u}_{\tau}$ in $L^{2}(0,T;{\cal H}^{1})$ in (4.9) and
the strict monotonicity of $\beta_{\varepsilon}$ in (4.1) implies that

$\displaystyle\hat{u}_{\tau},\,\bar{u}_{\tau}\rightharpoonup u\text{ weakly in
}L^{2}(0,T;{\cal H}^{1}),$ (4.18c)
$\displaystyle\hat{u}_{\tau},\,\bar{u}_{\tau}\to u\text{ strongly in
}L^{2}(Q).$ (4.18d)

Observe from (4.3) and (4.7) that $\hat{u}_{\tau}$ and $\bar{u}_{\tau}$
satisfy

$\displaystyle\int_{0}^{T}[\langle\zeta,\partial_{t}\beta_{\varepsilon}(\bar{u}_{\tau})\rangle+(\nabla\hat{u}_{\tau},\nabla\zeta)]=\int_{0}^{T}(f_{0}(\beta_{\varepsilon}(\hat{u}_{\tau}),\vec{s}),\zeta),$
(4.18e)

for all $\zeta\in C([0,T];{\cal H}^{1})$ which is dense in $L^{2}(0,T;{\cal
H}^{1})$. Passing to the limit $\tau\to 0$ we conclude from (4.18) that $u$
solves the system

$\int_{0}^{T}[\langle\zeta,\partial_{t}\beta_{\varepsilon}(u)\rangle+(\nabla
u,\nabla\zeta)]=\int_{0}^{T}(f_{0}(\beta_{\varepsilon}(u),\vec{s}),\zeta).$

Defining $M_{s,\varepsilon}=\beta_{\varepsilon}(u)$ we obtain the desired
solution. We conclude that $M_{s,\varepsilon}\in{\bf{\cal X}}\hookrightarrow
C([0,T];L^{2}(\Omega))$, see [11, Section 5.9] for the continuous embedding
result.

(Step 2) A-priori bounds: The a-priori estimate (3.4) follows by inserting
$\varphi=M_{s,\varepsilon}$ and
$\varphi=\Phi_{\varepsilon}(M_{s,\varepsilon})$ in (3.3) and proceeding
similar to the steps of the time-discrete case in Lemma 4.3. Lemma 4.4 shows
that if in addition $M_{0}\in H^{1}(\Omega)$ then
$\partial_{t}M_{s,\varepsilon}\in L^{2}(Q)$ which implies by the definition of
weak derivatives that

$\Delta\Phi_{\varepsilon}(M_{s,\varepsilon})=(\partial_{t}M_{s,\varepsilon}-f_{0})\in
L^{2}(Q).$

Multiplying the above equation with
$\partial_{t}\Phi_{\varepsilon}(M_{s,\varepsilon})=\partial_{t}u\in L^{2}(Q)$,
integrating in $Q$, and using integration by parts we conclude that

$\int_{Q}\partial_{t}u\,\Delta
u=-\int_{0}^{T}\partial_{t}(\tfrac{1}{2}\|\nabla
u\|^{2})=\tfrac{1}{2}\|\nabla\Phi_{\varepsilon}(M_{0})\|^{2}-\tfrac{1}{2}\|\nabla
u(T)\|^{2},$

which proves (3.1). The detailed steps mimic its discrete counterpart in Lemma
4.4.

###  4.2 Finite time blow-up

The model (1.1) breaks down when $M$ reaches $1$. Henceforth, we will refer to
this as _blow-up_. Unlike in the case of Dirichlet or mixed boundary
conditions, this situation cannot in general be excluded for homogeneous
Neumann conditions. Whether a solution will blow-up in finite time or not
depends on the initial values $M_{0},\,\vec{S}_{0}$. One can construct cases
when the solution will definitely blow-up in finite time. We give a simple
example below.

######  Example 4.1 (Constant initial states).

Let us focus on the cellulolytic biofilm model with a single substrate [6],
i.e., we look at the system

$\displaystyle\partial_{t}M$ $\displaystyle=\Delta\Phi(M)+f_{0}(M,S_{1})\text{
in }Q,$ $\displaystyle M(0)$ $\displaystyle=M_{0}\text{ in }\Omega,$
$\displaystyle[\nabla M\cdot\bm{\hat{n}}]|_{\partial\Omega}$
$\displaystyle=0,$ (4.19a) $\displaystyle\partial_{t}S_{1}$
$\displaystyle=f_{1}(M,S_{1})\text{ in }Q,$ $\displaystyle S_{1}(0)$
$\displaystyle=S_{0,1}\text{ in }\Omega.$ (4.19b)

Moreover, the reaction terms are given by a non-dimensionalized version of
(1.2),

$\displaystyle f_{0}(m,s)=\left(\frac{s}{1+s}-\lambda\right)m,\quad
f_{1}(m,s)=-\frac{s\,m}{1+s}.$ (4.20)

For the initial and boundary values we assume that

$\displaystyle M_{0}\equiv\bar{M}\in(0,1),\quad S_{0,1}\equiv\bar{S},$ (4.21)

where $\bar{M},\,\bar{S}>0$ are given constants. Then it is clear that the
solution $(M,S_{1})$ of (1.1) remains constant in space for a given time.
Hence, the system evolves according to the system of ODEs

$\displaystyle\partial_{t}M=\frac{MS_{1}}{1+S_{1}}-\lambda\,M,\quad\partial_{t}S_{1}=-\frac{MS_{1}}{1+S_{1}}\
\text{ for }t>0\ \text{ with }(M(0),S_{1}(0))=(\bar{M},\bar{S}).$ (4.22)

Clearly, for $\bar{M}$ close to 1, $\bar{S}$ large, and $\lambda$ small, the
biomass density $M$ reaches 1 in finite time.

It is possible to generalize Lemma 3.2 to provide a necessary condition for
blow-up in finite time, or a sufficient condition for $M$ to stay bounded away
from $1$. This is stated in the following proposition for the single substrate
($k=1$) case.

######  Proposition 4.1 (Upper and lower bounds of $(M,S_{1})$).

Let (P1)–(P6) and (P4enumi) be satisfied, $k=1$ (single substrate) and
$\Gamma_{1}=\emptyset$ (homogeneous Neumann condition) in (1.1). Recall (P5),
and let $f_{0}(m,s)$ be increasing with respect to $s\geq 0$ for fixed $m$,
and let $f_{1}(m,s)$ be decreasing with respect to $m\geq 0$ for fixed $s$.
Let $(\check{M},\check{S},\hat{M},\hat{S})\in C^{1}({\mathbb{R}}^{+})^{4}$ be
the solution of the ODE system

$\displaystyle\begin{cases}\partial_{t}\check{M}=f_{0}(\check{M},\check{S}),\qquad\partial_{t}\check{S}=f_{1}(\hat{M},\check{S}),\\\
\partial_{t}\hat{M}=f_{0}(\hat{M},\hat{S}),\qquad\partial_{t}\hat{S}=f_{1}(\check{M},\hat{S}),\end{cases}$
(4.23)

with
$(\check{M},\check{S},\hat{M},\hat{S})=(\underline{M},\underline{S},\overline{M},\overline{S})$
at $t=0$. Further, assume that if $\nu_{1}>0$, then for $h_{1}$ defined in
(P6), $\check{S}(t)\leq h_{1}\leq\hat{S}(t)$ a.e. in $\partial\Omega$ for all
$t\in[0,T]$. Let $(M,S)$ be the weak solution of (1.1) in the sense of
Equation 2.9. Then for all $t\in[0,T]$,

$\displaystyle\check{M}(t)\leq M(t)\leq\hat{M}(t)\ \text{ and }\
\check{S}(t)\leq S_{1}(t)\leq\hat{S}(t)\ \text{ a.e. in }\Omega.$ (4.24)

######  Remark 4.2 (Assumptions in Proposition 4.1).

Observe that the assumptions of Proposition 4.1 are satisfied by the reaction
terms in (1.2) and (4.20) which were considered, e.g., in [7, 6]. Moreover,
the assumption $h_{1}\in[\check{S}(t),\hat{S}(t)]$ a.e. in $\partial\Omega$ is
a consistency condition that can be omitted in the case of immobile substrates
($\nu_{1}=0$) which occurs in the models for cellulolytic biofilms [6], or
when homogeneous Neumann conditions are assumed for $S$.

###### Proof.

The proof generalizes the arguments in Lemma 3.2 and follows the proof of
Proposition 1 of [18]. The existence and uniqueness of the solution
$(\check{M},\check{S},\hat{M},\hat{S})$ is evident from the Picard-Lindelöf
Theorem. Moreover, $f_{0}(m,s)$ is increasing with $s$, $f_{1}(m,s)$ is
decreasing with $m$, along with
$\hat{M}(0)=\bar{M}\geq\underline{M}=\check{M}(0)$ and
$\hat{S}(0)=\bar{S}\geq\underline{S}=\check{S}(0)$, together imply for all
$t>0$,

$\displaystyle\hat{M}(t)\geq\check{M}(t),\quad\text{ and
}\quad\hat{S}(t)\geq\check{S}(t).$ (4.25)

This follows by writing from (4.23),

$\displaystyle\frac{1}{2}[\check{M}(t)-\hat{M}(t)]_{+}^{2}=\int_{0}^{t}[\check{M}-\hat{M}]_{+}(f_{0}(\check{M},\check{S})-f_{0}(\hat{M},\hat{S})),$
$\displaystyle\frac{1}{2}[\check{S}(t)-\hat{S}(t)]_{+}^{2}=\int_{0}^{t}[\check{S}-\hat{S}]_{+}(f_{1}(\hat{M},\check{S})-f_{1}(\check{M},\hat{S})),$

for $t>0$. Then, following the manipulations in Lemma 3.2 (also repeated
below), Gronwall’s Lemma yields
$[\check{M}(t)-\hat{M}(t)]_{+}=[\check{S}(t)-\hat{S}(t)]_{+}=0$. We omit the
detailed proof for brevity.

Insert the test functions $\varphi=[M-\hat{M}]_{+}$ and
$\zeta_{1}=[S_{1}-\hat{S}]_{+}$ in (2.9). Observe that, $\zeta_{1}\in
L^{2}(0,T;H^{1}_{0}(\Omega))$ is a valid test function for $\nu_{1}>0$ since
$S_{1}-\hat{S}=h_{1}-\hat{S}\leq 0$ on $\partial\Omega$. Then following the
manipulations in Lemma 3.2, one obtains from the first equation that

$\displaystyle\int_{0}^{T}\partial_{t}\left(\frac{1}{2}\|[M-\hat{M}]_{+}\|^{2}\right)\leq\int_{0}^{T}(f_{0}(M,S_{1})-f_{0}(\hat{M},\hat{S}),[M-\hat{M}]_{+})$
$\displaystyle=\int_{0}^{T}(f_{0}(M,S_{1})-f_{0}(\hat{M},S_{1}),[M-\hat{M}]_{+})+\int_{0}^{T}(f_{0}(\hat{M},S_{1})-f_{0}(\hat{M},\hat{S}),[M-\hat{M}]_{+})$
$\displaystyle\overset{\ref{prop:fl}}{\leq}C_{L}\int_{0}^{T}\|[M-\hat{M}]_{+}\|^{2}+C_{L}\int_{0}^{T}([S_{1}-\hat{S}]_{+},[M-\hat{M}]_{+})$
$\displaystyle\leq
C\int_{0}^{T}[\|[M-\hat{M}]_{+}\|^{2}+\|[S_{1}-\hat{S}]_{+}\|^{2}].$ (4.26a)
Here, we used that $f_{0}(\hat{M},\cdot)$ is increasing to conclude that
$(f_{0}(\hat{M},S_{1})-f_{0}(\hat{M},\hat{S}))[M-\hat{M}]_{+}\leq
C_{L}[S_{1}-\hat{S}]_{+}[M-\hat{M}]_{+}.$ Similarly, from the second equation,
noting that $f_{1}(\cdot,s)$ is decreasing for a given $s$, one obtains
$\displaystyle\int_{0}^{T}\partial_{t}\left(\frac{1}{2}\|[S_{1}-\hat{S}]_{+}\|^{2}\right)\leq\int_{0}^{T}(f_{1}(M,S_{1})-f_{1}(\check{M},\hat{S}),[S_{1}-\hat{S}]_{+})$
$\displaystyle=\int_{0}^{T}(f_{1}(M,S_{1})-f_{1}(\check{M},S_{1}),[S_{1}-\hat{S}]_{+})+\int_{0}^{T}(f_{1}(\check{M},S_{1})-f_{1}(\check{M},\hat{S}),[S_{1}-\hat{S}]_{+})$
$\displaystyle\leq
C_{L}\int_{0}^{T}([M-\check{M}]_{-},[S_{1}-\hat{S}]_{+})+C_{L}\int_{0}^{T}\|[S_{1}-\hat{S}]_{+}\|^{2}$
$\displaystyle\leq
C\int_{0}^{T}[\|[M-\check{M}]_{-}\|^{2}+\|[S_{1}-\hat{S}]_{+}\|^{2}].$ (4.26b)

Finally, inserting the test functions $\varphi=[M-\check{M}]_{-}$ and
$\zeta_{1}=[S_{1}-\check{S}]_{-}$ in (2.9) we get analogous estimates to
(4.26). Adding these inequalities and using Gronwall’s Lemma completes the
proof. ∎

######  Remark 4.3 (Guaranteed finite time blow-up/ boundedness).

If the solution $\check{M}$ in Proposition 4.1 reaches $1$ in finite time,
then it implies that the solution of the original system $(M,S_{1})$ blows up
in finite time. On the other hand, if the solution $\hat{M}$ remains bounded
by a constant strictly less than $1$, then $M$ does not blow up and hence, the
solution $(M,S_{1})$ is global-in-time. The bounds are sharp if
$|\bar{M}-\underline{M}|$ and $|\bar{S}-\underline{S}|$ are small.

##  5 Spatial regularity of the biomass density

In this section, we analyze the spatial regularity of solutions of the
degenerate diffusion equation (1.1a), i.e., we focus on the scalar equation

$\displaystyle\partial_{t}M=\nabla\cdot[D(M)\nabla M]+f(M,\cdot)\qquad\text{in
}Q,$ (5.1)

where $D:[0,1)\to[0,\infty)$ and $f:[0,\infty)\times Q\to{\mathbb{R}}$. The
regularity results we derive apply to a broad class of degenerate diffusion
problems, see Remark 5.1, including the biofilm growth models [7, 6].

It is well known that the degeneracy of the diffusion coefficient $D(0)=0$
causes a finite speed of propagation and sharp fronts at the interface between
the regions $\\{M>0\\}$ and $\\{M=0\\}$ corresponding to steep gradients of
$M$. Despite this fact, the solution $M$ is locally Hölder continuous. This
was shown for porous medium type equations in [5] and for equations with
degenerate and singular diffusion in [14]. The global space-time regularity of
solutions of the porous medium equation in ${\mathbb{R}}^{d}$ has also been
studied extensively using optimal regularity theory, see [12] and the
references therein. Assuming homogeneous Neumann boundary conditions, in this
section we show that $M$ can further inherit global spatial regularity in the
more general case (5.1), i.e. $M\in L^{2}(0,T;H^{r}(\Omega))$, where $r=1$ for
$a<2$, and $r<1$ otherwise. This fact is not only mathematically intriguing
but has important consequences in designing numerical tools and test functions
for such problems. We now specify the assumptions on the functions $D$ and
$f$.

######  Assumption 5.1 (Assumptions on $D$ and $f$).

The diffusion coefficient satisfies (P1). In addition, there exists
$a\in{\mathbb{R}}^{+}$ and a constant $C>0$ such that

$D(m)\geq Cm^{a}$

for all $m\in[0,1)$. The function $f:[0,\infty)\times Q\to{\mathbb{R}}$ is
Lipschitz continuous with respect to the first variable, and there exists a
non-negative function $f_{\max}\in\mathrm{Lip}({\mathbb{R}})$ such that
$f(\cdot,(\bm{x},t))\leq f_{\max}(\cdot)$. Moreover, we assume that
$f(0,(\bm{x},t))\geq 0$ for all $(\bm{x},t)\in Q$.

######  Theorem 5.1 (Global spatial regularity of $M$).

Let $\Gamma_{1}=\emptyset$ (homogeneous Neumann condition). Let $M\in{\bf{\cal
W}}$ with $\Phi(M)=\int_{0}^{M}D\in L^{2}(0,T;H^{1}(\Omega))$, and
$M(0)=M_{0}$ (see (P5)) be the weak solution of (5.1), i.e.,

$\displaystyle\int_{0}^{T}\langle\varphi,\partial_{t}M\rangle+\int_{0}^{T}(\nabla\Phi(M),\nabla\varphi)=\int_{0}^{T}(f(M,\cdot),\varphi),$
(5.2)

for all $\varphi\in L^{2}(0,T;H^{1}(\Omega))$. Then under the Assumption 5.1,
$M\in L^{2}(0,T;H^{r}(\Omega))$ for

  1. (a)

$r=1$ if either $a<2$ or $\underline{M}=\mathrm{ess}\inf\\{M_{0}\\}>0$.

  2. (b)

all $r<2/a$, if $a\geq 2$ and $\underline{M}=0$.

######  Remark 5.1 (Generality of Theorem 5.1).

Theorem 5.1 applies to the solution $M$ of the coupled system (1.1) under the
conditions (P1)–(P6) and (P4enumi). Since the spatial irregularity of $M$
stems from the degeneracy at $M=0$, our regularity results also cover
diffusion coefficients $D$ that are degenerate but non-singular, for instance,
porous medium type equations. In this case, the additional assumption that $f$
is bounded by $f_{\max}$ can be omitted as solutions are not required to take
values in $[0,1)$.

######  Remark 5.2 (Assumptions on the boundary conditions in Theorem 5.1).

To simplify notations Theorem 5.1 is stated for homogeneous Neumann boundary
conditions. However, the result remains valid for Dirichlet or mixed boundary
conditions provided that $\Phi(M)=\Phi(h^{e}_{0})$ at $\Gamma_{1}$ and the
functions $\Psi_{\varepsilon}(h^{e}_{0})\in H^{1}(\Omega)$ are uniformly
bounded with respect to $\varepsilon\in(0,1)$, where $\Psi_{\varepsilon}$ is
introduced in (5.3) and $h^{e}_{0}$ in (P6).

The rest of this section is dedicated to the proof of Theorem 5.1. The main
idea behind the proof is to use a test function $\varphi$ of the form
$M^{-\alpha}$ ($\alpha>0$) in (5.2). However, $\varphi$ might not be a valid
test function due to $M$ not being sufficiently regular, and $M^{-\alpha}$
having a singularity at $0$. To resolve this, we will construct a modified
function that is admissible.

###  5.1 Some auxiliary functions

As in the proof of Theorem 5.1 we consider the regularized problem introduced
in Lemma 3.1. The function $\Phi_{\varepsilon}$ is taken as in (3.2). For a
given constant $\alpha>0$, we further introduce the $C^{1}({\mathbb{R}})$
function

$\displaystyle\Psi_{\varepsilon}(m):=\int_{1}^{m}\frac{\mathrm{d}\varrho}{\min\\{\max\\{\varepsilon,\varrho^{\alpha}\\},1\\}}.$
(5.3)

Note that

$\displaystyle\Psi^{\prime}_{\varepsilon}\geq 0\quad\text{ and
}\quad\Psi_{\varepsilon}(m)<0\quad\text{ for }m<1.$ (5.4)

######  Lemma 5.1 (Growth of $\Psi_{\varepsilon}$).

For a given $\alpha>0$ and $\varepsilon\in(0,1)$, let $\Psi_{\varepsilon}$ be
defined as in (5.3). Then, the following estimate holds,

$\displaystyle|m\Psi_{\varepsilon}(m)|\lesssim
1+\int_{1}^{m}\Psi_{\varepsilon}\qquad\text{ for all }m\geq 0.$ (5.5)

###### Proof.

Case 1 ($1<m$): For $m>1$, $\Psi_{\varepsilon}(m)=m-1$, and the inequality can
be verified directly.

Case 2 ($\varepsilon^{\frac{1}{\alpha}}\leq m\leq 1$): If
$\varepsilon^{\frac{1}{\alpha}}\leq m\leq 1$ and $\alpha\not=1$, we have

$\displaystyle|m\,\Psi_{\varepsilon}(m)|\leq\left|m\int^{m}_{1}\frac{1}{\varrho^{\alpha}}\right|d\rho|=\left|m\left(\frac{m^{1-\alpha}-1}{1-\alpha}\right)\right|\lesssim|m^{2-\alpha}-m|.$
(5.6a) Observe that, if $\alpha\leq 2$ then the right hand side is bounded
since $m\leq 1,$ and (5.5) definitely holds. The case when $\alpha=1$ can also
be handled rather easily since it yields $\Psi_{\varepsilon}(m)=\log(m)$. The
interesting case is when $\alpha>2$. Then, we can estimate the right-hand side
of (5.5) as follows, $\displaystyle 1+\int_{1}^{m}\Psi_{\varepsilon}$
$\displaystyle=1+\int_{1}^{m}\int_{1}^{s}\frac{1}{\varrho^{\alpha}}d\rho
ds=1+\int_{1}^{m}\frac{s^{1-\alpha}-1}{\alpha-1}\,\mathrm{d}s\gtrsim\int_{1}^{m}s^{1-\alpha}\gtrsim
m^{2-\alpha}-1.$ (5.6b)

Combining (5.6) and noting that $m<1$, we have (5.5) for this case.

Case 3 ($0\leq m<\varepsilon^{\frac{1}{\alpha}}$): We only focus on $\alpha>2$
since the case $\alpha\leq 2$ can be shown exactly as in Case 2. We observe
that

$\displaystyle|m\,\Psi_{\varepsilon}(m)|=$
$\displaystyle\left|m\Psi_{\varepsilon}(\varepsilon^{\frac{1}{\alpha}})+m\smallint^{m}_{\varepsilon^{\frac{1}{\alpha}}}\tfrac{\mathrm{d}\varrho}{\max(\varepsilon,\varrho^{\alpha})}\right|=\left|m\Psi_{\varepsilon}(\varepsilon^{\frac{1}{\alpha}})+m\smallint^{m}_{\varepsilon^{\frac{1}{\alpha}}}\tfrac{1}{\varepsilon}\right|$
$\displaystyle\leq$
$\displaystyle|m\Psi_{\varepsilon}(\varepsilon^{\frac{1}{\alpha}})|+\frac{m(\varepsilon^{\frac{1}{\alpha}}-m)}{\varepsilon}\leq|\varepsilon^{\frac{1}{\alpha}}\Psi_{\varepsilon}(\varepsilon^{\frac{1}{\alpha}})|+\varepsilon^{\frac{1}{\alpha}-1}(\varepsilon^{\frac{1}{\alpha}}-m).$
(5.7a) From Case 2 we conclude that
$|\varepsilon^{\frac{1}{\alpha}}\Psi_{\varepsilon}(\varepsilon^{\frac{1}{\alpha}})|\lesssim
1+\int^{\varepsilon^{\frac{1}{\alpha}}}_{1}\Psi_{\varepsilon}$, and we obtain
$\displaystyle 1+\int_{1}^{m}\Psi_{\varepsilon}$
$\displaystyle=1+\int_{1}^{\varepsilon^{\frac{1}{\alpha}}}\Psi_{\varepsilon}+\int^{m}_{\varepsilon^{\frac{1}{\alpha}}}\int_{1}\tfrac{\mathrm{d}\varrho}{\max(\varepsilon,\varrho^{\alpha})}=1+\int_{1}^{\varepsilon^{\frac{1}{\alpha}}}\Psi_{\varepsilon}+\int^{\varepsilon^{\frac{1}{\alpha}}}_{m}\int^{1}\tfrac{\mathrm{d}\varrho}{\max(\varepsilon,\varrho^{\alpha})}$
$\displaystyle\geq
1+\int_{1}^{\varepsilon^{\frac{1}{\alpha}}}\Psi_{\varepsilon}+\int_{m}^{\varepsilon^{\frac{1}{\alpha}}}\int_{\varepsilon^{\frac{1}{\alpha}}}^{1}\tfrac{\mathrm{d}\varrho}{\max(\varepsilon,\varrho^{\alpha})}\gtrsim|\varepsilon^{\frac{1}{\alpha}}\Psi_{\varepsilon}(\varepsilon^{\frac{1}{\alpha}})|+(\varepsilon^{\frac{1}{\alpha}}-m)\int_{\varepsilon^{\frac{1}{\alpha}}}^{1}\tfrac{\mathrm{d}\varrho}{\varrho^{\alpha}}$
$\displaystyle\gtrsim|\varepsilon^{\frac{1}{\alpha}}\Psi_{\varepsilon}(\varepsilon^{\frac{1}{\alpha}})|+\varepsilon^{\frac{1}{\alpha}-1}(\varepsilon^{\frac{1}{\alpha}}-m)-(\varepsilon^{\frac{1}{\alpha}}-m).$
(5.7b) Hence, combining again (5.7) we have (5.5).

∎

###  5.2 Boundedness of $M$ in $L^{2}(0,T;H^{r}(\Omega))$

To prove Theorem 5.1 we first show the following lemma.

######  Lemma 5.2 (An estimate for the regularized solutions).

Let Assumption 5.1 hold and $\Gamma_{1}=\emptyset$. For $\varepsilon\in(0,1)$
and $\alpha>0$, let $\Phi_{\varepsilon}$ and $\Psi_{\varepsilon}$ be defined
by (3.2) and (5.3) respectively. Let $M_{\varepsilon}\in{\bf{\cal X}}$ satisfy
$M_{\varepsilon}(0)=M_{0}$ and

$\displaystyle\int_{0}^{T}\langle\varphi,\partial_{t}M_{\varepsilon}\rangle+\int_{0}^{T}(\nabla\Phi_{\varepsilon}(M_{\varepsilon}),\nabla\varphi)=\int_{0}^{T}(f(M_{\varepsilon},\cdot),\varphi),$

for all $\varphi\in L^{2}(0,T;H^{1}(\Omega))$. Then, we have

$\displaystyle\int_{0}^{T}\int_{\Omega}\min\\{M_{\varepsilon}^{a-\alpha},1\\}|\nabla
M_{\varepsilon}|^{2}\lesssim
1+\int_{\Omega}\int_{1}^{M_{0}}\Psi_{\varepsilon}.$ (5.8)

###### Proof.

First, we show that the following estimate holds,

$\displaystyle\Phi^{\prime}_{\varepsilon}(m)\,\Psi^{\prime}_{\varepsilon}(m)\gtrsim\min\\{1,m^{a-\alpha}\\}\qquad\text{
for all }m\geq 0.$ (5.9)

We distinguish several cases. Case 1: ${\Phi_{\varepsilon}}^{\prime}(m)=D(m)$.
This implies that
$\varepsilon\leq{\Phi_{\varepsilon}}^{\prime}(m)=D(m)\leq\frac{1}{\varepsilon}$
and hence, $m<1$. If ${\Psi_{\varepsilon}}^{\prime}(m)=m^{-\alpha}$, then the
result follows from Assumption 5.1. If
${\Psi_{\varepsilon}}^{\prime}(m)=\frac{1}{\varepsilon}$, then
$\Phi^{\prime}_{\varepsilon}(m)\,\Psi^{\prime}_{\varepsilon}(m)=D(m)/\varepsilon>1$.
Finally, if ${\Psi_{\varepsilon}}^{\prime}(m)=1$, then $m\geq 1$ which is
excluded.

Case 2: ${\Phi_{\varepsilon}}^{\prime}(m)=\varepsilon$. Consequently, $m<1$.
The definition of $\Phi_{\varepsilon}$ in (3.2) implies that $\varepsilon\geq
D(m)\gtrsim m^{a}$, where the last inequality holds by Assumption 5.1. Hence,
for $\Psi^{\prime}_{\varepsilon}=\varepsilon^{-1}$ the product
$\Phi^{\prime}_{\varepsilon}\,\Psi^{\prime}_{\varepsilon}=1$ and for
$\Psi^{\prime}_{\varepsilon}=m^{-\alpha}$ we have
$\Phi^{\prime}_{\varepsilon}\,\Psi^{\prime}_{\varepsilon}\gtrsim
m^{a-\alpha}$.

Case 3: ${\Phi_{\varepsilon}}^{\prime}=\frac{1}{\varepsilon}$. This case
follows similarly.

Inserting the test function $\varphi=\Psi_{\varepsilon}(M_{\varepsilon})$ in
(3.3), the first term becomes

$\displaystyle\int_{0}^{T}\langle\partial_{t}M_{\varepsilon},\Psi_{\varepsilon}(M_{\varepsilon})\rangle=\int_{\Omega}\int_{1}^{M_{\varepsilon}(T)}\Psi_{\varepsilon}-\int_{\Omega}\int_{1}^{M_{0}}\Psi_{\varepsilon}.$
(5.10a) The second term of (3.3) gives
$\displaystyle\int_{0}^{T}(\nabla\Phi_{\varepsilon}(M_{\varepsilon}),\nabla\Psi_{\varepsilon}(M_{\varepsilon}))$
$\displaystyle=\int_{0}^{T}\int_{\Omega}\Phi^{\prime}_{\varepsilon}(M_{\varepsilon})\,\Psi^{\prime}_{\varepsilon}(M_{\varepsilon})|\nabla
M_{\varepsilon}|^{2}$
$\displaystyle\overset{\eqref{eq:PhiPsi}}{\gtrsim}\int_{0}^{T}\int_{\Omega}\min\\{M_{\varepsilon}^{a-\alpha},1\\}|\nabla
M_{\varepsilon}|^{2}.$ (5.10b) Finally, the third term of (3.3) yields using
$f(0,\cdot)\geq 0$ and (5.4) that
$\displaystyle\int_{0}^{T}(f(M_{\varepsilon},\cdot),\Psi_{\varepsilon}(M_{\varepsilon}))=\int_{0}^{T}(f(M_{\varepsilon},\cdot)-f(0,\cdot),\Psi_{\varepsilon}(M_{\varepsilon}))+\int_{0}^{T}(f(0,\cdot),\Psi_{\varepsilon}(M_{\varepsilon}))$
$\displaystyle\quad\overset{\ref{prop:fl}}{\leq}C_{L}\int_{0}^{T}\int_{\Omega}|M_{\varepsilon}\Psi_{\varepsilon}(M_{\varepsilon})|\overset{f(0,\cdot)\geq
0}{+}\int_{0}^{T}(f(0,\cdot),[\Psi_{\varepsilon}(M_{\varepsilon})]_{+})$
$\displaystyle\quad\lesssim\int_{0}^{T}\int_{\Omega}|M_{\varepsilon}\Psi_{\varepsilon}(M_{\varepsilon})|\overset{\eqref{eq:PsiIneq}}{\lesssim}\int_{0}^{T}\int_{\Omega}[1+\int_{1}^{M_{\varepsilon}}\Psi_{\varepsilon}].$
(5.10c)

In the above, noting that $\Psi_{\varepsilon}(m)>0$ only when $m>1$, we
estimated $f(0,\cdot)[\Psi_{\varepsilon}(M_{\varepsilon})]_{+}\leq
f_{\max}(0)|M_{\varepsilon}\Psi_{\varepsilon}(M_{\varepsilon})|$. Combining
the inequalities (5.10) we have

$\displaystyle\int_{\Omega}\int_{1}^{M_{\varepsilon}(T)}\Psi_{\varepsilon}+\int_{0}^{T}\int_{\Omega}\min\\{M_{\varepsilon}^{a-\alpha},1\\}|\nabla
M_{\varepsilon}|^{2}\lesssim
1+\int_{\Omega}\int_{1}^{M_{0}}\Psi_{\varepsilon}+\int_{0}^{T}\int_{\Omega}\int_{1}^{M_{\varepsilon}}\Psi_{\varepsilon}.$
(5.11)

Using Gronwall’s Lemma (2.5a) the estimate (5.8) follows. ∎

To conclude the proof of Theorem 5.1 from (5.8), we need the following lemma.
For its proof we refer to Lemma 1.3 and Lemma B.1 of [24].

######  Lemma 5.3 (Property of $H^{r}(\Omega)$).

If $u^{\gamma}\in H^{1}(\Omega)$ for some $\gamma>1$ then $u\in H^{r}(\Omega)$
for all $r\in(0,\gamma^{-1}]$.

###### Proof of Theorem 5.1.

Case 1 ($\underline{M}=\mathrm{ess}\inf\\{M_{0}\\}>0$): In this case, taking
$\varepsilon_{1}^{\frac{1}{\alpha}}<\underline{M}$ we conclude by (5.3) that

$\int_{\Omega}\int_{1}^{M_{0}}\Psi_{\varepsilon}\quad\text{ is uniformly
bounded for all }\varepsilon\leq\varepsilon_{1}.$

Hence, taking $\alpha=a$ in Lemma 5.2 provides a uniform bound on
$\int_{0}^{T}\|\nabla M_{\varepsilon}\|^{2}$. Moreover,
$\|M_{\varepsilon}\|_{L^{\infty}(0,T;L^{\infty}(\Omega))}$ is bounded by Lemma
3.2. Hence, $M_{\varepsilon}$ is uniformly bounded in
$L^{2}(0,T;H^{1}(\Omega))$. Passing to the limit $\varepsilon\to 0$, the
convergence of $M_{\varepsilon}$ to a unique $M\in{\bf{\cal W}}$ follows from
Lemma 3.3 (see (3.12)). Consequently, the uniform bound implies that $M\in
L^{2}(0,T;H^{1}(\Omega))$.

Case 2 ($a<2$): In this case, put $\alpha=a$. Then, passing the limit
$\varepsilon\to 0$ on the right hand side of (5.8) one has

$\lim\limits_{\varepsilon\searrow
0}\int_{\Omega}\int_{1}^{M_{0}}\Psi_{\varepsilon}\lesssim\int_{\Omega}\int_{M_{0}}^{1}(1-\varrho^{1-a})\,\mathrm{d}\varrho\lesssim
1.$

Here we used that $M^{2-a}_{0}\leq 1$ a.e. in $\Omega$ by assumption (P5).
Hence, we again obtain a uniform bound on $\int_{0}^{T}\|\nabla
M_{\varepsilon}\|^{2}$ and consequently, $M\in L^{2}(0,T;H^{1}(\Omega))$.

Case 3 ($a\geq 2$): We set $\alpha=2-\delta$ for sufficiently small
$\delta>0$. Then the previous case and (5.8) gives that
$(M_{\varepsilon})^{1+\frac{a-\alpha}{2}}\in L^{2}(0,T;H^{1}(\Omega))$ and it
is uniformly bounded. From Lemma 5.3 it follows that $M_{\varepsilon}\in
L^{2}(0,T;H^{r}(\Omega))$ for

$r\leq\frac{2}{2+a-\alpha}=\frac{2}{a+\delta}\qquad\text{and sufficiently
small }\delta>0.$

This concludes the proof. ∎

### Acknowledgements

K. Mitra and S. Sonner would like to thank the Nederlandse Organisatie voor
Wetenschappelijk Onderzoek (NWO) for their support through the Grant
OCENW.KLEIN.358. K. Mitra was additionally supported by Fonds voor
Wetenschappelijk Onderzoek (FWO) through the Junior Postdoctoral Fellowship
during the completion of this work.

## References

  * [1] H.W. Alt and S. Luckhaus.  Quasilinear elliptic-parabolic differential equations.  Mathematische Zeitschrift, 183(3):311–341, 1983.
  * [2] H.W. Alt, S. Luckhaus, and A. Visintin.  On nonstationary flow through porous media.  Annali di Matematica Pura ed Applicata, 136(1):303–316, 1984.
  * [3] B.P. Andreianov and F. Bouhsiss.  Uniqueness for an elliptic-parabolic problem with Neumann boundary condition.  Journal of Evolution Equations, 4(2):273–295, 2004.
  * [4] H. Brezis.  Functional Analysis, Sobolev Spaces and Partial Differential Equations, vol. 2.  Springer, 2011.
  * [5] E. DiBenedetto and A. Friedman.  Hölder estimates for nonlinear degenerate parabolic systems.  Journal für die reine und angewandte Mathematik, 1985(357): 1–22, 1985.
  * [6] 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:141–151, 2017.
  * [7] H.J. Eberl, D.F. Parker, and M.C.M. van Loosdrecht.  A new deterministic spatio-temporal continuum model for biofilm development.  Computational and Mathematical Methods in Medicine, 3(3):161–175, 2001.
  * [8] M.A. Efendiev, M. Otani, and H.J. Eberl.  Mathematical analysis of a PDE-ODE coupled model of mitochondrial swelling with degenerate Calcium ion diffusion.  SIAM Journal on Mathematical Analysis, 52(1):543–569, 2020.
  * [9] M.A. Efendiev, S. Zelik, and H.J. Eberl.  Existence and longtime behavior of a biofilm model.  Communications on Pure & Applied Analysis, 8(2):509–531, 2009\.
  * [10] 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 and Engineering, 14(3):625–653, 2017.
  * [11] L.C. Evans.  Partial differential equations.  Wiley Online Library, 1988.
  * [12] B. Gess, J. Sauer, and E. Tadmor.  Optimal regularity in time and space for the porous medium equation.  Analysis and PDE, 13(8): 2441-2480, 2020.
  * [13] P. Hartman and G. Stampacchia.  On some non-linear elliptic differential-functional equations.  Acta Mathematica, 115:271–310, 1966.
  * [14] 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: 92, 2022.
  * [15] 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(1):125894, 2022.
  * [16] J. Kačur.  Method of Rothe in evolution equations, vol. 1192.  Springer, 1986.
  * [17] K. Kawasaki, A. Mochizuki, M. Matsushita, T. Umeda, and N. Shigesada.  Modeling spatio-temporal patterns generated bybacillus subtilis.  Journal of Theoretical Biology, 188(2):177–185, 1997.
  * [18] K. Mitra.  Existence and properties of solutions of extended play-type hysteresis model.  Journal of Differential Equations, 288: 118–140, 2020.
  * [19] K. Mitra, J.M. Hughes, S. Sonner, H.J. Eberl, and J.D. Dockery.  Travelling Waves in a PDE-ODE Coupled Model of Cellulolytic Biofilms with Nonlinear Diffusion.  Journal of Dynamics and Differential Equations, 2023.
  * [20] J.D. Murray.  Mathematical Biology I. An Introduction.  Springer, 2002.
  * [21] F. Otto.  $L^{1}$-contraction and uniqueness for quasilinear elliptic–parabolic equations.  Journal of Differential Equations, 131(1):20–38, 1996.
  * [22] S. Salsa.  Partial Differential Equations in Action. From Modelling to Theory. Third edition. vol. 99.  Springer, 2016.
  * [23] J. Simon.  Compact sets in the space ${L}^{p}(0,{T};{B})$.  Annali di Matematica pura ed applicata, 146(1):65–96, 1986.
  * [24] S. Sonner.  Systems of quasi-linear PDEs arising in the modelling of biofilms and related dynamical questions.  PhD thesis, Technische Universität München, 2012.
  * [25] S. Sonner, M.A. Efendiev, and H.J. Eberl.  On the well-posedness of mathematical models for multicomponent biofilms.  Mathematical Methods in the Applied Sciences, 38(17):3753–3775, 2015.
  * [26] E. Zeidler.  Applied Functional Analysis: Applications to Mathematical Physics (Applied Mathematical Sciences) vol. 108.  Springer, 1995.