\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \usepackage{subfigure} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2017 (2017), No. 194, pp. 1--26.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu} \thanks{\copyright 2017 Texas State University.} \vspace{8mm}} \begin{document} \title[\hfilneg EJDE-2017/194\hfil Stability in bioreactor processes] {Asymptotic stability of a coupled advection-diffusion-reaction system arising in bioreactor processes} \author[M. Crespo, B. Ivorra, A. M. Ramos \hfil EJDE-2017/194\hfilneg] {Mar\'ia Crespo, Benjamin Ivorra, \'Angel Manuel Ramos} \address{Mar\'ia Crespo (corresponding author) \newline UMR MISTEA - Math\'ematiques, Informatique et Statistique pour l\'Environnement et l\'Agronomie, (INRA/SupAgro). 2, Place P.Viala, 34060 Montpellier, France} \email{maria.crespo-moya@umontpellier.fr} \address{Benjamin Ivorra \newline Departamento de Matem\'atica Aplicada \& Instituto de Matem\'atica Interisciplinar, Universidad Complutense de Madrid, Plaza de Ciencias, 3, 28040 Madrid, Spain} \email{ivorra@ucm.es} \address{\'Angel Manuel Ramos \newline Departamento de Matem\'atica Aplicada \& Instituto de Matem\'atica Interisciplinar, Universidad Complutense de Madrid, Plaza de Ciencias, 3, 28040 Madrid, Spain} \email{angel@mat.ucm.es} \dedicatory{Communicated by Jesus Ildefonso Diaz} \thanks{Submitted March 1, 2017. Published August 7, 2017.} \subjclass[2010]{35B30, 35B35, 35B40, 35K51, 35Q35} \keywords{Asymptotic stability; advection-diffusion-reaction system; \hfill\break\indent dimensional analysis; continuous bioreactor} \begin{abstract} In this work, we present an asymptotic analysis of a coupled system of two advection-diffusion-reaction equations with Danckwerts boundary conditions, which models the interaction between a microbial population (e.g., bacteria), called biomass, and a diluted organic contaminant (e.g., nitrates), called substrate, in a continuous flow bioreactor. This system exhibits, under suitable conditions, two stable equilibrium states: one steady state in which the biomass becomes extinct and no reaction is produced, called washout, and another steady state, which corresponds to the partial elimination of the substrate. We use the linearization method to give sufficient conditions for the linear asymptotic stability of the two stable equilibrium configurations. Finally, we compare our asymptotic analysis with the usual asymptotic analysis associated to the continuous bioreactor when it is modeled with ordinary differential equations. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{remark}[theorem]{Remark} \newtheorem{definition}[theorem]{Definition} \newtheorem{conjecture}[theorem]{Conjecture} \allowdisplaybreaks \section{Introduction}\label{sec:1} A bioreactor is a vessel in which a microorganism (e.g., bacteria), called biomass, is used to degrade a considered diluted organic contaminant, called substrate. There exist various modes of operation in chemical reactor execution \cite{bailey1986,dochain2001}, among which continuous flow bioreactors are commonly used in the bioremediation of water resources (see, for instance, \cite{Rapaport2011,herbert1956,moreno1999}). These biological reactors are filled from a polluted resource with a flow rate $Q$ (m$^3$/s), and their output returns the treated water with the same flow rate $Q$, producing a desired quality effluent for a reasonable operating and maintenance cost. A simplified model for this process could be given by the equations \cite{SmithWaltman} \begin{equation}\label{eq:1_1} \begin{gathered} \frac{ {\rm d} S}{{\rm d} t} = -\mu(S)B + \frac{Q(t)}{V}(S_{\rm e}(t)- S) \quad t>0,\\ \frac{ {\rm d} B}{{\rm d} t} = \mu(S)B - \frac{Q(t)}{V}B \quad t>0, \end{gathered} \end{equation} where $S$ (kg/m$^3$) and $B$ (kg/m$^3$) are the concentrations of substrate and biomass, respectively; $S_{\rm e}(t)$ (kg/m$^3$) is the concentration of substrate that enters the reactor at time $t$; $V$ (m$^3$) is the reactor volume; $Q(t)$ (m$^3$/s) is the volumetric flow rate at time $t$; and $\mu(\cdot)$ (1/s) refers to the growth rate of the biomass in function of the substrate concentration. From a general point of view, due to experimental observations, we consider growth rate functions, that satisfy the following assumptions (see \cite{dochain2010,dochain2001,Harmand2015}): The function $\mu:[0,+\infty) \to [0,+\infty)$ satisfies $\mu(0)=0$, $\bar{\mu}\geq \mu(z)>0$ for $z>0$ and one of the following two properties: \begin{itemize} \item[(A1)] $\mu$ is increasing and concave. \item[(A2)] There exists $s>0$ such that $\mu$ is increasing on $(0,s)$ and decreasing on $(s,+\infty)$. \end{itemize} The Monod function \cite{SmithWaltman}, defined by $$ \mu(S)= \mu_{\rm max}\frac{S}{K_S+S}, $$ satisfies (A1), and the Haldane function \cite{Andrews}, described by $$\mu(S)= \mu^{\ast}\frac{S}{K_S+S+S^{2}/K_{\rm I}},$$ satisfies (A2). Both functions are extensively used in the literature. In the particular case when $S_{\rm e}$ and $Q$ are constant, system \eqref{eq:1_1} can be non-dimensionalized by setting $\hat{S}=S/S_{\rm e}$, $\hat{B}=B/S_{\rm e}$, $\hat{\mu}(\hat{S})= \mu(S_{\rm e}\hat{S})/\bar{\mu}$ and $\hat{t}=\bar{\mu}t$. For simplicity, we drop the $\widehat{\ }$ notation, and so $S$, $B$, $\mu$ and $t$ denote the non-dimensional variables. System \eqref{eq:1_1} in non-dimensional form is therefore given by \begin{equation}\label{eq:1_2} \begin{gathered} \frac{ {\rm d} S}{{\rm d} t} = -\mu(S)B + d(1- S) \quad \text{in }(0,+\infty),\\ \frac{ {\rm d} B}{{\rm d} t} = \mu(S)B - dB \quad \text{in }(0,+\infty),\\ \end{gathered} \end{equation} where $d= Q / V \bar{\mu}$ is the dimensionless dilution rate. If $\mu$ satisfies (A1), system \eqref{eq:1_2} has two equilibrium configurations $(S^{\ast}_1,B^{\ast}_1)=(1,0)$, usually called \textit{washout}, and $(S^{\ast}_2,B^{\ast}_2)=(S^{\ast}_2,1-S^{\ast}_2)$, where $S^{\ast}_2$ is such that $\mu(S^{\ast}_2)=d$ (see \cite{SmithWaltman}). In \cite{grivet,Rao2009,SmithWaltman} the authors conclude that the steady state $(1,0)$ is asymptotically stable if $d\geq\mu(1)$, while the steady state $(S^{\ast}_2,1-S^{\ast}_2)$ is asymptotically stable if $d<\mu(1)$. Similarly (see \cite{Andrews}), if $\mu$ satisfies (A2), system \eqref{eq:1_2} has three equilibrium configurations $(S^{\ast}_1,B^{\ast}_1)=(1,0)$, $(S^{\ast}_2,B^{\ast}_2)=(S^{\ast}_2,1-S^{\ast}_2)$ and $(S^{\ast}_{3},B^{\ast}_{3})=(S^{\ast}_{3},1-S^{\ast}_{3})$, where $\mu(S^{\ast}_2)=\mu(S^{\ast}_{3})=d$ and $S_2^{\ast}\mu(1)$, the steady state $(S^{\ast}_2,1-S^{\ast}_2)$ is asymptotically stable if $d< 1$ and the steady state $(S_{3}^{\ast},1-S_{3}^{\ast})$ is unstable. Thereby, if $\mu(1)< 1$, there is bistability when $\mu(1)< d< 1$. System \eqref{eq:1_2} describes the bioreactor dynamics under the assumption that both substrate and biomass concentrations are spatially uniform through the tank. It is of interest to consider more realistic models, for instance those based on partial differential equations, to study the influence of spatial inhomogeneities in the bioreactor (see the comparison between ordinary differential equation (ODE) and partial differential equation (PDE) bioreactor model approaches performed in \cite{aris1975,Crespo2016}). Particularly, system \eqref{eq:1_2} can be improved by considering a coupled system of spatio-temporal parabolic equations of the form \begin{equation}\label{eq:1_3} \begin{gathered} \frac{ {\rm d} S}{{\rm d} t} = L_S(S) -\mu(S)B \quad \text{in }\Omega \times (0,+\infty),\\ \frac{ {\rm d} B}{{\rm d} t} = L_B(B) + \mu(S)B \quad \text{in }\Omega \times (0,+\infty), \end{gathered} \end{equation} where $\Omega$ is the bioreactor domain and $L_S$, $L_B$ are linear second order elliptic partial differential operators on $\Omega$. Notice that system \eqref{eq:1_3} must be complemented with suitable boundary conditions that take into account the inflow--outflow balance of substrate and biomass in the bioreactor, which in the ODE system \eqref{eq:1_1} is modeled by the terms $d(1-S)$ and $-dB$, respectively. Particularly, models \eqref{eq:2_1} and \eqref{eq:2_4} presented in Section \ref{sec:2} include the so called Danckwerts boundary conditions, typically used for continuous flow systems (see \cite{Danckwerts1953,Drame2004,Drame2008}), which preserve the continuity of the substrante and biomass concentrations both at the inlet and outlet boundaries of the reactor. The asymptotic analysis of system \eqref{eq:1_3} should provide more accurate results, compared with the asymptotic ones detailed above for system \eqref{eq:1_2}, about the behavior of the substances in the bioreactor. The asymptotic behavior of parabolic equations has received a considerable attention in the literature \cite{Ildefonso2015,hale2010,Kuiper1980,liu2009,Martin1978,matano1978,matano1979}. Most theoretical studies focusing on bioreactor processes consider the assumption that both $L_S$ and $L_B$ are diffusion operators (see, e.g. \cite{Ishihara2007,Kuiper1980,morita2010,Pao1981}). For instance, in Yosida and Morita \cite{morita2010}, the authors show the existence of two different steady states (one constant, and another one spatially distributed) and develop bifurcation diagrams of the equilibrium solutions for specific model parameters. Nevertheless, such Diffusion-Reaction systems describe the behavior of batch type bioreactors, which are different to continuous flow type bioreactors, for which the addition of an advective term in operators $L_S$ and $L_B$ is required. Indeed, during batch operation no substrate is added to the initial charge and the product is not removed until the end of the process; whereas in continuous operation the substrate is continually added and the product is continually removed. The asymptotic behavior and stability analysis of advection-diffusion-reaction systems is mainly dedicated to the one-dimensional case \cite{Drame2004,Drame2008,McGowin1971,morphogenesis1995,Siebert2014,Dochain2000}. In \cite{Drame2004,Drame2008,McGowin1971,Dochain2000}, the authors study system \eqref{eq:1_3} together with Danckwerts boundary conditions under the assumption that $L_S= L_B$. Presuming that the diffusion rates of both substrate and biomass are the same, the authors discuss the asymptotic stability of the different steady states of the system. The case $L_S \neq L_B$ has been tackled in \cite{morphogenesis1995,Siebert2014}, where the authors consider periodic boundary conditions and analyze the influence of the model parameters on the stability of the different equilibrium configurations of the system. In this work, we carry out the asymptotic stability of a coupled system of two advection-diffusion-reaction equations completed with boundary conditions of mixed type, which models the interaction between substrate and biomass in a continuous flow bioreactor. We use the method of linearization to give sufficient conditions for the asymptotic stability of the two stable equilibrium configurations that the system may exhibit. In contrast to the works presented in \cite{Drame2004,McGowin1971,morphogenesis1995,Siebert2014,Dochain2000}, we consider cylindrical reactors with two spatial variables (height and radius), to study radial inhomogeneities of concentrations in the tank. We impose Danckwerts boundary conditions and allow the differential operators $L_S$ and $L_B$ to have different substrate and biomass diffusion rates. This article is organized as follows: In Section \ref{sec:2}, we introduce a PDE model describing the dynamics of the bioreactor by using a coupled system of parabolic semilinear equations together with Danckwerts boundary conditions. Additionally, we perform its dimensional analysis. In Section \ref{sec:3}, we present the steady states of the system and analyze the asymptotic stability using linearization methods. Then, Section \ref{sec:4} presents numerical experiments to analyze the validity and robustness of the stability results obtained in Section \ref{sec:3}. Finally, we perform a comparison with the asymptotic results related to system \eqref{eq:1_2}. \section{Mathematical modeling}\label{sec:2} In this section, we introduce an advection-diffusion-reaction system to model a continuous flow bioreactor and perform a dimensional analysis of this model. The bioreactor in consideration is a cylinder denoted by $\Omega^{\ast}$ (see Figure \ref{fig:1}(a)). Since this device's geometry is a solid of revolution, it can be simplified, in cylindrical coordinates, by a rectangular 2D domain, denoted by $\Omega$ and represented in Figure \ref{fig:1}(b). At the beginning of the process, there is an initial concentration of biomass in $\Omega$ that is reacting with the polluted water entering the device through the inlet $\Gamma_{\rm in}$ (i.e., the upper boundary of the rectangle $\Omega$). Treated water leaves the reactor through the outlet $\Gamma_{\rm out}$ (i.e., the lower boundary of the rectangle $\Omega$). Moreover, $\Gamma_{\rm sym} = \{0\} \times (0,H)$ is the axis of symmetry and $\Gamma_{\rm wall}= \delta \Omega \setminus (\Gamma_{\rm in} \cup \Gamma_{\rm out} \cup \Gamma_{\rm sym})$ is the bioreactor wall for which no flux passes through. \begin{figure}[ht] \begin{center} \subfigure[3D reactor]{\includegraphics[width=0.45\textwidth]{fig1a}} \subfigure[2D simplification]{\includegraphics[width=0.45\textwidth]{fig1b}} \end{center} \caption{Typical representation of the domain geometry.} \label{fig:1} \end{figure} By using cylindrical coordinates $(r,z)$, where $r$ is the distance to the symmetrical cylinder axis, we consider the following system describing the behavior of the continuous bioreactor \cite{crespo2015}: \begin{equation}\label{eq:2_1} \begin{gathered} \frac{\partial S}{\partial t}= \frac{1}{r} \frac\partial{\partial r} ( r D_S \frac{\partial S}{\partial r}) + \frac\partial{\partial z} (D_S \frac{\partial S}{\partial z}) + u \frac{\partial S}{\partial z}- \mu(S) B \quad \text{in } \Omega \times (0,T), \\ \frac{\partial B}{\partial t}= \frac{1}{r} \frac\partial{\partial r} ( r D_B \frac{\partial B}{\partial r}) + \frac\partial{\partial z} (D_B \frac{\partial B}{\partial z}) + u \frac{\partial B}{\partial z} + \mu(S) B \quad \text{in } \Omega \times (0,T), \\ D_S \frac{\partial S}{\partial z} + u S= u S_{\rm e} \quad \text{in } \Gamma_{\rm in} \times (0,T), \\ D_B \frac{ \partial B}{\partial z} + u B = 0 \quad \text{in } \Gamma_{\rm in} \times (0,T), \\ D_S\frac{\partial S}{\partial r} = 0 \quad \text{in } \Gamma_{\rm sym} \times (0,T), \\ D_B\frac{\partial B}{\partial r}=0 \quad \text{in } \Gamma_{\rm sym} \times (0,T), \\ D_S\frac{\partial S}{\partial r} = 0 \quad \text{in } \Gamma_{\rm wall} \times (0,T), \\ D_B\frac{\partial B}{\partial r}=0 \quad \text{in } \Gamma_{\rm wall} \times (0,T), \\ D_S\frac{\partial S}{\partial z} = 0 \quad \text{in } \Gamma_{\rm out} \times (0,T), \\ D_B\frac{\partial B}{\partial z}=0 \quad \text{in } \Gamma_{\rm out} \times (0,T),\\ S(\cdot,\cdot,0)=S_0 \quad \text{in }\Omega, \\ B(\cdot,\cdot,0)=B_0 \quad \text{in } \Omega, \end{gathered} \end{equation} where $T>0$ (s) is the length of the time interval for which we want to model the process; $S$ (kg/m$^{3}$) and $B$ (kg/m$^{3}$) are the substrate and biomass concentrations inside the bioreactor, which diffuse throughout the water in the vessel with diffusion coefficients $D_S$ (m$^{2}$/s) and $D_B$ (m$^{2}$/s), respectively; the fluid velocity is taken as $\mathbf{u}=(0,0,-u)$, where $u$ (m/s) is the flow speed; $S_{\rm e}$ (kg/m$^{3}$) is the concentration of substrate that enters into the bioreactor; $S_0$ (kg/m$^3$) and $B_0$ (kg/m$^{3}$) are the initial concentrations of substrate and biomass inside the bioreactor, respectively. Furthermore, as in system \eqref{eq:1_1}, we consider a term corresponding to the reaction between biomass and substrate, governed by the growth rate function $\mu$ (s$^{-1}$). \begin{remark}\label{rmk:2_1}\rm According to \cite{crespo2015}, if $\mu \in L^{\infty}(0,+\infty)$ is continuous and Lipschitz, $u \in L^{\infty}(\bar{\Omega} \times (0,T))$, $S_{\rm e} \in L^{\infty}(0,T)$, $S_{\rm e}\geq 0$ in $(0,T)$, $S_0 \in L^{\infty}(\Omega)$, $S_0\geq 0$ in $\Omega$, $B_0 \in L^{\infty}(\Omega)$ and $B_0\geq 0$ in $\Omega$, there exists a unique solution $(S,B) \in L^{2}(0,T,H^{1}(\Omega))^2 \cap \mathcal{C}([0,T],L^{2}(\Omega))^{2} \cap L^{\infty}(\Omega\times(0,T))^{2}$ of system \eqref{eq:2_1}. \end{remark} We perform a dimensional analysis of the model, following similar methodology as the one presented in \cite{Smith2014}. System \eqref{eq:2_1} is non-dimensionalized by setting $\hat{B}= B/b$, $\hat{S}= S/s$, $\hat{t}=t/\tau$, $\hat{u}= u/\gamma$, $\hat{S}_{\rm e}= S_{\rm e}/e$, $\hat{\mu}(\hat{S})= \mu(s\hat{S})/\nu$, $\hat{z}= z/Z$ and $\hat{r}= r/R$, where $b$, $s$, $\tau$, $\gamma$, $e$, $\nu$, $Z$ and $R$ are suitable scales. Thus, for $0\leq \hat{t} \leq \hat{T}=T /\tau$ and $(\hat{r},\hat{z})\in \hat{\Omega}$ (the nondimensional domain obtained from $\Omega$ with the change of variables $(\hat{r},\hat{z})=(r/R,z/Z)$ the first and second equations in system \eqref{eq:2_1} become \begin{equation}\label{eq:2_2} \frac{\partial \hat{S}}{\partial \hat{t}} = \frac{ \tau D_S}{R^{2}\hat{r}} \frac\partial{\partial \hat{r}} (\hat{r} \frac{\partial \hat{S}}{\partial \hat{r}}) + \frac{\tau D_S}{Z^{2}} \frac{\partial^{2} \hat{S}}{\partial \hat{z}^{2}} + \frac{\gamma \tau}{Z} \hat{u} \frac{\partial \hat{S}}{\partial \hat{z}} - \frac{b\tau\nu}{s} \hat{\mu}(\hat{S}) \hat{B} \end{equation} and \begin{equation}\label{eq:2_3} \frac{\partial \hat{B}}{\partial \hat{t}} = \frac{ \tau D_B}{R^{2}\hat{r}} \frac\partial{\partial \hat{r}} ( \hat{r} \frac{\partial \hat{B}}{\partial \hat{r}}) + \frac{\tau D_B}{Z^{2}} \frac{\partial^{2} \hat{B}}{\partial \hat{z}^{2}} + \frac{\gamma \tau}{Z} \hat{u} \frac{\partial \hat{B}}{\partial \hat{z}} + \tau\nu \hat{\mu}(\hat{S}) \hat{B}. \end{equation} The dimensionless groups of parameters in equations \eqref{eq:2_2} and \eqref{eq:2_3} are $\alpha_1= \tau D_S/R^{2}$, $\alpha_2= \tau D_S/Z^{2}$, $\alpha_{3}= \tau D_B/R^{2}$, $\alpha_{4}= \tau D_B/Z^{2}$, $\alpha_{5}= \tau \gamma / Z$, $\alpha_{6}= \tau \nu$ and $\alpha_{7}= \tau \nu b /s$. The radius and the height scales proposed here come from the dimensions of the bioreactor, giving $R=L$ and $Z=H$. We set $\nu=\bar{\mu}$ and $\gamma= \|u\|_{L^{\infty}(\bar{\Omega} \times (0,T))}$ for the reaction and velocity scales, respectively. Finally, for the entering substrate scale we set $e= \|S_{\rm e}\|_{L^{\infty}(0,T)}$ and, for the sake of simplicity, we choose $s=b=\|S_{\rm e}\|_{L^{\infty}(0,T)}$. The time scale $\tau$ is chosen from equations \eqref{eq:2_2} and \eqref{eq:2_3} depending on the process (diffusion, advection or reaction) we want to focus on. In particular, we can choose $$ \tau \in \{ L^{2}/D_S, L^{2}/D_B, H^{2}/D_S, H^{2}/D_B, H/\|u\|_{L^{\infty}(\Omega \times (0,T))}, 1/\bar{\mu} \}, $$ where $\tau=L^{2}/D_S$ (resp., $\tau=H^{2}/D_S$) corresponds to the case focusing on the substrate diffusion rate on the horizontal (resp., vertical) axis; $\tau= L^{2}/D_B$ (resp., $\tau=H^{2}/D_B$) focuses on the biomass diffusion rate on the horizontal (resp., vertical) axis; $\tau=H/\|u\|_{L^{\infty}(\bar{\Omega} \times (0,T))}$ focuses on the advection transport rate; and $\tau=1/\bar{\mu}$ focuses on the reaction rate. Since in next sections we perform a comparison with system \eqref{eq:1_2}, we center our study on the reaction process and take $\tau= 1/\bar{\mu}$. Two well-known dimensionless numbers (see \cite{marin2005advances}) appear now in the non-dimensional form of system \eqref{eq:2_1}: \begin{itemize} \item Damkh\"oler Number: ${\rm Da}=\frac{\text{reaction rate}}{\text{advective transport rate}} = \frac{\tau_{\rm a}}{\tau_{\rm r}}$, \item Thiele Modulus: ${\rm Th}= \frac{\text{reaction rate}} {\text{diffusive transport rate}}= \frac{\tau_{\rm d}}{\tau_{\rm r}}$, \end{itemize} where $\tau_{\rm d}$, $\tau_{\rm a}$ and $\tau_{\rm r}$ are diffusion, advection and reaction times scales, respectively. For ease of notation, we drop the $\widehat{\ }$ symbol, and so $B$, $S$, $t$, $u$, $S_{\rm e}$, $\mu$, $z$, $r$ and $T$ denote now the non-dimensional variables. Particularly, if $S_{\rm e}$ and $u$ are constants, system \eqref{eq:2_1} in its non-dimensional form is given by \begin{equation}\label{eq:2_4} \begin{gathered} \frac{\partial S}{\partial t}= \frac{\sigma^{2}}{{\rm Th}_S} \frac{1}{r} \frac\partial{\partial r} ( r \frac{\partial S}{\partial r}) + \frac{1}{{\rm Th}_S} \frac{\partial^{2} S}{\partial z^{2}} + \frac{1}{\rm Da} \frac{\partial S}{\partial z} - \mu(S) B \quad \text{in } \Omega \times (0,T), \\ \frac{\partial B}{\partial t} = \frac{\sigma^{2}}{{\rm Th}_B} \frac{1}{r} \frac\partial{\partial r} ( r \frac{\partial B}{\partial r}) + \frac{1}{{\rm Th}_B} \frac{\partial^{2} B}{\partial z^{2}} + \frac{1}{\rm Da} \frac{\partial B}{\partial z} + \mu(S) B \quad \text{in } \Omega \times (0,T), \\ \frac{1}{{\rm Th}_S} \frac{\partial S}{\partial z} + \frac{1}{\rm Da} S= \frac{1}{\rm Da} \quad \text{in } \Gamma_{\rm in} \times (0,T), \\ \frac{1}{{\rm Th}_B} \frac{ \partial B}{\partial z} + \frac{1}{\rm Da} B = 0 \quad \text{in } \Gamma_{\rm in} \times (0,T), \\ \frac{\partial S}{\partial r} = \frac{\partial B}{\partial r} = 0 \quad \text{in } \Gamma_{\rm sym} \times (0,T), \\ \frac{\partial S}{\partial r} = \frac{\partial B}{\partial r} = 0 \quad \text{in } \Gamma_{\rm wall} \times (0,T), \\ \frac{\partial S}{\partial z} = \frac{\partial B}{\partial z} = 0 \quad \text{in } \Gamma_{\rm out} \times (0,T), \end{gathered} \end{equation} complemented by the initial conditions \begin{equation}\label{eq:2_5} S(\cdot,\cdot,0)=S_{\rm init} \quad\text{and} \quad B(\cdot,\cdot,0)=B_{\rm init} \quad \text{in }\Omega, \end{equation} where $\Omega= (0,1)\times(0,1)$ is the nondimensional domain, $\Gamma_{\rm in}=(0,1)\times \{1\}$, $\Gamma_{\rm out}=(0,1)\times \{0\}$, $\Gamma_{\rm wall}=\{1\}\times (0,1)$ and $\Gamma_{\rm sym}=\{0\}\times (0,1)$ are the non-dimensional boundary edges. The final dimensionless parameters are ${\rm Da}= H\bar{\mu}/u$, ${\rm Th}_S= H^{2} \bar{\mu}/D_S$, ${\rm Th}_B= D_S {\rm Th}_S/D_B$ and $\sigma=H/L$, and the dimensionless initial conditions are $S_{\rm init}=S_0/S_{\rm e}$ and $B_{\rm init}=B_0/S_{\rm e}$. \begin{remark}\label{rmk:2_2} \rm Since the bioreactor in consideration is a cylinder of height $H$ and radius $L$, the reactor volume is $\pi H L^{2}$ and the volumetric flow rate in system \eqref{eq:1_1} can be written as $Q= \pi L^{2} u$, where $u$ (m/s) is the vertical inflow. Thus, the nondimensional dilution rate $d$ in system \eqref{eq:1_2} corresponds to the nondimensional flow rate $(1/{\rm Da})$ in system \eqref{eq:2_4}. \end{remark} \section{Steady states and stability analysis}\label{sec:3} In this section, we study the asymptotic behavior of system \eqref{eq:2_4}-\eqref{eq:2_5}. Firstly, we study the particular case for which diffusion terms in system \eqref{eq:2_4} are neglected. Then, we perform the stability analysis of system \eqref{eq:2_4} for the general case. The asymptotic stability of an equilibrium solution of system \eqref{eq:2_4} is defined as follows (see \cite{pao1992}). \begin{definition}\label{defn:3_1}\rm An equilibrium solution $(S^{\ast},B^{\ast})$ of system \eqref{eq:2_4} is said to be asymptotically stable if there exists $\epsilon >0$ such that if given $(S_{\rm init},B_{\rm init})\in (L^{\infty}(\Omega))^{2}$ satisfying $$ \| (S_{\rm init},B_{\rm init}) - (S^{\ast},B^{\ast})\|_{(L^{2}(\Omega))^2} < \epsilon, $$ then the corresponding unique solution $(S,B)$ of system \eqref{eq:2_4}-\eqref{eq:2_5} satisfies $$ \lim_{t\to \infty} \|(S(t),B(t)) - (S^{\ast},B^{\ast}))\|_{(L^{2}(\Omega))^2}= 0. $$ We point out that we require $(S_{\rm init},B_{\rm init}) \in (L^{\infty}(\Omega))^{2}$ to ensure the existence and uniqueness of solution of system \eqref{eq:2_4}-\eqref{eq:2_5} (see Remark \ref{rmk:2_1}). \end{definition} \subsection{Case {$(1/{\rm Th}_S)$, $(1/{\rm Th}_B)$, $(\sigma^2/{\rm Th}_S)$, $(\sigma^2/{\rm Th}_B)$ $\ll 1$}{Lg}.} \label{sec:3_1} We consider the particular case where the nondimensional diffusion coefficients are negligible with respect to the advection and reaction coefficients in system \eqref{eq:2_4}. For each fixed value of $r \in (0,1)$, the solution $S(r,\cdot)$, $B(r,\cdot)$ can be approximated by the solution of the following 1-dimensional advection reaction system: \begin{equation}\label{eq:3_1} \begin{gathered} \frac{\partial S}{\partial t}= \frac{1}{\rm Da} \frac{\partial S}{\partial z} - \mu(S) B \quad \text{in }(0,1)\times (0,T), \\ \frac{\partial B}{\partial t}= \frac{1}{\rm Da} \frac{\partial B}{\partial z} + \mu(S) B \quad \text{in }(0,1)\times (0,T), \\ S(r,1,t)= 1 \quad \forall t \in (0,T), \\ B(r,1,t) = 0 \quad \forall t \in (0,T), \\ S(r,\cdot,0)=S_{\rm init}(r,\cdot) \quad \text{in } (0,1),\\ B(r,\cdot,0)=B_{\rm init}(r,\cdot) \quad \text{in } (0,1). \end{gathered} \end{equation} Let us prove that $(1,0)$ (which is called the \textit{washout} state), is an asymptotically stable equilibrium. The following theorem shows, in fact, a property for $(1,0)$ stronger than asymptotic stability. \begin{theorem}\label{thm:3_2} For any arbitrary initial condition $(S_{\rm init},B_{\rm init}) \in (L^{\infty}(\Omega))^2$, the solution of \eqref{eq:3_1} satisfies that $S(r,z,t)$= 1 and $B(r,z,t)=0$, for all $(r,z)\in \Omega$ and $t\geq {\rm Da}$. \end{theorem} \begin{proof} We detail here this easy proof for the clarity of the reading. For any fixed value of $r \in (0,1)$, we apply the Euler-Lagrange transformation from $(r,z,t)$ to $(r,\tilde{z}(t,z),t)$, where $\tilde{z}(t,z)= z- (t/{\rm Da})$, so that for every fixed value of $(r,z)\in \Omega$, the second equation of system \eqref{eq:3_1} is rewritten as \begin{align*} \frac{{\rm d} B}{{\rm d} t}(r,\tilde{z}(t,z),t) & = \frac{\partial B}{\partial t}(r,\tilde{z}(t,z),t) - \frac{1}{{\rm Da}} \frac{\partial B}{\partial \tilde{z}} (r,\tilde{z}(t,z),t) \\ &= \mu(S (r,\tilde{z}(t,z),t))B(r,\tilde{z}(t,z),t). \end{align*} Thus, for any $(r,z)\in \Omega$, one has that $$ B(r,\tilde{z}(t,z),t)= B(r,\tilde{z}(0,z),0) + \int_0^{t} \mu(S (r,\tilde{z}(\tau,z),\tau))B(\tilde{z}(\tau,z),\tau){\rm d} \tau. $$ Particularly, for $z=1$, we obtain $$ B(r,\tilde{z}(t,1),t)= \int_0^{t} \mu(S (r,\tilde{z}(\tau,1),\tau)) B(r,\tilde{z}(\tau,1),\tau){\rm d} \tau $$ and, by applying the Gronwall's inequality, we have that $B(r,\tilde{z}(t,1),t)=0$ for all $t>0$. Using the same reasoning for the first equation of system \eqref{eq:3_1}, it follows that for $z=1$ $$ S(r,\tilde{z}(t,1),t)= 1+ \int_0^{t} \mu(S (r,\tilde{z}(\tau,1),\tau)) B(r,\tilde{z}(\tau,1),\tau){\rm d} \tau. $$ Since $B(r,\tilde{z}(t,1),t)=0$ for all $t>0$, we deduce that $S(r,\tilde{z}(t,1),t)=1$ for all $t>0$. Coming back to Eulerian coordinates, one has that $$ B(r,1-t/{\rm Da},t)=0 \text{ and } S(r,1-t/{\rm Da},t)=1 \text{ for all } t>0. $$ Consequently, if $t\geq {\rm Da}$, $B(r,z,t)=0$ and $S(r,z,t)=1$ for all $(r,z)\in \Omega$. \end{proof} \subsection{General Case}\label{sec:3_2} To obtain a parallelism with the asymptotic analysis of system \eqref{eq:1_2}, shown in Section \ref{sec:1}, we assume that $\mu$ satisfies properties (A1) or (A2). In both cases, the constant (washout) solution $(S^{\ast}_1,B^{\ast}_1)= (1,0)$ is a steady state of system \eqref{eq:2_4}. By analogy with system \eqref{eq:1_2}, we conjecture, supported by numerical experiments, that system \eqref{eq:2_4} has, under suitable conditions, another asymptotically stable steady state (different from the washout) denoted by $(S_2^{\ast},B_2^{\ast})$. In this section, we use linearization techniques to study the stability of the steady states. First, we introduce the concept of \textit{linear asymptotic stability}. Then, we give a sufficient condition for the linear asymptotic stability of the washout equilibrium. Finally, we use this result to infer a sufficient condition for the linear asymptotic stability of the other equilibrium solution. Let $(S^{\ast},B^{\ast})$ be an equilibrium solution of system \eqref{eq:2_4}. We define the concept of \textit{linear asymptotic stability} by following the reasoning below. Let us consider system \eqref{eq:2_4} with initial conditions in $(L^{\infty}(\Omega))^{2}$ and close to $(S^{\ast},B^{\ast})$ given by $S(r,z,0)=S^{\ast}+\delta S_{\rm init}\geq 0$ and $B(r,z,0)=B^{\ast}+\delta B_{\rm init}\geq 0$, with $\| \delta S_{\rm init}\|_{L^{2}(\Omega)}\ll 1$ and $\| \delta B_{\rm init}\|_{L^{2}(\Omega)}\ll 1$. The solution of system \eqref{eq:2_4} can be seen as \begin{equation}\label{eq:3_2} \begin{pmatrix} S(r,z,t) \\ B(r,z,t) \end{pmatrix} = \begin{pmatrix} S^{\ast}(r,z) \\ B^{\ast}(r,z) \end{pmatrix} + \begin{pmatrix} \delta S(r,z,t) \\ \delta B(r,z,t) \end{pmatrix} + \text{Higher order terms}, \end{equation} (see \cite[page 220]{Fischer1975}) where \begin{equation}\label{eq:3_3} \begin{gathered} \begin{aligned} \frac{{\rm d} \delta S}{{\rm d} t} &= \frac{\sigma^{2}}{{\rm Th}_S} \frac{1}{r} \frac{\rm d}{{\rm d} r} ( r \frac{{\rm d} \delta S}{{\rm d} r}) + \frac{1}{{\rm Th}_S} \frac{{\rm d}^{2} \delta S}{{\rm d} z^{2}} + \frac{1}{\rm Da} \frac{{\rm d} \delta S}{{\rm d} z} \\ &\quad - \mu(S^{\ast}) \delta B-\mu'(S^{\ast})B^{\ast}\delta S \quad \text{in }\Omega \times (0,T), \end{aligned} \\ \begin{aligned} \frac{{\rm d} \delta B}{{\rm d} t} &= \frac{\sigma^{2}}{{\rm Th}_B} \frac{1}{r} \frac{\rm d}{{\rm d} r} ( r \frac{{\rm d} \delta B}{{\rm d} r}) + \frac{1}{{\rm Th}_B} \frac{{\rm d}^{2} \delta B}{{\rm d} z^{2}} + \frac{1}{\rm Da} \frac{{\rm d} \delta B}{{\rm d} z} \\ &\quad + \mu(S^{\ast}) \delta B+\mu'(S^{\ast})B^{\ast}\delta S \quad \text{in }\Omega \times (0,T), \end{aligned} \\ \frac{1}{{\rm Th}_S}\frac{{\rm d} \delta S}{{\rm d} z} + \frac{1}{\rm Da}\delta S = \frac{1}{{\rm Th}_B} \frac{ {\rm d} \delta B}{{\rm d} z} + \frac{1}{\rm Da}\delta B = 0 \quad \text{in }\Gamma_{\rm in} \times (0,T), \\ \frac{{\rm d} \delta S}{{\rm d} r} = \frac{{\rm d} \delta B}{{\rm d} r} =0 \quad \text{in } \Gamma_{\rm sym} \times (0,T), \\ \frac{{\rm d} \delta S}{{\rm d} r} = \frac{{\rm d} \delta B}{{\rm d} r} = 0 \quad \text{in } \Gamma_{\rm wall} \times (0,T), \\ \frac{{\rm d} \delta S}{{\rm d} z} = \frac{{\rm d} \delta B}{{\rm d} z} = 0 \quad \text{in }\Gamma_{\rm out} \times (0,T), \\ \delta S(\cdot,\cdot,0)=\delta S_{\rm init} \quad \text{in }\Omega, \\ \delta B(\cdot,\cdot,0)=\delta B_{\rm init} \quad \text{in }\Omega. \end{gathered} \end{equation} \begin{definition}\label{defn:3_3}\rm An equilibrium solution $(S^{\ast},B^{\ast})$ of system \eqref{eq:2_4} is said to be linearly asymptotically stable if there exists $\epsilon >0$ such that if given $(\delta S_{\rm init},\delta B_{\rm init})\in (L^{\infty}(\Omega))^{2}$ satisfying $\| (\delta S_{\rm init},\delta B_{\rm init})\|_{(L^{2}(\Omega))^2} < \epsilon$, then the corresponding unique solution $(\delta S,\delta B)$ of system \eqref{eq:3_3} satisfies $\lim_{t\to \infty} \|(\delta S(t),\delta B(t))\|_{(L^{2}(\Omega))^2}= 0$. \end{definition} We now define the following functions, which will be used through the rest of the manuscript. \begin{definition}\label{defn:3_4} \rm $\bullet$ In terms of the dimensionless variables appearing in system \eqref{eq:2_4}, we define $\beta_1({\rm Da},{\rm Th}_B)$ as the smallest positive solution of the transcendental equation $\tan (\beta) = \frac{{\rm Th_{B}}\beta}{{\rm Da}\big(\beta^{2} - (\frac{{\rm Th}_B}{2{\rm Da}}\big)^2)}$ if ${\rm Th_{B}}\neq \pi {\rm Da}$. If ${\rm Th_{B}}=\pi {\rm Da}$, we define $\beta_1({\rm Da},{\rm Th}_B)=\pi/2$. $\bullet$ In terms of the variables with dimensions appearing in system \eqref{eq:2_1}, we define $\tilde{\beta}_1(H,u,D_B)$ as the smallest positive solution of the transcendental equation \[ \tan (\beta) = \frac{Hu\beta}{D_B\big(\beta^{2} - (\frac{Hu}{2D_B})^2\big)} \] if $Hu \neq \pi D_B$. If $Hu=\pi D_B$ we define $\tilde{\beta}_1(H,u,D_B)=\pi/2$. \end{definition} \begin{theorem}\label{thm:3_5} A sufficient condition for $(S^{\ast}_1,B^{\ast}_1)=(1,0)$ to be a linearly asymptotically stable steady state of system \eqref{eq:2_4} is that \begin{equation}\label{eq:3_4} \mu(1) < \frac{{\rm Th}_B}{(2{\rm Da})^{2}} + \frac{(\beta_1({\rm Da},{\rm Th_{B}}))^{2}}{{\rm Th_{B}}}. \end{equation} \end{theorem} \begin{remark}\rm In terms of the variables with dimensions appearing in system \eqref{eq:2_1}, the steady state is $(S_{\rm e},0)$ and inequality \eqref{eq:3_4} is reformulated as \begin{equation*} \mu(S_{\rm e}) < \frac{u^{2}}{4D_B} + \frac{D_B}{H^{2}}(\tilde{\beta}_1(H,u,D_B))^{2}. \end{equation*} \end{remark} \begin{proof}[Proof of Theorem \ref{thm:3_5}] To check the stability of the washout equilibrium solution, we replace $(S^{\ast},B^{\ast})$ by $(1,0)$ in System \eqref{eq:3_3} and, so, functions $\delta S$ and $\delta B$ fulfill \begin{equation}\label{eq:3_5} \begin{gathered} \frac{{\rm d} \delta S}{{\rm d} t} = \frac{\sigma^{2}}{{\rm Th}_S} \frac{1}{r} \frac{\rm d}{{\rm d} r} ( r \frac{{\rm d} \delta S}{{\rm d} r}) + \frac{1}{{\rm Th}_S} \frac{{\rm d}^{2} \delta S}{{\rm d} z^{2}} + \frac{1}{\rm Da} \frac{{\rm d} \delta S}{{\rm d} z} - \mu(1) \delta B \\ \text{in }\Omega \times (0,T), \\ \frac{{\rm d} \delta B}{{\rm d} t} = \frac{\sigma^{2}}{{\rm Th}_B} \frac{1}{r} \frac{\rm d}{{\rm d} r} ( r \frac{{\rm d} \delta B}{{\rm d} r}) + \frac{1}{{\rm Th}_B} \frac{{\rm d}^{2} \delta B}{{\rm d} z^{2}} + \frac{1}{\rm Da} \frac{{\rm d} \delta B}{{\rm d} z} + \mu(1) \delta B \\ \text{in }\Omega \times (0,T), \\ \frac{1}{{\rm Th}_S}\frac{{\rm d} \delta S}{{\rm d} z} + \frac{1}{\rm Da}\delta S= \frac{1}{{\rm Th}_B} \frac{ {\rm d} \delta B}{{\rm d} z} + \frac{1}{\rm Da}\delta B= 0 \quad \text{in }\Gamma_{\rm in} \times (0,T), \\ \frac{{\rm d} \delta S}{{\rm d} r} = \frac{{\rm d} \delta B}{{\rm d} r} = 0 \quad \text{in } \Gamma_{\rm sym} \times (0,T), \\ \frac{{\rm d} \delta S}{{\rm d} r} = \frac{{\rm d} \delta B}{{\rm d} r} = 0 \quad \text{in } \Gamma_{\rm wall} \times (0,T), \\ \frac{{\rm d} \delta S}{{\rm d} z} = \frac{{\rm d} \delta B}{{\rm d} z} = 0 \quad \text{in }\Gamma_{\rm out} \times (0,T), \\ \delta S(\cdot,\cdot,0)=\delta S_{\rm init} \quad \text{in }\Omega,\\ \delta B(\cdot,\cdot,0)=\delta B_{\rm init} \quad \text{in } \Omega, \end{gathered} \end{equation} with $\| \delta S_{\rm init}\|_{L^{2}(\Omega)}\ll 1$ and $\| \delta B_{\rm init}\|_{L^{2}(\Omega)}\ll 1$. We are going to prove that the steady state $(S_1^{\ast},B_1^{\ast})=(1,0)$ is linearly asymptotically stable by showing that (see Definition \ref{defn:3_3}) \[ \|\delta S(t)\|_{L^{2}(\Omega)}\to 0 \quad \text{and}\quad \|\delta B(t)\|_{L^{2}(\Omega)} \to 0 \text{ as }t \to \infty. \] \smallskip \noindent\textbf{Step 1.} Let us prove that $\|\delta B(t)\|_{L^{2}(\Omega)} \to 0 \text{ as }t \to \infty$: Notice that the equations involving the biomass in system \eqref{eq:3_5} are decoupled from those involving the substrate, and may be solved by separation of variables by imposing $$ \delta B(r,z,t)= R(r)Z(z)T(t). $$ \smallskip \noindent\textbf{Step 1.1.} Separation of variables. From the second equation in system \eqref{eq:3_5} one has $$ \frac{T'(t)}{T(t)} = \frac{\sigma^{2}}{{\rm Th}_B}\Big( \frac{R''(r)}{R(r)} + \frac{1}{r}\frac{R'(r)}{R(r)}\Big) + \frac{1}{{\rm Th}_B}\frac{Z''(z)}{Z(z)} + \frac{1}{{\rm Da}}\frac{Z'(z)}{Z(z)} + \mu(1). $$ If we equate this expression to a constant $\lambda$, it follows that \begin{gather*} T'(t)-\lambda T(t)=0 ,\\ \frac{\sigma^{2}}{{\rm Th}_B}\Big( \frac{R''(r)}{R(r)} + \frac{1}{r}\frac{R'(r)}{R(r)}\Big) = - \frac{1}{{\rm Th}_B}\frac{Z''(z)}{Z(z)} - \frac{1}{{\rm Da}}\frac{Z'(z)}{Z(z)} +\lambda- \mu(1). \end{gather*} Equating this expression to an arbitrary constant $\eta$, one obtains \begin{gather*} R''(r) + \frac{1}{r}R'(r)- \frac{{\rm Th}_B}{\sigma^{2}}\eta R(r)=0 , \\ \frac{1}{{\rm Th}_B}Z''(z) + \frac{1}{{\rm Da}}Z'(z) -\big(\lambda- \mu(1)-\eta\big)Z(z)=0. \end{gather*} Proceeding as in the proof of \cite[Theorem 3]{crespo2015}, it is easy to see that $$ \delta B(r,z,t)=|\delta B(r,z,t)| \leq \|\delta B_{\rm init}\|_{L^{\infty}(\Omega)}{\rm e}^{\mu(1)t} \quad \text{for a.e. } (r,z,t)\in \Omega \times (0,T). $$ Particularly, the function $R:[0,1]\to \mathbb{R}$ must be bounded in $(0,1)$ (this fact will be used in the step 1.2 of this proof). \smallskip \noindent\textbf{Step 1.2.} Calculation of $R(r)$. Using the boundary conditions of system \eqref{eq:3_5} on $\Gamma_{\rm wall}$ and $\Gamma_{\rm sym}$, it is clear that $R(r)$ is a solution of system \begin{equation}\label{eq:3_6} \begin{gathered} R''(r) + \frac{1}{r}R'(r)- \frac{{\rm Th}_B}{\sigma^{2}}\eta R(r)=0 \quad r\in(0,1),\\ R'(0)=R'(1)=0. \end{gathered} \end{equation} Using the change of variables $s=ar$, with $a=\sqrt{|\eta|\frac{{\rm Th}_B}{\sigma^{2}}}$, the differential equation for $R$ can be rewritten in one of the following forms \begin{enumerate} \item $s^{2}R''(s)+sR'(s) + s^{2}R(s)=0$ if $\eta<0$, \item $s^{2}R''(s)+sR'(s) - s^{2}R(s)=0$ if $\eta>0$, \item $s R''(s)+R'(s)=0$ if $\eta=0$. \end{enumerate} \smallskip \noindent\textbf{Case 1: $\eta<0$}. In this case the equation for $R(s)$ is known as the \textit{Bessel equation of order zero}, with general solution $$ R(s)= C_1 J_0(s) + C_2 Y_0(s), $$ where $C_1$, $C_2 \in \mathbb{R}$ and $J_n$ and $Y_n$ are, respectively, the Bessel functions of first and second kind of order $n$. Since $Y_0$ has a singularity at $s=0$, to ensure that function $R(s)$ is bounded, $C_2$ must be zero, and consequently, $R(s)=C_1J_0(s)$. It is well known that $J_0'(s)=-J_1(s)$ and $0\in \{ s\in [0,+\infty)\text{: }J_1(s)=0\}$, which is a countable set $\{ T_n\}_{n \in \mathbb{N}}$ with an infinite number of elements (see, e.g., \cite{asmar2005}). Therefore, $R'(0)=0$ is always satisfied and from the boundary condition at $s=a$ ($r=1$), one has that the eigenvalues $\eta$ are such that $J_0'\big(\sqrt{-\frac{\eta {\rm Th}_B}{\sigma^{2}}}\big)=0$. Consequently, $\eta \in \{\eta_n\}_{n\in \mathbb{N}}$, with \begin{equation}\label{eq:3_7} \eta_n= - \frac{(\sigma T_n)^{2}}{{\rm Th}_B}, \end{equation} and the solution $R(r)$ is given by $$ R(r)=\sum_{n\in \mathbb{N}} C_n J_0 \big(\frac{\sqrt{-\eta_n{\rm Th}_B}}{\sigma}r \big). $$ \smallskip \noindent\textbf{Case 2: $\eta>0$}. In this case the equation for $R(s)$ is known as the \textit{modified Bessel equation of order zero}, with general solution $$ R(r)= C_1 I_0(s) + C_2 K_0(s), $$ where $C_1$, $C_2 \in \mathbb{R}$ and $I_n$ and $K_n$ are, respectively, the modified Bessel functions of first and second kind of order $n$. Again, since $K_n$ has a singularity at $s=0$, we have that $R(s)=C_1I_0(s)$. It is well known that $I_0'(s)=I_1(s)$ and the boundary condition at $s=a$ implies that that the eigenvalues $\eta$ satisfy that \[ C_1I_0'\big(\frac{\sqrt{\eta {\rm Th}_B}}{\sigma}\big)=0. \] Nevertheless, $I_0'(s)=I_1(s)>0$, so that $C_1$ must be zero and the corresponding solution $R(s)$ is the trivial one. \smallskip \noindent\textbf{Case 3: $\eta=0$}. Denoting $Q(s)=R'(s)$, the second order differential equation in $R$ can be rewritten as $sQ'(s)+Q(s)=0$. Easy calculations lead to $$ R(s)=-C_1{\rm e}^{-s}+C_2, $$ where $C_1$ and $C_2$ are constants to be determined with the boundary conditions. Thus, since $R'(0)=0$, it follows that $C_1=0$ and one concludes that $R(s)=C_2$. Consequently, one has that the countable set of admissible eigenvalues $\eta$ is \begin{equation}\label{eq:3_8} E=\{0\}\cup \big\{- \frac{(\sigma T_n)^{2}}{{\rm Th}_B} \big\}_{n\in \mathbb{N}}, \end{equation} where $T_n$ is such that $J_1(T_n)=0$, $J_1$ being the Bessel function of first kind and order one. The general solution for the second order differential equation for $R$ is $$ R(r)= C_0 + \sum_{n \in \mathbb{N}} C_n J_0 \big(\frac{\sqrt{-{\rm Th}_B\eta_n}}{\sigma}r\big). $$ \smallskip \noindent\textbf{Step 1.3.} Calculation of $Z(z)$. Using the boundary conditions of system \eqref{eq:3_5} on $\Gamma_{\rm in}$ and $\Gamma_{\rm out}$, it is clear that function $Z(z)$ is solution of system \begin{equation}\label{eq:3_9} \begin{gathered} \frac{1}{{\rm Th}_B} Z''(z) + \frac{1}{\rm Da} Z'(z) - (\lambda-\mu(1) - \eta)Z(z)=0, \quad z \in (0,1), \\ \frac{1}{{\rm Th}_B} Z'(1) + \frac{1}{\rm Da}Z(1) = 0, \\ Z'(0) = 0, \end{gathered} \end{equation} which corresponds to a regular Sturm-Liouville eigenvalue problem (see \cite[Theorem 1.3]{Knapp2005}). The corresponding characteristic equation is $$ \frac{1}{{\rm Th}_B} \rho^{2} + \frac{1}{{\rm Da}} \rho - (\lambda-\mu(1) - \eta)=0, $$ with roots $$ \rho= \frac{-{\rm Th}_B}{2 {\rm Da}} \pm \frac{ {\rm Th}_B}{2} \sqrt{ (\frac{1}{{\rm Da}})^2 + \frac{4(\lambda-\mu(1) - \eta)}{{\rm Th}_B}}. $$ Now, depending on the value of $\Delta= (\frac{1}{{\rm Da}})^2 + \frac{4(\lambda-\mu(1) - \eta)}{{\rm Th}_B}$, three possible solutions appear. \smallskip \noindent\textbf{Case 1: $\Delta=0 \Leftrightarrow \lambda = \eta + \mu(1) - {\rm Th}_B(\frac{1}{2{\rm Da}})^{2}$.} In this case, the solution of system \eqref{eq:3_9} is $$ Z(z)= D_1 {\rm e}^{\alpha z} + D_2z {\rm e}^{\alpha z}, $$ where $\alpha= \frac{-{\rm Th}_B}{2 {\rm Da}}$ and $D_1$, $D_2$ are constants which are determined by the boundary conditions of the system. Since $$ Z'(z)= \alpha {\rm e}^{\alpha z}(D_1 + z D_2) + D_2 {\rm e}^{\alpha z}, $$ then $Z'(0)= \alpha D_1 + D_2 =0$ if and only if $D_2 = -\alpha D_1$. Thus, the solution and its derivative can be rewritten as $$ Z(z)= D_1 {\rm e}^{\alpha z}\big(1- \alpha z \big) \quad\text{and} \quad Z'(z)= -D_1 \alpha^{2} z {\rm e}^{\alpha z}. $$ From the boundary condition at $z=1$ it follows that \begin{equation}\label{eq:3_10} D_1 {\rm e}^{\alpha} \big ( \frac{1}{{\rm Da}}(1-\alpha) - \frac{\alpha^{2}}{{\rm Th}_B}\big)=0. \end{equation} By replacing $\alpha$ by its value into equation \eqref{eq:3_10}, we conclude that this equation is true either if ${\rm Th}_B=-4{\rm Da}$ or if $D_1=0$. The first option is not possible since constants ${\rm Da}$ and ${\rm Th}_B$ are assumed strictly positive. Thus, the only solution in this case is $Z(z)=0$. \smallskip \noindent\textbf{Case 2: $\Delta<0 \Leftrightarrow \lambda < \eta + \mu(1) - {\rm Th}_B(\frac{1}{2{\rm Da}})^{2}$.} In this case, we have two complex conjugate roots $\rho= \alpha \pm i \beta$, where $\alpha\in (-\infty,0)$ and $\beta\in (0,+\infty)$. Then, the solution of system \eqref{eq:3_9} is of the form $$ Z(z)= {\rm e}^{\alpha z}\big(D_1\cos(\beta z) + D_2 \sin(\beta z)\big), $$ where $D_1$ and $D_2$ are constants which will be determined by the boundary conditions. Since $$ Z'(z)= \alpha Z(z) + \beta {\rm e}^{\alpha z}\big(-D_1\sin(\beta z) + D_2 \cos(\beta z)\big), $$ then $Z'(0)= \alpha D_1 + \beta D_2 =0$ if and only if $D_2 = -\frac{\alpha}{\beta} D_1$. Thus, the solution and its derivative can be rewritten as $$ Z(z)= D_1 {\rm e}^{\alpha z}\big(\cos(\beta z) - \frac{\alpha}{\beta} \sin(\beta z)\big) \quad\text{and} \quad Z'(z)= -D_1 {\rm e}^{\alpha z}\sin(\beta z) (\frac{\alpha^{2}}{\beta} + \beta). $$ From the boundary condition at $z=1$ it follows that $$ D_1 {\rm e}^{\alpha} \Big ( \frac{1}{{\rm Da}}\cos(\beta) - \sin(\beta) \big( \frac{1}{{\rm Da}} \frac{\alpha}{\beta} + \frac{1}{{\rm Th}_B} (\frac{\alpha^{2}}{\beta} + \beta)\big)\Big)=0, $$ which solutions are $D_1=0$ or \begin{equation}\label{eq:3_11} \tan(\beta) = \frac{ \frac{1}{{\rm Da}}}{ \frac{\alpha}{\beta} \frac{1}{{\rm Da}} + (\frac{\alpha^{2}}{\beta}+\beta) \frac{1}{{\rm Th}_B}} = \frac{\beta}{ \frac{{\rm Da}}{{\rm Th}_B} \beta^{2} + \frac{\alpha}{2}} = \frac{2\alpha \beta}{-\beta^{2} + \alpha^{2}}. \end{equation} As $F(\beta)= \frac{2\alpha \beta}{\alpha^{2}-\beta^{2}}$ is a decreasing function and has an asymptote at $\beta=-\alpha$, there exists a countable set $\{ \beta_n \}_{n\in \mathbb{N}}$ with $\beta_n \in ((n-1)\pi, n \pi)$ satisfying $F(\beta_n)=\tan(\beta_n)$. Consequently, $$ Z(z)= \sum_{n \in \mathbb{N}} D_n {\rm e}^{- \frac{{\rm Th}_B} {2{\rm Da}} z}\big(\cos(\beta_n z) + \frac{{\rm Th}_B}{2{\rm Da}\beta_n} \sin(\beta_n z)\big), $$ where $\beta_n \in (0,+\infty)$ satisfies equation \eqref{eq:3_11}. \smallskip \noindent\textbf{Case 3: $\Delta>0 \Leftrightarrow \lambda > \eta + \mu(1) - {\rm Th}_B(\frac{1}{2{\rm Da}})^{2}$.} In this case, we have two different real roots $\rho_{1,2}= \alpha \pm \beta$, with $\alpha= \frac{-{\rm Th}_B}{2 {\rm Da}}$, $\beta= \frac{ {\rm Th}_B}{2} \sqrt{ (\frac{1}{{\rm Da}})^2 + \frac{4(\lambda-\mu(1) - \eta)}{{\rm Th}_B}}$, and the solution of equation \eqref{eq:3_9} is of the form $$ Z(z)= D_1 {\rm e}^{(\alpha+\beta) z} + D_2 {\rm e}^{(\alpha-\beta)z}, $$ where $D_1$ and $D_2$ are constants which will be determined by the boundary conditions. Since $Z'(z)= (\alpha+\beta) D_1 {\rm e}^{(\alpha+\beta)z} + (\alpha-\beta)D_2{\rm e}^{(\alpha-\beta) z},$ $\alpha<0$ and $\beta>0$, then $Z'(0)= (\alpha+\beta) D_1 + (\alpha-\beta) D_2 =0$ if and only if $D_2 = -\frac{(\alpha+\beta)}{(\alpha-\beta)} D_1$. Thus, the solution and its derivative can be rewritten as $$ Z(z)= D_1 \big( {\rm e}^{(\alpha + \beta) z} - \frac{(\alpha+\beta)}{(\alpha-\beta)} {\rm e}^{(\alpha-\beta) z}\big), \quad Z'(z)= D_1 (\alpha+\beta)\big( {\rm e}^{(\alpha+\beta) z} -{\rm e}^{(\alpha-\beta) z}\big). $$ From the boundary condition at $z=1$, it follows that $$ D_1{\rm e}^{\alpha} \frac{(\alpha+\beta)}{{\rm Th}_B} \big({\rm e}^{\beta}-{\rm e}^{-\beta}\big) + D_1{\rm e}^{\alpha} \frac{1}{{\rm Da}}\big( {\rm e}^{\beta} - \frac{(\alpha +\beta)}{(\alpha-\beta)} {\rm e}^{-\beta}\big)=0, $$ which implies $D_1=0$ or \begin{equation}\label{eq:3_12} \begin{aligned} & {\rm e}^{\beta}\big(\frac{(\alpha+\beta)}{{\rm Th}_B} + \frac{1}{{\rm Da}}\big)= {\rm e}^{-\beta} \big( \frac{(\alpha+\beta)}{{\rm Th}_B} + \frac{(\alpha+\beta)}{(\alpha-\beta)}\frac{1}{{\rm Da}}\big) \\ &\Leftrightarrow \; {\rm e}^{2\beta}= \frac{ \frac{(\alpha+\beta)}{{\rm Th}_B} + \frac{1}{{\rm Da}} \frac{(\alpha+\beta)}{(\alpha-\beta)}}{ \frac{(\alpha+\beta)} {{\rm Th}_B} + \frac{1}{\rm Da}} = \frac{(\alpha+\beta)}{(\alpha-\beta)} \big(\frac{-(\beta+\alpha){\rm Da}}{(\beta-\alpha){\rm Da}}\big) \\ &\Leftrightarrow \; {\rm e}^{2\beta} = (\frac{\alpha+\beta}{\alpha-\beta})^{2}. \end{aligned} \end{equation} Again, as $\beta > 0$ and $\alpha<0$, then $(\beta+\alpha)^{2} < (\alpha-\beta)^{2}$ and thus $(\frac{\alpha+\beta}{\alpha-\beta})^{2} < 1$. This implies that $D_1=0$ is the unique admissible solution and $Z(z)=0$. \smallskip \noindent\textbf{Step 1.4.} General expression of $\delta B(r,z,t)$. Given $\eta_n\in E$ (see equation \eqref{eq:3_8}), there exists a countable set of admissible eigenvalues $\lambda$ \begin{equation}\label{eq:3_13} \Lambda_n= \{\lambda_{n m}\}_{m \in \mathbb{N}} = \{ \mu(1) +\eta_n - \frac{{\rm Th}_B}{(2{\rm Da})^{2}} -\frac{\beta_{m}^{2}}{{\rm Th}_B}\}_{m\in \mathbb{N}}, \end{equation} where $\beta_{m}$ satisfies system \eqref{eq:3_11}. Consequently, \begin{align*} \delta B(r,z,t) &= \sum_{ n\in \mathbb{N} \cup \{0\} }\sum_{m\in \mathbb{N}} A_{nm} {\rm e}^{ \lambda_{nm} t} J_0\big(\frac{\sqrt{-{\rm Th}_B\eta_n}} {\sigma}r\big)\\ &\quad\times {\rm e}^{- \frac{{\rm Th}_B}{2{\rm Da}} z} \big(\cos(\beta_{m} z) + \frac{{\rm Th}_B}{2{\rm Da}\beta_{m}} \sin(\beta_{m} z)\big), \end{align*} where $\eta_n\in E$, $\beta_{m}$ satisfies \eqref{eq:3_11}, $\lambda_{nm}\in \Lambda_n$ and the constants $A_{n m}$ are chosen such that $\delta B(r,z,0)=\delta B_{\rm init}(r,z)$. Notice that the constants $A_{nm}$ are well defined since the two systems \eqref{eq:3_6} and \eqref{eq:3_9} are regular Sturm-Liouville eigenvalue problems (see, e.g. \cite[Theorem 1.3]{Knapp2005}). Using Parseval's equation (see, for instance, \cite{Weinberger}) one has that $$ \|\delta B(t)\|_{L^{2}(\Omega)}^{2} = \sum_{n\in \mathbb{N}\cup\{0\}} \sum_{m \in \mathbb{N}} A_{nm}^{2} {\rm e}^{2\lambda_{n m}t}. $$ Furthermore, it is straightforward to see that $$ \lambda_{n m}\leq \lambda_{01} = \mu(1) - \frac{{\rm Th}_B}{(2{\rm Da})^{2}} - \frac{\beta_1^{2}}{{\rm Th}_B} \quad \forall\text{ }(n,m)\in (\{0\}\cup\mathbb{N})\times \mathbb{N}. $$ Therefore, if \begin{equation}\label{eq:3_14} \lambda_{01}= \mu(1) - \frac{{\rm Th}_B}{(2{\rm Da})^{2}} - \frac{\beta_1^{2}}{{\rm Th}_B}< 0, \end{equation} (which is the same condition as \eqref{eq:3_4}) it follows that \begin{equation*} \|\delta B(t)\|_{L^{2}(\Omega)}^{2} \leq {\rm e}^{2\lambda_{01} t} \sum_{n\in \mathbb{N}\cup\{0\}} \sum_{m \in \mathbb{N}} A_{nm}^{2} = {\rm e}^{2\lambda_{01} t} \|\delta B(0)\|_{L^{2}(\Omega)}^{2} \xrightarrow{t \to \infty} 0. \end{equation*} Note that, if $\lambda_{01}<0$, one can also deduce inequality (that will be used at the end of this proof) \begin{equation}\label{eq:3_15} \|\delta B(t)\|_{L^{2}(\Omega)}^{2} \leq \|\delta B(0)\|_{L^{2}(\Omega)}^{2}\leq K^{2} \|\delta B_{\rm init}\|_{L^{\infty}(\Omega)}^{2}, \end{equation} where $K$ is a constant relating the norms $\|\cdot\|_{L^{2}(\Omega)}$ and $\|\cdot\|_{L^{\infty}(\Omega)}$. \smallskip \noindent\textbf{Step 2.} Let us prove that $\|\delta S(t)\|_{L^{2}(\Omega)} \to 0 \text{ as }t \to \infty$. Regarding $\delta S$, the main equation involving the substrate in system \eqref{eq:3_5} is an Advection-Diffusion equation with non-homogeneous term $-\mu(1)\delta B(r,z,t)$, which makes complex the use of separation of variables. Here, we prove that $\|\delta S(\cdot,\cdot,t)\|_{L^{2}(\Omega)} \xrightarrow{t \to \infty} 0$ by using variational techniques. To this aim, we multiply the first equation in system \eqref{eq:3_5} by $r\delta S$ and integrate as follows \begin{align*} &\int_0^{t}\int_{\Omega} r\frac{{\rm d}\delta S}{{\rm d} \tau} \delta S{\rm d} r{\rm d} z{\rm d} \tau\\ &= \frac{1}{{\rm Da}}\int_0^{t}\int_{\Omega} r\frac{{\rm d} \delta S}{{\rm d} z} \delta S{\rm d}r {\rm d}z {\rm d}\tau + \frac{1}{{\rm Th}_S}\int_0^{t}\int_{\Omega} r\frac{{\rm d}^{2} \delta S}{{\rm d} z^{2}}\delta S{\rm d}r {\rm d}z{\rm d}\tau \\ &\quad + \frac{\sigma^{2}}{{\rm Th}_S}\int_0^{t} \int_{\Omega }\frac{{\rm d}}{{\rm d} r}(r\frac{{\rm d} \delta S}{{\rm d}r}) \delta S {\rm d} r{\rm d}z{\rm d} \tau -\mu(1)\int_0^{t}\int_{\Omega}r\delta B\delta S{\rm d} r{\rm d}z {\rm d}\tau \\ &= - \frac{1}{\rm Da}\int_0^{t}\int_{\Omega} r\frac{{\rm d}\delta S}{{\rm d}z} \delta S {\rm d} r {\rm d}z {\rm d} \tau - \frac{\sigma^{2}}{{\rm Th}_S}\int_0^{t}\int_{\Omega} r(\frac{{\rm d} \delta S}{{\rm d}r})^{2} {\rm d} r {\rm d}z{\rm d} \tau \\ &\quad -\frac{1}{{\rm Th}_S} \int_0^{t} \int_{\Omega} r \big( \frac{{\rm d} \delta S}{{\rm d} z})^{2} {\rm d}r {\rm d}z{\rm d}\tau -\mu(1)\int_0^{t}\int_{\Omega}r\delta B\delta S{\rm d} r {\rm d} z {\rm d}\tau \\ &\quad + \frac{\sigma^{2}}{{\rm Th}_S}\int_0^{t}\int_{\Gamma_{\rm sym} \cup \Gamma_{\rm wall}}r\frac{{\rm d} \delta S}{{\rm d}r} \delta S {\rm d}z{\rm d} \tau +\int_0^{t}\int_{\Gamma_{\rm in}} r \big( \frac{1}{{\rm Th}_S} \frac{{\rm d} \delta S}{{\rm d} z} + \frac{1}{{\rm Da}} \delta S\big) \delta S {\rm d}r {\rm d}\tau \\ &\quad - \int_0^{t}\int_{\Gamma_{\rm out}} r \big( \frac{1}{{\rm Th}_S} \frac{{\rm d} \delta S}{{\rm d} z} + \frac{1}{{\rm Da}} \delta S\big) \delta S {\rm d}r {\rm d}\tau. \end{align*} However, applying the boundary conditions for $\delta S$ in system \eqref{eq:3_5} one has \begin{equation}\label{eq:3_16} \begin{aligned} &\int_0^{t}\int_{\Omega} r\frac{{\rm d}\delta S}{{\rm d} \tau}\delta S{\rm d} r{\rm d} z{\rm d} \tau \\ &= \underbrace{-\frac{1}{\rm Da}\int_0^{t}\int_{\Omega} r \frac{{\rm d} \delta S}{{\rm d}z} \delta S {\rm d} r {\rm d}z {\rm d} \tau}_{(I)} - \frac{1}{\rm Da} \int_0^{t}\int_{\Gamma_{\rm out}} r \delta S^{2}{\rm d} r {\rm d} z {\rm d} \tau \\ &\quad - \frac{\sigma^{2}}{{\rm Th}_S}\int_0^{t}\int_{\Omega} r(\frac{{\rm d} \delta S}{{\rm d}r})^{2} {\rm d} r {\rm d}z{\rm d} \tau - \frac{1}{{\rm Th}_S}\int_0^{t}\int_{\Omega} r (\frac{{\rm d} \delta S}{{\rm d}z})^{2} {\rm d} r {\rm d}z{\rm d} \tau \\ &\quad - \mu(1)\int_0^{t}\int_{\Omega} r\delta S\delta B{\rm d}r {\rm d} z {\rm d} \tau. \end{aligned} \end{equation} The integral denoted by $(I)$ in equation \eqref{eq:3_16} can be rewritten as $$ (I) = -\frac{1}{2\rm Da}\int_0^{t}\int_{\Gamma_{\rm in}} r \delta S^{2}{\rm d} r {\rm d} \tau + \frac{1}{2\rm Da} \int_0^{t}\int_{\Gamma_{\rm out}} r \delta S^{2} {\rm d} r {\rm d} \tau. $$ Thus, equation \eqref{eq:3_16} leads to \begin{equation}\label{eq:3_17} \begin{aligned} &\frac{1}{2}\int_0^{t}\int_{\Omega} r\frac{{\rm d} (\delta S^{2})}{{\rm d}\tau} {\rm d} r{\rm d} z{\rm d} \tau + \frac{\sigma^{2}}{{\rm Th}_S}\int_0^{t}\int_{\Omega} r(\frac{{\rm d} \delta S}{{\rm d}r})^{2} {\rm d} r {\rm d}z{\rm d} \tau \\ & + \frac{1}{{\rm Th}_S}\int_0^{t}\int_{\Omega} r (\frac{{\rm d} \delta S}{{\rm d}z})^{2} {\rm d} r {\rm d}z{\rm d} \tau +\frac{1}{2\rm Da}\int_0^{t}\int_{\Gamma_{\rm in}} r \delta S^{2} {\rm d} r {\rm d} \tau \\ & \frac{1}{2\rm Da}\int_0^{t}\int_{\Gamma_{\rm out}} r \delta S^{2} {\rm d} r {\rm d} \tau \\ &= - \mu(1)\int_0^{t}\int_{\Omega} r\delta S\delta B{\rm d}r {\rm d} z {\rm d} \tau. \end{aligned} \end{equation} By multiplying \eqref{eq:3_17} by $2\pi$ and applying Young's inequality (with $\epsilon>0$ to be chosen afterward), we obtain \begin{align*} & \frac{1}{2}\int_0^{t}\frac{{\rm d}}{{\rm d} \tau} \big(\|\delta S(\tau)\|_{L^{2}(\Omega^{\ast})}^{2}\big) {\rm d}\tau + \frac{\min(1,\sigma^{2})}{{\rm Th}_S} \int_0^{t}\|\nabla \delta S(\tau)\|_{L^{2}(\Omega^{\ast})}^{2}{\rm d} \tau \\ & + \frac{1}{2{\rm Da}}\int_0^{t} \|\delta S(\tau)\|^{2}_{L^{2} (\Gamma_{\rm out}^{\ast})}{\rm d}\tau \\ &\leq \epsilon\mu(1)\int_0^{t} \|\delta S(\tau)\|_{L^{2}(\Omega^{\ast})}^{2}{\rm d} \tau + \frac{\mu(1)}{4\epsilon} \int_0^{t}\|\delta B(\tau) \|_{L^{2}(\Omega^{\ast})}^{2}{\rm d}\tau. \end{align*} Considering $A= \min\{ \frac{1}{{\rm Th}_S }, \frac{\sigma^{2}}{{\rm Th}_S}, \frac{1}{2{\rm Da}}\}$, it follows that \begin{equation}\label{eq:3_18} \begin{aligned} & \frac{1}{2}\int_0^{t}\frac{{\rm d}}{{\rm d}\tau} \big(\|\delta S(\tau)\|_{L^{2}(\Omega^{\ast})}^{2}\big) {\rm d}\tau + A\int_0^{t}(\|\nabla \delta S(\tau)\|_{L^{2}(\Omega^{\ast})}^{2} +\|\delta S(\tau)\|^{2}_{L^{2}(\Gamma_{\rm out}^{\ast})}){\rm d} \tau\\ & \leq \epsilon\mu(1)\int_0^{t} \|\delta S(\tau)\|_{L^{2}(\Omega^{\ast})}^{2}{\rm d} \tau + \frac{\mu(1)}{4\epsilon} \|\delta B\|_{L^{2}(\Omega^{\ast} \times (0,t))}^{2}. \end{aligned} \end{equation} Now, applying Friedrich's inequality (see e.g., \cite[Theorem 6.1]{Lanzani2004}) to inequality \eqref{eq:3_18} with $E=\Gamma_{\rm out}^{\ast}$, there exits a constant $C$ depending on $\Omega^{\ast}$ and $\Gamma_{\rm out}^{\ast}$ such that \begin{equation}\label{eq:3_19} \begin{aligned} \frac{1}{2}\int_0^{t}\frac{{\rm d}}{{\rm d}\tau} \big(\|\delta S(\tau)\|_{L^{2}(\Omega^{\ast})}^{2}\big){\rm d}\tau &\leq \big(\epsilon\mu(1)-\frac{A}{C}\big)\int_0^{t} \|\delta S(\tau)\|_{L^{2}(\Omega^{\ast})}^{2}{\rm d} \tau \\ &\quad + \frac{\mu(1)}{4\epsilon} \|\delta B\|_{L^{2}(\Omega^{\ast} \times (0,t))}^{2}. \end{aligned} \end{equation} Next, applying the Gronwall's inequality in its integral form, it follows that $$ \|\delta S(t)\|_{L^{2}(\Omega^{\ast})}^{2} \leq \underbrace{\big(\|\delta S_{\rm init}\|_{L^{2}(\Omega^{\ast})}^{2} +\frac{\mu(1)}{2\epsilon}\|\delta B\|_{L^{2}(\Omega^{\ast} \times(0,t))}^{2}\big)}_{(:=m(t))} \exp\big(\underbrace{2(\epsilon\mu(1)- \frac{A}{C})}_{(:=\alpha)}t\big). $$ Since $\|\delta B(t)\|_{L^{2}(\Omega)}^{2} \leq K^{2} \|\delta B_{\rm init}\|_{L^{\infty}(\Omega)}^{2}$ for all $t>0$ (see equation \eqref{eq:3_15}), taking $\epsilon< \frac{A}{\mu(1)C}$ it follows that $\alpha<0$. Thus, $$ \|\delta S(t)\|_{L^{2}(\Omega^{\ast})}^{2} \leq \big(\|\delta S_{\rm init}\|_{L^{2}(\Omega^{\ast})}^{2} +\frac{\mu(1)}{2\epsilon}t K^{2} \|\delta B_{\rm init}\|_{L^{\infty}(\Omega^{\ast})}^{2} \big) {\rm e}^{\alpha t} \xrightarrow{ t \to \infty} 0. $$ \end{proof} \begin{remark} \label{rmk3.7} \rm To the best of our knowledge, if the inequality ``$<$'' is replaced by the equality ``$=$'' in condition \eqref{eq:3_4}, the stability analysis of the steady state $(1,0)$ requires different techniques from those used in the proof of Theorem \ref{thm:3_5}. This case has not been tackled here. \end{remark} Taking into account Theorem \ref{thm:3_5}, we conjecture (supported by the numerical experiments presented in Section \ref{sec:4}) that the following result holds: \begin{conjecture}\label{conj:3_8} \rm If $\mu$ satisfies (A1) (respectively, (A2)), a sufficient condition for $(S_2^{\ast},B_2^{\ast})$ to be a linearly asymptotically stable steady state of system \eqref{eq:2_4} is that \begin{equation}\label{eq:3_20} \mu(1) > \frac{\rm Th_{B}}{(2{\rm Da})^{2}} + \frac{(\beta_1({\rm Da},{\rm Th_{B}}))^{2}}{\rm Th_{B}} \end{equation} (respectively, \begin{equation}\label{eq:3_21} 1 > \frac{\rm Th_{B}}{(2{\rm Da})^{2}} + \frac{(\beta_1({\rm Da},{\rm Th_{B}}))^{2}}{\rm Th_{B}}). \end{equation} \end{conjecture} \begin{remark}\rm In terms of the variables with dimensions appearing in system \eqref{eq:2_1}, conditions \eqref{eq:3_20} and \eqref{eq:3_21} are reformulated, respectively, as \begin{equation*} \mu(S_{\rm e}) > \frac{u^{2}}{4D_B} + \frac{D_B}{H^{2}}(\tilde{\beta}_1(H,u,D_B))^{2}, \end{equation*} and \begin{equation*} \bar{\mu} > \frac{u^{2}}{4D_B} + \frac{D_B}{H^{2}}(\tilde{\beta}_1(H,u,D_B))^{2}. \end{equation*} \end{remark} \begin{remark}\label{rmk:3_10}\rm From Theorem \ref{thm:3_5} and Conjecture \ref{conj:3_8}, it follows that if $\mu$ satisfies (A2) and $\mu(1)<1$, there is bistability in system \eqref{eq:2_4} when $$ \mu(1) <\frac{\rm Th_{B}}{(2{\rm Da})^{2}} + \frac{(\beta_1({\rm Da},{\rm Th_{B}}))^{2}}{\rm Th_{B}} <1. $$ \end{remark} \subsubsection{Bounds for the flow rate assuring asymptotic stability of the steady states}\label{sec:3_2_1} Conditions \eqref{eq:3_4} and \eqref{eq:3_21} include in their analytical expression the model parameters ${\rm Da}$, ${\rm Th}_B$ and $\mu(1)$, among which the flow rate ${\rm Da}$ can be seen as a bioreactor control parameter. In this section, we present bounds for the parameter ${\rm Da}$ assuring the asymptotic stability of the steady states $(1,0)$ and $(S_2^{\ast},B_2^{\ast})$. To do so, we first define the following function. \begin{definition}\label{defn:3_11}\rm For a fixed value ${\rm Th}_B$, we define the function $f_{\rm Th_{B}}: [0,+\infty) \to [0,+\infty)$ \[ f_{\rm Th_{B}}({\rm Da})= \frac{\rm Th_{B}}{(2{\rm Da})^{2}} + \frac{(\beta_1({\rm Da},{\rm Th_{B}}))^{2}}{\rm Th_{B}}. \] \end{definition} In Figure \ref{fig:2} we plot the value of functions $\beta_1({\rm Da},{\rm Th_{B}})$ and $f_{\rm Th_{B}}({\rm Da})$ for ${\rm Th_{B}}\in \{ \frac{1}{5},1,5\}$ and ${\rm Da}\in [0,2]$. For a fixed value ${\rm Th}_B$, function $\beta_1(\cdot,{\rm Th_{B}})$ is decreasing, bounded by $\pi$ (see the proof of Theorem \ref{thm:3_5} for a detailed explanation of this feature) and $\beta_1({\rm Da},{\rm Th_{B}}) \xrightarrow{{\rm Da}\to +\infty} 0$. One can also conclude that, for a fixed value ${\rm Th}_B$, function $f_{\rm Th_{B}}$ is decreasing, $f_{\rm Th_{B}}({\rm Da})\xrightarrow{{\rm Da}\to 0} +\infty$ and $f_{\rm Th_{B}}({\rm Da})\xrightarrow{{\rm Da}\to +\infty} 0$. Taking into account these properties of $f_{\rm Th_{B}}$, we define the following variables. \begin{figure}[ht] \begin{center} \subfigure[$\beta_1({\rm Da},{\rm Th}_B)$] {\includegraphics[width=0.48\textwidth]{fig2a}} \subfigure[$f_{\rm Th_{B}}({\rm Da})$] {\includegraphics[width=0.48\textwidth]{fig2b}} \end{center} \caption{Graphical plots of functions $\beta_1({\rm Da},{\rm Th_{B}})$ and $f_{\rm Th_{B}}({\rm Da})$ (described in Definitions \ref{defn:3_4} and \ref{defn:3_11}, respectively) for ${\rm Th_{B}}\in \{\frac{1}{5},1,5\}$ and ${\rm Da}\in [0,2]$.} \label{fig:2} \end{figure} \begin{definition}\label{defn:3_12} \rm We define \begin{itemize} \item ${\rm Da}_{\eqref{eq:2_4}}^{{\rm W}}({\rm Th}_B,\mu(1)) := (f_{\rm Th_{B}})^{-1}(\mu(1))$. \item ${\rm Da}_{(A1),\eqref{eq:2_4}}^{\rm NW}({\rm Th}_B,\mu(1)) := {\rm Da}_{\eqref{eq:2_4}}^{{\rm W}}({\rm Th}_B,\mu(1))$. \item ${\rm Da}^{\rm NW}_{(A2),\eqref{eq:2_4}}({\rm Th}_B) := (f_{\rm Th_{B}})^{-1}(1)$. \end{itemize} \end{definition} \begin{remark}\label{rmk:3_13}\rm Following Theorem \ref{thm:3_5}, if Conjecture \ref{conj:3_8} is true, it follows that \begin{itemize} \item If ${\rm Da} < {\rm Da}_{\eqref{eq:2_4}}^{{\rm W}}({\rm Th}_B,\mu(1))$, then the equilibrium state $(1,0)$ of system \eqref{eq:2_4} is linearly asymptotically stable. \item If $\mu$ satisfies (A1) and ${\rm Da} > {\rm Da}_{(A1), \eqref{eq:2_4}}^{{\rm NW}}({\rm Th}_B,\mu(1))$, then the equilibrium state $(S_2^{\ast},B_2^{\ast})$ of system \eqref{eq:2_4} is linearly asymptotically stable. \item If $\mu$ satisfies (A2) and ${\rm Da} > {\rm Da}_{(A2), \eqref{eq:2_4}}^{{\rm NW}}({\rm Th}_B)$, then the equilibrium state $(S_2^{\ast},B_2^{\ast})$ of system \eqref{eq:2_4} is linearly asymptotically stable. \end{itemize} \end{remark} \section{Numerical experiments}\label{sec:4} In this section, we describe the results of the numerical experiments performed to analyze the validity and robustness of the stability analysis done in Section \ref{sec:3}. In Section \ref{sec:4_1}, we study the sensitivity of variables ${\rm Da}^{\rm W}_{\eqref{eq:2_4}}({\rm Th_{B}},\mu(1))$ and ${\rm Da}^{\rm NW}_{(A2),\eqref{eq:2_4}}({\rm Th_{B}})$, defined in Section \ref{sec:3_2_1}, regarding the model parameters. Then, In Section \ref{sec:4_2}, we carry out the numerical implementation of system \eqref{eq:2_4}-\eqref{eq:2_5} in order to check the interest of these functions. Finally, in Section \ref{sec:4_3}, we compare the results of the stability analysis of systems \eqref{eq:1_2} and \eqref{eq:2_4}. In this section, the value of functions ${\rm Da}_{\eqref{eq:2_4}}^{{\rm W}}({\rm Th}_B,\mu(1))$ and ${\rm Da}^{\rm NW}_{(A2),\eqref{eq:2_4}}({\rm Th}_B)$ is approximated numerically using a self-implemented \textit{Dichotomy method} (see, e.g. \cite{Infante2015}). Moreover, for each pair $({\rm Th}_B,{\rm Da})$, the value of $\beta_1({\rm Th}_B,{\rm Da})$ (see Definition \ref{defn:3_4}) was computed by using the MATLAB function (see \\ www.mathworks.com/help/symbolic/vpasolve.html). \subsection{Sensitivity to model parameters}\label{sec:4_1} In this section, we perform the sensitivity analysis of ${\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,\mu(1))$ with respect to the nondimensional parameters ${\rm Th}_B$ and $\mu(1)$ (the sensitivity analysis of ${\rm Da}_{(A2),\eqref{eq:2_4}}^{\rm NW}({\rm Th}_B)$ can be obtained with a similar methodology). \subsubsection{Sensitivity with respect to {$\mu(1)$}{Lg}} Taking into account that \[ {\rm Da}^{\rm W}_{\eqref{eq:2_4}}({\rm Th}_B,\mu(1)) = (f_{\rm Th_{B}})^{-1}(\mu(1)) \] and $f_{\rm Th_{B}}$ is decreasing, one concludes that, for any fixed value ${\rm Th}_B$, the function ${\rm Da}^{\rm W}_{\eqref{eq:2_4}}({\rm Th}_B,\mu(1))$ decreases as $\mu(1)$ increases. This is physically reasonable since, as parameter $\mu(1)$ increases, the range of flow rates $(1/{\rm Da})$ suitable to avoid washout also increases (see, e.g, \cite{Crespo2016,Drame2008,Vetjasa1970}). \subsubsection{Sensitivity with respect to {${\rm Th}_B$}{Lg}} To easily analyze the of ${\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,\mu(1))$ with respect to ${\rm Th}_B$, we aim to approximate ${\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,\mu(1))$ by using the following variables: $\bullet$ $\overline{\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,\mu(1)) :=\frac{1}{2}\sqrt{\frac{{\rm Th}_B}{\mu(1)}}$. This should be a good approximation of ${\rm Da}^{\rm W}_{\eqref{eq:2_4}}({\rm Th}_B,\mu(1))$ assuming that the second term of the right hand side of condition \eqref{eq:3_4} is negligible. $\bullet$ $\widehat{\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,\mu(1)) := (g_{{\rm Th}_B})^{-1}(\mu(1))$, where $g_{\rm Th_{B}}: [0,+\infty) \to [0,+\infty)$, \[ g_{\rm Th_{B}}({\rm Da})= \to \frac{(\beta_1({\rm Th}_B, {\rm Da}))^{2}}{{\rm Th}_B}. \] This should be a good approximation of ${\rm Da}^{\rm W}_{\eqref{eq:2_4}}({\rm Th}_B,\mu(1))$ assuming that the first term of the right hand side of condition \eqref{eq:3_4} is negligible. Since $\beta_1({\rm Th}_B,{\rm Da})<\pi$ (see the proof of Theorem \ref{thm:3_5} for a detailed explanation of this fact), if ${\rm Th}_B\mu(1) > \pi^{2}$, then the function $\widehat{\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,\mu(1))$ is not defined. We approximate numerically $\widehat{\rm Da}_{\eqref{eq:2_4}}^{{\rm W}}({\rm Th}_B,\mu(1))$ applying the same methodology that the one used to approximate numerically ${\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,\mu(1))$, described above. Figure \ref{fig:3} illustrates the difference between ${\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,\mu(1))$, $\overline{\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,\mu(1))$ and $\widehat{\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,\mu(1))$ when $\mu(1)=0.5$ and ${\rm Th}_B\in [5\cdot 10^{-3}, 5\cdot 10^{3}]$. We observe that $\widehat{\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,0.5)$ approximates ${\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,0.5)$ for values smaller than $\log({\rm Th}_B) = -2$ (${\rm Th}_B \approx 0.1$) while $\overline{\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,0.5)$ approximates ${\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,0.5)$ for values larger than $\log({\rm Th}_B)=6$ (${\rm Th}_B \approx 400$). \begin{figure}[ht] \begin{center} \includegraphics[width=0.7\textwidth]{fig3} \end{center} \caption{Comparison between the functions ${\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,0.5)$, $\overline{\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,0.5)$ and $\widehat{\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,0.5)$ (depicted with solid, dotted and dashed lines, respectively), when ${\rm Th}_B\in [5\cdot 10^{-3},5\cdot 10^{3}]$.} \label{fig:3} \end{figure} The comparison between the functions ${\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,\mu(1))$, $\overline{\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,\mu(1))$ and $\widehat{\rm Da}_{\eqref{eq:1_2}}^{\rm W}({\rm Th}_B,\mu(1))$, shown in Figure \ref{fig:3} for $\mu(1)=0.5$, has been reproduced for reaction values $\mu(1)\in \{i/20\}_{i=1}^{20}$ and the results seems to indicate that in general: if ${\rm Th}_B \geq 10^{4}$, the function $\overline{\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,\mu(1))$ can be used as an approximation of ${\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,\mu(1))$; and if ${\rm Th}_B \leq 0.1$, the function $\widehat{\rm Da}_{\eqref{eq:2_4}}^{\rm W}(\mu(1))$ can be used as an approximation of ${\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,\mu(1))$. Taking into account the approximations of ${\rm Da}^{\rm W}_{\eqref{eq:2_4}}({\rm Th}_B,\mu(1))$ presented above and Figure \ref{fig:3}, the sensitivity of ${\rm Da}^{\rm W}_{\eqref{eq:2_4}}({\rm Th}_B,\mu(1))$ with respect to ${\rm Th}_B$ reads as follows: $\bullet$ If ${\rm Th}_B \leq 0.1$, the variable ${\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,\mu(1))$ is not sensible to parameter ${\rm Th}_B$. Indeed, small values of ${\rm Th}_B$ correspond, for instance, to high diffusion coefficients implying almost spatial homogeneous biomass concentration. In this case, there would be no differences when considering even higher diffusion coefficients. As we will see in Section \ref{sec:4_3}, if ${\rm Th}_B \leq 0.1$, the dynamics of the bioreactor can be modeled with ordinary differential equations. $\bullet$ If ${\rm Th}_B>0.1$, the variable ${\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,\mu(1))$ seems to increase with parameter ${\rm Th}_B$. This outcome is physically reasonable, since as parameter ${\rm Th}_B$ increases (equivalently, the diffusion coefficient decreases) the flow rate $(1/{\rm Da})$ should be chosen smaller to favor the reaction between the substrate and the biomass (see \cite{Crespo2016}). $\bullet$ If ${\rm Th}_B \geq 10^{4}$, the variable ${\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,\mu(1))$ is quadratically proportional to ${\rm Th}_B$. \subsection{Numerical validation of the results}\label{sec:4_2} In this section, we check the properties given in Remark \ref{rmk:3_13} for the threshold values ${\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,\mu(1))$ and \\ ${\rm Da}_{(A2),\eqref{eq:2_4}}^{\rm NW}({\rm Th}_B,\mu(1))$ by using the numerical solution of system \eqref{eq:2_4}-\eqref{eq:2_5}. To do that computation, we use the software COMSOL Multiphysics 5.0 (see the web page www.comsol.com), based on the Finite Element Method (see \cite{vazquez2014,angel}). The numerical experiments were carried out in a 2.8Ghz Intel i7-930 64bits computer with 12Gb of RAM. We used a triangular mesh with around 1000 elements and final nondimensional time $T= 300$. To validate the properties of the threshold values ${\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,\mu(1))$ and \\ ${\rm Da}_{(A2),\eqref{eq:2_4}}^{\rm NW}({\rm Th}_B,\mu(1))$, we define the following variables: $\bullet$ $\widetilde{\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,\mu(1)) := \sup\{ {\rm Da}$: the numerical solution of system \eqref{eq:2_4}-\eqref{eq:2_5} (with parameters ${\rm Th}_B$, ${\rm Da}$, ${\rm Th}_S={\rm Th}_B$, $\sigma=1$, $\mu$ the nondimensional Monod function with $K_S=\frac{1-\mu(1)}{\mu(1)}$, $S_{\rm init}=0.1$ and $B_{\rm init}=0.9$) approaches asymptotically the steady state $(1,0)\}$. $\bullet$ $\widetilde{\rm Da}_{(A2),\eqref{eq:2_4}}^{\rm NW}({\rm Th}_B) :=\inf\{ {\rm Da}$: the numerical solution of system \eqref{eq:2_4}-\eqref{eq:2_5} (with parameters ${\rm Th}_B$, ${\rm Da}$, ${\rm Th}_S={\rm Th}_B$, $\sigma=1$, $\mu$ the nondimensional Haldane function with $\mu^{\ast}/\bar{\mu}=1.7071$, $K_S= 0.3536$ and $K_{\rm I}= 2.8284$) approaches asymptotically a steady state different from $(1,0)\}$. We approximate numerically $\widetilde{\rm Da}_{\eqref{eq:2_4}}^{{\rm W}}({\rm Th}_B,\mu(1))$ and $\widetilde{\rm Da}^{\rm NW}_{(A2),\eqref{eq:2_4}}({\rm Th}_B)$ by using again a self-implemented \textit{Dichotomy method}. Figure \ref{fig:4}(a) illustrates the difference between ${\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,\mu(1))$ and $\widetilde{\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,\mu(1))$ when ${\rm Th}_B\in [5\cdot 10^{-3}, 1.5\cdot 10^{2}]$ and $\mu(1)=0.5$. Similarly, Figure \ref{fig:4}-(b) shows the difference between ${\rm Da}_{(A2),\eqref{eq:2_4}}^{\rm NW}({\rm Th}_B)$ and \\ $\widetilde{\rm Da}_{(A2),\eqref{eq:2_4}}^{\rm NW}({\rm Th}_B)$ when ${\rm Th}_B\in [5\cdot 10^{-3}, 1.5\cdot10^{2}]$. We point out that these comparisons were also performed with $\widetilde{\rm Da}^{\rm W}_{\eqref{eq:2_4}}({\rm Th_{B}},\mu(1))$ and $\widetilde{\rm Da}^{\rm NW}_{(A2),\eqref{eq:2_4}}({\rm Th_{B}})$ defined using other model parameters $\sigma$, ${\rm Th}_S$ and $\mu$ and similar results were obtained. \begin{figure}[ht] \begin{center} \subfigure[Comparing $\log({\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,0.5))$ (solid line) and $\log(\widetilde{\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,0.5))$ (dashed lines).] {\includegraphics[width=0.47\textwidth]{fig4a} }\; \subfigure[Comparing $\log({\rm Da}_{(A2), \eqref{eq:2_4}}^{\rm NW} ({\rm Th}_B)$ (solid line) and $\log(\widetilde{\rm Da}_{(A2),\eqref{eq:2_4}}^{\rm NW}({\rm Th}_B))$ (dashed line).] { \includegraphics[width=0.47\textwidth]{fig4b} } \end{center} \caption{Numerical validation of results.} \label{fig:4} \end{figure} In Figure \ref{fig:5}, we plot the steady-state solution $(S_2^{\ast},B_2^{\ast})$ of system \eqref{eq:2_4}, computed numerically when ${\rm Th}_B={\rm Th}_S={\rm e}^{4}$, ${\rm Da}={\rm e}^{2}$, $\sigma=1$, $S_{\rm init}=0.1$, $B_{\rm init}=0.9$ and $\mu$ being the nondimensional Monod function with $K_S=1$ (so that $\mu(1)=0.5$). With these parameters, e.g. when $\log({\rm Th}_B)=4$ and $\log({\rm Da})=2$, the equilibrium solution $(S_2^{\ast},B_2^{\ast})$ is linearly asymptotically stable (see Figure \ref{fig:4}-(a)). Notice that the same steady-state solution can be obtained with nonhomogeneous initial conditions (for instance, $S_{\rm init}(r,z)=r z$ and $B_{\rm init}(r,z)=r (1-z)$). \begin{figure}[ht] \begin{center} \subfigure[$S_2^{\ast}(r,z)$] {\includegraphics[width=0.47\textwidth]{fig5a}} \subfigure[$B_2^{\ast}(r,z)$] {\includegraphics[width=0.47\textwidth]{fig5b}} \end{center} \caption{Representation of the steady-state solution $(S_2^{\ast},B_2^{\ast})$ of \eqref{eq:2_4} computed numerically when ${\rm Th}_B={\rm Th}_S={\rm e}^{4}$, ${\rm Da}={\rm e}^{2}$, $\sigma=1$, $S_{\rm init}=0.1$, $B_{\rm init}=0.9$ and $\mu$ being the nondimensional Monod Function with $K_S=1$ (so that $\mu(1)=0.5$).} \label{fig:5} \end{figure} The bistability of system \eqref{eq:2_4}, stated in Remark \ref{rmk:3_10}, is perceivable when numerically solving system \eqref{eq:2_4}. For instance, if ${\rm Th_{B}}=0.01$, ${\rm Da}=1.5$ and $\mu(1)=0.5$, we observe that the solution of system \eqref{eq:2_4} (computed with parameters $\sigma=1$, ${\rm Th}_S=0.01$ and $\mu$ the nondimensional Haldane function with $\mu^{\ast}/\bar{\mu}=1.7071$, $K_S=0.0529$ and $K_{\rm I}=0.4235$) approaches $(1,0)$ if we choose $S_{\rm init}=0.9$ and $B_{\rm init}=0.1$, while it approaches a different equilibrium (similar to the one represented in Figure \ref{fig:5}) solution if we set $S_{\rm init}=0.1$ and $B_{\rm init}=0.9$. \subsection{Comparison with the stability analysis of system {\eqref{eq:1_2}}{Lg}}\label{sec:4_3} In this section, we compare the stability analysis conditions associated to the ODE and PDE systems \eqref{eq:1_2} and \eqref{eq:2_4}, respectively. As done in Section \ref{sec:3_2_1} for system \eqref{eq:2_4}, we define the following variables: \begin{definition}\label{defn:4_1}\rm \text{ } \begin{itemize} \item ${\rm Da}_{\eqref{eq:1_2}}^{\rm W}(\mu(1)):= 1/\mu(1)$. \item ${\rm Da}_{(A1),\eqref{eq:1_2}}^{\rm NW}(\mu(1)):= {\rm Da}_{\eqref{eq:1_2}}^{\rm W}(\mu(1))$. \item ${\rm Da}_{(A2),\eqref{eq:1_2}}^{\rm NW}:= 1$. \end{itemize} \end{definition} \begin{remark}\rm According to Remark \ref{rmk:2_2} and Definition \ref{defn:4_1}, the stability analysis of system \eqref{eq:1_2} (shown in Section \ref{sec:1}) can be rewritten as \begin{itemize} \item If ${\rm Da}<{\rm Da}^{\rm W}_{\eqref{eq:1_2}}(\mu(1))$, then the equilibrium solution $(1,0)$ of system \eqref{eq:1_2} is asymptotically stable. \item If $\mu$ satisfies (A1) and ${\rm Da}>{\rm Da}^{\rm NW}_{(A1),\eqref{eq:1_2}} (\mu(1))$, then the equilibrium solution $(S_2^{\ast},B_2^{\ast})$ of system \eqref{eq:1_2} is asymptotically stable. \item If $\mu$ satisfies (A2) and ${\rm Da}>{\rm Da}^{\rm NW}_{(A2), \eqref{eq:1_2}}(1)$, then the equilibrium solution $(S_2^{\ast},B_2^{\ast})$ of system \eqref{eq:1_2} is asymptotically stable. \end{itemize} \end{remark} Figure \ref{fig:6} illustrates the difference between the variable ${\rm Da}_{(A2),\eqref{eq:2_4}}^{\rm NW}({\rm Th}_B)$ and the constant ${\rm Da}_{(A2),\eqref{eq:1_2}}^{\rm NW}=1$ (and the difference, when $\mu(1)=0.5$, between the variable ${\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,\mu(1))$ and the constant ${\rm Da}_{\eqref{eq:1_2}}^{\rm W}(\mu(1))=2$). In both cases ${\rm Th}_B\in [5\cdot 10^{-3}, 1.5\cdot 10^{2}]$. Notice that the area limited between the curves ${\rm Da}_{(A2),\eqref{eq:2_4}}^{\rm NW}({\rm Th}_B)$ and ${\rm Da}_{\eqref{eq:2_4}}^{\rm W}({\rm Th}_B,0.5)$ is the region of bistability of system \eqref{eq:2_4} (see Remark \ref{rmk:3_10}). \begin{figure}[ht] \begin{center} \includegraphics[width=0.7\textwidth]{fig6} \end{center} \caption{Comparisong between $\log({\rm Da}_{\eqref{eq:2_4}}^{\rm W} ({\rm Th}_B,0.5))$, $\log({\rm Da}_{(A2), \eqref{eq:2_4}}^{\rm NW}({\rm Th}_B))$ (solid lines) and constant values $\log({\rm Da}_{\eqref{eq:1_2}}^{\rm W}(0.5))=\log(2)$, $\log({\rm Da}_{(A2),\eqref{eq:1_2}}^{\rm NW})=0$ (dashed line) when ${\rm Th}_B\in [5\times 10^{-3}, 1.5\times 10^{2}]$.} \label{fig:6} \end{figure} We observe that $\log({\rm Da}^{\rm NW}_{(A2),\eqref{eq:2_4}}({\rm Th}_B)) \approx 0$ for values smaller than $\log({\rm Th}_B) \approx -2$ (${\rm Th}_B \approx 0.1$). Similarly, for the particular case when $\mu(1)=0.5$, we observe that $\log({\rm Da}^{\rm W}_{\eqref{eq:2_4}}({\rm Th}_B,0.5))\approx \log(2)$ also for values smaller than $\log({\rm Th}_B) \approx -2$ (${\rm Th}_B \approx 0.1$). This comparison, performed with other reaction values $\mu(1) \in \{ i/20\}_{i=1}^{20}$, lead to the same conclusion, and consequently, we can deduce that if ${\rm Th}_B < 0.1$, the stability results obtained for the ODE and PDE systems \eqref{eq:1_2} and \eqref{eq:2_4} are similar. This result is consistent with the physics of the problem. Indeed, small values of ${\rm Th}_B$ correspond, for instance, to high diffusion coefficients implying almost spatial homogeneous biomass concentration. In this case, the dynamics in the reactor can be modeled with an ordinary differential equation cheaper to implement numerically (see \cite{Crespo2016}). \section{Conclusions}\label{sec:conclusions} In this work, we have performed an asymptotic analysis of a coupled system of two advection-diffusion-reaction equations with Danckwerts boundary conditions, which models the interaction between a microbial population and a diluted substrate in a continuous flow bioreactor. First, we have showed that for the particular case where the diffusion coefficients are negligible, after some finite time, the biomass becomes extinct and no reaction is produced (this state is usually called \textit{washout}). Next, we have studied the case when the diffusion coefficients are not negligible, and in this case the system exhibits, under suitable conditions, two stable equilibrium states: the \textit{washout} state and another steady state, which corresponds to the partial elimination of substrate. We have also taken into account that, depending on the reaction function, the system may exhibit either single stability or bistability. We have used the method of linearization to give a sufficient condition for the linear asymptotic stability of the washout equilibrium, and used this result, together with numerical experiments, to conjecture a sufficient condition for the linear asymptotic stability of the other stable equilibrium solution. These conditions were written in terms of nondimensional parameters ${\rm Da}$ (Damkh\"oler Number, relating reaction and advective rates), ${\rm Th}_B$ (Thi\'ele Modulus, relating reaction and biomass diffusion rates) and $\mu(1)$ (nondimensional reaction rate). Finally, our asymptotic stability results have been validated numerically and compared to the stability analysis results associated to the continuous bioreactor when it is modeled with ordinary differential equations. Results seem to indicate that the stability analysis results for the ODE are also valid for values of Thi\'ele Modulus (${\rm Th}_B$) lower than $0.1$, but not valid for values of Thi\'ele Modulus above this value. \begin{remark}\rm The Editor in charge of the review process of this work pointed out that the stability results presented here (i.e., the stability in the sense of $L^{2}$ norm) could be improved to some sharper functional spaces. A brief description of the arguments needed to obtain such a stability property is proposed in \cite{Ildefonso2017} and published in the same issue of this article. \end{remark} \subsection*{Acknowledgements} This work was carried out thanks to the financial support of the Spanish ``Ministry of Economy and Competitiveness'' under projects MTM2011-22658 and MTM2015-64865-P; the research group MOMAT (Ref. 910480) supported by ``Banco Santander'' and ``Universidad Complutense de Madrid''; and the ``Junta de Andaluc\'ia'' and the European Regional Development Fund through project P12-TIC301. The authors thank the French Labex Numev (convention ANR-10-LABX-20) for the postdoctoral grant of the first author at UMR MISTEA, Montpellier, France. \begin{thebibliography}{99} \bibitem{Andrews} J. F. Andrews; \newblock A mathematical model for the continuous culture of microorganisms utilizing inhibitory substrates. \newblock {\em Biotechnology and Bioengineering}, 10(6):707--723, 1968. \bibitem{Aris1974} R. Aris; \newblock Phenomena of multiplicity, stability and symmetry. \newblock {\em Annals of the New York Academy of Sciences}, 231(1):86--98, 1974. \bibitem{aris1975} R. Aris; \newblock {\em The Mathematical Theory of Diffusion and Reaction in Permeable Catalysts: Vol. 1: The theory of the steady state}. \newblock Oxford University Press, USA, 1975. \bibitem{asmar2005} N. H. Asmar; \newblock {\em Partial Differential Equations with Fourier Series and Boundary Value Problems}. \newblock Pearson Prentice Hall, 2005. \bibitem{bailey1986} J. E. Bailey and D. F. Ollis; \newblock {\em Biochemical Engineering Fundamentals}. \newblock Chemical engineering. McGraw-Hill, 1986. \bibitem{Bush1976} A. W. Bush, A. E. Cool; \newblock The effect of time delay and growth rate inhibition in the bacterial treatment of wastewater. \newblock {\em Journal of Theoretical Biology}, 63(2):385--395, 1976. \bibitem{crespo2015} M. Crespo, B. Ivorra, A.M. Ramos; \newblock Existence and uniqueness of solution of a continuous flow bioreactor model with two species. \newblock {\em Real Academia de Ciencias Exactas, F\'isicas y Naturales. Serie A: Matem\'aticas}, 110(2):357--377, 2016. \bibitem{Crespo2016} M. Crespo, A. M. Ramos, B. Ivorra, A. Rapaport; \newblock Modeling and optimization of activated sludge bioreactors for wastewater treatment taking into account spatial inhomogeneities. \newblock {\em Journal of Process Control}, 54:118--128, 2017. \bibitem{Danckwerts1953} P. V. Danckwerts; \newblock Continuous flow systems. \newblock {\em Chemical Engineering Science}, 2(1):1--13, 1953. \bibitem{Ildefonso2017} J. I. D\'iaz; \newblock Remarks on a precedent paper by Crespo, Ivorra and Ramos on the stability of bioreactor processes. \newblock This issue of EJDE-2017. \bibitem{Ildefonso2015} J. I. D\'iaz, A. C. Casal; \newblock On the principle of pseudo-linearized stability: {A}pplications to some delayed nonlinear parabolic equations. \newblock {\em Nonlinear Analysis}, 63(5):e997--e1007, 2005. \bibitem{dochain2010} D. Dochain; \newblock {\em Automatic Control of Bioprocesses}. \newblock ISTE. Wiley, 2010. \bibitem{dochain2001} D. Dochain and P. Vanrolleghem. \newblock {\em Dynamical Modelling and Estimation in Wastewater Treatment Processes}. \newblock IWA Publishing, 2001. \bibitem{Drame2004} A. K. Dram\'e; \newblock A semilinear parabolic boundary-value problem in bioreactors theory. \newblock {\em Electronic Journal of Differential Equations}, 2004. \bibitem{Drame2008} A. K. Dram\'e, C. Lobry, J. Harmand, A. Rapaport, F. Mazenc; \newblock Multiple stable equilibrium profiles in tubular bioreactors. \newblock {\em Mathematical and Computer Modelling}, 48(11–12):1840--1853, 2008. \bibitem{Fischer1975} A. E. Fischer, J. E. Marsden; \newblock Linearization stability of nonlinear partial differential equations. \newblock {\em Proceedings of Symposia in Pure Mathematics}, 27:219--263, 1975. \bibitem{Rapaport2011} P. Gajardo, J. Harmand, H.C. Ram\'irez, A. Rapaport; \newblock Minimal time bioremediation of natural water resources. \newblock {\em Automatica}, 47(8):1764--1769, 2011. \bibitem{grivet} J. P. Grivet; \newblock Nonlinear population dynamics in the chemostat. \newblock {\em Computing in Science Engineering}, 3, 2001. \bibitem{hale2010} J. K. Hale; \newblock {\em Asymptotic Behavior of Dissipative Systems}. \newblock Mathematical Surveys and Monographs. American Mathematical Society, 2010. \bibitem{herbert1956} D. Herbert, R. Elsworth, R. C. Telling; \newblock The continuous culture of bacteria; a theoretical and experimental study. \newblock {\em Journal of general microbiology}, 1956. \bibitem{Infante2015} J. A. Infante, J. M. Rey; \newblock {\em M\'etodos Num\'ericos. Teor\'ia, problemas y pr\'acticas con {MATLAB}.} \newblock Ediciones Pir\'amide. Grupo Anaya, 2015. \bibitem{Ishihara2007} S. Ishihara, M. Otsuji, A. Mochizuki; \newblock Transient and steady state of mass-conserved reaction-diffusion systems. \newblock {\em Phys. Rev. E}, 75:015203, 2007. \bibitem{Knapp2005} A. W. Knapp; \newblock {\em Advanced Real Analysis}. \newblock Cornerstones. Birkh{\"a}user Boston, 2005. \bibitem{Kuiper1980} H. J. Kuiper; \newblock Invariant sets for nonlinear elliptic and parabolic systems. \newblock {\em Journal on Mathematical Analysis}, 11(6):1075–--1103, 1980. \bibitem{Lanzani2004} L. Lanzani, Z. W. Shen; \newblock On the {R}obin boundary condition for {L}aplace's equation in {L}ipschitz domains. \newblock {\em Communications in Partial Differential Equations}, 29(1 \& 2): 91--109, 2004. \bibitem{liu2009} W. Liu; \newblock{\em Elementary Feedback Stabilization of the Linear Reaction-Convection-Diffusion Equation and the Wave Equation}. \newblock Math{\'e}matiques et Applications. Springer Berlin Heidelberg, 2009. \bibitem{marin2005advances} G. B. Marin; \newblock Advances in {C}hemical {E}ngineering: Multiscale analysis. \newblock {\em Advances in Chemical Engineering}, 2005. \bibitem{Martin1978} R. H. Martin; \newblock Asymptotic behavior for semilinear differential equations in {B}anach spaces. \newblock {\em SIAM Journal on Mathematical Analysis}, 9(6):1105--1119, 1978. \bibitem{matano1978} H. Matano; \newblock Convergence of solutions of one-dimensional semilinear parabolic equations. \newblock {\em J. Math. Kyoto Univ.}, 18(2):221--227, 1978. \bibitem{matano1979} H. Matano; \newblock {Asymptotic behavior and stability of solutions of semilinear diffusion equations.} \newblock {\em {Publ. Res. Inst. Math. Sci.}}, 15:401--454, 1979. \bibitem{McGowin1971} C. R. McGowin, D. D. Perlmutter; \newblock Tubular reactor steady state and stability characteristics. {I}. {E}ffect of axial mixing. \newblock {\em AIChE Journal}, 17(4):831--837, 1971. \bibitem{moreno1999} J. Moreno; \newblock Optimal time control of bioreactors for the wastewater treatment. \newblock {\em Optimal Control Applications and Methods}, 20(3):145--164, 1999. \bibitem{morita2010} Y. Morita, T. Ogawa; \newblock Stability and bifurcation of nonconstant solutions to a reaction – diffusion system with conservation of mass. \newblock {\em Nonlinearity}, 23(6):1387, 2010. \bibitem{Pao1981} C. V. Pao; \newblock Asymptotic stability of reaction-diffusion systems in chemical reactor and combustion theory. \newblock {\em Journal of Mathematical Analysis and Applications}, 82(2):503--526, 1981. \bibitem{pao1992} C. V. Pao; \newblock {\em Nonlinear Parabolic and Elliptic Equations}. \newblock Fems Symposium. Springer US, 1992. \bibitem{vazquez2014} C. Par{\'e}s, C. V{\'a}zquez, F. Coquel; \newblock {\em Advances in Numerical Simulation in Physics and Engineering: Lecture Notes of the XV 'Jacques-Louis Lions' Spanish-French School}. \newblock SEMA SIMAI Springer Series. Springer International Publishing, 2014. \bibitem{morphogenesis1995} A. J. {Perumpanani}, J. A. {Sherratt}, P. K. {Maini}; \newblock {Phase differences in reaction-diffusion-advection systems and applications to morphogenesis}. \newblock {\em IMA Journal of Applied Mathematics}, 55: 19--33, 1995. \bibitem{angel} A. M. Ramos; \newblock {\em Introducci\'on al An\'alisis Matem\'atico del M\'etodo de Elementos Finitos}. \newblock Editorial Complutense, 2012. \bibitem{Rao2009} V. S. H. Rao, P. {Raja Sekhara Rao}; \newblock Basic chemostat model revisted. \newblock {\em Differential Equations and Dynamical Systems}, 17(1): 3--16, 2009. \bibitem{Harmand2015} A. Rapaport, I. Haidar, J. Harmand; \newblock {Global dynamics of the buffered chemostat for a general class of response functions}. \newblock {\em {Journal of Mathematical Biology}}, 71(1): 69--98, 2015. \bibitem{Siebert2014} J. Siebert, S. Alonso, M. B\"ar, E. Sch\"oll; \newblock Dynamics of reaction-diffusion patterns controlled by asymmetric nonlocal coupling as a limiting case of differential advection. \newblock {\em Phys. Rev. E}, 89: 052909, 2014. \bibitem{SmithWaltman} H. Smith, P. Waltman; \newblock {\em The theory of the Chemostat. In Cambridge studies in Mathematical biology}, volume 13. \newblock Cambridge: Cambridge University Press, 1995. \bibitem{Smith2014} N. A. S. Smith, S. L. Mitchell, A. M. Ramos; \newblock Analysis and simplification of a mathematical model for high-pressure food processes. \newblock {\em Applied Mathematics and Computation}, 226: 20--37, 2014. \bibitem{Vetjasa1970} S. A. Vejtasa, R. A. Schmitz; \newblock An experimental study of steady state multiplicity and stability in an adiabatic stirred reactor. \newblock {\em AIChE Journal}, 16(3): 410--419, 1970. \bibitem{Weinberger} H. F. Weinberger; \newblock {\em A First Course in Partial Differential Equations with Complex Variables and Transform Methods}. \newblock Dover Publications, 1995. \bibitem{Dochain2000} J. J. Winkin, D. Dochain, P. Ligarius; \newblock Dynamical analysis of distributed parameter tubular reactors. \newblock {\em Automatica}, 36(3): 349--361, 2000. \end{thebibliography} \end{document}