text:
\section*{Introduction}
Catalysis is a key technology in Chemical Engineering, employed not only to increase the speed of chemical
reactions by up to several orders of magnitude, but also to change the selectivity in favor of a desired product
against other possible output components of a chemical reaction network.
In heterogeneous catalysis, the catalytic substance forms a separate phase which is advantageous
concerning the separation of the products from the catalytic material. A prototypical setting,
which also underlies the mathematical model below, consists of a solid phase catalyst brought into contact with a gas or liquid which carries the educts as well as the product species inside the chemical reactor.
In this case, the overall chemical conversion consists of the following steps:
\begin{enumerate}
\item
the educt species are transported to the surface of the catalytic substance;
\item
molecules of at least one educt species adsorb at the catalyst surface;
\item
adsorbed molecules react, either with other adsorbed molecules or
with molecules in the bulk phase directly adjacent to the surface;
\item
the product molecules are desorbed.
\end{enumerate}
Of course, further processes will usually appear as well. For instance, adsorbed educt molecules may desorb
back into the bulk before a chemical reaction occurs, or they can be transported along the surface by means
of surface diffusion processes.
For a recent view on the complexity of heterogeneous catalysis modeling see \cite{Keil:Complexities}.

In the present paper, we only consider the case of pure surface chemistry, i.e.\
chemical reactions are only allowed between adsorbed species. This is actually no restriction, since one may
otherwise introduce an artificial adsorbed form of the reaction partner which is in the bulk adjacent to the surface and assign to it an infinite adsorption rate such that all arriving bulk molecules immediately adsorb
and, hence, are available for surface reaction.

In order for a heterogenous catalytic process to be efficient, a large surface area is required.
Therefore, in classical heterogeneous catalysis with solid phase catalyst, the latter is often provided as
a porous structure, e.g.\ in so-called packed-bed reactors. In this case, the smallest unit is a single
pore, into which the educts have to be transported in order to reach the pore wall, i.e.\ the catalytic surface.
More information on this classical reactor concept can be found, e.g., in \cite{Aris:Permeable-Catalysts},
\cite{Levenspiel:Chemical-Engineering} or \cite{White:HetCat}.
In recent years, with the advent of microreactor engineering technology, new reactor designs became feasible.
Due to the large area-to-volume ratio at the micro scale, multichannel microreactors with catalytic wall
coatings can replace classical porous structures and still provide fast and intense diffusive transport
to the channel walls in order to facilitate the reaction speed or selectivity enhancement; see, e.g.,
\cite{Wirth:Microreactors}, \cite{Renken:Catalytic-Microreactors}.
Since the given and precise structure of microreactors together with modern control and measurement techniques
allows for defined and reproducible operating conditions, this approach is much better accessible for
detailed quantitative modeling and simulation; cf.\ \cite{Bothe-Lojewski-Instantaneous},
\cite{Bothe-Lojewski-Parabolized}.
Structured catalytic microreactors are also employed for efficient screening of potential catalysts for
new reaction pathways; see, e.g., \cite{Jensen:HetCat}. To tap the full potential of such microsystems approaches and to intensify also more classical heterogeneous catalysis processes, realistic and sound mathematical models
are required as the basis for any numerical simulation. The most fundamental question then is whether a given
model is well-posed, a necessary requirement to enable any reasonable numerical treatment.

In what follows, we consider a single pore as a prototypical element, where we allow for convection
through the pore with a solenoidal velocity field which is assumed to be known and to satisfy the no-slip
boundary condition at fixed walls. We focus on pores having smoothly bounded cross shapes.
Let therefore $\Omega:=A\times(-h,h)\subset \mathbb{R}^3$ denote a finite three-dimensional cylinder of height $2h>0$ with cross section $A\subset \mathbb{R}^2$ being a bounded simply connected $C^2$-domain, such that $\partial A$ is a closed regular $C^2$-curve.
The boundary of $\Omega$ decomposes into bottom $\Gamma_\text{in}$, top $\Gamma_\text{out}$ and lateral surface $\Sigma$, standing for inflow area, outflow area and active surface.
The mathematical model consists of the partial mass balances for all involved chemical components,
both within the bulk phase $\Omega$ (representing the interior of the pore) and on the active surface $\Sigma$
(representing the catalytic surface). Inside the bulk phase, the species mass fluxes are due to
advection and diffusion, where we assume the latter to be governed by Fick's law.
On the active surface, we only consider diffusive fluxes along the surface, again assuming Fick's law
to be a reasonable constitutive relation. We allow for different diffusivities but the model ignores
cross-diffusion effects. Let us note in passing that for high surface coverage, cross-effects
between the transport of different constituents will appear which are not accounted for by our model.
The mass exchange between bulk and active surface is due to ad- and desorption phenomena, which are
usually modeled via kinetic relations in analogy to chemical reaction kinetics. Examples will be discussed below.

Insertion of the flux relations into the partial mass balances for continua yields the following mathematical
model for the unknown concentrations $(c_i,c^{\scriptscriptstyle\Sigma}_i)$ with $i=1,...,N$:
\begin{equation}\label{eq:cat}
\left\{\begin{array}{rclcl}
\partial_t c_i+(u\cdot \nabla) c_i- d_i\Delta c_i &=&0& \text{in}& (0,T)\times\Omega, \\
\partial_t c^{\scriptscriptstyle\Sigma}_i - d^{\scriptscriptstyle\Sigma}_i \Delta_\Sigma c^{\scriptscriptstyle\Sigma}_i &=&r^\text{sorp}_i(c_i,c^{\scriptscriptstyle\Sigma}_i)+r^\text{ch}_i(c^{\scriptscriptstyle\Sigma}) & \text{on}& (0,T)\times \Sigma, \\
(u\cdot\nu)c_i-d_i\partial_\nu c_i&=&g^\text{in}_i& \text{on}& (0,T)\times \Gamma_\text{in}, \\
-d_i\partial_\nu c_i &=&r^\text{sorp}_i(c_i,c^{\scriptscriptstyle\Sigma}_i)& \text{on}& (0,T)\times \Sigma,\\
-d_i\partial_\nu c_i &=& 0 & \text{on}& (0,T)\times \Gamma_\text{out}, \\
-d^{\scriptscriptstyle\Sigma}_i\partial_{\nu_\Sigma} c^{\scriptscriptstyle\Sigma}_i &=&0 &\text{on}& (0,T)\times \partial \Sigma, \\
{c_i}|_{t=0}&=&c_{0,i}&\text{in}&\Omega,\\
{c^{\scriptscriptstyle\Sigma}_i}|_{t=0}&=&c^{\scriptscriptstyle\Sigma}_{0,i} & \text{on}& \Sigma.
\end{array}\right.
\end{equation}
In \eqref{eq:cat}, the $d_i,d^{\scriptscriptstyle\Sigma}_i>0$ are given constant diffusivities, $u=u(t,x)$ denotes the velocity,
$r^\text{sorp}_i$ are the sorption rate functions, $r^\text{ch}_i$ the rates of molar mass production
due to surface chemistry, $g^\text{in}_i$ the inflow rates of molar mass, $c_{0,i}$ and $c^{\scriptscriptstyle\Sigma}_{0,i}$ the
initial concentrations in the bulk and on the active surface
and $\nu$ is the outer normal to $\Omega$.

{\bfseries ($\text{A}^\text{vel}$)}
Throughout this paper we assume that the velocity field satisfies
\begin{align}\label{eq:regularity-calss-velocity}
u \in \mathbb{U}^\Omega_p(T):=W^{1,p}((0,T),L^p(\Omega,\mathbb{R}^3))\cap L^p((0,T),W^{2,p}(\Omega,\mathbb{R}^3)),
\end{align} for given time $T>0$. Moreover, we assume
\begin{align*} u\cdot \nu\leq 0 \quad \text{on}\; \Gamma_\text{in}, \qquad u\cdot \nu =0\quad \text{on}\; \Sigma, \qquad u\cdot \nu\geq 0\quad \text{on}\;\Gamma_\text{out},
\end{align*} and $\div u=0$ in the distributional sense.

\subsection*{Examples for sorption and reaction rates}
We give a few examples for sorption and reaction rate functions.

\begin{enumerate}
\item[(S1)] Let $k^\text{ad}_i,k^\text{de}_i>0$ denote adsorption and
desorption rate constants. The the simplest sorption rate is given by the linear {\it Henry law}, i.e.\
\begin{align*}
r^\text{sorp}_{H,i}(c_i,c^{\scriptscriptstyle\Sigma}_i)=k^\text{ad}_i c_i - k^\text{de}_i c^{\scriptscriptstyle\Sigma}_i.
\end{align*}
This law only applies for dilute systems.
\item[(S2)] For moderate concentrations, {\it Langmuir's law} given by
\begin{align*}
r^\text{sorp}_{L,i}(c_i,c^{\scriptscriptstyle\Sigma}_i)=k^\text{ad}_i c_i\left(1-\frac{c^{\scriptscriptstyle\Sigma}_i}{c^{\scriptscriptstyle\Sigma}_{\infty,i}}\right)-k^\text{de}_ic^{\scriptscriptstyle\Sigma}_i
\end{align*} may be employed. Here $c^{\scriptscriptstyle\Sigma}_{\infty,i}>0$ denotes the maximum capacity for adsorption of species $i$.
In an application of our main results, we actually consider a modified version; see Remark \ref{rk:revisiting-examples}, which satisfies all of our assumptions on the sorption rate stated in Section \ref{sec:loc}.
\end{enumerate}

\begin{enumerate}
\item[(R1)] A standard example considers a reversible chemical reaction of type $A+B\rightleftharpoons P$
with $N=3$ components. If mass action kinetics is employed, the mass productions are governed by the rate function
\begin{align*}
r^\text{ch}(c^{\scriptscriptstyle\Sigma})=\begin{pmatrix}
-k^\text{re}(c^{\scriptscriptstyle\Sigma}_1c^{\scriptscriptstyle\Sigma}_2-\kappa c^{\scriptscriptstyle\Sigma}_3)\\
-k^\text{re}(c^{\scriptscriptstyle\Sigma}_1c^{\scriptscriptstyle\Sigma}_2-\kappa c^{\scriptscriptstyle\Sigma}_3)\\
+k^\text{re}(c^{\scriptscriptstyle\Sigma}_1c^{\scriptscriptstyle\Sigma}_2-\kappa c^{\scriptscriptstyle\Sigma}_3) \end{pmatrix}.
\end{align*} Here $k^\text{re}>0$ denotes the rate constant of the
forward reaction, while $\kappa$ is the equilibrium constant for this reaction, determined as the ratio between forward and backward reaction rates.
\end{enumerate}

Due to the nonlinear coupling between bulk and surface in \eqref{eq:cat}, the extension of local and global
existence results from classical bulk reaction-diffusion systems (see \cite{Pie10} for a recent survey)
to the considered advection-diffusion-sorption-reaction system is not straightforward and there are only
few papers dealing with related models.
In \cite{Knabner-Otto}, a similar system but without chemical reactions has been studied. The authors have
shown that an $L^1$-contraction principle holds for the evolution operator and thereby weak solutions are unique. In \cite{Marpeau-Saad}, the case of fast sorption is considered in which local equilibrium between the
adsorbed concentration and the adjacent bulk concentration yields an algebraic relation between these quantities. Surface chemical reactions are not included.
As mentioned above, heterogeneous catalysis processes are often performed in porous media in which case
homogenization is a useful technique to obtain scale-reduced models. The mathematical details of such
a homogenization for periodic porous media have been worked out in \cite{Hornung}.
Another scale-reduced model of heterogeneous catalysis itself has been analyzed in \cite{Bothe:HetCat}
concerning the existence of time-periodic solutions.

\section{Main Results}\seclabel{results}

The main results of this paper are the local-in-time existence of a unique nonnegative strong $L^p$-solution
and the global-in-time existence of a unique nonnegative strong $L^2$-solution.
In all sections, $(0,T)$ denotes a finite time interval. For a bounded simply connected $C^2$-domain $A\subset \mathbb{R}^2$, such that $\partial A$ is a closed regular $C^2$-curve, let $\Omega:=A\times(-h,h)$ be a finite cylinder in $\mathbb{R}^3$.
The boundary of $\Omega$ decomposes into three parts, which are
the lateral surface $\Sigma=\partial A\times(-h,h)$, an inflow area
$\Gamma_\text{in}=A\times\{-h\}$ and an outflow area $\Gamma_\text{out}=A\times\{h\}$, where
$\Gamma_\text{in},\Gamma_\text{out}$ are the bottom and the top of the cylinder $\Omega$, respectively.
In particular, we have $\partial\Omega=\Gamma_\text{in}\cup\overline{\Sigma}\cup\Gamma_\text{out}$.
The local existence result reads as follows. For the assumptions imposed on the sorption and reaction rates see the beginning of Section \ref{sec:loc}.

\begin{theorem}\thmlabel{Local-WP}(Local existence) Let $J=(0,T')\subset \mathbb{R}$ and $5/3<p<\infty$ with $p\neq 3$.
Suppose $u$ satisfies ($\text{A}^\text{vel}$), $r^\text{sorp}$ satisfies ($A_\text{F}^\text{sorp}$), ($A_\text{M}^\text{sorp}$), ($A_\text{B}^\text{sorp}$)
and $r^\text{ch}$ fulfills ($A^\text{ch}_\text{F}$), ($A^\text{ch}_\text{N}$), ($A^\text{ch}_\text{P}$).
Then for every set of data
\begin{align*}
g^\text{in}_i &\in W^{1/2-1/2p}_p(J,L^p(\Gamma_\text{in}))\cap L^p(J,W^{1-1/p}_p(\Gamma_\text{in})),\\
c_{0,i} &\in W^{2-2/p}_p(\Omega),\\
c^{\scriptscriptstyle\Sigma}_{0,i} &\in W^{2-2/p}_p(\Sigma),
\end{align*}
which, if $p>3$, satisfies the compatibility conditions
\begin{align*}
\begin{array}{rcll}
c_{0,i} u(0) \cdot \nu -d_i\partial_\nu c_{0,i}&=&g^\text{in}_i (0)&\text{on}\;\Gamma_{\textrm{in}},\\
-d_i\partial_\nu c_{0,i}&=&r^\text{sorp}_i(c_{0,i},c^{\scriptscriptstyle\Sigma}_{0,i})& \text{on}\;\Sigma,\\
-d_i\partial_\nu c_{0,i}&=&0&\text{on}\;\Gamma_{\textrm{out}},\\
-d^{\scriptscriptstyle\Sigma}_i\partial_{\nu_\Sigma} c^{\scriptscriptstyle\Sigma}_{0,i} &=&0 &\text{on}\; \partial \Sigma,
\end{array}
\end{align*}
there exists a $T^\ast\in(0,T')$ and a unique strong solution $(c_i,c^{\scriptscriptstyle\Sigma}_i)$ of (\ref{eq:cat})
satisfying
\begin{align*}
c_i&\in W^{1,p}((0,T),L^p(\Omega))\cap L^p((0,T),W^{2,p}(\Omega)),\\
c^{\scriptscriptstyle\Sigma}_i &\in W^{1,p}((0,T),L^p(\Sigma))\cap L^p((0,T),W^{2,p}(\Sigma)),
\end{align*} for all $T\in (0,T^*)$. If, in addition, $g^\text{in}_i\leq 0$ on $\Gamma_\text{in}$, $c_{0,i}\geq
0$ in $\Omega$ and $c^{\scriptscriptstyle\Sigma}_{0,i}\geq 0$ on $\Sigma$, then $c_i$ and $c^{\scriptscriptstyle\Sigma}_i $ are nonnegative.
\end{theorem}

\begin{remark} If $c_i,c^{\scriptscriptstyle\Sigma}_i$ are not continuous functions, then $c_i\geq 0$ has to be understood in the a.e.\ sense with respect to Lebesgue measure on $\Omega$, and $c^{\scriptscriptstyle\Sigma}_i\geq 0$ with respect to the surface measure on $\Sigma$.
\end{remark}

The global existence is given by
\begin{theorem}\thmlabel{Global-WP} (Global existence)
	Let the assumptions of \Thmref{Local-WP} be satisfied for $p=2$
	and $T'=\infty$. Additionally, assume that $r^\text{ch}$ fulfills
	($\textrm{A}^{\textrm{ch}}_{\textrm{S}}$)
	(see (\ref{eq:A^ch_S-structure-of-reaction}))
	and that $g^\text{in}_i\in BC((0,T)\times\Gamma_{\text{in}})$, $c_{0,i}\in BC(\Omega)$, $c^{\scriptscriptstyle\Sigma}_{0,i}\in
	BC(\Sigma)$,
	and that $-g^\text{in}_i$,
	$c_{0,i}$, and $c^{\scriptscriptstyle\Sigma}_{0,i}$ are nonnegative.
	Then the local solution $(c_i,c^{\scriptscriptstyle\Sigma}_i)$
	extends to a global solution of (\ref{eq:cat}), i.e., for $p=2$ the
	assertions of \Thmref{Local-WP} hold for every finite $T>0$.
\end{theorem}

\section{Notation}\seclabel{notation}
Let $X$ be a Banach space and let $A:D(A)\subset X\rightarrow X$ denote a linear and densely defined operator. Let $Y$ be another Banach space. We denote the space of bounded linear operators mapping from $X$ to $Y$ by $\mathscr{L}(X,Y)$.

All appearing constants, e.g.\ $C,M>0$ denote generic constants which may
vary from line to line, as long as it is not explicitly stated otherwise.

When working in time-space sets we make use of the notation $\Omega_T:=(0,T)\times\Omega$ and $\Sigma_T:=(0,T)\times\Sigma$ for a finite $T>0$.
When dealing with half infinite cylinders we write,
$\Omega_{(-\infty,h)}:=A\times(-\infty,h)\subset \mathbb{R}^3$ or $\Sigma_{(-\infty,h)}:=\partial A\times (-\infty,h)$ for its lateral surface.

For a Banach space $X$, a domain $G\subset \mathbb{R}^n$, $m\in\mathbb{N}_0:=\mathbb{N}\cup\{0\}$, and
$s\in(0,\infty)\setminus\mathbb{N}$, $W^{m,p}(G,X)$ denotes the usual $X$-valued Sobolev space and
$W^s_p(G,X)$ denotes the $X$-valued Sobolev-Slobodeckij space.
The norm of $X$ will be denoted by
$\Vert\cdot\Vert_X$. We also set $H^k:=W^{k,2}$.
In the same manner we employ $C^m$ for $m$-times continuously differentiable functions, and $BC^m$ for those with bounded derivatives up to order $m\in \mathbb{N}_0$.

We denote by $f^+:=\max\{0,f\}$, $f^-:=-\min\{0,f\}$ the positive and negative part of a function $f$.
Moreover, we use the superscripts $\pm$ to denote sets of
functions, whose elements are nonnegative or nonpositive; e.g., we write
$ L^\infty(\Omega)^+$ for functions which admit a bounded essential
supremum on $\Omega$ and which are nonnegative a.e.\ in $\Omega$. With
corresponding meaning we employ, e.g., $L^\infty(\Omega)^-$.

Throughout this work, let
$\nabla_\Sigma u:=(\nabla u)|_{\partial\Omega}
-\nu(\nu\cdot (\nabla u)|_{\partial\Omega})$
denote the surface gradient and let $\Delta_\Sigma u=\nabla_\Sigma \cdot \nabla_\Sigma u$ denote the Laplace-Beltrami operator on $\Sigma$.

\subsection{Maximal Regularity Spaces}
For $1<p<\infty$, we employ the following maximal regularity
spaces.
The solution spaces for the unknown functions $c_i,c^\Sigma_i$ are given by
\begin{align*}
\mathbb{E}^\Omega_p(T)&:=W^{1,p}((0,T),L^p(\Omega))\cap L^p((0,T),W^{2,p}(\Omega)),\\
\mathbb{E}^\Sigma_p(T)&:= W^{1,p}((0,T),L^p(\Sigma))\cap L^p((0,T),W^{2,p}(\Sigma)).
\end{align*}
For the data spaces we first establish appropriate regularity classes.
Then we give necessary compatibility conditions in order to guarantee well-posedness. We set
\begin{align*}
\mathbb{F}^\Omega_p(T)&:=L^p((0,T)\times \Omega),\\
\mathbb{F}^\Sigma_p(T)&:=L^p((0,T)\times \Sigma),\\
\mathbb{G}^\text{in}_p(T)&:=W^{1/2-1/2p}_p((0,T),L^p(\Gamma_\text{in})) \cap L^p((0,T),W^{1-1/p}_p(\Gamma_\text{in})),\\
\mathbb{G}^\Sigma_p(T)&:= W^{1/2-1/2p}_p((0,T),L^p(\Sigma)) \cap L^p((0,T),W^{1-1/p}_p(\Sigma)),\\
\mathbb{G}^\text{out}_p(T)&:= W^{1/2-1/2p}_p((0,T),L^p(\Gamma_\text{out})) \cap L^p((0,T),W^{1-1/p}_p(\Gamma_\text{out})),\\
\mathbb{I}^\Omega_p&:=W^{2-2/p}_p(\Omega), \\
\mathbb{I}^\Sigma_p&:= W^{2-2/p}_p(\Sigma).
\end{align*} We define the tupel data space for the heterogeneous catalysis equations without initial data through
\begin{align*}
\mathbb{F}^{\Omega,\Sigma}_p(T):=\mathbb{F}^\Omega_p(T)\times\mathbb{F}^\Sigma_p(T)\times\mathbb{G}^\text{in}_p(T)\times\mathbb{G}^\Sigma_p(T)\times \mathbb{G}^\text{out}_p(T)\times\{0\
\end{align*} and the tupel data space with initial spaces through
\begin{align*}
\mathbb{F}^{\Omega,\Sigma}_{p,I}(T):= \mathbb{F}^{\Omega,\Sigma}_p(T) \times\mathbb{I}^\Omega_p\times\mathbb{I}^\Sigma_p.
\end{align*} In some statements we also employ the Dirichlet trace space on $\Sigma$, which is given by
\begin{align*}
\mathbb H^\Sigma_p(T)&:= W^{1-1/2p}_p((0,T),L^p(\Sigma)) \cap L^p((0,T),W^{2-1/p}_p(\Sigma)).
\end{align*}

Moreover, we need subspaces of functions having zero time trace. For instance, the $X$-valued
Sobolev-Slobodeckij space with zero time trace is defined as
\begin{equation}\label{eq:zerotime} \,_0W^s_p((0,T),X):=\{c_i\in W^s_p((0,T),X):\; {c_i}|_{t=0}=0\}
\end{equation} for $s\in(0,\infty)\setminus\mathbb{N}$ and $sp>1$. We use this notation for all
appearing spaces whenever zero time trace makes sense and write, e.g., $\,_0\mathbb{E}^\Omega_p(T)$, $_0\mathbb{E}^\Sigma_p(T)$ etc.

\section{Linear Equations}\label{sec:linear_equations}

We discuss a suitable linearization of the heterogeneous catalysis equations as it is given below in (\ref{eq:withPerturbation}) and show maximal regularity by
means of cylindrical $L^p$-theory, the surjectivity of the Neumann trace
operator, and a perturbation argument. For given data
\begin{equation*} (f_i,f^{\scriptscriptstyle\Sigma}_i,g^\text{in}_i,g^{\scriptscriptstyle\Sigma}_i,g^\text{out}_i,0,c_{0,i},c^{\scriptscriptstyle\Sigma}_{0,i})\in\mathbb{F}^{\Omega,\Sigma}_{p,I}(T),
\end{equation*} we consider the linear system:
\begin{equation}\label{eq:withPerturbation}
\left\{\begin{array}{rclcl}
\partial_t c_i+ u\cdot\nabla c_i- d_i\Delta c_i &=&f_i& \text{in}& (0,T)\times\Omega, \\
\partial_t c^{\scriptscriptstyle\Sigma}_i - d^{\scriptscriptstyle\Sigma}_i \Delta_\Sigma c^{\scriptscriptstyle\Sigma}_i &=&f^{\scriptscriptstyle\Sigma}_i & \text{on}& (0,T)\times \Sigma, \\
c_i u\cdot\nu-d_i\partial_\nu c_i&=&g^\text{in}_i& \text{on}& (0,T)\times \Gamma_\text{in}, \\
-d_i\partial_\nu c_i &=&g^{\scriptscriptstyle\Sigma}_i& \text{on}& (0,T)\times \Sigma,\\
-d_i\partial_\nu c_i &=& g^\text{out}_i & \text{on}& (0,T)\times \Gamma_\text{out}, \\
-d^{\scriptscriptstyle\Sigma}_i\partial_{\nu_\Sigma} c^{\scriptscriptstyle\Sigma}_i &=&0 &\text{on}& (0,T)\times \partial \Sigma, \\
c_i(0)&=&c_{0,i}   &\text{in}&\Omega,   \\
c^{\scriptscriptstyle\Sigma}_i(0)&=&c^{\scriptscriptstyle\Sigma}_{0,i} & \text{in}& \Sigma.
\end{array}\right.
\end{equation} Our purpose is to solve (\ref{eq:withPerturbation}) for the unknown concentrations $c_i, c^{\scriptscriptstyle\Sigma}_i$.
The main result of this section is given by

\begin{proposition}{\label{prop:maxRegPerturbedSystem}} Let $5/3<p<\infty$ with  $p\neq 3$ and let $T>0$ be finite. Suppose the velocity field $u$ satisfies assumption ($A^\text{vel}$). Then (\ref{eq:withPerturbation}) admits a unique solution
\begin{equation*} (c_i,c^{\scriptscriptstyle\Sigma}_i)\in \mathbb{E}^\Omega_p(T)\times \mathbb{E}^\Sigma_p(T)
\end{equation*} if and only if the data satisfy the regularity condition
\begin{equation*}
(f_i,f^{\scriptscriptstyle\Sigma}_i,g^\text{in}_i, g^{\scriptscriptstyle\Sigma}_i, g^\text{out}_i,0,c_{0,i},c^{\scriptscriptstyle\Sigma}_{0,i})\in\mathbb{F}^{\Omega,\Sigma}_{p,I}(T),
\end{equation*} and in case of $p>3$ the compatibility conditions
\begin{align*}\begin{array}{rcll}
(u|_{t=0}\cdot\nu)c_{0,i}-d_i\partial_\nu c_{0,i}&=&g^\text{in}_i|_{t=0} &\text{on}\;\Gamma_\text{in},\\
-d_i\partial_\nu c_{0,i}&=&g^{\scriptscriptstyle\Sigma}_i|_{t=0}& \text{on}\;\Sigma,\\
-d_i\partial_\nu c_{0,i}&=&g^\text{out}_i|_{t=0}& \text{on}\;\Gamma_\text{out},\\
-d^{\scriptscriptstyle\Sigma}_i\partial_{\nu_\Sigma} c^{\scriptscriptstyle\Sigma}_{0,i}&=&0 &\text{on}\;\partial\Sigma.\end{array}
\end{align*} Additionally, the corresponding solution operator
${_0{\mathcal S}}_T$ with respect to zero time trace spaces satisfies
\begin{equation*}
\Vert
{_0{\mathcal S}}_\tau\Vert_{\mathscr{L}\left({{_0\mathbb{F}}^{\Omega,\Sigma}_{p}(\tau)}^N,\,{_0\mathbb{E}^\Omega_p(\tau)}^N\times{_0\mathbb{E}^\Sigma_p(\tau)}^N\right)}
\leq M
\end{equation*} with a constant $M>0$ independent of $\tau<T$.
\end{proposition}

Proposition~\ref{prop:maxRegPerturbedSystem} shows that the map ${\mathcal L} : \mathbb{E} \to \mathbb{F}$,
where $\mathbb{E}= \mathbb{E}^\Omega_p(T)^N\times \mathbb{E}^\Sigma_p(T)^N$ and $\mathbb{F}$ is the subspace of $\mathbb{F}^{\Omega,\Sigma}_{p,I}(T)^N$
containing all elements which satisfy the required compatibility conditions,
is an isomorphism between the Banach spaces $\mathbb{E}$ and $\mathbb{F}$.
In this context, one speaks of {\em maximal regularity} of problem (\ref{eq:withPerturbation}),
or of the operator ${\mathcal L}:\mathbb{E}\to\mathbb{F}$, acting on $(c_i,c^{\scriptscriptstyle\Sigma}_i)$.

{\it Plan of the proof:} System (\ref{eq:withPerturbation}) decomposes into two systems: One for the bulk concentrations $c_i$ in $\Omega$ and one for the surface concentrations $c^{\scriptscriptstyle\Sigma}_i$ on $\Sigma$. In the first step we
neglect the velocity terms $(u\cdot \nabla ) c_i$ and $(u\cdot\nu)c_i$ -- playing the role of perturbation terms -- and consider only homogeneous boundary data, i.e.\ we start with
\begin{equation}\label{eq1.2}
\left\{\begin{array}{rclcl}
\partial_t c_i- d_i\Delta c_i &=&f_i& \text{in}& (0,T)\times\Omega, \\
-d_i\partial_\nu c_i &=&0& \text{on}&  (0,T)\times\partial\Omega, \\
c_i(0)&=&c_{0,i}& \text{in}& \Omega, \end{array} \right. \quad (i=1,...,N)
\end{equation} and
\begin{equation}\label{eq1.3}
\left\{\begin{array}{rclcl}
\partial_t c^{\scriptscriptstyle\Sigma}_i - d^{\scriptscriptstyle\Sigma}_i \Delta_\Sigma c^{\scriptscriptstyle\Sigma}_i &=&f^{\scriptscriptstyle\Sigma}_i & \text{on}& (0,T)\times \Sigma, \\
-d^{\scriptscriptstyle\Sigma}_i\partial_{\nu_\Sigma} c^{\scriptscriptstyle\Sigma}_i &=&0 &\text{on}& (0,T)\times \partial \Sigma, \\
c^{\scriptscriptstyle\Sigma}_i(0)&=&c^{\scriptscriptstyle\Sigma}_{0,i}&\text{on}& \Sigma. \end{array}\right. \quad (i=1,...,N)
\end{equation} We proceed as follows: We solve (\ref{eq1.2}) and (\ref{eq1.3}) separately via cylindrical $L^p$-theory. For more information on this topic, see \cite{NS2011}, \cite{NS2012}, \cite{Nau2012} and \cite{Nau2013}. Then a symmetric extension in axial direction of $\Omega$ yields the surjectivity of the Neumann trace operator and,
consequently, the solvability of the inhomogeneous initial boundary value problem. By  perturbation arguments the obtained result carries over to (\ref{eq:withPerturbation}).

\subsection{Maximal Regularity of the Laplacian}\label{sec1}
We define
\begin{align*}
A_i:=-d_i \Delta, \quad A_i:D(A_i)\subset L^p(\Omega)\rightarrow L^p(\Omega), \\
D(A_i):=\{c_i\in W^{2,p}(\Omega): \; -d_i\partial_\nu c_i=0\;\text{on}\; \partial\Omega\},\quad(i=1,...,N)
\end{align*} and, analogously,
\begin{align*}
&A^{\scriptscriptstyle\Sigma}_i:=-d^{\scriptscriptstyle\Sigma}_i\Delta_\Sigma, \quad A^{\scriptscriptstyle\Sigma}_i:D(A^{\scriptscriptstyle\Sigma}_i)\subset L^p(\Sigma)\rightarrow L^p(\Sigma), \\
&D(A^{\scriptscriptstyle\Sigma}_i):=\{c^{\scriptscriptstyle\Sigma}_i\in W^{2,p}(\Sigma): \; -d^{\scriptscriptstyle\Sigma}_i\partial_{\nu_\Sigma}c^{\scriptscriptstyle\Sigma}_i=0\;\text{on}\;\partial\Sigma\}\quad (i=1,...,N).
\end{align*}
In this section we will show that $A_i$ admits
a bounded ${\mathcal H}^\infty$-calculus and that $A^{\scriptscriptstyle\Sigma}_i$ is ${\mathcal R}$-sectorial with angles strictly less than $\pi/2$ which implies the
desired maximal regularity, cf.\ \cite{Denk-Hieber-Pruess:Maximal-Regularity}, \cite{Kunstmann-Weiss}, \cite{Weis-Maximal-Regularity-R-Sectoriality}. For $A_i$
we directly apply \cite[Theorem 4.1]{Nau2013}.
To this end, we employ the following cylindrical decomposition.
We set $V_1:=A\subset\mathbb{R}^2$ and $V_2:=(-h,h)\subset\mathbb{R}$ for the cylinder $\Omega$ with height $2h>0$.
\begin{itemize}
\item
We first consider the problem which results on the sections for fixed $x_3$, yielding the following problems on
$V_1=A\subset \mathbb{R}^2$:
\begin{align*}
A_{i,1}:D(A_{i,1})\subset L^p(V_1)\rightarrow L^p(V_1), \quad A_{i,1}c_i:=-d_i(\partial_{x_1}^2+\partial_{x_2}^2)c_i,\\
D(A_{i,1}):=\{c_i\in W^{2,p}(V_1):\, B_{i,1} c_i:=-d_i\partial_\nu c_i
=0\; \text{on}\; \partial V_1\}.
\end{align*}
\item The direction along the cylinder leads to an operator on an interval $V_2$:
\begin{align*}
A_{i,2}:D(A_{i,2})\subset L^p(V_2)\rightarrow L^p(V_2), \quad A_{i,2}c_i:=-d_i \partial_{x_3}^2 c_i,\\
D(A_{i,2}):=\{ c_i\in W^{2,p}(V_2):\; B_{i,2} c_i
:=-d_i\partial_\nu c_i=0\;\text{on}\; \partial V_2\}.
\end{align*}
\end{itemize}
It is a well-known fact that both operators, $A_{i,1}$ and $A_{i,2}$
admit a bounded ${\mathcal H}^\infty$-calculus with zero ${\mathcal H}^\infty$-angle. For these results and an introduction of the ${\mathcal H}^\infty$-calculus we refer to e.g.\ \cite{Denk-Hieber-Pruess:Maximal-Regularity} and \cite{Kunstmann-Weiss}.
Thus, we are exactly in the setting of \cite{Nau2013} in the case of the
strong Neumann-Laplacian given on both intersections  $V_1,V_2$.
Therefore \cite[Theorem 4.1 a)]{Nau2013} yields that $A_i+\delta$
for some $\delta>0$ admits a bounded ${\mathcal H}^\infty$-calculus on
$L^p(\Omega)$ with ${\mathcal H}^\infty$-angle
$\phi^{\infty}_{A_i+\delta}<\frac{\pi}{2}$.
This implies maximal regularity for (\ref{eq1.2}) on finite intervals
$(0,T)$, see e.g.\ \cite{Denk-Hieber-Pruess:Maximal-Regularity}.

We continue with the discussion of $A^{\scriptscriptstyle\Sigma}_i$.
Here we first employ a parametrization of the lateral surface and afterwards apply
a result from \cite{Nau2012}.
By the assumption on the cross section of $\Omega$, $\partial A$ constitutes a closed regular $C^2$-curve. Hence we may choose a parametrization $\psi^\Sigma:(0,2\pi)\rightarrow\mathbb{R}^2$ of $\partial A$. Let
\begin{equation*}\label{eq1.9}
\Psi^\Sigma:L^p(\partial V_1)\rightarrow L^p((0,2\pi)), \quad \Psi^\Sigma c^{\scriptscriptstyle\Sigma}_i:=c^{\scriptscriptstyle\Sigma}_i\circ \psi^\Sigma
\end{equation*} denote the induced pull-back. Then the $\Psi^\Sigma$-transformed Laplace-Beltrami operator $\Delta_{\Sigma,1}$ on $\partial V_1$ is a second-order elliptic differential operator on $L^p((0,2\pi))$ subject to periodic boundary conditions:
\begin{align*}
A^{\scriptscriptstyle\Sigma}_{i,1}&:D(A^{\scriptscriptstyle\Sigma}_{i,1})\subset L^p((0,2\pi))\rightarrow L^p((0,2\pi)), \quad
A^{\scriptscriptstyle\Sigma}_{i,1} c^{\scriptscriptstyle\Sigma}_i:= \Psi^\Sigma(-d^{\scriptscriptstyle\Sigma}_i\Delta_{\Sigma,1}c^{\scriptscriptstyle\Sigma}_i),\\
D(A^{\scriptscriptstyle\Sigma}_{i,1})&=\{c^{\scriptscriptstyle\Sigma}_i\in W^{2,p}((0,2\pi)): \;
\partial_\varphi^j c^{\scriptscriptstyle\Sigma}_i|_{\varphi=0}=
\partial_\varphi^jc^{\scriptscriptstyle\Sigma}_i|_{\varphi=2\pi}\;\;(j=0,1) \}.
\end{align*}
We have
\begin{lemma}\label{lm1.3} Let the pull-back $\Psi^\Sigma$ be given as above, then
\begin{enumerate}
\item[(1)] $\Psi^\Sigma$ is an isomorphism.
\item[(2)] $\Psi^\Sigma: D(\Delta_{\Sigma,1})\rightarrow D(A^{\scriptscriptstyle\Sigma}_{i,1})$
is an isomorphism.
\end{enumerate}
\end{lemma}

Analogously to $A_i$, we resolve the transformation of $A^{\scriptscriptstyle\Sigma}_i$ into two
cylindrical parts: $A^{\scriptscriptstyle\Sigma}_{i,1}$ defined as above and $A^{\scriptscriptstyle\Sigma}_{i,2}:=A_{i,2}$.

Applying this time \cite[Theorem 8.10]{Nau2012}
we infer that there is a shift $\delta>0$ such that
$\delta+ A^{\scriptscriptstyle\Sigma}_i$ is ${\mathcal R}$-sectorial with $\phi^{\mathcal R}_{A^{\scriptscriptstyle\Sigma}_i+\delta}<\frac{\pi}{2}$. This implies maximal regularity of (\ref{eq1.3}) on finite intervals $(0,T)$ by \cite[Theorem 4.2]{Weis-Maximal-Regularity-R-Sectoriality}.

We summarize the results of Subsection~\ref{sec1} in

\begin{lemma}\label{lm:Homogeneous-Equations} Let $1<p<\infty$ with $p\neq 3$ and $T>0$ be given.
\begin{enumerate}
\item[a)] System (\ref{eq1.2}) admits a unique solution $c_i\in \mathbb{E}^\Omega_p(T)$ if and only if
the data satisfies the regularity condition $f_i\in\mathbb{F}^\Omega_p(T)$, $c_{0,i}\in\mathbb{I}^\Omega_p$ and, in case of $p>3$, the compatibility condition $-d_i\partial_\nu c_{0,i}=0$ on $\partial\Omega$.
\item[b)] System (\ref{eq1.3}) admits a unique solution $c^{\scriptscriptstyle\Sigma}_i\in \mathbb{E}^\Sigma_p(T)$ if and only if
the data satisfies the regularity condition $f^{\scriptscriptstyle\Sigma}_i\in\mathbb{F}^\Sigma_p(T)$, $c^{\scriptscriptstyle\Sigma}_{0,i}\in\mathbb{I}^\Sigma_p$ and, in case of $p>3$, the compatibility condition $-d^{\scriptscriptstyle\Sigma}_i\partial_{\nu_{\scriptscriptstyle\Sigma}} c_{0,i}=0$ on $\partial\Sigma$.
\end{enumerate}
\end{lemma}

\subsection{Inhomogeneous Neumann Boundary Conditions}\label{sec4}

We turn to the discussion of inhomogeneous boundary values.
We show surjectivity of the trace operator
which leads to
the solvability of the corresponding inhomogeneous boundary value problem.

\begin{lemma}\label{lm:Surjectivity-of-Trace} Let $1<p<\infty$ with
$p\neq 3$ and let $T>0$. Then the Neumann trace operator
\begin{align*}
\gamma_1:\mathbb{E}^\Omega_p(T)&\rightarrow \mathbb{G}^\text{in}_p(T)\times\mathbb{G}^\Sigma_p(T)\times\mathbb{G}^\text{out}_p(T)\\
&c\mapsto (\partial_\nu c|_{\Gamma_\text{out}}, \partial_\nu c|_\Sigma,\partial_\nu c|_{\Gamma_\text{in}})
\end{align*} is a retraction.
\end{lemma}

\begin{proof} Let $g^\text{in}\in \mathbb{G}^\text{in}_p(T),\; g^{\scriptscriptstyle\Sigma}\in\mathbb{G}^\Sigma_p(T)$ and $g^\text{out}\in\mathbb{G}^\text{out}_p(T)$.
We show that there exists a $c\in\mathbb{E}^\Omega_p(T)$ such that
\begin{align}
\partial_\nu c&=g^\text{in}\quad \text{on}\;\Gamma_\text{in},\label{eq:neum1}\\
\partial_\nu c&=g^{\scriptscriptstyle\Sigma}\quad\; \text{on}\;\Sigma,\label{eq:neum2}\\
\partial_\nu c&=g^\text{out}\;\;\; \text{on}\;\Gamma_\text{out}. \label{eq:neum3}
\end{align} Note that we skip the index $i$ for better readability.
As before we work with $\Omega=A\times(-h,h)$, such that $\Gamma_\text{in}=A\times\{-h\}$ and $\Gamma_\text{out}=A\times\{h\}$ in the present section. Let us define the halfspaces
\begin{equation*}
H_{-h}:=\mathbb{R}^2\times(-h,\infty),\quad H_{h}:=\mathbb{R}^2\times(-\infty,h).
\end{equation*} We proceed in two steps. \\
{\bfseries Step 1.} Due to \cite[Chapter 5]{Adams-Fournier:Sobolev-Spaces} there is an extension of $g^\text{in}$ to
\begin{equation*}
\tilde{g}^\text{in} \in \mathbb{G}^{\partial H_{-h}}_p(T)
\end{equation*} and of $g^\text{out}$ to
\begin{equation*}
\tilde{g}^\text{out} \in \mathbb{G}^{\partial H_{h}}_p(T).
\end{equation*} We have
\begin{equation*}
\tilde{g}^\text{in}(0) \in W^{1-3/p}_p(\partial H_{-h}), \quad \tilde{g}^\text{out}(0) \in W^{1-3/p}_p(\partial H_{h}).
\end{equation*} In case that $p>3$ we may choose a $\tilde{c}^\text{in}_0\in W^{2-2/p}_p(H_{-h})$ and $\tilde{c}^\text{out}_0\in W^{2-2/p}_p(H_{h})$ such that
\begin{align*}
\partial_\nu \tilde{c}^\text{in}_0&=\tilde{g}^\text{in}(0)\quad \text{on}\;\partial H_{-h}, \\
\partial_\nu \tilde{c}^\text{out}_0&=\tilde{g}^\text{out}(0)\quad \text{on}\;\partial H_{h},
\end{align*}
since $\partial_\nu : W^{2-2/p}_p(H_{-h}) \to W^{1-3/p}_p(\partial H_{-h})$ is
a retraction due to \cite[Theorem 2]{Marschall87}.
In case $p<3$ we let
$\tilde{c}^\text{in}_0\in W^{2-2/p}_p(H_{-h})$ and $\tilde{c}^\text{out}_0\in W^{2-2/p}_p(H_{h})$ be arbitrary. Due to \cite[Proposition 4.6]{Denk-Hieber-Pruess:Maximal-Regularity-Inhomogeneous} we can solve the parabolic problem
\begin{equation*}
\left\{\begin{array}{rllll}\partial_t v^\text{in}-\Delta v^\text{in}&=&0&\text{in}&(0,T)\times H_{-h} \\
\partial_\nu v^\text{in} &=&\tilde{g}^\text{in}&\text{on}&(0,T)\times\partial H_{-h} \\
v^\text{in}(0)&=&\tilde{c}^\text{in}_0&\text{in}& H_{-h} \end{array}\right.
\end{equation*} with a unique
\begin{equation*}
v^\text{in} \in \mathbb{E}^{H_{-h}}_p(T)
\end{equation*} and analogously we solve
\begin{equation*}
\left\{\begin{array}{rllll}\partial_t v^\text{out}-\Delta v^\text{out}&=&0&\text{in}&(0,T)\times H_{h} \\
\partial_\nu v^\text{out} &=&\tilde{g}^\text{out}&\text{on}&(0,T)\times\partial H_{h} \\
v^\text{out}(0)&=&\tilde{c}^\text{out}_0&\text{in}& H_{h} \end{array}\right.
\end{equation*} with a unique
\begin{equation*}
v^\text{out} \in \mathbb{E}^{H_h}_p(T).
\end{equation*}
Let the cut-off function $\zeta \in C^\infty(\mathbb{R},[0,1])$ satisfy
\begin{equation}\label{eq:retraction-cut-off-zeta}
\zeta|_{(-\infty,-h/3)}=1, \quad \zeta|_{(h/3,\infty)}=0.
\end{equation} Then the convex combination
\begin{equation*}
v:=\zeta v^\text{in}|_{\Omega} + (1-\zeta) v^\text{out}|_\Omega \in \mathbb{E}^\Omega_p(T)
\end{equation*} fulfills the boundary conditions on top and bottom of $\Omega$ by construction. \\
{\bfseries Step 2.} It remains to show that there exists a $w\in\mathbb{E}^\Omega_p(T)$ such that
\begin{align*}
\partial_\nu w&=0\qquad\qquad \text{on}\;\Gamma_\text{in},\\
\partial_\nu w&=g^{\scriptscriptstyle\Sigma}-\partial_\nu v\;\; \text{on}\;\Sigma,\\
\partial_\nu w&=0\qquad\qquad \text{on}\;\Gamma_\text{out}.
\end{align*}
This problem can be reduced to an equation on a bounded $C^2$-domain, which works as follows.
To this end, first define $\Omega_{-h}$ as the domain resulting from extending $\Omega$
in some way boundedly and smoothly (at least in the $C^2$-sense) on the top.
We also set $\Sigma_{-h}:=\partial\Omega_{-h}\setminus\overline{\Gamma}_{\text{in}}$.
In a similar manner we define $\Omega_{+h}$ and $\Sigma_{+h}$ by extending
$\Omega$ suitably at the bottom. Then, let $G_\pm$ denote the domains resulting from reflecting $\Omega_{\pm h}$ at $h_\pm$ and set $\Gamma_\pm:=\partial G_\pm$.
For instance, if the cross-section $A$ of $\Omega$ is a circle, we connect
half of a ball to $\Omega$ at $\Gamma_{\text{out}}$ resp.\
$\Gamma_{\text{in}}$. Then $G_\pm$ has the form of a `pill'.
It is clear that this way we always can find a suitable extension such that
$G_\pm$ is of class  $C^2$.

Let $\zeta$ be the cut-off function with the properties given in (\ref{eq:retraction-cut-off-zeta}). We extend $\zeta(g^{\scriptscriptstyle\Sigma}-\partial_\nu v)$ to
\begin{align*}
\tilde{g}^{\scriptscriptstyle\Sigma}_- &\in \mathbb{G}^{\Sigma_{-h}}_p(T)\\
&=W^{1/2-1/2p}_p((0,T),L^p(\Sigma_{-h}))\cap L^p((0,T),W^{1-1/p}_p(\Sigma_{-h}))
\end{align*}
by $0$ and $\tilde{g}^{\scriptscriptstyle\Sigma}_-$ to
\begin{equation*}
\hat{g}^{\scriptscriptstyle\Sigma}_- \in \mathbb{G}^{\Gamma_-}_p(T)
\end{equation*} by even reflection at $-h$.
Note that the extension by even reflection conserves the regularity
$W^{1-1/p}_p$ in axial direction.
Next, note that $\hat{g}^{\scriptscriptstyle\Sigma}_-(0)\in W^{1-1/3p}_p(\Gamma_-)$, if $p>3$. Hence we can use \cite[Theorem 2]{Marschall87} to obtain $\hat{c}^{\scriptscriptstyle\Sigma}_\text{in} \in W^{2-2/p}_p(G_-)$ with
\begin{equation*}
\partial_\nu \hat{c}^{\scriptscriptstyle\Sigma}_\text{in} =\hat{g}^{\scriptscriptstyle\Sigma}_-(0)\quad  \text{on}\;\Gamma_-,
\end{equation*} if $p>3$.
Here we can assume w.l.o.g.\ that $\hat{c}^{\scriptscriptstyle\Sigma}_\text{in}$ is even with respect to $\Gamma_\text{in}$.
We solve the problem
\begin{equation}\label{eq:reflection-w^in}
\left\{\begin{array}{rllll}\partial_t w^\text{in}-\Delta w^\text{in}&=&0&\text{in}&(0,T)\times G_- \\
\partial_\nu w^\text{in} &=&\hat{g}^{\scriptscriptstyle\Sigma}_-&\text{on}&(0,T)\times\Gamma_- \\
w^\text{in}(0)&=&\hat{c}^{\scriptscriptstyle\Sigma}_\text{in}&\text{in}& G_- \end{array}\right.
\end{equation} by \cite[Theorem 2.1]{Denk-Hieber-Pruess:Maximal-Regularity-Inhomogeneous} to obtain
\begin{equation*}
w^\text{in} \in \mathbb{E}^{G_-}_p(T).
\end{equation*}
Defining
\begin{equation*} \tilde{w}^\text{in}:=w^\text{in}|_{\Omega_{-h}}
\in  \mathbb{E}^{\Omega_{-h}}_p(T),
\end{equation*} we have $\partial_\nu \tilde{w}^\text{in}=\tilde{g}^{\scriptscriptstyle\Sigma}_-$ on $\Sigma_{-h}$. Since $\hat{g}^{\scriptscriptstyle\Sigma}_-$ and $\hat{c}^{\scriptscriptstyle\Sigma}_\text{in}$ are even in axial direction we have $\partial_\nu \tilde{w}^\text{in}=0$ on $\Gamma_\text{in}$.
Analogously we proceed with $\Gamma_\text{out}$. Here we extend
$(1-\zeta)(g^{\scriptscriptstyle\Sigma}-\partial_\nu v)$ to obtain a
\begin{equation*}
\tilde{w}^\text{out}\in \mathbb{E}^{\Omega_{+h}}_p(T)
\end{equation*} with
$\partial_\nu \tilde{w}^\text{out}=(1-\zeta)(g^{\scriptscriptstyle\Sigma}-\partial_\nu v)$ on $\Sigma$ and $\partial_\nu \tilde{w}^\text{out}=0$ on $\Gamma_\text{out}$.
Let the cut-off functions \\
$\tilde{\zeta}_1,\tilde{\zeta}_2\in C^\infty(\mathbb{R},[0,1])$ satisfy
\begin{equation*}
{\tilde{\zeta}_1}=
\left\{\begin{array}{lll}
	1&: &(-\infty,h/2)\\
	0&:&(2h/3,\infty)
\end{array}\right., \qquad {\tilde{\zeta}_2}=
\left\{\begin{array}{lll}
	0&:&(-\infty,-2h/3)\\
	1&:&(-h/2,\infty)
\end{array}\right..
\end{equation*} Then the combination
\begin{equation*}
w:=\tilde{\zeta}_1 \tilde{w}^\text{in}|_{\Omega} + \tilde{\zeta}_2 \tilde{w}^\text{out}|_{\Omega} \in \mathbb{E}^\Omega_p(T)
\end{equation*} satisfies by construction $\partial_\nu w=0$ on $(0,T)\times\Gamma_\text{in}$ and $\partial_\nu w=0$ on $(0,T)\times\Gamma_\text{out}$. The remaining inhomogeneous boundary condition on $\Sigma$ is satisfied as well, since
\begin{align*}
\partial_\nu w&=\tilde{\zeta}_1\partial_\nu \tilde{w}^\text{in} + \tilde{\zeta}_2\partial_\nu \tilde{w}^\text{out}\\
&=\tilde{\zeta}_1\zeta(g^{\scriptscriptstyle\Sigma}-\partial_\nu v) +\tilde{\zeta}_2(1-\zeta )(g^{\scriptscriptstyle\Sigma}-\partial_\nu v)\\
&=\zeta (g^{\scriptscriptstyle\Sigma}-\partial_\nu v) + (1-\zeta)(g^{\scriptscriptstyle\Sigma}-\partial_\nu v) =g^{\scriptscriptstyle\Sigma}-\partial_\nu v \qquad \text{on}\;(0,T)\times\Sigma.
\end{align*}
Putting together step (i) and (ii) we define $c:=v+w \in \mathbb{E}^\Omega_p(T)$
and obtain that $c$ satisfies (\ref{eq:neum1})-(\ref{eq:neum3}).
Thus, we have proved that
there exists a bounded linear right-inverse to $\gamma_1$ which yields
that the trace operator
$\gamma_1 : \mathbb{E}^\Omega_p(T) \rightarrow \mathbb{G}^\text{in}_p(T)\times\mathbb{G}^\Sigma_p(T)\times\mathbb{G}^\text{out}_p(T)$ in fact is a retraction.
\end{proof}

We turn to the fully inhomogeneous Neumann system
\begin{equation}\label{eq:withoutPerturbation}
\left\{\begin{array}{rclcl}
\partial_t c_i- d_i\Delta c_i &=&f_i& \text{in}& (0,T)\times\Omega, \\
\partial_t c^{\scriptscriptstyle\Sigma}_i - d^{\scriptscriptstyle\Sigma}_i \Delta_\Sigma c^{\scriptscriptstyle\Sigma}_i &=&f^{\scriptscriptstyle\Sigma}_i & \text{on}& (0,T)\times \Sigma, \\
-d_i\partial_\nu c_i&=&g^\text{in}_i& \text{on}& (0,T)\times \Gamma_\text{in}, \\
-d_i\partial_\nu c_i &=&g^{\scriptscriptstyle\Sigma}_i& \text{on}& (0,T)\times \Sigma,\\
-d_i\partial_\nu c_i &=& g^\text{out}_i & \text{on}& (0,T)\times \Gamma_\text{out}, \\
-d^{\scriptscriptstyle\Sigma}_i\partial_{\nu_\Sigma} c^{\scriptscriptstyle\Sigma}_i &=&0 &\text{on}& (0,T)\times \partial \Sigma, \\
c_i(0)&=&c_{0,i}   &\text{in}&\Omega,   \\
c^{\scriptscriptstyle\Sigma}_i(0)&=&c^{\scriptscriptstyle\Sigma}_{0,i} & \text{in}& \Sigma.
\end{array}\right. \quad (i=1,...,N)
\end{equation}

\begin{lemma}\label{lm:withoutPerturbation}
Let $1<p<\infty$ with $p\neq 3$ and let $T>0$ be given. Then (\ref{eq:withoutPerturbation}) admits a unique solution
\begin{equation*} (c_i,c^{\scriptscriptstyle\Sigma}_i)\in \mathbb{E}^\Omega_p(T)\times \mathbb{E}^\Sigma_p(T)
\end{equation*} if and only if the data satisfies the regularity conditions
\begin{align*}
(f_i,f^{\scriptscriptstyle\Sigma}_i,g^\text{in}_i,g^{\scriptscriptstyle\Sigma}_i,g^\text{out}_i,0,c_{0,i},c^{\scriptscriptstyle\Sigma}_{0,i})\in \mathbb{F}^{\Omega,\Sigma}_{p,I}(T)
\end{align*} and, in case of $p>3$, the compatibility conditions
\begin{align*}\begin{array}{rcll}
-d_i\partial_\nu c_{0,i}&=&g^\text{in}_i|_{t=0}& \text{on}\;\Gamma_\text{in},\\
-d_i\partial_\nu c_{0,i}&=&g^{\scriptscriptstyle\Sigma}_i|_{t=0}& \text{on}\;\Sigma,\\
-d_i\partial_\nu c_{0,i}&=&g^\text{out}_i|_{t=0}&  \text{on}\;\Gamma_\text{out},\\
-d^{\scriptscriptstyle\Sigma}_i\partial_{\nu_\Sigma} c^{\scriptscriptstyle\Sigma}_{0,i}&=&0&  \text{on}\;\partial\Sigma.\end{array}
\end{align*} Additionally, the corresponding solution operator ${_0{\mathcal S}}_T$ with respect to homogeneous initial values satisfies
\begin{equation}\label{eq:solution-operator-estimate}
\Vert {_0{\mathcal S}}_\tau\Vert_{\mathscr{L}({{_0\mathbb{F}}^{\Omega,\Sigma}_{p}(\tau)}^N,{_0\mathbb{E}^\Omega_p(\tau)}^N\times{_0\mathbb{E}^\Sigma_p(\tau)}^N)}
\leq M\qquad (0<\tau<T)
\end{equation} for a constant $M>0$ independent of $\tau$.
\end{lemma}

\begin{proof}
By Lemma \ref{lm:Surjectivity-of-Trace} for given $g^\text{in}_i\in\mathbb{G}^\text{in}_p(T)$, $g^{\scriptscriptstyle\Sigma}_i\in\mathbb{G}^\Sigma_p(T)$, $g^\text{out}_i\in\mathbb{G}^\text{out}_p(T)$ there exists a $c^1_i\in\mathbb{E}^\Omega_p(T)$ with
\begin{align*}
-d_i\partial_\nu c_i^1&=g^\text{in}_i \quad \text{on}\;(0,T)\times\Gamma_\text{in}, \\
-d_i\partial_\nu c_i^1&=g^{\scriptscriptstyle\Sigma}_i \quad \text{on}\;(0,T)\times\Sigma, \\
-d_i\partial_\nu c_i^1&=g^\text{out}_i \quad \text{on}\;(0,T)\times\Gamma_\text{out}.
\end{align*}
Secondly, due to Lemma \ref{lm:Homogeneous-Equations} for $f_i\in \mathbb{F}^\Omega_p(T)$, $c_{0,i}\in\mathbb{I}^\Omega_p(T)$ we find a unique $c^2_i\in\mathbb{E}^\Omega_p(T)$ such that
\begin{equation*}
\left\{\begin{array}{rclcl}
\partial_t c_i^2- d_i\Delta c_i^2 &=&f_i-(\partial_t-d_i\Delta)c_i^1& \text{in}& (0,T)\times\Omega, \\
-d_i\partial_\nu c_i^2&=&0& \text{on}& (0,T)\times \Gamma_\text{in}, \\
-d_i\partial_\nu c_i^2 &=&0& \text{on}&  (0,T)\times\Sigma, \\
-d_i\partial_\nu c_i^2 &=&0& \text{on}&  (0,T)\times\Gamma_\text{out}, \\
c_i^2|_{t=0}&=&c_{0,i}-c_i^1|_{t=0}& \text{in}& \Omega. \end{array} \right.
\end{equation*} By construction $c_i:=c_i^1+c_i^2\in\mathbb{E}^\Omega_p(T)$
satisfies (\ref{eq:withoutPerturbation}).
By employing the extension operator in zero time trace spaces from
\cite[Proposition 6.1]{Pruess-Saal-Simonett:Stefan-Analytic-Classical}
and that its norm is independent of $\tau<T$ the estimate for
the solution operator readily follows.
\end{proof}

\subsection{Advection terms}\label{subsec:perturbation}
Let $5/3<p<\infty$ with $p\neq 3$ and $T>0$ be given. Assume that $u$ satisfies $(A^\text{vel})$.
We prove Proposition \ref{prop:maxRegPerturbedSystem} by a perturbation argument.
To this end, we have to show that the results obtained in Lemma \ref{lm:withoutPerturbation} carry over when adding the two perturbation terms $(u\cdot \nabla)c_i$ and $(u\cdot\nu)c_i$.
Let $\mathbb{U}^\Omega_p(T)$ be given as in
(\ref{eq:regularity-calss-velocity}) and set
\begin{align*}
\mathbb{U}^\text{in}_p(T):=W^{1-1/2p}_p((0,T),L^p(\Gamma_\text{in},\mathbb{R}^3))\cap L^p((0,T),W^{2-1/p}_p(\Gamma_\text{in},\mathbb{R}^3))
\end{align*} for the Dirichlet trace space.

{\it Proof of Proposition \ref{prop:maxRegPerturbedSystem}.} The proof is carried out in three steps.

{\bfseries Step 1.} We estimate both perturbation terms occuring in (\ref{eq:withPerturbation}). Let $\varepsilon\in(0,1)$ be sufficiently small. Then the following algebra properties hold for all $5/3<p<\infty$:
\begin{align*}
\mathbb{U}^\Omega_p(T)&\cdot W^{1/2-\varepsilon/2}_p((0,T),L^p(\Omega,\mathbb{R}^3))\cap L^p((0,T),W^{1-\varepsilon}_p(\Omega,\mathbb{R}^3)) \hookrightarrow
\mathbb{F}^\Omega_p(T),\\
\mathbb{U}^\text{in}_p(T)&\cdot
W^{1-1/2p-\varepsilon/2}_p((0,T),L^p(\Gamma_\text{in}))\cap L^p((0,T),
W^{2-1/p-\varepsilon}_p(\Gamma_\text{in}))\hookrightarrow \mathbb{G}^\text{in}_p(T).
\end{align*}
The first embedding follows by a direct calculation and the second
by taking trace of the embedding
\begin{align*}
\mathbb{U}^\Omega_p(T)&\cdot W^{1-\varepsilon/2}_p((0,T),L^p(\Omega,\mathbb{R}^3))\cap
L^p((0,T),W^{2-\varepsilon}_p(\Omega,\mathbb{R}^3)) \\
&\hookrightarrow
W^{1/2}_p((0,T),L^p(\Omega))\cap L^p((0,T),W^{1,p}(\Omega)),
\end{align*}
which follows by a straight forward calculation, too.
For $0<\tau< T$ we infer the following estimates:
\begin{align*}
\Vert (u\cdot \nabla)c_i\Vert_{_0\mathbb{F}^{\Omega}_p(\tau)} &\leq C \Vert u\Vert_{\mathbb{U}^\Omega_p(T)} \Vert \nabla c_i\Vert_{{_0W}^{1/2-\varepsilon/2}_p((0,\tau),L^p(\Omega,\mathbb{R}^3))\cap L^p((0,\tau),W^{1-\varepsilon}_p(\Omega,\mathbb{R}^3))} \\
&\leq C \Vert u\Vert_{\mathbb{U}^p_p(T)} \Vert c_i\Vert_{{_0W}^{1-\varepsilon/2}_p((0,\tau),L^p(\Omega))\cap L^p((0,\tau),W^{2-\varepsilon}_p(\Omega))} \\
&\leq C\tau^\eta \Vert c_i\Vert_{_0\mathbb{E}^\Omega_p(\tau)} \qquad\qquad (c_i\in{_0\mathbb{E}}^\Omega_p(\tau))
\end{align*} with a constant $C>0$ and an exponent $\eta>0$ both being independent of $\tau$. Analogously
\begin{align*}
\Vert (u\cdot \nu)c_i\Vert_{_0\mathbb{G}^\text{in}_p(\tau)} &\leq C \Vert u\Vert_{\mathbb{U}^\text{in}_p(T)} \Vert c_i\Vert_{{_0W}^{1-1/2p-\varepsilon/2}_p((0,\tau),L^p(\Gamma_\text{in}))\cap L^p((0,\tau),W^{2-1/p-\varepsilon}_p(\Gamma_\text{in}))} \\
&\leq C \Vert u\Vert_{\mathbb{U}^\Omega_p(T)} \Vert c_i\Vert_{{_0W}^{1-\varepsilon/2}_p((0,\tau),L^p(\Omega))\cap L^p((0,\tau),W^{2-\varepsilon}_p(\Omega))} \\
&\leq C\tau^\eta \Vert c_i\Vert_{_0\mathbb{E}^\Omega_p(\tau)} \qquad\qquad (c_i\in{_0\mathbb{E}}^\Omega_p(\tau))
\end{align*} with a constant $C>0$ and an exponent $\eta>0$ both being independent of $\tau$. It follows that the linear operator
\begin{align*}
{\mathcal B}:{_0\mathbb{E}^\Omega_p(\tau)}^N\!\!\!\times{_0\mathbb{E}^\Sigma_p(\tau)}^N\rightarrow{_0\mathbb{F}^{\Omega,\Sigma}_{p}(\tau)}^N\\
{\mathcal B}(c,c^{\scriptscriptstyle\Sigma})=\left((u\cdot \nabla) c_i,\,0,\,(u\cdot\nu)c_i,\,0,\,0,\,0\right)_{i=1,...,N}
\end{align*} may be estimated by
\begin{align}\label{eq:perturbation-operator-estimate}
&\Vert {\mathcal B}(c,c^{\scriptscriptstyle\Sigma})\Vert_{{_0\mathbb{F}^{\Omega,\Sigma}_{p}(\tau)}^N}  \nonumber \\
&\leq C\tau^\eta \Vert(c,c^{\scriptscriptstyle\Sigma})\Vert_{{_0\mathbb{E}^\Omega_p(\tau)}^N\times{_0\mathbb{E}^\Sigma_p(\tau)}^N}, \;\;((c,c^{\scriptscriptstyle\Sigma})\in
{_0\mathbb{E}^\Omega_p(\tau)}^N\!\!\!\times{_0\mathbb{E}^\Sigma_p(\tau)}^N)
\end{align} with a constant $C>0$ and an exponent $\eta>0$ both independent of $\tau<T$.

{\bf Step 2.} We give the construction of the solution of (\ref{eq:withPerturbation}) as a sum
$c_i=\hat{c}_i+\bar{c}_i$, $c^{\scriptscriptstyle\Sigma}_i=\hat{c}^{\scriptscriptstyle\Sigma}_i+\bar{c}^{\scriptscriptstyle\Sigma}_i$.
Let $(\hat{c}_i,\hat{c}^{\scriptscriptstyle\Sigma}_i)\in \mathbb{E}^\Omega_p(\tau)\times\mathbb{E}^\Sigma_p(\tau)$
be the solution to (\ref{eq:withoutPerturbation}) with $g^\text{in}_i$
replaced by some $\hat g^\text{in}_i$ satisfying the compatibility
condition $-d_i\partial_\nu c_{0,i}=\hat g^\text{in}_i(0)$ in case $p>3$
and which exists according to Lemma~\ref{lm:withoutPerturbation}.
Next, we set
\begin{align*}
\bar{f}_i=-(u\cdot\nabla)\hat{c}_i,\qquad
\bar{g}^\text{in}_i = g^\text{in}_i-\hat{g}^\text{in}_i
-(u\cdot\nu)\hat{c}_i.
\end{align*}
Note that then for $p>3$ the compatibility condition
\begin{align*}
\bar{g}^\text{in}_i|_{t=0}=g^\text{in}_i|_{t=0}-\hat{g}^\text{in}_i|_{t=0}-(u|_{t=0}\cdot\nu)c_{0,i} =0
\end{align*} is satisfied by construction.
Thus, the task is reduced to prove that for $0<\tau\leq T$
there exists a unique
solution $(\bar{c}_i,\bar{c}^{\scriptscriptstyle\Sigma}_i)\in\mathbb{E}^\Omega_p(\tau)\times\mathbb{E}^\Sigma_p(\tau)$ of
\begin{equation}\label{eq:cBar}
\left\{\begin{array}{rclcl}
\partial_t \bar{c}_i+(u\cdot\nabla)\bar{c}_i- d_i\Delta \bar{c}_i &=&\bar{f}_i& \text{in}& (0,\tau)\times\Omega, \\
\partial_t \bar{c}^{\scriptscriptstyle\Sigma}_i- d^{\scriptscriptstyle\Sigma}_i \Delta_\Sigma \bar{c}^{\scriptscriptstyle\Sigma}_i &=&0& \text{on}& (0,\tau)\times \Sigma, \\
(u\cdot\nu)\bar{c}_i-d_i\partial_\nu \bar{c}_i&=&\bar{g}^\text{in}_i& \text{on}& (0,\tau)\times \Gamma_\text{in}, \\
-d_i\partial_\nu \bar{c}_i &=&0& \text{on}& (0,\tau)\times \Sigma,\\
-d_i\partial_\nu \bar{c}_i &=&0 & \text{on}& (0,\tau)\times \Gamma_\text{out}, \\
-d^{\scriptscriptstyle\Sigma}_i\partial_{\nu_\Sigma} \bar{c}^{\scriptscriptstyle\Sigma}_i &=&0 &\text{on}& (0,\tau)\times \partial \Sigma, \\
\bar{c}_i|_{t=0}&=&0   &\text{in}&\Omega,   \\
\bar{c}^{\scriptscriptstyle\Sigma}_i|_{t=0}&=&0 & \text{in}& \Sigma.
\end{array}\right.
\end{equation}
This will be done in the final step.

{\bf Step 3.} We show the unique solvability of (\ref{eq:cBar}) on
some interval $(0,\tau)$. The proof will show that $\tau$ is independent
of the data
$f_i,f^{\scriptscriptstyle\Sigma}_i,g^\text{in}_i,g^{\scriptscriptstyle\Sigma}_i,g^\text{out}_i,c_{0,i},c^{\scriptscriptstyle\Sigma}_{0,i}$.
Due to the linearity of the system solvability then carries over
to the whole time interval $(0,T)$.

We apply a Neumann series argument to (\ref{eq:cBar}). To this end, let us reformulate (\ref{eq:cBar}) by means of the operators ${_0{\mathcal L}}_\tau$ induced by the left-hand side of (\ref{eq:cBar}) and ${\mathcal B}$. Let
\begin{equation*}
\bar{F}_i=(\bar{f}_i,0,\bar{g}^\text{in}_i,0,0,0)\in {_0\mathbb{F}}^{\Omega,\Sigma}_p(\tau)^N, \qquad (i=1,...,N)
\end{equation*} such that (\ref{eq:cBar}) is equivalent to
\begin{equation*}
{_0{\mathcal L}_\tau}(\bar{c},\bar{c}^{\scriptscriptstyle\Sigma})+{\mathcal B}(\bar{c},\bar{c}^{\scriptscriptstyle\Sigma})=\bar{F} \qquad ((\bar{c},\bar{c}^{\scriptscriptstyle\Sigma})\in
{_0\mathbb{E}^\Omega_p(\tau)}^N\!\!\!\times{_0\mathbb{E}^\Sigma_p(\tau)}^N).
\end{equation*} Due to
\begin{equation*}
{_0{\mathcal L}_\tau} +{\mathcal B} =(I+{\mathcal B}\,{_0{\mathcal S}}_\tau)\,{_0{\mathcal L}}_\tau
\end{equation*} with the solution operator $_0{\mathcal S}_\tau={_0{\mathcal L}}_\tau^{-1}$ from Lemma \ref{lm:withoutPerturbation}, the invertibility of
$(I+{\mathcal B}\,{_0{\mathcal S}_\tau})$ from ${_0\mathbb{F}}^{\Omega,\Sigma}_p(\tau)$ to
${_0\mathbb{E}^\Omega_p(\tau)}^N\!\!\!\times{_0\mathbb{E}^\Sigma_p(\tau)}^N$ and, in turn, of ${_0{\mathcal L}_\tau} +{\mathcal B} $ from ${_0\mathbb{E}^\Omega_p(\tau)}^N\!\!\!\times{_0\mathbb{E}^\Sigma_p(\tau)}^N$ to ${_0\mathbb{F}}^{\Omega,\Sigma}_p(\tau)$ readily follows from (\ref{eq:perturbation-operator-estimate}) if we choose $\tau$ so small that $C\tau^\eta\,M<1$ with $M$ from  (\ref{eq:solution-operator-estimate}). Note that this is possible since $M$ is independent of $\tau<T$.\hfill $\square$

\section{Local well-posedness}\label{sec:loc}

In this section we derive unique solvability of
\begin{equation}\label{eq:locExSystem}
\left\{\begin{array}{rclcl}
\partial_t c_i+(u\cdot \nabla)c_i- d_i\Delta c_i &=&f_i& \text{in}& (0,T)\times\Omega, \\
\partial_t c^{\scriptscriptstyle\Sigma}_i - d^{\scriptscriptstyle\Sigma}_i \Delta_\Sigma c^{\scriptscriptstyle\Sigma}_i &=&r^\text{sorp}_i(c_i,c^{\scriptscriptstyle\Sigma}_i)+r^\text{ch}_i(c^{\scriptscriptstyle\Sigma}) & \text{on}& (0,T)\times \Sigma, \\
(u\cdot\nu)c_i-d_i\partial_\nu c_i&=&g^\text{in}_i& \text{on}& (0,T)\times \Gamma_\text{in}, \\
-d_i\partial_\nu c_i &=&r^\text{sorp}_i(c_i,c^{\scriptscriptstyle\Sigma}_i)& \text{on}& (0,T)\times \Sigma,\\
-d_i\partial_\nu c_i &=& 0 & \text{on}& (0,T)\times \Gamma_\text{out}, \\
-d^\Sigma_i\partial_{\nu_\Sigma} c^{\scriptscriptstyle\Sigma}_i &=&0 &\text{on}& (0,T)\times \partial \Sigma, \\
c_i|_{t=0}&=&c_{0,i}&\text{in}&\Omega,\\
c^{\scriptscriptstyle\Sigma}_i|_{t=0}&=&c^{\scriptscriptstyle\Sigma}_{0,i} & \text{on}& \Sigma,
\end{array}\right.
\end{equation} in the strong $L^p$-sense, locally in time. Recall when the index $i$ is used we mean $i=1,...,N$, and have e.g.
$c=(c_i)_{i=1,...,N}$ and $c^{\scriptscriptstyle\Sigma}=(c^{\scriptscriptstyle\Sigma}_i)_{i=1,...,N}$.

Sorption and reaction terms are required to satisfy the following assumptions:

\vspace{0,1cm}
{\bfseries ($\text{A}_\text{F}^{\text{sorp}}$)} The sorption rate acts as a function
satisfying
\begin{align*}
&r^\text{sorp}_i=r^\text{sorp}_i(c_i,c^{\scriptscriptstyle\Sigma}_i), \quad\!\! r^\text{sorp}_i\in C^2(\mathbb{R}^2), \quad\!\!
\nabla r^\text{sorp}_i\in BC^1(\mathbb{R}^2,\mathbb{R}^2).
\end{align*}

\vspace{0,1cm}
{\bfseries ($\text{A}_\text{M}^{\text{sorp}}$)}
The sorption rate increases monotonically in $c_i$ and decreases monotonically in $c^{\scriptscriptstyle\Sigma}_i$.

\vspace{0,1cm}
{\bfseries ($\text{A}_\text{B}^{\text{sorp}}$)}
The sorption rate admits linear bounds
\begin{align*}
-k^\text{de}_i c^{\scriptscriptstyle\Sigma}_i \leq r^\text{sorp}_i(c_i,c^{\scriptscriptstyle\Sigma}_i) \leq k^\text{ad}_i c_i \qquad (c_i,c^{\scriptscriptstyle\Sigma}_i\geq 0)
\end{align*} for given adsorption and desorption constants $k^\text{ad}_i,k^\text{de}_i>0$.

\vspace{0,25cm}
{\bfseries ($\text{A}^{\text{ch}}_\text{F}$)}
We assume that the chemical reactions fulfill
\begin{align*}
& r^\text{ch}_i=r^\text{ch}_i(c^{\scriptscriptstyle\Sigma}), \quad r^\text{ch} \in C^1([0,\infty)^N, \mathbb{R}^N), \qquad i=1,..,N.
\end{align*}

\vspace{0,1cm}
{\bfseries ($\text{A}^{\text{ch}}_\text{N}$)} The reaction is supposed to be quasi-positive, i.e.
\begin{align*}
r^\text{ch}_i(y)\geq 0, \qquad (y\in [0,\infty)^N, \;y_i=0).
\end{align*}

\vspace{0,1cm}
{\bfseries ($\text{A}^{\text{ch}}_\text{P}$)} The reaction admits polynomial growth, i.e.\ there exist a constant $M>0$ and an exponent $\gamma\in[1,\infty)$ if $p\in[2,\infty)$ and $\gamma=\frac{1}{1-p/2}\in [1,6]$ if $p\in(5/3,2)$, such that
\begin{align*}
|r^\text{ch}(y)|\leq M \left(1+|y|^\gamma\right) \qquad (y\in [0,\infty)^N).
\end{align*} Additionally, suppose the Jacobian fulfills
\begin{align*}
|({r^\text{ch}})'(y)|\leq M \left(1+|y|^{\gamma-1}\right) \qquad (y\in [0,\infty)^N ).
\end{align*}

\vspace{0,25cm}
\begin{remark}\label{rk:polynomial-growth} Let us comment on the polynomial growth conditions of $r^\text{ch}$ and $({r^\text{ch}})'$.
\begin{itemize}
\item[a)] In ($\text{A}^{\text{ch}}_\text{P}$) we restrict to $\gamma\in[1,6]$ if $p\in(5/3, 2)$. This is due to the embedding
\begin{equation*}
{_0\mathbb{E}}^\Sigma_p(T)\hookrightarrow L^{p\gamma}(\Sigma_T)
\end{equation*} for $p>5/3$, cf.\ \cite{Amann-Anisotropic-Function-Spaces}, which we employ in the proof of the local existence result. In case $p\geq 2$ only an arbitrary polynomial growth is required.
\item[b)] The growth rate $\gamma$ of $r^\text{ch}$ in ($\text{A}^{\text{ch}}_\text{P}$) yields that $r^\text{ch}$ acts as a Nemytskij operator
\begin{align*}
r^\text{ch}:{L^{p\gamma}(\Sigma_T)}^N\rightarrow {L^p(\Sigma_T)}^N,
\end{align*} cf.\ \cite[Theorem 3.1]{AZ90}. Analogously the growth rate $\gamma-1$ of $(r^\text{ch})'$ in ($\text{A}^{\text{ch}}_\text{P}$) yields
\begin{align*}
({r^\text{ch}})':L^{p\gamma}(\Sigma_T)^N\rightarrow L^{p\gamma/(\gamma-1)}(\Sigma_T)^{N\times N}.
\end{align*} In particular, $({r^\text{ch}})'$ maps a ball $\bar{B}_\delta$ of radius $\delta>0$ in $L^{p\gamma}(\Sigma_T)^N$ into a ball of radius $k(\delta)>0$ in $L^{p\gamma/(\gamma-1)}(\Sigma_T)^{N\times N}$. Hence an  application of the mean value theorem to the function $r^\text{ch}$ and H\"older's inequality with
\begin{align*}
\frac{1}{p}=\frac{1}{p\gamma/(\gamma-1)}+\frac{1}{p\gamma}
\end{align*} yields
\begin{align*}
\Vert r^\text{ch}&(c^{\scriptscriptstyle\Sigma})-r^\text{ch}(z^{\scriptscriptstyle\Sigma})\Vert_{L^p(\Sigma_T)^N} \\
&\leq
\sup_{v^{\scriptscriptstyle\Sigma} \in \bar{B}_\delta}\Vert (r^\text{ch})'(v^{\scriptscriptstyle\Sigma})\Vert_{L^{p\gamma/(\gamma-1)}(\Sigma_T)^{N\times N}}\Vert c^{\scriptscriptstyle\Sigma}-z^{\scriptscriptstyle\Sigma} \Vert_{L^{p\gamma}(\Sigma_T)^N}\\
&\leq k(\delta) \Vert c^{\scriptscriptstyle\Sigma}-z^{\scriptscriptstyle\Sigma} \Vert_{L^{p\gamma}(\Sigma_T)^N} \qquad (c^{\scriptscriptstyle\Sigma},z^{\scriptscriptstyle\Sigma}\in \bar{B}_\delta),
\end{align*} i.e.\ $r^\text{ch}$ acts as a locally Lipschitz continuous Nemytskij operator, cf.\ \cite[Theorem 3.10]{AZ90}.
\end{itemize}
\end{remark}

\vspace{0,25cm}
For $y\in\mathbb{R}^N$ let us denote $y^+=(y_i^+)_{i=1,...,N}$ where as before $y^+_i=\max\{0,y_i\}$. Since we do not know a priori whether a corresponding solution $(c,c^{\scriptscriptstyle\Sigma})$ is nonnegative, we extend $r^\text{ch}$ as
\begin{align*}
r^\text{ch}_{i,+}:\mathbb{R}^N\rightarrow \mathbb{R}, \qquad r^\text{ch}_{i,+}(y):= r^\text{ch}_i({y}^+) \qquad (y\in \mathbb{R}^N).
\end{align*} Then (\ref{eq:locExSystem}) remains meaningful even if $c$ or $c^{\scriptscriptstyle\Sigma}$ take negative values.

The outcome of this section is the proof of \Thmref{Local-WP}. It
is based on maximal regularity of the linear system
(\ref{eq:withPerturbation}) proved in Section~\ref{sec:linear_equations}
and the contraction mapping principle. We start by proving the nonnegativity
of $c_i$ and of $c^{\scriptscriptstyle\Sigma}_i$.

\subsection{Nonnegativity of Concentrations}

For nonnegative initial
concentrations we have the following result.

\begin{lemma}\label{lm:nonneg}\lemlabel{Nonneg} (Nonnegativity)
Let $1<p<\infty$. Let $T>0$ and let $c_{0,i}\in \mathbb{I}^{\Omega,+}_p$, $c^{\scriptscriptstyle\Sigma}_{0,i}\in \mathbb{I}^{\Sigma,+}_p$ and $g^\text{in}_i \in \mathbb{G}^\text{in}_p(T)^-$ be given. Suppose $u$ satisfies ($\text{A}^\text{vel}$), $r^\text{sorp}$ satisfies ($A_\text{F}^\text{sorp}$), ($A_\text{M}^\text{sorp}$), ($A_\text{B}^\text{sorp}$)
and $r^\text{ch}$ fulfills ($A^\text{ch}_\text{F}$), ($A^\text{ch}_\text{N}$), ($A^\text{ch}_\text{P}$).
Moreover, suppose $(c_i,c^{\scriptscriptstyle\Sigma}_i)\in\mathbb{E}^\Omega_p(T)\times\mathbb{E}^\Sigma_p(T)$ is a strong $L^p$-solution of (\ref{eq:locExSystem}).  Then
\begin{align*}
c_i\geq 0 \quad \text{a.e.}\;\text{in}\;\Omega_T, \qquad c^{\scriptscriptstyle\Sigma}_i \geq 0 \quad \text{a.e.}\; \text{on}\; \Sigma_T
\end{align*} hold true.
\end{lemma}

\begin{proof}
Let $\phi_\varepsilon \in C^\infty(\mathbb{R})$ be a pointwise approximation of
\begin{equation*}
\phi(r)=\left\{\begin{array}{cr} -r & :r\leq 0 \\ 0 & :r>0 \end{array}\right.,
\end{equation*} as $\varepsilon \rightarrow 0+$, which satisfies $\phi_\varepsilon\geq 0$, $\phi_\varepsilon'\leq 0$ and $\phi_\varepsilon''\geq0$, e.g.
\begin{align*}
\phi_\varepsilon(r):=\left\{
	\begin{array}{rl}
	-re^{\varepsilon/r}&:r\leq 0\\
	0&:r>0
	\end{array}\right..
\end{align*}
Then we have for $c_i^-=\max\{0,-c_i\}$ that
\begin{align}\label{eq:properties_of_phi}
\phi_\varepsilon (c_i) \rightarrow c_i^-, \qquad c_i \phi_\varepsilon'(c_i) \rightarrow c_i
\left\{\begin{array}{cr} -1&: c_i< 0\\ 0 &: c_i\geq 0 \end{array}\right\} =c_i^-
\end{align} as $\varepsilon\rightarrow 0+$. We show
\begin{equation*}
\lim_{\varepsilon\rightarrow 0}\left[\intop\limits_\Omega \phi_\varepsilon(c_i)dx+\intop\limits_\Sigma \phi_\varepsilon(c^{\scriptscriptstyle\Sigma}_i) d\sigma\right]\leq 0.
\end{equation*} Applying $\phi_\varepsilon$ to $c_i$ we obtain by partial integration
\begin{align}
\frac{d}{dt}\intop\limits_\Omega \phi_\varepsilon (c_i) dx &= \intop\limits_\Omega \phi_\varepsilon'(c_i) \partial_t c_i dx = \intop\limits_\Omega \phi_\varepsilon'(c_i) \div(d_i\nabla c_i -uc_i) dx\nonumber\\
&=\intop\limits_{\partial\Omega} \phi_\varepsilon'(c_i)(d_i \partial_\nu c_i -(u\cdot\nu)c_i) d\sigma - \intop\limits_\Omega \phi_\varepsilon''(c_i) \nabla c_i \cdot (d_i \nabla c_i - uc_i) dx\nonumber\\
&=-\intop\limits_{\Gamma_\text{in}} \phi_\varepsilon'(c_i) g^\text{in}_i d\sigma - \intop\limits_\Sigma \phi_\varepsilon'(c_i) r^\text{sorp}_i d\sigma - \intop\limits_{\Gamma_\text{out}} \phi_\varepsilon'(c_i)(u\cdot\nu) c_i d\sigma\nonumber\\
&- d_i \intop\limits_\Omega \phi_\varepsilon''(c_i) |\nabla c_i|^2 dx + \intop\limits_\Omega \phi_\varepsilon''(c_i) \nabla c_i u c_i dx \label{eq:nonneg_bulk}
\end{align}due to the boundary conditions. In the same way we have
\begin{align}
\frac{d}{dt}\intop\limits_\Sigma \phi_\varepsilon(c^{\scriptscriptstyle\Sigma}_i) d\sigma &= \intop\limits_\Sigma \phi_\varepsilon'(c^{\scriptscriptstyle\Sigma}_i) \partial_tc^{\scriptscriptstyle\Sigma}_i d\sigma =\intop\limits_\Sigma \phi_\varepsilon'(c^{\scriptscriptstyle\Sigma}_i)(d^{\scriptscriptstyle\Sigma}_i\Delta_\Sigma c^{\scriptscriptstyle\Sigma}_i +r^\text{sorp}_i(c_i,c^{\scriptscriptstyle\Sigma}_i)+r^\text{ch}_{i,+}(c^{\scriptscriptstyle\Sigma}))d\sigma\nonumber\\
&=-d^{\scriptscriptstyle\Sigma}_i\intop\limits_\Sigma \phi_\varepsilon''(c^{\scriptscriptstyle\Sigma}_i) |\nabla_\Sigmac^{\scriptscriptstyle\Sigma}_i|^2 d\sigma+\intop\limits_\Sigma \phi_\varepsilon'(c^{\scriptscriptstyle\Sigma}_i)(r^\text{sorp}_i(c_i,c^{\scriptscriptstyle\Sigma}_i)+r^\text{ch}_{i,+}(c^{\scriptscriptstyle\Sigma}))d\sigma.\label{eq:nonneg_surface}
\end{align} Let us go through all the integrals appearing on the right-hand side of (\ref{eq:nonneg_bulk}) and (\ref{eq:nonneg_surface}). The first and the fourth integrals on the right-hand side of (\ref{eq:nonneg_bulk}) and the first integral on the right-hand side of (\ref{eq:nonneg_surface}) are negative or zero such that we may drop them. The remaining four integrals are treated as follows: We combine the sorption boundary integrals to
\begin{equation*} \intop\limits_\Sigma (\phi_\varepsilon'(c^{\scriptscriptstyle\Sigma}_i)-\phi_\varepsilon'(c_i))r^\text{sorp}_i(c_i,c^{\scriptscriptstyle\Sigma}_i) d\sigma
\end{equation*} and split up this integral into three integrals on
\begin{equation*}
\{\text{sign}(c_i)=\text{sign}(c^{\scriptscriptstyle\Sigma}_i)\}, \qquad \{c_i\leq 0\leq c^{\scriptscriptstyle\Sigma}_i\},\qquad \{c^{\scriptscriptstyle\Sigma}_i\leq 0\leq c_i\}.
\end{equation*}
When $c_i$ and $c^{\scriptscriptstyle\Sigma}_i$ have the same sign this integral tends to $0$ as $\varepsilon\rightarrow 0+$.
On $\{c_i\leq 0\leq c^{\scriptscriptstyle\Sigma}_i\}$ we have
\begin{align*} \intop\limits_\Sigma \mathds{1}_{\{c_i\leq 0\leq c^{\scriptscriptstyle\Sigma}_i\}}(\phi_\varepsilon'(c^{\scriptscriptstyle\Sigma}_i)-\phi_\varepsilon'(c_i))r^\text{sorp}_i(c_i,c^{\scriptscriptstyle\Sigma}_i) d\sigma\\ \rightarrow \intop\limits_\Sigma r^\text{sorp}_i(c_i,c^{\scriptscriptstyle\Sigma}_i) d\sigma \leq \intop\limits_\Sigma r^\text{sorp}_i(0,c^{\scriptscriptstyle\Sigma}_i) d\sigma \leq 0
\end{align*} as $\varepsilon\rightarrow 0+$ by monotonicity of $r^\text{sorp}_i$ and $r^\text{sorp}_i (0,c^{\scriptscriptstyle\Sigma}_i)\leq 0$. In the same way on $\{c^{\scriptscriptstyle\Sigma}_i\leq 0\leq c_i\}$ we obtain
\begin{align*} \intop\limits_\Sigma\mathds{1}_{\{c^{\scriptscriptstyle\Sigma}_i\leq 0\leq c_i\}} (\phi_\varepsilon'(c^{\scriptscriptstyle\Sigma}_i)-\phi_\varepsilon'(c_i))r^\text{sorp}_i(c_i,c^{\scriptscriptstyle\Sigma}_i) d\sigma\\ \rightarrow -\intop\limits_\Sigma r^\text{sorp}_i(c_i,c^{\scriptscriptstyle\Sigma}_i) d\sigma \leq -\intop\limits_\Sigma r^\text{sorp}_i(c_i,0) d\sigma \leq 0
\end{align*} as $\varepsilon\rightarrow 0+$ by monotonicity of $r^\text{sorp}_i$ and $r^\text{sorp}_i(c_i,0)\geq0$.
Because of $u\cdot \nu \geq 0$ on $\Gamma_\text{out}$ we see that
\begin{align*}
-\intop\limits_{\Gamma_\text{out}} \phi_\varepsilon'(c_i) (u\cdot \nu) c_i d\sigma \rightarrow -\intop\limits_{\Gamma_\text{out}} c_i^- (u\cdot\nu) d\sigma \leq 0.
\end{align*} We treat the reaction boundary integral by the quasi-positivity of $r^\text{ch}$ as follows.
We show
\begin{equation*}
\intop\limits_\Sigma \phi_\varepsilon'(c^{\scriptscriptstyle\Sigma}_i) r^\text{ch}_{i,+}(c^{\scriptscriptstyle\Sigma}) d\sigma\leq 0
\end{equation*}
through
\begin{align*}
\intop\limits_\Sigma \phi_\varepsilon'(c^{\scriptscriptstyle\Sigma}_i) r^\text{ch}_{i,+}(c^{\scriptscriptstyle\Sigma}) d\sigma
&=\intop\limits_\Sigma \phi_\varepsilon'(c^{\scriptscriptstyle\Sigma}_i)\mathds{1}_{\{c^{\scriptscriptstyle\Sigma}_i>0\}} r^\text{ch}_{i,+}(c^{\scriptscriptstyle\Sigma}) d\sigma \\
&+\intop\limits_\Sigma \phi_\varepsilon'(c^{\scriptscriptstyle\Sigma}_i) \mathds{1}_{\{c^{\scriptscriptstyle\Sigma}_i=0\}}r^\text{ch}_{i,+}(c^{\scriptscriptstyle\Sigma}) d\sigma
+ \intop\limits_\Sigma \phi_\varepsilon'(c^{\scriptscriptstyle\Sigma}_i)\mathds{1}_{\{c^{\scriptscriptstyle\Sigma}_i<0\}} r^\text{ch}_{i,+}(c^{\scriptscriptstyle\Sigma}) d\sigma.
\end{align*} The first integral vanishes by the properties of $\phi_\epsilon'$, the second one is less than or equal to zero
by quasi-positivity and $\phi_\varepsilon'(0)\leq 0$ as $\varepsilon\rightarrow 0+$. The third one is less than or equal to zero by definition of the extension of $r^\text{ch}$ to $\mathbb{R}^N$, i.e.\ because of $c^{\scriptscriptstyle\Sigma}_i<0$ implies ${c^{\scriptscriptstyle\Sigma}_i}^+=0$ and $r^\text{ch}_i({c^{\scriptscriptstyle\Sigma}}^+)\geq 0$ by quasi-positivity.
We turn to the remaining integral $\intop\limits\nolimits_\Omega \phi_\varepsilon''(c_i) \nabla c_i u c_i dx$. Here we make use of $\nabla(\phi_\varepsilon'(c_i))=\phi_\varepsilon''(c_i)\nabla c_i$ and in the same manner of $\nabla(\phi_\varepsilon(c_i))=\phi_\varepsilon'(c_i)\nabla c_i$ and integrate by parts twice, such that
\begin{align*}
\intop\limits_\Omega \phi_\varepsilon''(c_i)\nabla c_i \cdot u c_i dx &=\intop\limits_\Omega \nabla (\phi_\varepsilon'(c_i))\cdot uc_i dx \\
&= \intop\limits_{\partial\Omega} \phi_\varepsilon' (c_i) (u\cdot \nu) c_i d\sigma - \intop\limits_\Omega \phi_\varepsilon'(c_i) \underbrace{\div(u c_i)}_{=u\cdot \nabla c_i} dx\\
&= \intop\limits_{\partial\Omega} \phi_\varepsilon' (c_i) (u\cdot \nu) c_i d\sigma - \intop\limits_\Omega\nabla(\phi_\varepsilon (c_i)) u dx\\
&= \intop\limits_{\partial\Omega} \phi_\varepsilon'(c_i)(u\cdot\nu) c_i d\sigma - \intop\limits_{\partial\Omega} \phi_\varepsilon(c_i)(u\cdot\nu) d\sigma,
\end{align*} where in the second and in the last step we made use of $\div u=0$. Employing
(\ref{eq:properties_of_phi}) we see
\begin{align*}
\intop\limits_\Omega \phi_\varepsilon''(c_i) \nabla c_i \cdot u c_i dx= \intop\limits_{\partial\Omega} (\phi_\varepsilon'(c_i)c_i - \phi_\varepsilon(c_i))(u\cdot\nu) d\sigma \rightarrow 0
\end{align*} as $\varepsilon\rightarrow 0+$. Therefore summing up (\ref{eq:nonneg_bulk}) and (\ref{eq:nonneg_surface}), integration in time over $[0,t]$ and taking the limit $\varepsilon\rightarrow 0+$ yields
\begin{align*}
\intop\limits_\Omega c_i^-(t) dx +\intop\limits_\Sigma {c^{\scriptscriptstyle\Sigma}_i}^-(t)d\sigma=\intop\limits_\Omega \phi(c_i(t)) dx +\intop\limits_\Sigma \phi(c^{\scriptscriptstyle\Sigma}_i(t))d\sigma\\
\leq \intop\limits_\Omega \phi (c_{0,i}) dx +\intop\limits_\Sigma \phi(c^{\scriptscriptstyle\Sigma}_{0,i})d\sigma= \intop\limits_\Omega (c_{0,i})^- dx +\intop\limits_\Sigma (c^{\scriptscriptstyle\Sigma}_{0,i})^-d\sigma=0
\end{align*} which in turn gives $c_i^-=0$ a.e.\ in $\Omega_T$, ${c^{\scriptscriptstyle\Sigma}_i}^-=0$ a.e.\ on $\Sigma_T$ and therefore $c_i\geq 0$ a.e.\ in $\Omega_T$, $c^{\scriptscriptstyle\Sigma}_i\geq 0$ a.e.\ on $\Sigma_T$.
Note that for $\varepsilon \rightarrow 0+$ we make use of Lebesgue's theorem on dominated convergence.
\end{proof}
\subsection{Existence of Solutions}
Let $T'>0$ be given and $T\leq T'$. Assume a set of (fixed) data
\begin{equation*} (f_i,0,g^\text{in}_i,0,0,0,c_{0,i},c^{\scriptscriptstyle\Sigma}_{0,i})\in\mathbb{F}^{\Omega,\Sigma}_{p,I}(T')
\end{equation*}
to be given. We denote by
\begin{equation*} {\mathcal L}_{T,i}:\mathbb{E}^\Omega_p(T)\times \mathbb{E}^\Sigma_p(T)\rightarrow \mathbb{F}^{\Omega,\Sigma}_{p}(T)
\end{equation*}
the isomorphism induced by Proposition~\ref{prop:maxRegPerturbedSystem},
that is, ${\mathcal L}_{T,i}$ is the full linear operator on the right-hand side
of (\ref{eq:withPerturbation}) (except for the time traces).
The full nonlinear problem (\ref{eq:cat}) then is reformulated as
\begin{align}
	{\mathcal L}_{T,i}(c_i,c^{\scriptscriptstyle\Sigma}_i)&=(f_i,0,g^\text{in}_i,0,0,0)
	+{\mathcal N}_{T,i}(c,c^{\scriptscriptstyle\Sigma}),\label{reformnp}\\
	c_i(0)&=c_{0,i},\quad c^{\scriptscriptstyle\Sigma}_i(0)=c^{\scriptscriptstyle\Sigma}_{i,0},
	\quad i=1,\ldots,N,\nonumber
\end{align}
where ${\mathcal N}_{T,i}$ includes the nonlinear sorption and reaction terms, i.e.,
\begin{align*}
{\mathcal N}_{T,i}(c,c^{\scriptscriptstyle\Sigma}):=\big(0,r^\text{sorp}_i(c_i,c^{\scriptscriptstyle\Sigma}_i)+r^\text{ch}_{i,+}(c^{\scriptscriptstyle\Sigma}),0,r^\text{sorp}_i(c_i,c^{\scriptscriptstyle\Sigma}_i),0,0\big).
\end{align*}
In order keep the constants resulting from the estimates below independent of
$T$, we employ a suitable zero time trace splitting as described in
the following.

First we take care of the compatibility condition
arising from the nonlinear boundary condition on $\Sigma$. Taking time
trace results in $r^\text{sorp}_i(c_{0,i},c^{\scriptscriptstyle\Sigma}_{0,i})\in
W^{1-3/p}_p(\Sigma)$, which will be extended to $\mathbb{G}^\Sigma_p(T)$ by setting
\begin{equation*}
r^*_i:={\mathcal R}_\Sigma
e^{t\Delta_{\Sigma_{(-\infty,\infty)}}}{\mathcal E}_{\Sigma_{(-\infty,\infty)}}
r^\text{sorp}_i(c_{0,i},c^{\scriptscriptstyle\Sigma}_{0,i}).
\end{equation*}
Here ${\mathcal E}_{\Sigma_{(-\infty,\infty)}}$ denotes the extension operator
from the lateral surface $\Sigma$ to the surface of the infinite cylinder $\Sigma_{(-\infty,\infty)}$
and ${\mathcal R}_\Sigma$ the corresponding restriction operator (note that both act
as bounded operators on the function classes considered here,
cf.\ \cite{Adams-Fournier:Sobolev-Spaces}).
Since $e^{t\Delta_{\Sigma_{(-\infty,\infty)}}}$ has the same regularizing
properties as the Laplacian on the whole space $\mathbb{R}^n$, for which the desired
regularity is well known \cite{Pruess-Saal-Simonett:Stefan-Analytic-Classical},
we see that $r^*_i\in \mathbb{G}^\Sigma_p(T)$.

Now we define the reference solution
$(c^*_i,{c^{\scriptscriptstyle\Sigma}_i}^*)\in\mathbb{E}^\Omega_p(T')\times\mathbb{E}^\Sigma_p(T')$,
existing according to Proposition~\ref{prop:maxRegPerturbedSystem}, via
\begin{equation}\label{eq:OpLAndreferenceSolution}
{\mathcal L}_{T,i}(c_i^*,{c^{\scriptscriptstyle\Sigma}_i}^*)=(f_i,0,g^\text{in}_i,r^*_i,0,0), \quad c_i(0)=c_{0,i} \quad \text{in}\; \Omega,   \quad c^{\scriptscriptstyle\Sigma}_i(0)=c^{\scriptscriptstyle\Sigma}_{0,i} \quad \text{on}\;\Sigma.
\end{equation}
Decomposing $(c_i,c^{\scriptscriptstyle\Sigma}_i)$ as
\begin{equation*} c_i=\bar{c}_i+c^\ast_i, \quad c^{\scriptscriptstyle\Sigma}_i=\bar{c}^{\scriptscriptstyle\Sigma}_i+{c^{\scriptscriptstyle\Sigma}_i}^\ast
\end{equation*}
and subtracting (\ref{eq:OpLAndreferenceSolution}) from (\ref{reformnp}),
we end up with the reduced zero time trace problem
\begin{equation*}
{_0{\mathcal L}}_{T,i}(\bar{c},\bar{c}^{\scriptscriptstyle\Sigma})={_0{\mathcal N}}_{T,i}(\bar{c},\bar{c}^{\scriptscriptstyle\Sigma}) \qquad (i=1,...,N).
\end{equation*}
Here ${_0{\mathcal L}}_{T,i}$ denotes
the restriction of ${\mathcal L}_{T,i}$ to
$\,_0\mathbb{E}^\Omega_p(T)\times{_0\mathbb{E}^\Sigma_p(T)}$ and
\[
	{_0{\mathcal N}}_{T,i}(\bar{c},\bar{c}^{\scriptscriptstyle\Sigma})
	:={\mathcal N}_{T,i}(\bar{c}_i+c^\ast_i,\bar{c}^{\scriptscriptstyle\Sigma}_i+{c^{\scriptscriptstyle\Sigma}_i}^\ast)
	-(0,0,0,r_i^\ast,0,0).
\]
Next, we define $\Phi_{T}:=(\Phi_{T,i})_{i=1,...,N}$ through
\begin{align*}
&{_0\Phi}_{T}:{\,_0\mathbb{E}^\Omega_p(T)}^N\!\!\times{_0\mathbb{E}^\Sigma_p(T)}^N\rightarrow{\,_0\mathbb{E}^\Omega_p(T)}^N\!\!\times{_0\mathbb{E}^\Sigma_p(T)}^N,\\
&{_0\Phi}_{T,i}(\bar{c},\bar{c}^{\scriptscriptstyle\Sigma}):={_0\mathcal{S}}_{T,i}\,{_0\mathcal{N}}_{T,i}(\bar{c},\bar{c}^{\scriptscriptstyle\Sigma}), \qquad (i=1,...,N),
\end{align*} with the bounded linear inverse ${_0\mathcal{S}}_{T,i}$ of ${_0{\mathcal L}}_{T,i}$ given in Proposition \ref{prop:maxRegPerturbedSystem}.

{\it Proof of \Thmref{Local-WP}.} We apply the contraction mapping principle to ${_0\Phi}_T$, i.e., we show that there exists a
$\delta>0$, such that the mapping ${_0\Phi_T}$ constitutes a contraction on the closed ball $\bar{B}_\delta(0)\subset {\,_0\mathbb{E}^\Omega_p(T)}^N\!\!\times{_0\mathbb{E}^\Sigma_p(T)}^N$ and fulfills ${_0\Phi}_T:\bar{B}_\delta(0)\rightarrow \bar{B}_\delta(0)$.

{\it (i) Contraction property:} Let $(\bar{c},\bar{c}^{\scriptscriptstyle\Sigma}), (\bar{z},\bar{z}^{\scriptscriptstyle\Sigma})\in
\bar{B}_\delta(0)$. Then we have
\begin{align}
&\Vert{_0\Phi}_T(\bar{c},\bar{c}^{\scriptscriptstyle\Sigma})-{_0\Phi}_T(\bar{z},\bar{z}^{\scriptscriptstyle\Sigma})
\Vert_{{\,_0\mathbb{E}^\Omega_p(T)}^N\!\!\times{_0\mathbb{E}^\Sigma_p(T)}^N}\nonumber\\
&\leq C \Vert{_0\mathcal{N}}_T(\bar{c},\bar{c}^{\scriptscriptstyle\Sigma})- {_0\mathcal{N}}_T(\bar{z},\bar{z}^{\scriptscriptstyle\Sigma}) \Vert_{{_0\mathbb{F}}^{\Omega,\Sigma}_p(T)^N}\nonumber\\
&=C\Vert r^\text{sorp}(\bar{c}+c^*,\bar{c}^{\scriptscriptstyle\Sigma}+{c^{\scriptscriptstyle\Sigma}}^*)-r^\text{sorp}(\bar{z}+c^*,\bar{z}^{\scriptscriptstyle\Sigma}+{c^{\scriptscriptstyle\Sigma}}^*) \Vert _{{(\mathbb{F}^\Sigma_p(T)\cap {_0\mathbb{G}}^\Sigma_p(T))}^N} \nonumber\\
&\qquad\strut+C\Vert r^\text{ch}(({\bar{c}^{\scriptscriptstyle\Sigma}+{c^{\scriptscriptstyle\Sigma}}^*})^+)- r^\text{ch}(({\bar{z}^{\scriptscriptstyle\Sigma}+{c^{\scriptscriptstyle\Sigma}}^*})^+)\Vert_{\mathbb{F}^\Sigma_p(T)^N} \label{eq:ContractionEstimate}
\end{align}
with
\begin{align*}
C:=\sup \{\Vert {_0\mathcal{S}}_T\Vert_{\mathscr{L}({\,_0\mathbb{F}^{\Omega,\Sigma}_{p}(T)^N},\,{\,_0\mathbb{E}^\Omega_p(T)}^N\!\!\times{_0\mathbb{E}^\Sigma_p(T)}^N)}:\; T\in (0,T']\}
\end{align*} independent of $T$, cf.\! Proposition \ref{prop:maxRegPerturbedSystem}.
>From Remark \ref{rk:polynomial-growth} we infer that
\begin{align*}
\Vert r^\text{ch}(({\bar{c}^{\scriptscriptstyle\Sigma}+{c^{\scriptscriptstyle\Sigma}}^*})^+)- r^\text{ch}(({\bar{z}^{\scriptscriptstyle\Sigma}+{c^{\scriptscriptstyle\Sigma}}^*})^+)\Vert_{L^p(\Sigma_T)^N}\leq L \Vert \bar{c}^{\scriptscriptstyle\Sigma}-\bar{z}^{\scriptscriptstyle\Sigma} \Vert_{L^{p\gamma}(\Sigma_T)^N}
\end{align*} for a constant $L>0$ depending on $\delta$ but not $T$
and $\gamma$. Note in passing that we also used that $h\mapsto h^+$ is globally
Lipschitz continuous from $L^{p\gamma}(\Sigma_T)$ to $L^{p\gamma}(\Sigma_T)$ with
Lipschitz constant $1$. By the fact that $p>5/3$ we can estimate as
\begin{align*}
\Vert \bar{c}^{\scriptscriptstyle\Sigma}-\bar{z}^{\scriptscriptstyle\Sigma}\Vert_{L^{p\gamma}(\Sigma_T)^N} \leq K T^\eta \Vert \bar{c}^{\scriptscriptstyle\Sigma}-\bar{z}^{\scriptscriptstyle\Sigma} \Vert_{{_0\mathbb{E}}^\Sigma_p(T)^N}
\end{align*} with a constant $K>0$ and an exponent $\eta>0$ independent of $T$. We arrive at
\begin{align}\label{eq:ReactionEstimate}
\Vert r^\text{ch}(({\bar{c}^{\scriptscriptstyle\Sigma}+{c^{\scriptscriptstyle\Sigma}}^*})^+)- r^\text{ch}(({\bar{z}^{\scriptscriptstyle\Sigma}+{c^{\scriptscriptstyle\Sigma}}^*})^+)\Vert_{\mathbb{F}^\Sigma_p(T)^N}\leq LK T^\eta \Vert \bar{c}^{\scriptscriptstyle\Sigma}-\bar{z}^{\scriptscriptstyle\Sigma}\Vert_{_0\mathbb{E}^\Sigma_p(T)^N}.
\end{align} We turn to the estimate of the sorption rate. By
($\text{A}_\text{F}^{\text{sorp}}$), ($\text{A}_\text{B}^{\text{sorp}}$)
the mapping $r^\text{sorp}_i$ acts as a Nemytskij operator from
$W^s_p(\Sigma_T)\times W^s_p(\Sigma_T)$ to $W^s_p(\Sigma_T)$ for $s\in(0,1)$,
cf.\ Section~3.1 in \cite{Sickel96}.
Note that ($\text{A}_\text{B}^{\text{sorp}}$) implies $r^\text{sorp}_i(0,0)=0$.
An application of the mean value theorem to the function $r^\text{sorp}_i$ and using $\nabla r^\text{sorp}_i\in BC^1(\mathbb{R}^2,\mathbb{R}^2)$ yields the global Lipschitz continuity of its induced Nemytskij operator. Hence we may employ
\begin{align*}
W^s_p(\Sigma_T)=W^s_p((0,T),L^p(\Sigma))\cap L^p((0,T),W^s_p(\Sigma))
\end{align*} for $s\in(0,1)$ and
\begin{align*}
{_0\mathbb H^\Sigma_p}(T)\times {_0\mathbb{E}^\Sigma_p}(T)\hookrightarrow \big({_0W^{1-1/p+\epsilon}_p}((0,T), L^p(\Sigma))\cap L^p((0,T),W^{1-1/p+\epsilon}_p(\Sigma))\big)^2 \\
 \overset{r^\text{sorp}_i}{\longrightarrow} {_0W^{1-1/p+\epsilon}_p}((0,T), L^p(\Sigma))\cap L^p((0,T), W^{1-1/p+\epsilon}_p(\Sigma))
 \hookrightarrow {_0\mathbb{G}^\Sigma_p}(T),
\end{align*}
for sufficiently small $\epsilon>0$, such that we obtain, similarly as for the estimate of the reaction term,
\begin{align}
\Vert r^\text{sorp}_i(\bar{c}+c^*,\bar{c}^{\scriptscriptstyle\Sigma}+{c^{\scriptscriptstyle\Sigma}}^*)-r^\text{sorp}_i(\bar{z}+c^*,\bar{z}^{\scriptscriptstyle\Sigma}+{c^{\scriptscriptstyle\Sigma}}^*) \Vert _{ {_0\mathbb{G}}^\Sigma_p(T)^N} \nonumber\\
\leq L'K'T^\eta \left(\Vert \bar{c}_i-\bar{z}_i\Vert_{_0\mathbb{E}^\Omega_p(T)} +\Vert \bar{c}^{\scriptscriptstyle\Sigma}_i-\bar{z}^{\scriptscriptstyle\Sigma}_i \Vert_{_0\mathbb{E}^\Sigma_p(T)} \right) \label{eq:SorptionEstimate}
\end{align} with constants $L',K'>0$ and an exponent $\eta>0$
independent of $T<T'$. Combining (\ref{eq:ReactionEstimate}) and
(\ref{eq:SorptionEstimate}) yields
\begin{align*}
&\Vert{_0\Phi}_T(\bar{c},\bar{c}^{\scriptscriptstyle\Sigma})-{_0\Phi}_T(\bar{z},\bar{z}^{\scriptscriptstyle\Sigma})
\Vert_{{\,_0\mathbb{E}^\Omega_p(T)}^N\!\!\times{_0\mathbb{E}^\Sigma_p(T)}^N}\\
&\leq C(LK+L'K')T^\eta\,\Vert (\bar{c},\bar{c}^{\scriptscriptstyle\Sigma})-(\bar{z},\bar{z}^{\scriptscriptstyle\Sigma})\Vert_{{\,_0\mathbb{E}^\Omega_p(T)}^N\!\!\times{_0\mathbb{E}^\Sigma_p(T)}^N}
\end{align*}
for $(\bar{c},\bar{c}^{\scriptscriptstyle\Sigma}), (\bar{z},\bar{z}^{\scriptscriptstyle\Sigma})\in \bar{B}_\delta(0)$.
We choose $T$ so small that
\begin{align}\label{eq:ChoosingT}
C(LK+L'K')T^\eta\leq \frac{1}{2},
\end{align} which is possible since all other constants appearing in (\ref{eq:ChoosingT}) are independent of $T<T'$.
Hence ${_0\Phi}_T$ is a contraction on $\bar{B}_\delta(0)$.

\vspace{0,25cm}
{\it (ii) Self mapping property:} Let $(\bar{c},\bar{c}^{\scriptscriptstyle\Sigma})\in \bar{B}_\delta(0)$.
Then we have
\begin{align*}
&\Vert {_0\Phi}_T(\bar{c},\bar{c}^{\scriptscriptstyle\Sigma})\Vert_{{\,_0\mathbb{E}^\Omega_p(T)}^N\!\!\times{_0\mathbb{E}^\Sigma_p(T)}^N} \leq C\Vert {_0\mathcal{N}}_T(\bar{c}+c^*,\bar{c}^{\scriptscriptstyle\Sigma}+{c^{\scriptscriptstyle\Sigma}}^*)\Vert_{_0\mathbb{F}^{\Omega,\Sigma}_p(T)^N} \\
&\leq C\Vert r^\text{sorp}(\bar{c}+c^*,\bar{c}^{\scriptscriptstyle\Sigma}+{c^{\scriptscriptstyle\Sigma}}^*)-r^* \Vert_{_0\mathbb{G}^\Sigma_p(T)^N}
+C\Vert r^\text{ch}({(\bar{c}^{\scriptscriptstyle\Sigma}+{c^{\scriptscriptstyle\Sigma}}^*)}^+) \Vert_{\mathbb{F}^\Sigma_p(T)^N}.
\end{align*}
Analogously to (i) we estimate the reaction term by
\begin{align}
&\Vert r^\text{ch}(({\bar{c}^{\scriptscriptstyle\Sigma}+{c^{\scriptscriptstyle\Sigma}}^*})^+)\Vert_{L^p(\Sigma_T)^N}\nonumber\\
&\leq \Vert r^\text{ch}(({\bar{c}^{\scriptscriptstyle\Sigma}+{c^{\scriptscriptstyle\Sigma}}^*})^+)-r^\text{ch}({({c^{\scriptscriptstyle\Sigma}}^*)}^+)\Vert_{L^p(\Sigma_T)^N}+\Vert r^\text{ch}(({c^{\scriptscriptstyle\Sigma}}^*)^+)\Vert_{L^p(\Sigma_T)^N}\nonumber\\
&\leq L \Vert \bar{c}^{\scriptscriptstyle\Sigma}\Vert_{L^{p\gamma}(\Sigma_T)^N}+\Vert r^\text{ch}(({c^{\scriptscriptstyle\Sigma}}^*)^+)\Vert_{L^p(\Sigma_T)^N}\nonumber\\
&\leq L K T^\eta \Vert \bar{c}^{\scriptscriptstyle\Sigma}\Vert_{{_0\mathbb{E}^\Sigma_p(T)}^N}+\Vert r^\text{ch}(({c^{\scriptscriptstyle\Sigma}}^*)^+)\Vert_{L^p(\Sigma_T)^N} \label{eq:SelfmappingReaction}
\end{align} with constants $L,K>0$ being independent of
$T$. In the same manner as in
(i) for the sorption term we obtain
\begin{align}
&\Vert r^\text{sorp}_i(\bar{c}_i+c_i^*,\bar{c}^{\scriptscriptstyle\Sigma}_i+{c^{\scriptscriptstyle\Sigma}_i}^*)-r^*_i\Vert_{_0\mathbb{G}^\Sigma_p(T)}
\nonumber\\
&\leq \Vert
r^\text{sorp}_i(\bar{c}_i+c_i^*,\bar{c}^{\scriptscriptstyle\Sigma}_i+{c^{\scriptscriptstyle\Sigma}_i}^*)-r^\text{sorp}_i(c_i^*,{c^{\scriptscriptstyle\Sigma}_i}^*)\Vert_{_0\mathbb{G}^\Sigma_p(T)}\nonumber\\
&\qquad\strut +\Vert r^\text{sorp}_i(c_i^*,{c^{\scriptscriptstyle\Sigma}_i}^*)\Vert_{\mathbb{G}^\Sigma_p(T)} +\Vert r_i^*\Vert_{_0\mathbb{G}^\Sigma_p(T)}\nonumber\\
& \leq L'K'T^\eta\Vert (\bar{c}_i,\bar{c}^{\scriptscriptstyle\Sigma}_i)\Vert_{{\,_0\mathbb{E}^\Omega_p(T)}^N\!\!\times{_0\mathbb{E}^\Sigma_p(T)}^N}
 +\Vert r^\text{sorp}_i(c_i^*,{c^{\scriptscriptstyle\Sigma}_i}^*)\Vert_{\mathbb{G}^\Sigma_p(T)} +\Vert r_i^*\Vert_{_0\mathbb{G}^\Sigma_p(T)} \label{eq:SelfmappingSorption}
\end{align}
with constants $L',K'>0$ and an exponent $\eta>0$ being all independent of $T$. Putting together (\ref{eq:SelfmappingReaction}) and (\ref{eq:SelfmappingSorption}) yields
\begin{align*}
\Vert
{_0\Phi}_T(\bar{c},\bar{c}^{\scriptscriptstyle\Sigma})\Vert_{{\,_0\mathbb{E}^\Omega_p(T)}^N\!\!\times{_0\mathbb{E}^\Sigma_p(T)}^N}
&\leq C(LK+L'K')T^\eta \delta +\Vert r^\text{ch}(({c^{\scriptscriptstyle\Sigma}}^*)^+)\Vert_{L^p(\Sigma_T)^N}\\
&\qquad\strut+\Vert r^\text{sorp}_i(c_i^*,{c^{\scriptscriptstyle\Sigma}_i}^*)\Vert_{\mathbb{G}^\Sigma_p(T)} +\Vert r_i^*\Vert_{_0\mathbb{G}^\Sigma_p(T)}.
\end{align*}
Since $c_i^*$, ${c^{\scriptscriptstyle\Sigma}_i}^*$, and $r_i^*$ are fixed functions,
notice that the latter three terms can be made small by choosing
$T>0$ small. Thus, by choosing $T$ so small that the sum of those
three terms is less than $\delta/2$ and such that (\ref{eq:ChoosingT})
is satisfied, we arrive at
\begin{align*}
\Vert {_0\Phi}_T(\bar{c},\bar{c}^{\scriptscriptstyle\Sigma})\Vert_{{\,_0\mathbb{E}^\Omega_p(T)}^N\!\!\times{_0\mathbb{E}^\Sigma_p(T)}^N} \leq\delta \qquad
((\bar{c},\bar{c}^{\scriptscriptstyle\Sigma}) \in \bar{B}_\delta(0)).
\end{align*}
The proof of \Thmref{Local-WP} is now complete.
\hfill $\square$

\section{Global Well-Posedness}
\seclabel{Global-WP}
In this section we show that the local-in-time strong solutions obtained in \Thmref{Local-WP}
in fact exist globally, provided the reaction rates satisfy some structural condition that allows for
the derivation of a priori estimates.
We suppose in addition to ($\textrm{A}^{\textrm{ch}}_{\textrm{F}}$), ($\textrm{A}^{\textrm{ch}}_{\textrm{N}}$) and ($\textrm{A}^{\textrm{ch}}_{\textrm{F}}$)
the following assumption on the structure of the reaction rates.

		\vspace*{0.5em}
		{\bfseries ($\textrm{A}^{\textrm{ch}}_{\textrm{S}}$)}
		There exists a lower triangular matrix $Q=(q_{ij})_{1\leq i,j\leq N} \in \mathbb R^{N \times N}$
		with strictly positive diagonal entries such that
		\begin{equation}\label{eq:A^ch_S-structure-of-reaction}
			Q r^{\textrm{ch}}(y) \leq C \left( 1 + \sum^N_{j = 1} y_j \right) v, \qquad y \in [0,\infty)^N
		\end{equation}
		for some constant $C > 0$ and $v = (1,\,\dots,\,1)$.
		\vspace*{0.5em}

\smallskip
Reaction-diffusion systems with this triangular condition have been
widely studied by several authors; see \cite{Pie10} and the references
cited therein. When proving global existence results, condition (\ref{eq:A^ch_S-structure-of-reaction}) allows for an iteration scheme which has been applied successfully in many situations. A major objective of this section is to generalize this iteration scheme for standard reaction-diffusion systems subject to ($\textrm{A}^{\textrm{ch}}_{\textrm{S}}$) to heterogeneous catalysis. The main difference compared to standard systems lies in the fact that the reaction takes place on the boundary instead inside the bulk and that we also have to deal with terms arising from sorption processes.

In this setting, we obtain the global-in-time well-posedness result, given by \Thmref{Global-WP}.
The proof of \Thmref{Global-WP} requires, besides the maximal regularity estimates obtained in \Secref{linear_equations},
also some comparison principles and some weak-type estimates, which are provided by the following results.
\begin{lemma}\lemlabel{Linear-Positivity}
	Let $T > 0$ and let $1 < p < \infty$.
	Let $\alpha,\,\beta > 0$. \vspace*{-0.5\baselineskip}
	\begin{enumerate}[a)]
		\item Assume $f \in \mathbb F^\Omega_p(T)$, $g^{\textrm{in}} \in \mathbb G^{\textrm{in}}_p(T)$, $g^{\scriptscriptstyle\Sigma} \in \mathbb G^\Sigma_p(T)$
			and $v_0 \in \mathbb I^\Omega_p$ with
			\begin{equation*}
				(u \cdot \nu) v_0 - \beta \partial_\nu v_0 = g^{\textrm{in}}(0) \quad \mbox{on}\ \Gamma_{\textrm{in}}, \qquad
				- \beta \partial_\nu v_0 = g^{\scriptscriptstyle\Sigma}(0) \quad \mbox{on}\ \Sigma,
			\end{equation*}
			if $p > 3$.
			Let $v \in \mathbb E^\Omega_p(T)$ be a strong solution to
			\begin{equation*}
				\left\{\begin{array}{rclll}
					\alpha \partial_t v + (u \cdot \nabla) v - \beta \Delta v & = & f               & \quad \mbox{in} & (0,\,T) \times \Omega,                \\[0.5em]
				                     (u \cdot \nu) v - \beta \partial_\nu v & = & g^{\textrm{in}} & \quad \mbox{on} & (0,\,T) \times \Gamma_{\textrm{in}},  \\[0.5em]
				                                     - \beta \partial_\nu v & = & g^{\scriptscriptstyle\Sigma}        & \quad \mbox{on} & (0,\,T) \times \Sigma,                \\[0.5em]
				                                     - \beta \partial_\nu v & = & 0               & \quad \mbox{on} & (0,\,T) \times \Gamma_{\textrm{out}}, \\[0.5em]
					                                                     v(0) & = & v_0             & \quad \mbox{in} & \Omega.
				\end{array}\right.
			\end{equation*}
			If $f,\,v_0 \geq 0$ and $g^{\textrm{in}},\,g^{\scriptscriptstyle\Sigma} \leq 0$, then $v \geq 0$. \vspace*{0.5\baselineskip}
		\item Assume $f \in \mathbb F^\Sigma_p(T)$ and $v_0 \in \mathbb I^\Sigma_p$ with $- \beta \partial_{\nu_\Sigma} v_0 = 0$ on $\partial \Sigma$, if $p > 3$.
			Let $v \in \mathbb E^\Sigma_p(T)$ be a strong solution to
			\begin{equation*}
				\left\{\begin{array}{rclll}
					\alpha \partial_t v - \beta \Delta_\Sigma v & = & f   & \quad \mbox{on} & (0,\,T) \times \Sigma,          \\[0.5em]
					     - \beta \partial_{\nu_\Sigma} v & = & 0   & \quad \mbox{on} & (0,\,T) \times \partial \Sigma, \\[0.5em]
			                                    v(0) & = & v_0 & \quad \mbox{on} & \Sigma.
				\end{array}\right.
			\end{equation*}
			If $f,\,v_0 \geq 0$, then $v \geq 0$.
	\end{enumerate}
\end{lemma}
\begin{proof}
	The proof follows the lines of the proof of \Lemref{Nonneg},
	except that here we deal with a linear problem, only.
\end{proof}
\begin{lemma}\lemlabel{Weak-Estimates-Domain}
	Let $T^\ast > 0$ and let $2 \leq p < \infty$ and $2 \leq q \leq \infty$.
	Let $\alpha,\,\beta > 0$.
	Assume $f \in \mathbb F^\Omega_p(T)$, $g^{\textrm{in}} \in \mathbb G^{\textrm{in}}_p(T)$, $g^{\scriptscriptstyle\Sigma} \in \mathbb G^\Sigma_p(T)$
	and $v_0 \in \mathbb I^\Omega_p \cap BC(\Omega)$ with
	\begin{equation*}
		(u \cdot \nu) v_0 - \beta \partial_\nu v_0 = g^{\textrm{in}}(0) \quad \mbox{on}\ \Gamma_{\textrm{in}}, \qquad
		- \beta \partial_\nu v_0 = g^{\scriptscriptstyle\Sigma}(0) \quad \mbox{on}\ \Sigma,
	\end{equation*}
	if $p > 3$.
	Let $v \in \mathbb E^\Omega_p(T)$ be a strong solution to
	\begin{equation*}
		\left\{\begin{array}{rclll}
			\alpha \partial_t v + (u \cdot \nabla) v - \beta \Delta v & = & f               & \quad \mbox{in} & (0,\,T) \times \Omega,                \\[0.5em]
		                     (u \cdot \nu) v - \beta \partial_\nu v & = & g^{\textrm{in}} & \quad \mbox{on} & (0,\,T) \times \Gamma_{\textrm{in}},  \\[0.5em]
		                                     - \beta \partial_\nu v & = & g^{\scriptscriptstyle\Sigma}        & \quad \mbox{on} & (0,\,T) \times \Sigma,                \\[0.5em]
		                                     - \beta \partial_\nu v & = & 0               & \quad \mbox{on} & (0,\,T) \times \Gamma_{\textrm{out}}, \\[0.5em]
			                                                     v(0) & = & v_0             & \quad \mbox{in} & \Omega.
		\end{array}\right.
	\end{equation*}
	Then
	\begin{equation*}
		\|v\|_{L^q(\Omega_T)} + \|v\|_{L^q(\Sigma_T)} \leq C \Big( \|f\|_{L^q(\Omega_T)} + \|g^{\textrm{in}}\|_{L^q(\Gamma_{\textrm{in},T})} + \|g^{\scriptscriptstyle\Sigma}\|_{L^q(\Sigma_T)} + \|v_0\|_{L^q(\Omega)} \Big)
	\end{equation*}
	for some constant $C> 0$ that is independent of $p$ and $0 < T < T^\ast$.
\end{lemma}
\begin{proof} If the right-hand side is infinite, nothing has to be proved. So, we now assume that the data admits finite $L^q$-norm. We distinguish two cases:\\
{\it Case 1:} $p>5/2$, such that $\mathbb{E}^\Omega_p(T)\hookrightarrow BC(\overline{\Omega_T})$, cf.\ \cite{Amann-Anisotropic-Function-Spaces}. \\
{\it $q=2$:} standard (multiply the equation by $v$ to obtain
$L^2$-estimates).\\
{\it $q=\infty$:} Let ${\mathcal L}_T$ denote the operator given by the left-hand side of the system under consideration, such that ${\mathcal L}_T v=(f,g^\text{in},g^{\scriptscriptstyle\Sigma},0,v_0)$. For fixed $5/2<r<3$ let $\phi\in \mathbb{E}^\Omega_r(T)$ denote the solution of
\begin{align*}
{\mathcal L}_T\phi =(1,-1,-1,0,1)=:F.
\end{align*} Note that $\mathbb{E}^\Omega_r(T)\hookrightarrow BC(\overline{\Omega_T})$ and that no compatibility conditions for $F$ occur. Let
\begin{align*}
\delta:=\Vert f\Vert_{L^\infty(\Omega_T)}+\Vert g^\text{in}\Vert_{L^\infty(\Gamma_{\text{in},T})}+\Vert g^{\scriptscriptstyle\Sigma}\Vert_{L^\infty(\Sigma_T)}+\Vert v_0\Vert_{L^\infty(\Omega)}
\end{align*} and let $\bar{v}\in \mathbb{E}^\Omega_p(T)$ be given through ${\mathcal L}_T\bar{v}=\delta F$. Then the comparison principle (\Lemref{Linear-Positivity} a)) applied to $\bar{v}-v$ yields $v\leq \bar{v}$ on $\overline{\Omega_T}$.
Since $\bar{v}=\delta \phi$ and $\Vert\bar{v}\Vert_{BC}=\delta \Vert \phi \Vert_{BC}$ it follows that
\begin{align*}
v\leq \delta \left(\Vert \phi \Vert_{L^\infty(\Omega_T)}+\Vert \phi \Vert_{L^\infty(\Sigma_T)}\right)=:\delta C
\end{align*} on $\overline{\Omega_T}$. Analogously $-\delta C\leq v$, such that $\Vert v \Vert_{L^\infty(\Omega_T)}\leq C\delta$.
Since in this case $v\in BC(\overline{\Omega_T})$, we trivially have
\begin{align*}
\Vert v\Vert_{L^\infty(\Sigma_T)}\leq \Vert v \Vert_{L^\infty(\Omega_T)}\leq C\delta,
\end{align*} hence that
\begin{align*}
\Vert v\Vert_{L^\infty(\Sigma_T)}+\Vert v \Vert_{L^\infty(\Omega_T)}\leq C\delta.
\end{align*}
For $2\leq q\leq\infty$ we set
\begin{align*}
X_q:=L^q(\Omega_T)\times L^q(\Sigma_T), \quad Y_q:=L^q(\Omega_T)\times L^q(\Gamma_{\text{in},T})\times L^q(\Sigma_T)\times \{0\}\times L^q(\Omega).
\end{align*} Let ${\mathcal S}_T$ denote the system's solution operator. By the $L^2$- and the $L^\infty$-estimates obtained above we have
\begin{align*}
{\mathcal S}_T \in \mathscr{L}(Y_2,X_2)\cap \mathscr{L}(Y_\infty,X_\infty).
\end{align*} By interpolation ${\mathcal S}_T\in \mathscr{L}(Y_q,X_q)$ which yields the assertion for Case 1.

\medskip
{\it Case 2:} $2\leq p\leq 5/2$.\\
For convenience set $F:=(f,0,g^\text{in},g^{\scriptscriptstyle\Sigma},0,0,v_0,0)\in \mathbb{F}^{\Omega,\Sigma}_{p,I}(T)$. Let $r>5/2$ and choose
\begin{align*}
F_n:=&(f_n,0,g^\text{in}_n,g^{\scriptscriptstyle\Sigma}_n,0,0,v_{0,n},0) \in \mathbb{Y}_{r,q}:= \\
\mathbb{F}^{\Omega,\Sigma}_{r,I}(T)\cap
(L^q(\Omega_T)\times\{0\}\times &L^q(\Gamma_{\text{in},T})\times L^q(\Sigma_T)\times\{0\}\times\{0\}\times L^q(\Omega)\times\{0\})
\end{align*} for $n\in \mathbb{N}$, such that $F_n\rightarrow F$ in $\mathbb{Y}_{p,q}$ (e.g.\ by extension of $F$ to $\mathbb{R}\times\mathbb{R}^3$ respectively $\mathbb{R}\times (\partial A \times \mathbb{R})$ and mollification). Let $v_n\in \mathbb{E}^\Omega_r(T)$ denote the corresponding solution of ${\mathcal L}_T v_n=F_n$. Then Case 1 applies to $F_n,v_n$ and there exists a $C>0$ independent of $n\in\mathbb{N}$, such that
\begin{align}\label{eq:Lq-estimate-approximation}
\Vert v_n \Vert_{L^q(\Omega_T)}+\Vert v_n\Vert_{L^q(\Sigma_T)} \leq C \Vert F_n\Vert_{L^q}.
\end{align} Obviously $(v_n)_n$ is a Cauchy sequence in
\begin{align*}
\mathbb{X}_{p,q}:=\mathbb{E}^\Omega_p(T)\cap \left(L^q(\Omega_T)\times
L^q(\Sigma_T)\right),
\end{align*}
such that we may pass to the limit $n\rightarrow \infty$  in (\ref{eq:Lq-estimate-approximation}). Hence we obtain $v_n\rightarrow v$ in $\mathbb{X}_{p,q}$ with $v$ being the solution of ${\mathcal L}_T v=F$ and
 \begin{align*}
\Vert v\Vert_{L^q(\Omega_T)}+\Vert v\Vert_{L^q(\Sigma_T)} \leq C \Vert F\Vert_{L^q}.
\end{align*}
\end{proof}

The next Lemma is standard for equations on standard domains $\Omega$, cf. \cite[Lemma 3.4]{Pie10}. Here we give a proof since we employ it on $\Sigma$.
\begin{lemma}\lemlabel{Weak-Estimates-Boundary}
	Let $T^\ast >0$ and let $1 <p,\,q<\infty$.
	Let $\mu > 0$ and let arbitrary coefficients $\alpha_1,\,\dots,\,\alpha_N,\,\beta_1,\,\dots,\,\beta_N \in \mathbb R$ be given.
	Assume the data satisfies $f,\,g_1,\,\dots,\,g_n \in \mathbb F^\Sigma_p(T)$ and $u_0,\,v^{1}_0,\,\dots,\,v^{N}_0 \in \mathbb I^\Sigma_p \cap BC(\Sigma)$
	with $\partial_{\nu_\Sigma} u_0 = \partial_{\nu_\Sigma} v^{j}_0 = 0$ on $\partial \Sigma$, if $p > 3$.
	Furthermore, let $u,\,v_1,\,\dots,\,v_N \in \mathbb E^\Sigma_p(T)$ be strong solutions to
	\begin{equation}\label{eq:duality-Lm-forward-heat-equation}
		\left\{\begin{array}{rclrclll}
			\partial_t u - \mu \Delta_\Sigma u & = & f,   & \qquad \alpha_j \partial_t v_j - \beta_j \Delta_\Sigma v_j & = & g_j       & \qquad \mbox{on} & (0,\,T) \times \Sigma,          \\[0.5em]
			  - \mu \partial_{\nu_\Sigma} u & = & 0,   & \qquad              - \beta_j \partial_{\nu_\Sigma} v_j & = & 0         & \qquad \mbox{on} & (0,\,T) \times \partial \Sigma, \\[0.5em]
			                           u(0) & = & u_0, & \qquad                                           v_j(0) & = & v^{j}_0 & \qquad \mbox{on} & \Sigma
		\end{array}\right.
	\end{equation}
	for some $0 < T < T^\ast$.
	If $f \leq \sum^N_{j = 1} g_j$, then
\begin{equation*} \|u^+\|_{L^q(\Sigma_T)} \leq C \left(1 + \sum^N_{j = 1} \|v_j\|_{L^q(\Sigma_T)}\right)
\end{equation*}
 for some constant $C = C(\|u_0\|_{BC(\Sigma)},\,\|v^{1}_0\|_{BC(\Sigma)},\,\dots,\,\|v^{N}_0\|_{BC(\Sigma)}) > 0$
	that is independent of $0 < T < T^\ast$.
\end{lemma}
\begin{proof}
Let $\theta \in L^{q^\prime}(\Sigma_T)^+$ with $1/q+1/q'=1$ and let $0<\tau<T$. By the transformation $t\mapsto \tau -t$ applied to the equation for $u$ in (\ref{eq:duality-Lm-forward-heat-equation}) with data $f=\theta$, $u_0=0$ we obtain the backward heat equation,
\begin{align}\label{eq:backward-heat-equation}
\left\{\begin{array}{rcll}
-[\partial_t \phi +\mu \Delta_\Sigma \phi]&=&\theta(\tau-t)&\text{on}\; (0,\tau)\times\Sigma,\\
-\mu \partial_{\nu_{\scriptscriptstyle\Sigma}}\phi&=&0&\text{on}\;(0,\tau)\times\partial\Sigma,\\
\phi(\tau)&=&0&\text{on}\;\Sigma,
 \end{array}\right.
\end{align}
which admits maximal $L^{q'}$-regularity, cf.\ Section \ref{sec:linear_equations}. In particular, we have $\phi \in BUC([0,\tau],L^{q'}(\Sigma))$ and
\begin{align*}
\sup_{s\in[0,\tau]}\Vert \phi(s)\Vert_{L^{q'}(\Sigma)}+\Vert \partial_t\phi\Vert_{L^{q'}(\Sigma_\tau)}+\mu \Vert \Delta_\Sigma\phi \Vert_{L^{q'}(\Sigma_\tau)} \leq C \Vert \theta\Vert_{L^{q'}(\Sigma_\tau)}
\end{align*} for a constant $C>0$ which is independent of $\tau$. Observe that $\phi\geq 0$, since $\theta\geq 0$. Hence by plugging in the first line of (\ref{eq:backward-heat-equation}) and multiple partial integrations we have
\begin{align*}
\intop\limits_{\Sigma_\tau} u\theta d(t,\sigma) &=\intop\limits_\Sigma u_0\phi(0) d\sigma + \intop\limits_{\Sigma_\tau}\phi(\partial_t u-\mu \Delta_\Sigma u) d(t,\sigma)  \\
&\leq \intop\limits_\Sigma u_0\phi(0) d\sigma +\sum_{j=1}^N \intop\limits_{\Sigma_\tau} \phi (\alpha_j\partial_t v_j -\beta_j \Delta_\Sigma v_j) d(t,\sigma) \\
&=\intop\limits_\Sigma (u_0-\sum_{j=1}^N\alpha_j v^{j}_0)\phi(0) d\sigma + \sum_{j=1}^N \intop\limits_{\Sigma_\tau}\left(-\alpha_j\partial_t\phi -\beta_j \Delta_\Sigma \phi\right)v_j d(t,\sigma) \\
&\leq C \Vert \theta\Vert_{L^{q'}(\Sigma_\tau)}\left(1+\sum_{j=1}^N \Vert v_j\Vert_{L^q(\Sigma_\tau)}\right).
\end{align*} Employing
\begin{equation*}
		\|u^+\|_{L^q(\Sigma_T)} = \sup \left\{\ \intop\limits^T_0 \intop\limits_\Sigma u \theta\,\mbox{d}\sigma(x)\,\mbox{d}t\,:\,\theta \in L^{q^\prime}(\Sigma_T)^+,\ \|\theta\|_{L^{q^\prime}(\Sigma_T)} < 1\ \right\}
	\end{equation*}
	the assertion follows.
\end{proof}

Now, we are in position to prove \Thmref{Global-WP}.
\begin{proof}[Proof of \Thmref{Global-WP}]
	Since according to \Thmref{Local-WP} the (local-in-time) solutions to (\ref{eq:cat}) generate a local semi-flow in the
	phase space $\mathbb I^\Omega_2 \times \mathbb I^\Sigma_2 = H^1(\Omega)
	\times H^1(\Sigma)$,
	we may assume $T^\ast < \infty$ and show that the solution stays bounded in $H^1(\Omega) \times H^1(\Sigma)$ on $(0,\,T^\ast)$
	in order to obtain a contradiction.
	It is sufficient to establish $L^\infty$-bounds for the solution in order to obtain boundedness in the phase space.
Then the $H^1$-boundedness of solutions follows, as it is shown in the last part of this proof.

	We will now derive $L^\infty$-bounds, which requires several steps.
	Note that we may use the fact that $c_i,\,c^{\scriptscriptstyle\Sigma}_i \geq 0$ on $(0,\,T^\ast)$ thanks to \Lemref{Nonneg}.

	{\bfseries Step 1.}
	We have $c_i \in \mathbb E^\Omega_2(T) = H^1((0,\,T),\,L^2(\Omega)) \cap L^2((0,\,T),\,H^2(\Omega))$ and
	\begin{equation}
		\eqnlabel{Global-WP-Base-Domain}
		\left\{\begin{array}{rclll}
			\partial_t c_i + (u \cdot \nabla) c_i - d_i \Delta c_i & = & f_i                                   & \quad \mbox{in} & (0,\,T) \times \Omega,                \\[0.5em]
		                (u \cdot \nu) c_i - d_i \partial_\nu c_i & = & g^{\textrm{in}}_i                     & \quad \mbox{on} & (0,\,T) \times \Gamma_{\textrm{in}},  \\[0.5em]
		                                  - d_i \partial_\nu c_i & = & r^{\textrm{sorp}}_i(c_i,\,c^{\scriptscriptstyle\Sigma}_i) & \quad \mbox{on} & (0,\,T) \times \Sigma,                \\[0.5em]
		                                  - d_i \partial_\nu c_i & = & 0                                     & \quad \mbox{on} & (0,\,T) \times \Gamma_{\textrm{out}}, \\[0.5em]
			                                                c_i(0) & = & c_{0,i}                               & \quad \mbox{in} & \Omega
		\end{array}\right.
	\end{equation}
	for $i = 1,\,\dots,\,N$ and all $0 < T < T^\ast$.
	Now, $r^{\textrm{sorp}}_i(c_i,\,c^{\scriptscriptstyle\Sigma}_i) \geq - k^{\textrm{de}}_i c^{\scriptscriptstyle\Sigma}_i$.
	Thus, by \Lemref{Linear-Positivity}~a) we have $0 \leq c_i \leq z_i$ for the unique maximal regular solution $z_i$ to
	\begin{equation*}
		\left\{\begin{array}{rclll}
			\partial_t z_i + (u \cdot \nabla) z_i - d_i \Delta z_i & = & f_i               & \quad \mbox{in} & (0,\,T) \times \Omega,                \\[0.5em]
		                (u \cdot \nu) z_i - d_i \partial_\nu z_i & = & g^{\textrm{in}}_i & \quad \mbox{on} & (0,\,T) \times \Gamma_{\textrm{in}},  \\[0.5em]
		                                  - d_i \partial_\nu z_i & = &-C c^{\scriptscriptstyle\Sigma}_i      & \quad \mbox{on} & (0,\,T) \times \Sigma,                \\[0.5em]
		                                  - d_i \partial_\nu z_i & = & 0                 & \quad \mbox{on} & (0,\,T) \times \Gamma_{\textrm{out}}, \\[0.5em]
			                                                z_i(0) & = & c_{0,i}           & \quad \mbox{in} & \Omega
		\end{array}\right.
	\end{equation*}
	with some appropriate constant $C = C((k^{\textrm{de}}_j)_{j = 1, \dots, N}) > 0$.
	Note that this problem allows for a strong solution in the $L^2$-setting without any compatibility conditions
	between the right hand sides of the boundary conditions and the initial value.
	Since $z_i$ is a solution to a linear problem, we may write $z_i = u_i + v_i$,
	where $u_i$ solves
	\begin{equation*}
		\left\{\begin{array}{rclll}
			\partial_t u_i + (u \cdot \nabla) u_i - d_i \Delta u_i & = & 0            & \quad \mbox{in} & (0,\,T) \times \Omega,                \\[0.5em]
		                (u \cdot \nu) u_i - d_i \partial_\nu u_i & = & 0            & \quad \mbox{on} & (0,\,T) \times \Gamma_{\textrm{in}},  \\[0.5em]
		                                  - d_i \partial_\nu u_i & = & -C c^{\scriptscriptstyle\Sigma}_i & \quad \mbox{on} & (0,\,T) \times \Sigma,                \\[0.5em]
		                                  - d_i \partial_\nu u_i & = & 0            & \quad \mbox{on} & (0,\,T) \times \Gamma_{\textrm{out}}, \\[0.5em]
			                                                u_i(0) & = & 0            & \quad \mbox{in} & \Omega
		\end{array}\right.
	\end{equation*}
	and $v_i$ solves
	\begin{equation*}
		\left\{\begin{array}{rclll}
			\partial_t v_i + (u \cdot \nabla) v_i - d_i \Delta v_i & = & f_i               & \quad \mbox{in} & (0,\,T) \times \Omega,                \\[0.5em]
		                (u \cdot \nu) v_i - d_i \partial_\nu v_i & = & g^{\textrm{in}}_i & \quad \mbox{on} & (0,\,T) \times \Gamma_{\textrm{in}},  \\[0.5em]
		                                  - d_i \partial_\nu v_i & = & 0                 & \quad \mbox{on} & (0,\,T) \times \Sigma,                \\[0.5em]
		                                  - d_i \partial_\nu v_i & = & 0                 & \quad \mbox{on} & (0,\,T) \times \Gamma_{\textrm{out}}, \\[0.5em]
			                                                v_i(0) & = & c_{0,i}            & \quad \mbox{in} & \Omega
		\end{array}\right.
	\end{equation*}
	For these solutions we have
	\begin{equation*}
		\|u_i\|_{L^q(\Omega_T)} + \|u_i\|_{L^q(\Sigma_T)} \leq C^\prime \|c^{\scriptscriptstyle\Sigma}_i\|_{L^q(\Sigma_T)}
	\end{equation*}
	as well as
	\begin{equation*}
		\|v_i\|_{L^q(\Omega_T)} + \|v_i\|_{L^q(\Sigma_T)} \leq A_i
	\end{equation*}
	provided that $2 \leq q \leq\infty$.
	Here, we employed \Lemref{Weak-Estimates-Domain} to obtain constants $C^\prime = C^\prime(C,\,q) > 0$
	and $A_i = A_i(\|f_i\|_{L^q(\Omega_T)},\,\|g^{\textrm{in}}_i\|_{L^q(\Gamma_{\textrm{in},T})},\,\|c_{0,i}\|_{L^q(\Omega)},\,q) > 0$
	that are independent of $0 < T < T^\ast$.
	Since $\|c_i\|_{L^q} \leq \|z_i\|_{L^q} \leq \|u_i\|_{L^q} + \|v_i\|_{L^q}$ we may sum up the above estimates to obtain
	\begin{equation}
		\eqnlabel{Global-WP-Estimates-Domain}
		\sum^N_{i = 1} \|c_i\|_{L^q(\Omega_T)} + \sum^N_{i = 1} \|c_i\|_{L^q(\Sigma_T)} \leq C^\ast \left( 1 + \sum^N_{j = 1} \|c^{\scriptscriptstyle\Sigma}_j\|_{L^q(\Sigma_T)} \right),
	\end{equation}
	where $C^\ast = C^\ast(C^\prime, (A_j)_{j = 1, \dots, N}) > 0$ denotes a constant that is independent of $0 < T < T^\ast$.
	Note that this estimate is available for all $2 \leq q \leq \infty$.

	{\bfseries Step 2.}
	We have $c^{\scriptscriptstyle\Sigma}_i \in \mathbb E^\Sigma_2(T) = H^1((0,\,T),\,L^2(\Sigma)) \cap L^2((0,\,T),\,H^2(\Sigma))$ and
	\begin{equation}
		\eqnlabel{Global-WP-Base-Boundary}
		\left\{\begin{array}{rclll}
			\partial_t c^{\scriptscriptstyle\Sigma}_i - d^{\scriptscriptstyle\Sigma}_i \Delta_\Sigma c^{\scriptscriptstyle\Sigma}_i & = & r^{\textrm{sorp}}_i(c_i,\,c^{\scriptscriptstyle\Sigma}_i) + r^{\textrm{ch}}_i(c^{\scriptscriptstyle\Sigma}) & \quad \mbox{on} & (0,\,T) \times \Sigma,          \\[0.5em]
			              - d^{\scriptscriptstyle\Sigma}_i \partial_{\nu_\Sigma} c^{\scriptscriptstyle\Sigma}_i & = & 0                                                                   & \quad \mbox{on} & (0,\,T) \times \partial \Sigma, \\[0.5em]
			                                              c^{\scriptscriptstyle\Sigma}_i(0) & = & c^{\scriptscriptstyle\Sigma}_{0,i}                                                       & \quad \mbox{on} & \Sigma
		\end{array}\right.
	\end{equation}
	for $i = 1,\,\dots,\,N$ and all $0 < T < T^\ast$.
	Now we use the triangular structure of the reaction rates that is guaranteed by ($\textrm{A}^{\textrm{ch}}_{\textrm{S}}$) to
	treat the cases $i = 1$ and $i = 2,\,\dots,\,N$ separately.

	{\bfseries Step 2.1.}
	According to assumption ($\textrm{A}^{\textrm{sorp}}_{\textrm{B}}$) we have $r^{\textrm{sorp}}_1(c_1,\,c^{\scriptscriptstyle\Sigma}_1) \leq k^{\textrm{ad}}_1 c_1$
	and according to ($\textrm{A}^{\textrm{ch}}_{\textrm{S}}$) we have $q_{11} r^{\textrm{ch}}_1(c^{\scriptscriptstyle\Sigma}) \leq C (1 + \sum^N_{j = 1} c^{\scriptscriptstyle\Sigma}_j)$
	for some fixed constant $C > 0$.
	Thus, by \Lemref{Linear-Positivity}~b) we have $0 \leq c^{\scriptscriptstyle\Sigma}_1 \leq z^{\scriptscriptstyle\Sigma}_1$ for the unique solution $z^{\scriptscriptstyle\Sigma}_1$ to
	\begin{equation*}
		\left\{\begin{array}{rclll}
			\partial_t z^{\scriptscriptstyle\Sigma}_1 - d^{\scriptscriptstyle\Sigma}_1 \Delta_\Sigma z^{\scriptscriptstyle\Sigma}_1 & = & C^\prime \Big( 1 + c_1 + \sum^N_{j = 1} c^{\scriptscriptstyle\Sigma}_j \Big) & \quad \mbox{on} & (0,\,T) \times \Sigma,          \\[0.5em]
			              - d^{\scriptscriptstyle\Sigma}_1 \partial_{\nu_\Sigma} z^{\scriptscriptstyle\Sigma}_1 & = & 0                                                        & \quad \mbox{on} & (0,\,T) \times \partial \Sigma, \\[0.5em]
			                                              z^{\scriptscriptstyle\Sigma}_1(0) & = & c^{\scriptscriptstyle\Sigma}_{0,1}                                            & \quad \mbox{on} & \Sigma
		\end{array}\right.
	\end{equation*}
	with some appropriate constant $C^\prime = C^\prime(C,\,k^{\textrm{ad}}_1,\,q_{11}) > 0$.
	Since $z^{\scriptscriptstyle\Sigma}_1$ is a solution to a linear problem, we may write $z^{\scriptscriptstyle\Sigma}_1 = u^{\scriptscriptstyle\Sigma}_1 + v^{\scriptscriptstyle\Sigma}_1$,
	where $u^{\scriptscriptstyle\Sigma}_1$ solves
	\begin{equation*}
		\left\{\begin{array}{rclll}
			\partial_t u^{\scriptscriptstyle\Sigma}_1 - d^{\scriptscriptstyle\Sigma}_1 \Delta_\Sigma u^{\scriptscriptstyle\Sigma}_1 & = & C^\prime \Big( 1 + c_1 + \sum^N_{j = 1} c^{\scriptscriptstyle\Sigma}_j \Big) & \quad \mbox{on} & (0,\,T) \times \Sigma,          \\[0.5em]
			              - d^{\scriptscriptstyle\Sigma}_1 \partial_{\nu_\Sigma} u^{\scriptscriptstyle\Sigma}_1 & = & 0                                                        & \quad \mbox{on} & (0,\,T) \times \partial \Sigma, \\[0.5em]
			                                              u^{\scriptscriptstyle\Sigma}_1(0) & = & 0                                                        & \quad \mbox{on} & \Sigma
		\end{array}\right.
	\end{equation*}
	and $v^{\scriptscriptstyle\Sigma}_1$ solves
	\begin{equation*}
		\left\{\begin{array}{rclll}
			\partial_t v^{\scriptscriptstyle\Sigma}_1 - d^{\scriptscriptstyle\Sigma}_1 \Delta_\Sigma v^{\scriptscriptstyle\Sigma}_1 & = & 0             & \quad \mbox{on} & (0,\,T) \times \Sigma,          \\[0.5em]
			              - d^{\scriptscriptstyle\Sigma}_1 \partial_{\nu_\Sigma} v^{\scriptscriptstyle\Sigma}_1 & = & 0             & \quad \mbox{on} & (0,\,T) \times \partial \Sigma, \\[0.5em]
			                                              v^{\scriptscriptstyle\Sigma}_1(0) & = & c^{\scriptscriptstyle\Sigma}_{0,1} & \quad \mbox{on} & \Sigma.
		\end{array}\right.
	\end{equation*}
	For these solutions we have
	\begin{equation*}
		\|u^{\scriptscriptstyle\Sigma}_1\|_{L^q(\Sigma_T)} \leq M \|u^{\scriptscriptstyle\Sigma}_1\|_{{}_0 \mathbb E^\Sigma_p(T)} \leq L M C^\prime \left( \Theta + \|c_1\|_{L^p(\Sigma_T)} + \sum^N_{j = 1} \|c^{\scriptscriptstyle\Sigma}_j\|_{L^p(\Sigma_T)} \right)
	\end{equation*}
	with $\Theta = \Theta(p) = \|1\|_{L^p(\Sigma_{T^\ast})}$ as well as
	\begin{equation*}
		\|(v^{\scriptscriptstyle\Sigma}_1)^+\|_{L^q(\Sigma_T)} \leq A_1,
	\end{equation*}
	provided that $2 \leq p,\,q < \infty$.
	Here, $M = M(p,\,q) > 0$ denotes the norm of the embedding ${}_0 \mathbb E^\Sigma_p(T) \hookrightarrow L^q(\Sigma_T)$, cf.\ \cite{Amann-Anisotropic-Function-Spaces}
	and $L = L(p) > 0$ denotes the norm of the solution operator in the $L^p$-setting for the time interval $(0,\,T)$,
	which are both independent of $0 < T < T^\ast$ thanks to the homogeneous initial condition.
	Furthermore, $A_1 = A_1(\|c^{\scriptscriptstyle\Sigma}_{0, 1}\|_{BC(\Sigma)},\,q) > 0$ denotes the constant delivered by \Lemref{Weak-Estimates-Boundary},
	which is also independent of $0 < T < T^\ast$. Observe, that a standard maximum principle could have been applied here, too.
	Note that $v^{\scriptscriptstyle\Sigma}_1\leq (v^{\scriptscriptstyle\Sigma}_1)^+$ implies $0 \leq c^{\scriptscriptstyle\Sigma}_1\leq z^{\scriptscriptstyle\Sigma}_1 \leq u^{\scriptscriptstyle\Sigma}_1 + (v^{\scriptscriptstyle\Sigma}_1)^+$.
	Therefore,
	\begin{equation*}
		\begin{array}{rcl}
			\|c^{\scriptscriptstyle\Sigma}_1\|_{L^q(\Sigma_T)} & \leq & \|z^{\scriptscriptstyle\Sigma}_1\|_{L^q(\Sigma_T)} \leq \|u^{\scriptscriptstyle\Sigma}_1\|_{L^q(\Sigma_T)} + \|(v^{\scriptscriptstyle\Sigma}_1)^+\|_{L^q(\Sigma_T)} \\[1.0em]
				& \leq & C^{\prime\prime} {\displaystyle \left( 1 + \|c_1\|_{L^p(\Sigma_T)} + \sum^N_{j = 1} \|c^{\scriptscriptstyle\Sigma}_j\|_{L^p(\Sigma_T)} \right),}
		\end{array}
	\end{equation*}
	provided that $2 \leq p,\,q < \infty$.
	Here $C^{\prime\prime} = C^{\prime\prime}(L M C^\prime,\,\Theta,\,A_1) > 0$ is independent of $0 < T < T^\ast$.

	{\bfseries Step 2.2.}
	Now fix $i \in \{\,2,\,\dots,\,N\,\}$.
	By \eqnrefbr{Global-WP-Base-Boundary} we obtain
	\begin{equation*}
		q_{ii} (\partial_t c^{\scriptscriptstyle\Sigma}_i - d^{\scriptscriptstyle\Sigma}_i \Delta_\Sigma c^{\scriptscriptstyle\Sigma}_i) + \sum_{j < i} q_{ij} \big( \partial_t c^{\scriptscriptstyle\Sigma}_j - d^{\scriptscriptstyle\Sigma}_j \Delta_\Sigma c^{\scriptscriptstyle\Sigma}_j \big) = [Q r^{\textrm{sorp}}(c,\,c^{\scriptscriptstyle\Sigma})]_i + [Q r^{\textrm{ch}}(\,c^{\scriptscriptstyle\Sigma})]_i
	\end{equation*}
	on $(0,\,T) \times \Sigma$.
	According to assumption ($\textrm{A}^{\textrm{sorp}}_{\textrm{B}}$) we have
	\begin{equation*}
		[Q r^{\textrm{sorp}}(c,\,c^{\scriptscriptstyle\Sigma})]_i = \sum_{j \leq i}q_{ij} r^{\textrm{sorp}}_j(c_j,\,c^{\scriptscriptstyle\Sigma}_j) \leq \sum_{j \leq i, q_{ij}>0}q_{ij}k^\text{ad}_jc_j - \sum_{j \leq i, q_{ij}<0}q_{ij}k^\text{de}_jc^{\scriptscriptstyle\Sigma}_j
	\end{equation*}
	and according to ($\textrm{A}^{\textrm{ch}}_{\textrm{S}}$) we have $[Q r^{\textrm{ch}}(\,c^{\scriptscriptstyle\Sigma})]_i \leq C (1 + \sum^N_{j = 1} c^{\scriptscriptstyle\Sigma}_j)$
	for some fixed constant $C > 0$.
	Thus, by \Lemref{Linear-Positivity}~b) we have $0 \leq c^{\scriptscriptstyle\Sigma}_i \leq z^{\scriptscriptstyle\Sigma}_i$ for the unique maximal regular solution $z^{\scriptscriptstyle\Sigma}_i$ to
	\begin{equation*}
		\left\{\begin{array}{rclll}
			\partial_t z^{\scriptscriptstyle\Sigma}_i - d^{\scriptscriptstyle\Sigma}_i \Delta_\Sigma z^{\scriptscriptstyle\Sigma}_i & = & C^\prime \Big( 1 + \sum_{j \leq i} c_j + \sum^N_{j = 1} c^{\scriptscriptstyle\Sigma}_j \Big) \\[1.0em]
			                                                            &   & \quad - \sum_{j < i} r_{ij} \big( \partial_t c^{\scriptscriptstyle\Sigma}_j - d^{\scriptscriptstyle\Sigma}_j \Delta_\Sigma c^{\scriptscriptstyle\Sigma}_j \big) & \quad \mbox{in} & (0,\,T) \times \Sigma,          \\[0.5em]
			              - d^{\scriptscriptstyle\Sigma}_i \partial_{\nu_\Sigma} z^{\scriptscriptstyle\Sigma}_i & = & 0                                                                                & \quad \mbox{on} & (0,\,T) \times \partial \Sigma, \\[0.5em]
			                                              z^{\scriptscriptstyle\Sigma}_i(0) & = & c^{\scriptscriptstyle\Sigma}_{0,i}                                                                                & \quad \mbox{on} & \Sigma
		\end{array}\right.
	\end{equation*}
	with some appropriate constant $C^\prime = C^\prime(C,\,k^{\textrm{ad}}_j,\,k^\text{de}_j,\,q_{ij}) > 0$ and $r_{ij} = q_{ij} / q_{ii}$.
	Since $z^{\scriptscriptstyle\Sigma}_i$ is a solution to a linear problem, we may write $z^{\scriptscriptstyle\Sigma}_i = u^{\scriptscriptstyle\Sigma}_i + v^{\scriptscriptstyle\Sigma}_i$,
	where $u^{\scriptscriptstyle\Sigma}_i$ solves
	\begin{equation*}
		\left\{\begin{array}{rclll}
			\partial_t u^{\scriptscriptstyle\Sigma}_i - d^{\scriptscriptstyle\Sigma}_i \Delta_\Sigma u^{\scriptscriptstyle\Sigma}_i & = & C^\prime \Big( 1 + \sum_{j \leq i} c_j + \sum^N_{j = 1} c^{\scriptscriptstyle\Sigma}_j \Big) & \quad \mbox{in} & (0,\,T) \times \Sigma,          \\[0.5em]
			              - d^{\scriptscriptstyle\Sigma}_i \partial_{\nu_\Sigma} u^{\scriptscriptstyle\Sigma}_i & = & 0                                                                        & \quad \mbox{on} & (0,\,T) \times \partial \Sigma, \\[0.5em]
			                                              u^{\scriptscriptstyle\Sigma}_i(0) & = & 0                                                                        & \quad \mbox{on} & \Sigma
		\end{array}\right.
	\end{equation*}
	and $v^{\scriptscriptstyle\Sigma}_i$ solves
	\begin{equation*}
		\left\{\begin{array}{rclll}
			\partial_t v^{\scriptscriptstyle\Sigma}_i - d^{\scriptscriptstyle\Sigma}_i \Delta_\Sigma v^{\scriptscriptstyle\Sigma}_i & = & - \sum_{j < i} r_{ij} \big( \partial_t c^{\scriptscriptstyle\Sigma}_j - d^{\scriptscriptstyle\Sigma}_j \Delta_\Sigma c^{\scriptscriptstyle\Sigma}_j \big) & \quad \mbox{in} & (0,\,T) \times \Sigma,          \\[0.5em]
			              - d^{\scriptscriptstyle\Sigma}_i \partial_{\nu_\Sigma} v^{\scriptscriptstyle\Sigma}_i & = & 0                                                                                & \quad \mbox{on} & (0,\,T) \times \partial \Sigma, \\[0.5em]
			                                              v^{\scriptscriptstyle\Sigma}_i(0) & = & c^{\scriptscriptstyle\Sigma}_{0,i}                                                                                & \quad \mbox{on} & \Sigma.
		\end{array}\right.
	\end{equation*}
	For these solutions we have
	\begin{equation*}
		\|u^{\scriptscriptstyle\Sigma}_i\|_{L^q(\Sigma_T)} \leq M \|u^{\scriptscriptstyle\Sigma}_i\|_{{}_0 \mathbb E^\Sigma_p(T)} \leq L M C^\prime \left( \Theta + \sum_{j \leq i} \|c_j\|_{L^p(\Sigma_T)} + \sum^N_{j = 1} \|c^{\scriptscriptstyle\Sigma}_j\|_{L^p(\Sigma_T)} \right)
	\end{equation*}
	with $\Theta = \Theta(p) = \|1\|_{L^p(\Sigma_{T^\ast})}$ as well as
	\begin{equation*}
		\|(v^{\scriptscriptstyle\Sigma}_i)^+\|_{L^q(\Sigma_T)} \leq A_i \left( 1 + \sum_{j < i} \|c^{\scriptscriptstyle\Sigma}_j\|_{L^q(\Sigma_T)} \right),
	\end{equation*}
	provided that $2 \leq p,\,q < \infty$.
	Here, $M = M(p,\,q) > 0$ denotes the norm of the embedding ${}_0 \mathbb E^\Sigma_p(T) \hookrightarrow L^q(\Sigma_T)$
	and $L = L(p) > 0$ denotes the norm of the solution operator in the $L^p$-setting for the time interval $(0,\,T)$,
	which are both independent of $0 < T < T^\ast$ thanks to the homogeneous initial condition as in Step 2.1.
	Furthermore, $A_i = A_i((\|c^{\scriptscriptstyle\Sigma}_{0, j}\|_{BC(\Sigma)})_{j = 1, \dots, i},\,q) > 0$ denotes the constant delivered by \Lemref{Weak-Estimates-Boundary},
	which is also independent of $0 < T < T^\ast$.
	We again have $c^{\scriptscriptstyle\Sigma}_i \leq z^{\scriptscriptstyle\Sigma}_i \leq u^{\scriptscriptstyle\Sigma}_i + (v^{\scriptscriptstyle\Sigma}_i)^+$.
	Therefore,
	\begin{equation*}
		\begin{array}{rcl}
			\|c^{\scriptscriptstyle\Sigma}_i\|_{L^q(\Sigma_T)} & \leq & \|z^{\scriptscriptstyle\Sigma}_i\|_{L^q(\Sigma_T)} \leq \|u^{\scriptscriptstyle\Sigma}_i\|_{L^q(\Sigma_T)} + \|(v^{\scriptscriptstyle\Sigma}_i)^+\|_{L^q(\Sigma_T)} \\[1.0em]
				& \leq & C^{\prime\prime} {\displaystyle \left( 1 + \sum_{j \leq i} \|c_j\|_{L^p(\Sigma_T)} + \sum^N_{j = 1} \|c^{\scriptscriptstyle\Sigma}_j\|_{L^p(\Sigma_T)} \right) + A_i \sum_{j < i} \|c^{\scriptscriptstyle\Sigma}_j\|_{L^q(\Sigma_T)},}
		\end{array}
	\end{equation*}
	provided that $2 \leq p,\,q < \infty$.
	Here $C^{\prime\prime} = C^{\prime\prime}(L M C^\prime,\,A_i,\,\Theta) > 0$ is independent of $0 < T < T^\ast$.

	{\bfseries Step 2.3.}
	Now we may combine the estimates obtained in Steps 2.1 and 2.2, recursively, and infer that
	\begin{equation}
		\eqnlabel{Global-WP-Estimates-Boundary}
		\sum^N_{i = 1} \|c^{\scriptscriptstyle\Sigma}_i\|_{L^q(\Sigma_T)} \leq C^\ast \left( 1 + \sum^N_{j = 1} \|c_j\|_{L^p(\Sigma_T)} + \sum^N_{j = 1} \|c^{\scriptscriptstyle\Sigma}_j\|_{L^p(\Sigma_T)} \right)
	\end{equation}
	provided that $2 \leq p,\,q < \infty$.
	Here, $C^\ast = C^\ast(C^{\prime\prime},\,(A_j)_{j = 1, \dots, N}) > 0$ is independent of $0 < T < T^\ast$.

	{\bfseries Step 3.}
	Now we combine estimates (\ref{eqn:Global-WP-Estimates-Domain}) and (\ref{eqn:Global-WP-Estimates-Boundary}) to obtain
	\begin{equation}
		\eqnlabel{Global-WP-Estimates-Final}
		\begin{array}{l}
			{\displaystyle \sum^N_{i = 1} \|c_i\|_{L^q(\Omega_T)} + \sum^N_{i = 1} \|c_i\|_{L^q(\Sigma_T)} + \sum^N_{i = 1} \|c^{\scriptscriptstyle\Sigma}_i\|_{L^q(\Sigma_T)}} \\[2.0em]
				\qquad \qquad \leq {\displaystyle C^\ast \left( 1 + \sum^N_{j = 1} \|c_j\|_{L^p(\Omega_T)} + \sum^N_{j = 1} \|c_j\|_{L^p(\Sigma_T)} + \sum^N_{j = 1} \|c^{\scriptscriptstyle\Sigma}_j\|_{L^p(\Sigma_T)} \right)},
		\end{array}
	\end{equation}
	provided that $2 \leq p,\,q < \infty$.
	Here, $C^\ast > 0$ is independent of $0 < T < T^\ast$.

	Using this inequality for $p = 2$ we may in particular obtain $L^q$-$L^2$-estimates for arbitrary $2 \leq q < \infty$.

	{\bfseries Step 4.1.}
	The surface concentrations satisfy
	\begin{equation*}
		\left\{\begin{array}{rclcl}
			\partial_t c^{\scriptscriptstyle\Sigma}_i - d^{\scriptscriptstyle\Sigma}_i \Delta_\Sigma c^{\scriptscriptstyle\Sigma}_i &=& f^{\scriptscriptstyle\Sigma}_i & \text{on}& (0,T)\times \Sigma, \\
			-d^{\scriptscriptstyle\Sigma}_i\partial_{\nu_\Sigma} c^{\scriptscriptstyle\Sigma}_i &=&0 &\text{on}& (0,T)\times \partial \Sigma, \\
			{c^{\scriptscriptstyle\Sigma}_i}|_{t=0}&=&c^{\scriptscriptstyle\Sigma}_{0,i} & \text{on}& \Sigma.
		\end{array}\right.
	\end{equation*}
	with $f^{\scriptscriptstyle\Sigma}_i := r^\text{sorp}_i(c_i,c^{\scriptscriptstyle\Sigma}_i)+r^\text{ch}_i(c^{\scriptscriptstyle\Sigma})$.
	Due to the polynomial growth of the nonlinearities we may estimate $f^{\scriptscriptstyle\Sigma}$ in terms of $c$ and $c^{\scriptscriptstyle\Sigma}$
	and employ (\ref{eqn:Global-WP-Estimates-Final}) to obtain
	\begin{equation*}
		{\displaystyle \sum^N_{i = 1} \|f^{\scriptscriptstyle\Sigma}_i\|_{L^r(\Sigma_T)}}
			\leq {\displaystyle C^\ast \left( 1 + \sum^N_{j = 1} \|c_j\|_{L^p(\Omega_T)} + \sum^N_{j = 1} \|c_j\|_{L^p(\Sigma_T)} + \sum^N_{j = 1} \|c^{\scriptscriptstyle\Sigma}_j\|_{L^p(\Sigma_T)} \right)},
	\end{equation*}
	where $2 \leq p,\,r < \infty$ and $C^\ast > 0$ is independent of $0 < T < T^\ast$.
	Thus, for given $2 \leq q \leq \infty$ we may use this estimate for sufficiently large $2 \leq r < \infty$
	together with a classical result from \cite{LSU}, which yields the estimate
	\begin{equation}
		\eqnlabel{Global-WP-Estimates-Boundary-Infinity}
		{\displaystyle \sum^N_{i = 1} \|c^{\scriptscriptstyle\Sigma}_i\|_{L^q(\Sigma_T)}}
			\leq {\displaystyle C^\ast \left( 1 + \sum^N_{j = 1} \|c_j\|_{L^p(\Omega_T)} + \sum^N_{j = 1} \|c_j\|_{L^p(\Sigma_T)} + \sum^N_{j = 1} \|c^{\scriptscriptstyle\Sigma}_j\|_{L^p(\Sigma_T)} \right)}.
	\end{equation}
	Let us note that \cite[Theorem~III.7.1]{LSU} is stated for Dirichlet boundary conditions,
	but the result remains true in the Neumann case; see \cite[Theorem~4]{Bothe-Rolland:Global-Existence}, whose proof carries over to smooth manifolds as $\Sigma$.
	Note that in contrast to (\ref{eqn:Global-WP-Estimates-Boundary}) obtained in the second step,
	the estimate (\ref{eqn:Global-WP-Estimates-Boundary-Infinity}) is available for all $2 \leq p < \infty$ and all $2 \leq q \leq \infty$,
	while $C^\ast > 0$ is still independent of $0 < T < T^\ast$.

	{\bfseries Step 4.2.}
	Now we combine estimates (\ref{eqn:Global-WP-Estimates-Domain}) and (\ref{eqn:Global-WP-Estimates-Boundary-Infinity}) to obtain
	\begin{equation}
		\eqnlabel{Global-WP-Estimates-Final-Infinity}
		\begin{array}{l}
			{\displaystyle \sum^N_{i = 1} \|c_i\|_{L^q(\Omega_T)} + \sum^N_{i = 1} \|c_i\|_{L^q(\Sigma_T)} + \sum^N_{i = 1} \|c^{\scriptscriptstyle\Sigma}_i\|_{L^q(\Sigma_T)}} \\[2.0em]
				\qquad \qquad \leq {\displaystyle C^\ast \left( 1 + \sum^N_{j = 1} \|c_j\|_{L^p(\Omega_T)} + \sum^N_{j = 1} \|c_j\|_{L^p(\Sigma_T)} + \sum^N_{j = 1} \|c^{\scriptscriptstyle\Sigma}_j\|_{L^p(\Sigma_T)} \right)},
		\end{array}
	\end{equation}
	provided that $2 \leq p < \infty$ and $2 \leq q \leq \infty$.
	Here, $C^\ast > 0$ is independent of $0 < T < T^\ast$.
	Using this inequality for $p = 2$ we may in particular obtain $L^q$-$L^2$-estimates for arbitrary $2 \leq q \leq \infty$.

	{\bfseries Step 5.1.}
	The estimate (\ref{eqn:Global-WP-Estimates-Final-Infinity}) applied, for $p = 2$ and $q = \infty$, yields
	\begin{equation*}
		\begin{array}{l}
			{\displaystyle \sum^N_{i = 1} \|c_i(t)\|_{L^2(\Omega)}^2 + \sum^N_{i = 1} \|c_i(t)\|_{L^2(\Sigma)}^2 + \sum^N_{i = 1} \|c^{\scriptscriptstyle\Sigma}_i(t)\|_{L^2(\Sigma)}^2} \\[2.0em]
				\qquad \qquad \leq {\displaystyle C^\ast \left( 1 + \sum^N_{j = 1} \intop\limits^t_0 \|c_j(s)\|_{L^2(\Omega)}^2\,\mbox{d}s + \sum^N_{j = 1} \intop\limits^t_0 \|c_j(s)\|_{L^2(\Sigma)}^2\,\mbox{d}s + \sum^N_{j = 1} \intop\limits^t_0 \|c^{\scriptscriptstyle\Sigma}_j(s)\|_{L^2(\Sigma)}^2\,\mbox{d}s \right)}
		\end{array}
	\end{equation*}
	for all $0 < t < T < T^\ast$, where $C^\ast > 0$ is independent of $0 < T < T^\ast$.
	Thus, a standard Gronwall argument implies
	\begin{equation}
		\eqnlabel{Global-WP-Estimates-Gronwall}
		\|c_i\|_{L^2(\Omega_T)},\ \|c_i\|_{L^2(\Sigma_T)},\ \|c^{\scriptscriptstyle\Sigma}_i\|_{L^2(\Sigma_T)} \leq M e^{\omega T}, \qquad \qquad 0 < T < T^\ast,
	\end{equation}
	with some constants $M,\,\omega > 0$, which are independent of $0 < T < T^\ast$.

	{\bfseries Step 5.2.}
	The estimate (\ref{eqn:Global-WP-Estimates-Final-Infinity}) again applied for $p = 2$ and $q = \infty$ together with (\ref{eqn:Global-WP-Estimates-Gronwall}) yields
	\begin{equation}
		\eqnlabel{Global-WP-Estimates-Final-Gronwall}
		\|c_i\|_{L^\infty(\Omega_T)},\ \|c_i\|_{L^\infty(\Sigma_T)},\ \|c^{\scriptscriptstyle\Sigma}_i\|_{L^\infty(\Sigma_T)} \leq M e^{\omega T}, \qquad \qquad 0 < T < T^\ast,
	\end{equation}
	with some constants $M,\,\omega > 0$, which are independent of $0 < T < T^\ast$.

{\bfseries Step 6.}
Now the obtained a priori estimates (\ref{eqn:Global-WP-Estimates-Final-Gronwall}) carry over from $L^\infty$ to $H^1(\Omega)$ and $H^1(\Sigma)$. This may be seen by the following argument: Due to the $L^\infty$-estimates, the $L^2$-solution of (1) is contained in $\mathbb{E}^\Omega_p(T)\times \mathbb{E}^\Sigma_p(T)$ for each $1<p<\infty$ with $p\neq 3$ by bootstrapping. Here the crucial
estimate is
\begin{equation*}
\Vert r^\text{sorp}_i(c_i,c^{\scriptscriptstyle\Sigma}_i)\Vert_{\mathbb{G}^\Sigma_p(T)}\leq C \Vert (c_i,c^{\scriptscriptstyle\Sigma}_i)\Vert_{\mathbb{E}^\Omega_q(T)\times \mathbb{E}^\Sigma_q(T)}
\end{equation*} with suitable $q<p$; see the proof of \Thmref{Local-WP}. This, in turn, yields
\begin{equation*}
c_i\in BC(\overline{\Omega_T}), \qquad c^{\scriptscriptstyle\Sigma}_i\in BC(\overline{\Sigma_T}), \qquad (T<T^\ast).
\end{equation*} Hence, by plugging in $c$, $c^{\scriptscriptstyle\Sigma}$ into $r^\text{sorp}$, $r^\text{ch}$, we may consider (\ref{eq:cat}) again, as a linear problem, this time for data being continuous in time.
More precisely, we consider (\ref{eq:withPerturbation}) for data
\begin{align*}
f_i&\in BC([0,T^*),L^2(\Omega)),\qquad f^{\scriptscriptstyle\Sigma}_i \in BC([0,T^*),L^2(\Sigma)),\\
g^\text{in}_i&\in BC([0,T^*),L^2(\Gamma_\text{in})), \quad g^{\scriptscriptstyle\Sigma}_i\in BC([0,T^*),L^2(\Sigma)), \\ g^\text{out}_i&\in BC([0,T^*),L^2(\Gamma_\text{out})),
\end{align*} and
\begin{align*}
c_{0,i}\in H^1(\Omega),\qquad c^{\scriptscriptstyle\Sigma}_{0,i}\in H^1(\Sigma).
\end{align*} Following the strategy of the proof of our linear result in Section \ref{sec:linear_equations},
in particular by transferring Lemma \ref{lm:Surjectivity-of-Trace} to handle the inhomogeneous boundary data $g^\text{in}_i$, $g^\text{out}_i$ and $g^{\scriptscriptstyle\Sigma}_i$
and by employing a semigroup representation, we obtain that the unique solution of (\ref{eq:cat}) satisfies
\begin{align*}
c_i\in BC([0,T^*),H^1(\Omega)), \qquad c^{\scriptscriptstyle\Sigma}_i\in BC([0,T^*),H^1(\Sigma)),
\end{align*} and the corresponding a priori estimates
\begin{align}\label{eq:final-BC-estimate}
\sup_{t\in[0,T]}(\Vert c_i(t)\Vert_{H^1(\Omega)}+ \Vert c^{\scriptscriptstyle\Sigma}_i(t)\Vert_{H^1(\Sigma)}) \leq M e^{\omega T} \sup_{t\in[0,T]}\big( \Vert f_i(t)\Vert_{L^2(\Omega)}
+ \Vert f^{\scriptscriptstyle\Sigma}_i(t)\Vert_{L^2(\Sigma)}\nonumber\\
+  \Vert g^\text{in}_i(t)\Vert_{L^2(\Gamma_\text{in})} +\Vert g^\Sigma_i(t)\Vert_{L^2(\Sigma)}
+ \Vert g^\text{out}_i(t)\Vert_{L^2(\Gamma_\text{out})}  \big), \qquad (T<T^*)
\end{align} for constants $M,\omega>0$ independent of $T$. Hence, we may pass to the limit $T\rightarrow T^\ast$
and see that both sides of (\ref{eq:final-BC-estimate}) stay finite.
The proof of \Thmref{Global-WP} is now complete.
\end{proof}

\begin{remark}[Sorption and reaction examples]\label{rk:revisiting-examples}
A few remarks on the examples given in the introduction are in order here concerning the assumptions
($A_\text{F}^\text{sorp}$), ($A_\text{M}^\text{sorp}$), ($A_\text{B}^\text{sorp}$),
($A_\text{F}^\text{ch}$), ($A_\text{N}^\text{ch}$), ($A_\text{P}^\text{ch}$), ($A_\text{S}^\text{ch}$)
stated in Section \ref{sec:loc} and \Secref{Global-WP}. Evidently, Henry's law (S1) satisfies all of our assumptions. However, Langmuir's law (S2) needs to be modified in order to meet all assumptions on the sorption.
To this end we introduce $\zeta^+$, a smooth cut-off function, which approximates $(\cdot)^+$ pointwise and $\zeta^B$ a smooth and bounded function with bounded derivatives up to order $2$ which is monotonically increasing. In addition suppose $\zeta^+(0)=0$ and $\zeta^B(0)=0$. Then we consider
\begin{align*}
\tilde{r}^\text{sorp}_{L,i}(c_i,c^{\scriptscriptstyle\Sigma}_i)=k^\text{ad}_i\zeta^\text{B}(c_i)\zeta^+\left(1-\frac{\zeta^+(c^{\scriptscriptstyle\Sigma}_i)}{c^{\scriptscriptstyle\Sigma}_{\infty,i}}\right)-k^\text{de}_ic^{\scriptscriptstyle\Sigma}_i,
\end{align*} which indeed satisfies ($A_\text{F}^\text{sorp}$),
($A_\text{M}^\text{sorp}$), ($A_\text{B}^\text{sorp}$) and therefore is
covered by our main results. This modification is only necessary due to
technical reasons, since
\begin{itemize}
\item ($A_\text{F}^\text{sorp}$) is violated due to $\nabla r^\text{sorp}_{L,i} \notin BC^1(\mathbb{R}^2,\mathbb{R}^2)$.
\item ($A_\text{M}^\text{sorp}$) is not guaranteed since $c^{\scriptscriptstyle\Sigma}_i\leq c^{\scriptscriptstyle\Sigma}_{\infty, i}$ is not postulated.
\item For the same reason ($A_\text{B}^\text{sorp}$) is not satisfied since $(1-c^{\scriptscriptstyle\Sigma}_i/c^{\scriptscriptstyle\Sigma}_{\infty,i})$ could be negative.
\end{itemize}
For the time being it is not clear to the authors whether there are still global solutions in case we omit the cut-off functions. Observe that in our model there is no maximal capacity on the active surface, which
is required in the original Langmuir law (S2) to gain nonnegativity of concentrations. Nonnegativity in turn is employed in the proof of the global existence result.

The reaction rate $r^\text{ch}$ given in (R1) satisfies ($A^\text{ch}_\text{F}$), admits quadratic growth ($A^\text{ch}_\text{P}$), is quasi-positive ($A^\text{ch}_\text{N}$) and respects the triangular structure
($\textrm{A}^{\textrm{ch}}_{\textrm{S}}$) with corresponding matrix
\begin{align*}
Q=\begin{pmatrix}
1&0&0\\
0&1&0\\
1&0&1 \end{pmatrix}.
\end{align*}
\end{remark}

\bibliographystyle{plain}