Asymptotic stability of a coupled
advection-diffusion-reaction system arising in bioreactor
processes

Maria Crespo Moya, Benjamin Ivorra, Angel Manuel Ramos

► To cite this version:

Maria Crespo Moya, Benjamin Ivorra, Angel Manuel Ramos. Asymptotic stability of a coupled advection-diffusion-reaction system arising in bioreactor processes. Electronic Journal of Differential Equations, 2017, 194, pp.1-26. hal-01607898

HAL Id: hal-01607898

https://hal.science/hal-01607898

Submitted on 25 May 2020

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
ASYMPTOTIC STABILITY OF A COUPLED
ADVECTION-DIFFUSION-REACTION SYSTEM ARISING IN
BIOREACTOR PROCESSES

MARÍA CRESPO, BENJAMIN IVORRA, ÁNGEL MANUEL RAMOS

Communicated by Jesus Ildefonso Diaz

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.

# 1. INTRODUCTION

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 [5, 13], among which continuous flow bioreactors are commonly used in the bioremediation of water resources (see, for instance, [17, 20, 32]). These biological reactors are filled from a polluted resource with a flow rate $Q$ (m³/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 [42]

$$ 
\begin{aligned}
\frac{dS}{dt} &= -\mu(S)B + \frac{Q(t)}{V}(S_e(t) - S) \quad t > 0, \\
\frac{dB}{dt} &= \mu(S)B - \frac{Q(t)}{V}B \quad t > 0,
\end{aligned}
\quad (1.1) $$

where $S$ (kg/m³) and $B$ (kg/m³) are the concentrations of substrate and biomass, respectively; $S_e(t)$ (kg/m³) is the concentration of substrate that enters the reactor

<footnote n="2010">2010 Mathematics Subject Classification. 35B30, 35B35, 35B40, 35K51, 35Q35.</footnote>

<footnote>Key words and phrases. Asymptotic stability; advection-diffusion-reaction system; dimensional analysis; continuous bioreactor.</footnote>

<footnote>©2017 Texas State University.</footnote>

<footnote>Submitted March 1, 2017. Published August 7, 2017.</footnote>
at time $t$; $V$ (m³) is the reactor volume; $Q(t)$ (m³/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 [12, 13, 40]):

The function $\mu: [0, +\infty) \to [0, +\infty)$ satisfies $\mu(0) = 0$, $\bar{\mu} \ge \mu(z) > 0$ for $z > 0$ and one of the following two properties:

(A1) $\mu$ is increasing and concave.

(A2) There exists $s > 0$ such that $\mu$ is increasing on $(0, s)$ and decreasing on $(s, +\infty)$.

The Monod function [42], defined by

$$\mu(S) = \mu_{\max} \frac{S}{K_S + S},$$

satisfies (A1), and the Haldane function [1], described by

$$\mu(S) = \mu^* \frac{S}{K_S + S + S^2/K_I},$$

satisfies (A2). Both functions are extensively used in the literature.

In the particular case when $S_e$ and $Q$ are constant, system (1.1) can be non-dimensionalized by setting $\hat{S} = S/S_e$, $\hat{B} = B/S_e$, $\hat{\mu}(\hat{S}) = \mu(S_e\hat{S})/\bar{\mu}$ and $\hat{t} = \bar{\mu}t$. For simplicity, we drop the $\hat{}$ notation, and so $S$, $B$, $\mu$ and $t$ denote the non-dimensional variables. System (1.1) in non-dimensional form is therefore given by

$$\frac{dS}{dt} = -\mu(S)B + d(1 - S) \quad \text{in } (0, +\infty), \qquad (1.2)$$

$$\frac{dB}{dt} = \mu(S)B - dB \quad \text{in } (0, +\infty),$$

where $d = Q/V\bar{\mu}$ is the dimensionless dilution rate.

If $\mu$ satisfies (A1), system (1.2) has two equilibrium configurations $(S_1^*, B_1^*) = (1, 0)$, usually called washout, and $(S_2^*, B_2^*) = (S_2^*, 1 - S_2^*)$, where $S_2^*$ is such that $\mu(S_2^*) = d$ (see [42]). In [18, 39, 42] the authors conclude that the steady state $(1, 0)$ is asymptotically stable if $d \ge \mu(1)$, while the steady state $(S_2^*, 1 - S_2^*)$ is asymptotically stable if $d < \mu(1)$.

Similarly (see [1]), if $\mu$ satisfies (A2), system (1.2) has three equilibrium configurations $(S_1^*, B_1^*) = (1, 0)$, $(S_2^*, B_2^*) = (S_2^*, 1 - S_2^*)$ and $(S_3^*, B_3^*) = (S_3^*, 1 - S_3^*)$, where $\mu(S_2^*) = \mu(S_3^*) = d$ and $S_2^* < S_3^*$. In [6], [12] and [40] the authors show that the steady state $(1, 0)$ is asymptotically stable if $d > \mu(1)$, the steady state $(S_2^*, 1 - S_2^*)$ is asymptotically stable if $d < 1$ and the steady state $(S_3^*, 1 - S_3^*)$ is unstable. Thereby, if $\mu(1) < 1$, there is bistability when $\mu(1) < d < 1$.

System (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 [3, 8]). Particularly, system (1.2) can be improved by considering a coupled system of spatio-temporal parabolic equations of the form

$$ \begin{aligned} \frac{dS}{dt} &= L_S(S) - \mu(S)B \quad \text{in } \Omega \times (0, +\infty), \\ \frac{dB}{dt} &= L_B(B) + \mu(S)B \quad \text{in } \Omega \times (0, +\infty), \end{aligned} \qquad (1.3) $$

where $\Omega$ is the bioreactor domain and $L_S, L_B$ are linear second order elliptic partial differential operators on $\Omega$. Notice that system (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 (1.1) is modeled by the terms $d(1-S)$ and $-dB$, respectively. Particularly, models (2.1) and (2.4) presented in Section 2 include the so called Danckwerts boundary conditions, typically used for continuous flow systems (see [9, 14, 15]), 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 (1.3) should provide more accurate results, compared with the asymptotic ones detailed above for system (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 [11, 19, 24, 26, 28, 29, 30]. Most theoretical studies focusing on bioreactor processes consider the assumption that both $L_S$ and $L_B$ are diffusion operators (see, e.g. [22, 24, 33, 34]). For instance, in Yosida and Morita [33], 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 [14, 15, 31, 37, 41, 46]. In [14, 15, 31, 46], the authors study system (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 [37, 41], 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 [14, 31, 37, 41, 46], 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 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 3, we present the steady states of the system and analyze the asymptotic stability using linearization methods. Then, Section 4 presents numerical experiments to analyze the validity and robustness of the stability results obtained in Section 3. Finally, we perform a comparison with the asymptotic results related to system (1.2).

## 2. MATHEMATICAL MODELING

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^*$ (see Figure 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 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_{in}$ (i.e., the upper boundary of the rectangle $\Omega$). Treated water leaves the reactor through the outlet $\Gamma_{out}$ (i.e., the lower boundary of the rectangle $\Omega$). Moreover, $\Gamma_{sym} = \{0\} \times (0, H)$ is the axis of symmetry and $\Gamma_{wall} = \delta\Omega \setminus (\Gamma_{in} \cup \Gamma_{out} \cup \Gamma_{sym})$ is the bioreactor wall for which no flux passes through.

<figure><img src="image_2.png" /><figcaption>FIGURE 1. Typical representation of the domain geometry.</figcaption></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 [7]:

$$
\begin{align*}
\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} + uS &= uS_e \quad \text{in } \Gamma_{\text{in}} \times (0, T), \\
D_B \frac{\partial B}{\partial z} + uB &= 0 \quad \text{in } \Gamma_{\text{in}} \times (0, T), \\
D_S \frac{\partial S}{\partial r} &= 0 \quad \text{in } \Gamma_{\text{sym}} \times (0, T), \\
D_B \frac{\partial B}{\partial r} &= 0 \quad \text{in } \Gamma_{\text{sym}} \times (0, T), \tag{2.1} \\
D_S \frac{\partial S}{\partial r} &= 0 \quad \text{in } \Gamma_{\text{wall}} \times (0, T), \\
D_B \frac{\partial B}{\partial r} &= 0 \quad \text{in } \Gamma_{\text{wall}} \times (0, T), \\
D_S \frac{\partial S}{\partial z} &= 0 \quad \text{in } \Gamma_{\text{out}} \times (0, T), \\
D_B \frac{\partial B}{\partial z} &= 0 \quad \text{in } \Gamma_{\text{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{align*}
$$

where $T > 0$ (s) is the length of the time interval for which we want to model the process; $S$ (kg/m³) and $B$ (kg/m³) are the substrate and biomass concentrations inside the bioreactor, which diffuse throughout the water in the vessel with diffusion coefficients $D_S$ (m²/s) and $D_B$ (m²/s), respectively; the fluid velocity is taken as $\mathbf{u} = (0, 0, -u)$, where $u$ (m/s) is the flow speed; $S_e$ (kg/m³) is the concentration of substrate that enters into the bioreactor; $S_0$ (kg/m³) and $B_0$ (kg/m³) are the initial concentrations of substrate and biomass inside the bioreactor, respectively. Furthermore, as in system (1.1), we consider a term corresponding to the reaction between biomass and substrate, governed by the growth rate function $\mu$ (s⁻¹).

**Remark 2.1.** According to [7], if $\mu \in L^\infty(0, +\infty)$ is continuous and Lipschitz, $u \in L^\infty(\hat{\Omega} \times (0, T))$, $S_e \in L^\infty(0, T)$, $S_e \ge 0$ in $(0, T)$, $S_0 \in L^\infty(\Omega)$, $S_0 \ge 0$ in $\Omega$, $B_0 \in L^\infty(\Omega)$ and $B_0 \ge 0$ in $\Omega$, there exists a unique solution $(S, B) \in L^2(0, T, H^1(\Omega))^2 \cap C([0, T], L^2(\Omega))^2 \cap L^\infty(\Omega \times (0, T))^2$ of system (2.1).

We perform a dimensional analysis of the model, following similar methodology as the one presented in [43]. System (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_e} = S_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 \le \hat{t} \le \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 (2.1) become

$$
\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} \quad (2.2)
$$
and

$$ \frac{\partial \hat{B}}{\partial 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}. \quad (2.3) $$

The dimensionless groups of parameters in equations (2.2) and (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_e\|_{L^\infty(0,T)}$ and, for the sake of simplicity, we choose $s = b = \|S_e\|_{L^\infty(0,T)}$. The time scale $\tau$ is chosen from equations (2.2) and (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(\bar{\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 (1.2), we center our study on the reaction process and take $\tau = 1/\bar{\mu}$. Two well-known dimensionless numbers (see [27]) appear now in the non-dimensional form of system (2.1):

* Damkhöler Number: $Da = \frac{\text{reaction rate}}{\text{advective transport rate}} = \frac{\tau_a}{\tau_r}$,

* Thiele Modulus: $Th = \frac{\text{reaction rate}}{\text{diffusive transport rate}} = \frac{\tau_d}{\tau_r}$,

where $\tau_d$, $\tau_a$ and $\tau_r$ are diffusion, advection and reaction times scales, respectively. For ease of notation, we drop the $\hat{}$ symbol, and so $B$, $S$, $t$, $u$, $S_e$, $\mu$, $z$, $r$ and $T$ denote now the non-dimensional variables. Particularly, if $S_e$ and $u$ are constants, system (2.1) in its non-dimensional form is given by

$$ \frac{\partial S}{\partial t} = \frac{\sigma^2}{\text{Th}_S} \frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial S}{\partial r}) + \frac{1}{\text{Th}_S} \frac{\partial^2 S}{\partial z^2} + \frac{1}{\text{Da}} \frac{\partial S}{\partial z} - \mu(S)B \quad \text{in } \Omega \times (0, T), $$

$$ \frac{\partial B}{\partial t} = \frac{\sigma^2}{\text{Th}_B} \frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial B}{\partial r}) + \frac{1}{\text{Th}_B} \frac{\partial^2 B}{\partial z^2} + \frac{1}{\text{Da}} \frac{\partial B}{\partial z} + \mu(S)B \quad \text{in } \Omega \times (0, T), $$

$$ \frac{1}{\text{Th}_S} \frac{\partial S}{\partial z} + \frac{1}{\text{Da}} S = \frac{1}{\text{Da}} \quad \text{in } \Gamma_{\text{in}} \times (0, T), $$ $$ \frac{1}{\text{Th}_B} \frac{\partial B}{\partial z} + \frac{1}{\text{Da}} B = 0 \quad \text{in } \Gamma_{\text{in}} \times (0, T), \quad (2.4) $$

$$ \frac{\partial S}{\partial r} = \frac{\partial B}{\partial r} = 0 \quad \text{in } \Gamma_{\text{sym}} \times (0, T), $$

$$ \frac{\partial S}{\partial r} = \frac{\partial B}{\partial r} = 0 \quad \text{in } \Gamma_{\text{wall}} \times (0, T), $$

$$ \frac{\partial S}{\partial z} = \frac{\partial B}{\partial z} = 0 \quad \text{in } \Gamma_{\text{out}} \times (0, T), $$

complemented by the initial conditions

$$ S(\cdot, \cdot, 0) = S_{\text{init}} \quad \text{and} \quad B(\cdot, \cdot, 0) = B_{\text{init}} \quad \text{in } \Omega, \quad (2.5) $$

where $\Omega = (0, 1) \times (0, 1)$ is the nondimensional domain, $\Gamma_{\text{in}} = (0, 1) \times \{1\}$, $\Gamma_{\text{out}} = (0, 1) \times \{0\}$, $\Gamma_{\text{wall}} = \{1\} \times (0, 1)$ and $\Gamma_{\text{sym}} = \{0\} \times (0, 1)$ are the non-dimensional
boundary edges. The final dimensionless parameters are Da = Hµ̄/u, Th_S = H²µ̄/D_S, Th_B = D_STh_S/D_B and σ = H/L, and the dimensionless initial conditions are S_init = S₀/S_e and B_init = B₀/S_e.

**Remark 2.2.** Since the bioreactor in consideration is a cylinder of height H and radius L, the reactor volume is πHL² and the volumetric flow rate in system (1.1) can be written as Q = πL²u, where u (m/s) is the vertical inflow. Thus, the nondimensional dilution rate d in system (1.2) corresponds to the nondimensional flow rate (1/Da) in system (2.4).

### 3. STEADY STATES AND STABILITY ANALYSIS

In this section, we study the asymptotic behavior of system (2.4)-(2.5). Firstly, we study the particular case for which diffusion terms in system (2.4) are neglected. Then, we perform the stability analysis of system (2.4) for the general case.

The asymptotic stability of an equilibrium solution of system (2.4) is defined as follows (see [35]).

**Definition 3.1.** An equilibrium solution $(S^*, B^*)$ of system (2.4) is said to be asymptotically stable if there exists $\epsilon > 0$ such that if given $(S_{\text{init}}, B_{\text{init}}) \in (L^\infty(\Omega))^2$ satisfying

$$ \|(S_{\text{init}}, B_{\text{init}}) - (S^*, B^*)\|(L^2(\Omega))^2 < \epsilon, $$

then the corresponding unique solution $(S, B)$ of system (2.4)-(2.5) satisfies

$$ \lim_{t \to \infty} \|(S(t), B(t)) - (S^*, B^*)\|(L^2(\Omega))^2 = 0. $$

We point out that we require $(S_{\text{init}}, B_{\text{init}}) \in (L^\infty(\Omega))^2$ to ensure the existence and uniqueness of solution of system (2.4)-(2.5) (see Remark 2.1).

3.1. **Case** $(1/\text{Th}_S)$, $(1/\text{Th}_B)$, $(\sigma^2/\text{Th}_S)$, $(\sigma^2/\text{Th}_B) \ll 1$ **Lg.** We consider the particular case where the nondimensional diffusion coefficients are negligible with respect to the advection and reaction coefficients in system (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{aligned} \frac{\partial S}{\partial t} &= \frac{1}{\text{Da}} \frac{\partial S}{\partial z} - \mu(S)B \quad \text{in } (0,1) \times (0,T), \\ \frac{\partial B}{\partial t} &= \frac{1}{\text{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_{\text{init}}(r, \cdot) \quad \text{in } (0,1), \\ B(r, \cdot, 0) &= B_{\text{init}}(r, \cdot) \quad \text{in } (0,1). \end{aligned} \quad (3.1) $$

Let us prove that (1,0) (which is called the washout state), is an asymptotically stable equilibrium. The following theorem shows, in fact, a property for (1,0) stronger than asymptotic stability.

**Theorem 3.2.** For any arbitrary initial condition $(S_{\text{init}}, B_{\text{init}}) \in (L^\infty(\Omega))^2$, the solution of (3.1) satisfies that $S(r, z, t)= 1$ and $B(r, z, t) = 0$, for all $(r, z) \in \Omega$ and $t \ge Da$.
**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/\text{Da})$, so that for every fixed value of $(r, z) \in \Omega$, the second equation of system (3.1) is rewritten as

$$ \begin{aligned} \frac{d\mathcal{B}}{dt}(r, \tilde{z}(t, z), t) &= \frac{\partial \mathcal{B}}{\partial t}(r, \tilde{z}(t, z), t) - \frac{1}{\text{Da}} \frac{\partial \mathcal{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{aligned} $$

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)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)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 (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)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/\text{Da}, t) = 0 \text{ and } S(r, 1 - t/\text{Da}, t) = 1 \text{ for all } t > 0. $$

Consequently, if $t \ge \text{Da}$, $B(r, z, t) = 0$ and $S(r, z, t) = 1$ for all $(r, z) \in \Omega$. $\square$

<figure><img src="image_3.png" /></figure>

## 3.2. General Case.
To obtain a parallelism with the asymptotic analysis of system (1.2), shown in Section 1, we assume that $\mu$ satisfies properties (A1) or (A2). In both cases, the constant (washout) solution $(S_1^*, B_1^*) = (1, 0)$ is a steady state of system (2.4). By analogy with system (1.2), we conjecture, supported by numerical experiments, that system (2.4) has, under suitable conditions, another asymptotically stable steady state (different from the washout) denoted by $(S_2^*, B_2^*)$. In this section, we use linearization techniques to study the stability of the steady states. First, we introduce the concept of 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^*, B^*)$ be an equilibrium solution of system (2.4). We define the concept of linear asymptotic stability by following the reasoning below.

Let us consider system (2.4) with initial conditions in $(L^\infty(\Omega))^2$ and close to $(S^*, B^*)$ given by $S(r, z, 0) = S^* + \delta S_{\text{init}} \ge 0$ and $B(r, z, 0) = B^* + \delta B_{\text{init}} \ge 0$, with $\|\delta S_{\text{init}}\|_{L^2(\Omega)} \ll 1$ and $\|\delta B_{\text{init}}\|_{L^2(\Omega)} \ll 1$. The solution of system (2.4) can be seen as

$$ \begin{pmatrix} S(r, z, t) \\ B(r, z, t) \end{pmatrix} = \begin{pmatrix} S^*(r, z) \\ B^*(r, z) \end{pmatrix} + \begin{pmatrix} \delta S(r, z, t) \\ \delta B(r, z, t) \end{pmatrix} + \text{Higher order terms}, \quad (3.2) $$
(see [16, page 220]) where

$$
\begin{align}
\frac{d\delta S}{dt} &= \frac{\sigma^2}{\text{Th}_S} \frac{1}{r} \frac{d}{dr} \left( r \frac{d\delta S}{dr} \right) + \frac{1}{\text{Th}_S} \frac{d^2\delta S}{dz^2} + \frac{1}{\text{Da}} \frac{d\delta S}{dz} \nonumber \\
&\quad - \mu(S^*)\delta B - \mu'(S^*)B^*\delta S \quad \text{in } \Omega \times (0, T), \\
\frac{d\delta B}{dt} &= \frac{\sigma^2}{\text{Th}_B} \frac{1}{r} \frac{d}{dr} \left( r \frac{d\delta B}{dr} \right) + \frac{1}{\text{Th}_B} \frac{d^2\delta B}{dz^2} + \frac{1}{\text{Da}} \frac{d\delta B}{dz} \nonumber \\
&\quad + \mu(S^*)\delta B + \mu'(S^*)B^*\delta S \quad \text{in } \Omega \times (0, T), \\
\frac{1}{\text{Th}_S} \frac{d\delta S}{dz} + \frac{1}{\text{Da}} \delta S &= \frac{1}{\text{Th}_B} \frac{d\delta B}{dz} + \frac{1}{\text{Da}} \delta B = 0 \quad \text{in } \Gamma_{\text{in}} \times (0, T), \tag{3.3} \\
\frac{d\delta S}{dr} &= \frac{d\delta B}{dr} = 0 \quad \text{in } \Gamma_{\text{sym}} \times (0, T), \\
\frac{d\delta S}{dr} &= \frac{d\delta B}{dr} = 0 \quad \text{in } \Gamma_{\text{wall}} \times (0, T), \\
\frac{d\delta S}{dz} &= \frac{d\delta B}{dz} = 0 \quad \text{in } \Gamma_{\text{out}} \times (0, T), \\
\delta S(\cdot, \cdot, 0) &= \delta S_{\text{init}} \quad \text{in } \Omega, \\
\delta B(\cdot, \cdot, 0) &= \delta B_{\text{init}} \quad \text{in } \Omega.
\end{align}
$$

**Definition 3.3.** An equilibrium solution $(S^*, B^*)$ of system (2.4) is said to be linearly asymptotically stable if there exists $\epsilon > 0$ such that if given $(\delta S_{\text{init}}, \delta B_{\text{init}}) \in (L^{\infty}(\Omega))^2$ satisfying $\|(\delta S_{\text{init}}, \delta B_{\text{init}})\|_{(L^2(\Omega))^2} < \epsilon$, then the corresponding unique solution $(\delta S, \delta B)$ of system (3.3) satisfies $\lim_{t \to \infty} \|(\delta S(t), \delta B(t))\|_{(L^2(\Omega))^2} = 0$.

We now define the following functions, which will be used through the rest of the manuscript.

**Definition 3.4.** • In terms of the dimensionless variables appearing in system (2.4), we define $\beta_1(\text{Da}, \text{Th}_B)$ as the smallest positive solution of the transcendental equation $\tan(\beta) = \frac{\text{Th}_B\beta}{\text{Da}(\beta^2 - (\frac{\text{Th}_B}{2\text{Da}})^2)}$ if $\text{Th}_B \neq \pi\text{Da}$. If $\text{Th}_B = \pi\text{Da}$, we define $\beta_1(\text{Da}, \text{Th}_B) = \pi/2$.

• In terms of the variables with dimensions appearing in system (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(\beta^2 - (\frac{Hu}{2D_B})^2)} $$

if $Hu \neq \pi D_B$. If $Hu = \pi D_B$ we define $\tilde{\beta}_1(H, u, D_B) = \pi/2$.

**Theorem 3.5.** A sufficient condition for $(S_1^*, B_1^*) = (1, 0)$ to be a linearly asymptotically stable steady state of system (2.4) is that

$$ \mu(1) < \frac{\text{Th}_B}{(2\text{Da})^2} + \frac{(\beta_1(\text{Da}, \text{Th}_B))^2}{\text{Th}_B}. \qquad (3.4) $$

**Remark 3.6.** In terms of the variables with dimensions appearing in system (2.1), the steady state is $(S_e, 0)$ and inequality (3.4) is reformulated as

$$ \mu(S_e) < \frac{u^2}{4D_B} + \frac{D_B}{H^2}(\tilde{\beta}_1(H, u, D_B))^2. $$
*Proof of Theorem 3.5.* To check the stability of the washout equilibrium solution, we replace $(S^*, B^*)$ by $(1,0)$ in System (3.3) and, so, functions $\delta S$ and $\delta B$ fulfill

$$ 
\begin{aligned}
\frac{d\delta S}{dt} &= \frac{\sigma^2}{\text{Th}_S} \frac{1}{r} \frac{d}{dr} (r \frac{d\delta S}{dr}) + \frac{1}{\text{Th}_S} \frac{d^2\delta S}{dz^2} + \frac{1}{\text{Da}} \frac{d\delta S}{dz} - \mu(1)\delta B \\
&\qquad \text{in } \Omega \times (0, T), \\
\frac{d\delta B}{dt} &= \frac{\sigma^2}{\text{Th}_B} \frac{1}{r} \frac{d}{dr} (r \frac{d\delta B}{dr}) + \frac{1}{\text{Th}_B} \frac{d^2\delta B}{dz^2} + \frac{1}{\text{Da}} \frac{d\delta B}{dz} + \mu(1)\delta B \\
&\qquad \text{in } \Omega \times (0, T), \\
\frac{1}{\text{Th}_S} \frac{d\delta S}{dz} + \frac{1}{\text{Da}}\delta S &= \frac{1}{\text{Th}_B} \frac{d\delta B}{dz} + \frac{1}{\text{Da}}\delta B = 0 \quad \text{in } \Gamma_{\text{in}} \times (0, T), \\
\frac{d\delta S}{dr} &= \frac{d\delta B}{dr} = 0 \quad \text{in } \Gamma_{\text{sym}} \times (0, T), \\
\frac{d\delta S}{dr} &= \frac{d\delta B}{dr} = 0 \quad \text{in } \Gamma_{\text{wall}} \times (0, T), \\
\frac{d\delta S}{dz} &= \frac{d\delta B}{dz} = 0 \quad \text{in } \Gamma_{\text{out}} \times (0, T), \\
\delta S(\cdot, \cdot, 0) &= \delta S_{\text{init}} \quad \text{in } \Omega, \\
\delta B(\cdot, \cdot, 0) &= \delta B_{\text{init}} \quad \text{in } \Omega,
\end{aligned}
\quad (3.5) $$

with $\|\delta S_{\text{init}}\|_{L^2(\Omega)} \ll 1$ and $\|\delta B_{\text{init}}\|_{L^2(\Omega)} \ll 1$. We are going to prove that the steady state $(S_1^*, B_1^*) = (1,0)$ is linearly asymptotically stable by showing that (see Definition 3.3)

$$ \|\delta S(t)\|_{L^2(\Omega)} \to 0 \quad \text{and} \quad \|\delta B(t)\|_{L^2(\Omega)} \to 0 \quad \text{as} \quad t \to \infty. $$

**Step 1.** Let us prove that $\|\delta B(t)\|_{L^2(\Omega)} \to 0$ as $t \to \infty$: Notice that the equations involving the biomass in system (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). $$

**Step 1.1.** Separation of variables. From the second equation in system (3.5) one has

$$ \frac{T'(t)}{T(t)} = \frac{\sigma^2}{\text{Th}_B} \left( \frac{R''(r)}{R(r)} + \frac{1}{r} \frac{R'(r)}{R(r)} \right) + \frac{1}{\text{Th}_B} \frac{Z''(z)}{Z(z)} + \frac{1}{\text{Da}} \frac{Z'(z)}{Z(z)} + \mu(1). $$

If we equate this expression to a constant $\lambda$, it follows that

$$ \frac{\sigma^2}{\text{Th}_B} \left( \frac{R''(r)}{R(r)} + \frac{1}{r} \frac{R'(r)}{R(r)} \right) = - \frac{1}{\text{Th}_B} \frac{Z''(z)}{Z(z)} - \frac{1}{\text{Da}} \frac{Z'(z)}{Z(z)} + \lambda - \mu(1). $$

Equating this expression to an arbitrary constant $\eta$, one obtains

$$ \begin{aligned} & R''(r) + \frac{1}{r} R'(r) - \frac{\text{Th}_B}{\sigma^2} \eta R(r) = 0, \\ & \frac{1}{\text{Th}_B} Z''(z) + \frac{1}{\text{Da}} Z'(z) - (\lambda - \mu(1) - \eta) Z(z) = 0. \end{aligned} $$
Proceeding as in the proof of [7, Theorem 3], it is easy to see that

$$ \delta B(r, z, t) = |\delta B(r, z, t)| \le \|\delta B_{\text{init}}\|_{L^{\infty}(\Omega)} 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).

**Step 1.2.** Calculation of $R(r)$. Using the boundary conditions of system (3.5) on $\Gamma_{\text{wall}}$ and $\Gamma_{\text{sym}}$, it is clear that $R(r)$ is a solution of system

$$ \begin{aligned} R''(r) + \frac{1}{r}R'(r) - \frac{\text{Th}_B}{\sigma^2}\eta R(r) = 0 \quad r \in (0, 1), \\ R'(0) = R'(1) = 0. \end{aligned} \quad (3.6) $$

Using the change of variables $s = ar$, with $a = \sqrt{|\eta|\frac{\text{Th}_B}{\sigma^2}}$, the differential equation for $R$ can be rewritten in one of the following forms

$$ (1) \quad s^2 R''(s) + sR'(s) + s^2 R(s) = 0 \quad \text{if } \eta < 0, $$

$$ (2) \quad s^2 R''(s) + sR'(s) - s^2 R(s) = 0 \quad \text{if } \eta > 0, $$

$$ (3) \quad sR''(s) + R'(s) = 0 \quad \text{if } \eta = 0. $$

**Case 1:** $\eta < 0$. In this case the equation for $R(s)$ is known as the *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_1 J_0(s)$. It is well known that $J'_0(s) = -J_1(s)$ and $0 \in \{s \in [0, +\infty): J_1(s) = 0\}$, which is a countable set $\{T_n\}_{n \in \mathbb{N}}$ with an infinite number of elements (see, e.g., [4]). 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(\sqrt{-\frac{\eta \text{Th}_B}{\sigma^2}}) = 0$. Consequently, $\eta \in \{\eta_n\}_{n \in \mathbb{N}}$, with

$$ \eta_n = - \frac{(\sigma T_n)^2}{\text{Th}_B}, \quad (3.7) $$

and the solution $R(r)$ is given by

$$ R(r) = \sum_{n \in \mathbb{N}} C_n J_0\left(\frac{\sqrt{-\eta_n \text{Th}_B}}{\sigma} r\right). $$

**Case 2:** $\eta > 0$. In this case the equation for $R(s)$ is known as the *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_1 I_0(s)$. It is well known that $I'_0(s) = I_1(s)$ and the boundary condition at $s = a$ implies that the eigenvalues $\eta$ satisfy that

$$ C_1 I'_0\left(\frac{\sqrt{\eta \text{Th}_B}}{\sigma}\right) = 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.
**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_1e^{-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

$$E = \{0\} \cup \left\{ -\frac{(\sigma T_n)^2}{\text{Th}_B} \right\}_{n \in \mathbb{N}}, \qquad (3.8)$$

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\left(\frac{\sqrt{-\text{Th}_B \eta_n}}{\sigma} r\right).$$

**Step 1.3.** Calculation of $Z(z)$. Using the boundary conditions of system (3.5) on $\Gamma_{in}$ and $\Gamma_{out}$, it is clear that function $Z(z)$ is solution of system

$$\begin{align}
\frac{1}{\text{Th}_B} Z''(z) + \frac{1}{\text{Da}} Z'(z) - (\lambda - \mu(1) - \eta)Z(z) &= 0, \quad z \in (0, 1), \nonumber \\
\frac{1}{\text{Th}_B} Z'(1) + \frac{1}{\text{Da}} Z(1) &= 0, \tag{3.9} \\
Z'(0) &= 0, \nonumber
\end{align}$$

which corresponds to a regular Sturm-Liouville eigenvalue problem (see [23, Theorem 1.3]). The corresponding characteristic equation is

$$\frac{1}{\text{Th}_B} \rho^2 + \frac{1}{\text{Da}} \rho - (\lambda - \mu(1) - \eta) = 0,$$

with roots

$$\rho = \frac{-\text{Th}_B}{2\text{Da}} \pm \frac{\text{Th}_B}{2} \sqrt{\left(\frac{1}{\text{Da}}\right)^2 + \frac{4(\lambda - \mu(1) - \eta)}{\text{Th}_B}}.$$

Now, depending on the value of $\Delta = (\frac{1}{\text{Da}})^2 + \frac{4(\lambda - \mu(1) - \eta)}{\text{Th}_B}$, three possible solutions appear.

**Case 1:** $\Delta = 0 \Leftrightarrow \lambda = \eta + \mu(1) - \text{Th}_B(\frac{1}{2\text{Da}})^2$. In this case, the solution of system (3.9) is

$$Z(z) = D_1 e^{\alpha z} + D_2 z e^{\alpha z},$$

where $\alpha = \frac{-\text{Th}_B}{2\text{Da}}$ and $D_1, D_2$ are constants which are determined by the boundary conditions of the system. Since

$$Z'(z) = \alpha e^{\alpha z}(D_1 + zD_2) + D_2 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 e^{\alpha z}(1 - \alpha z) \quad \text{and} \quad Z'(z) = -D_1 \alpha^2 z e^{\alpha z}.$$

From the boundary condition at $z = 1$ it follows that

$$D_1 e^{\alpha} \left( \frac{1}{\text{Da}} (1 - \alpha) - \frac{\alpha^2}{\text{Th}_B} \right) = 0. \qquad (3.10)$$
By replacing $\alpha$ by its value into equation (3.10), we conclude that this equation is true either if $\text{Th}_B = -4\text{Da}$ or if $D_1 = 0$. The first option is not possible since constants Da and $\text{Th}_B$ are assumed strictly positive. Thus, the only solution in this case is $Z(z) = 0$.

**Case 2:** $\Delta < 0 \Leftrightarrow \lambda < \eta + \mu(1) - \text{Th}_B(\frac{1}{2\text{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 (3.9) is of the form

$$Z(z) = e^{\alpha z}(D_1 \cos(\beta z) + D_2 \sin(\beta z)),$$

where $D_1$ and $D_2$ are constants which will be determined by the boundary conditions. Since

$$Z'(z) = \alpha Z(z) + \beta e^{\alpha z}(-D_1 \sin(\beta z) + D_2 \cos(\beta z)),$$

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 e^{\alpha z} \left( \cos(\beta z) - \frac{\alpha}{\beta} \sin(\beta z) \right) \quad \text{and} \quad Z'(z) = -D_1 e^{\alpha z} \sin(\beta z) \left( \frac{\alpha^2}{\beta} + \beta \right).$$

From the boundary condition at $z = 1$ it follows that

$$D_1 e^{\alpha} \left( \frac{1}{\text{Da}} \cos(\beta) - \sin(\beta) \left( \frac{1}{\text{Da}} \frac{\alpha}{\beta} + \frac{1}{\text{Th}_B} \left( \frac{\alpha^2}{\beta} + \beta \right) \right) \right) = 0,$$

which solutions are $D_1 = 0$ or

$$\tan(\beta) = \frac{\frac{1}{\text{Da}}}{\frac{\alpha}{\beta} \frac{1}{\text{Da}} + \left(\frac{\alpha^2}{\beta} + \beta\right) \frac{1}{\text{Th}_B}} = \frac{\beta}{\frac{\text{Da}}{\text{Th}_B} \beta^2 + \frac{\alpha}{2}} = \frac{2\alpha\beta}{-\beta^2 + \alpha^2}. \quad (3.11)$$

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 e^{-\frac{\text{Th}_B}{2\text{Da}} z} \left( \cos(\beta_n z) + \frac{\text{Th}_B}{2\text{Da}\beta_n} \sin(\beta_n z) \right),$$

where $\beta_n \in (0, +\infty)$ satisfies equation (3.11).

**Case 3:** $\Delta > 0 \Leftrightarrow \lambda > \eta + \mu(1) - \text{Th}_B(\frac{1}{2\text{Da}})^2$. In this case, we have two different real roots $\rho_{1,2} = \alpha \pm \beta$, with $\alpha = -\frac{\text{Th}_B}{2\text{Da}}$, $\beta = \frac{\text{Th}_B}{2} \sqrt{\left(\frac{1}{\text{Da}}\right)^2 + \frac{4(\lambda - \mu(1) - \eta)}{\text{Th}_B}}$, and the solution of equation (3.9) is of the form

$$Z(z) = D_1 e^{(\alpha+\beta)z} + D_2 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 e^{(\alpha+\beta)z} + (\alpha - \beta)D_2 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 \left( e^{(\alpha+\beta)z} - \frac{(\alpha+\beta)}{(\alpha-\beta)} e^{(\alpha-\beta)z} \right), \quad Z'(z) = D_1(\alpha+\beta) \left( e^{(\alpha+\beta)z} - e^{(\alpha-\beta)z} \right).$$

From the boundary condition at $z = 1$, it follows that

$$D_1 e^{\alpha} \frac{(\alpha + \beta)}{\text{Th}_B} (e^{\beta} - e^{-\beta}) + D_1 e^{\alpha} \frac{1}{\text{Da}} (e^{\beta} - \frac{(\alpha + \beta)}{(\alpha - \beta)} e^{-\beta}) = 0,$$
which implies $D_1 = 0$ or

$$ 
\begin{aligned}
e^{\beta} \left( \frac{(\alpha + \beta)}{\mathrm{Th}_B} + \frac{1}{\mathrm{Da}} \right) &= e^{-\beta} \left( \frac{(\alpha + \beta)}{\mathrm{Th}_B} + \frac{(\alpha + \beta)}{(\alpha - \beta)} \frac{1}{\mathrm{Da}} \right) \\
\Leftrightarrow e^{2\beta} &= \frac{\frac{(\alpha+\beta)}{\mathrm{Th}_B} + \frac{1}{\mathrm{Da}} \frac{(\alpha+\beta)}{(\alpha-\beta)}}{\frac{(\alpha+\beta)}{\mathrm{Th}_B} + \frac{1}{\mathrm{Da}}} = \frac{(\alpha + \beta)}{(\alpha - \beta)} \left( \frac{-(\beta + \alpha)\mathrm{Da}}{(\beta - \alpha)\mathrm{Da}} \right) \\
\Leftrightarrow e^{2\beta} &= \left( \frac{\alpha + \beta}{\alpha - \beta} \right)^2.
\end{aligned} 
\quad (3.12) $$

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$.

**Step 1.4.** General expression of $\delta B(r, z, t)$. Given $\eta_n \in E$ (see equation (3.8)), there exists a countable set of admissible eigenvalues $\lambda$

$$ \Lambda_n = \{\lambda_{nm}\}_{m \in \mathbb{N}} = \{\mu(1) + \eta_n - \frac{\mathrm{Th}_B}{(2\mathrm{Da})^2} - \frac{\beta_m^2}{\mathrm{Th}_B}\}_{m \in \mathbb{N}}, \quad (3.13) $$

where $\beta_m$ satisfies system (3.11). Consequently,

$$ \delta B(r, z, t) = \sum_{n \in \mathbb{N} \cup \{0\}} \sum_{m \in \mathbb{N}} A_{nm} e^{\lambda_{nm} t} J_0 \left( \frac{\sqrt{-\mathrm{Th}_B \eta_n}}{\sigma} r \right) \\ \times e^{-\frac{\mathrm{Th}_B}{2\mathrm{Da}} z} \left( \cos(\beta_m z) + \frac{\mathrm{Th}_B}{2\mathrm{Da} \beta_m} \sin(\beta_m z) \right), $$

where $\eta_n \in E$, $\beta_m$ satisfies (3.11), $\lambda_{nm} \in \Lambda_n$ and the constants $A_{nm}$ are chosen such that $\delta B(r, z, 0) = \delta B_{\mathrm{init}}(r, z)$. Notice that the constants $A_{nm}$ are well defined since the two systems (3.6) and (3.9) are regular Sturm-Liouville eigenvalue problems (see, e.g. [23, Theorem 1.3]).

Using Parseval's equation (see, for instance, [45]) 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 e^{2\lambda_{nm}t}. $$

Furthermore, it is straightforward to see that

$$ \lambda_{nm} \le \lambda_{01} = \mu(1) - \frac{\mathrm{Th}_B}{(2\mathrm{Da})^2} - \frac{\beta_1^2}{\mathrm{Th}_B} \quad \forall (n, m) \in (\{0\} \cup \mathbb{N}) \times \mathbb{N}. $$

Therefore, if

$$ \lambda_{01} = \mu(1) - \frac{\mathrm{Th}_B}{(2\mathrm{Da})^2} - \frac{\beta_1^2}{\mathrm{Th}_B} < 0, \quad (3.14) $$

(which is the same condition as (3.4)) it follows that

$$ \|\delta B(t)\|_{L^2(\Omega)}^2 \le e^{2\lambda_{01}t} \sum_{n \in \mathbb{N} \cup \{0\}} \sum_{m \in \mathbb{N}} A_{nm}^2 = e^{2\lambda_{01}t} \|\delta B(0)\|_{L^2(\Omega)}^2 \xrightarrow{t \to \infty} 0. $$

Note that, if $\lambda_{01} < 0$, one can also deduce inequality (that will be used at the end of this proof)

$$ \|\delta B(t)\|_{L^2(\Omega)}^2 \le \|\delta B(0)\|_{L^2(\Omega)}^2 \le K^2 \|\delta B_{\mathrm{init}}\|_{L^\infty(\Omega)}^2, \quad (3.15) $$

where $K$ is a constant relating the norms $\| \cdot \|_{L^2(\Omega)}$ and $\| \cdot \|_{L^\infty(\Omega)}$.
**Step 2.** Let us prove that $\|\delta S(t)\|_{L^2(\Omega)} \to 0$ as $t \to \infty$. Regarding $\delta S$, the main equation involving the substrate in system (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 (3.5) by $r\delta S$ and integrate as follows

$$ 
\begin{aligned}
& \int_0^t \int_\Omega r \frac{d\delta S}{d\tau} \delta S dr dz d\tau \\
&= \frac{1}{\text{Da}} \int_0^t \int_\Omega r \frac{d\delta S}{dz} \delta S dr dz d\tau + \frac{1}{\text{Th}_S} \int_0^t \int_\Omega r \frac{d^2\delta S}{dz^2} \delta S dr dz d\tau \\
&\quad + \frac{\sigma^2}{\text{Th}_S} \int_0^t \int_\Omega \frac{d}{dr} \left( r \frac{d\delta S}{dr} \right) \delta S dr dz d\tau - \mu(1) \int_0^t \int_\Omega r \delta B \delta S dr dz d\tau \\
&= -\frac{1}{\text{Da}} \int_0^t \int_\Omega r \frac{d\delta S}{dz} \delta S dr dz d\tau - \frac{\sigma^2}{\text{Th}_S} \int_0^t \int_\Omega r \left( \frac{d\delta S}{dr} \right)^2 dr dz d\tau \\
&\quad - \frac{1}{\text{Th}_S} \int_0^t \int_\Omega r \left( \frac{d\delta S}{dz} \right)^2 dr dz d\tau - \mu(1) \int_0^t \int_\Omega r \delta B \delta S dr dz d\tau \\
&\quad + \frac{\sigma^2}{\text{Th}_S} \int_0^t \int_{\Gamma_{\text{sym} \cup \Gamma_{\text{wall}}}} r \frac{d\delta S}{dr} \delta S dz d\tau + \int_0^t \int_{\Gamma_{\text{in}}} r \left( \frac{1}{\text{Th}_S} \frac{d\delta S}{dz} + \frac{1}{\text{Da}} \delta S \right) \delta S dr d\tau \\
&\quad - \int_0^t \int_{\Gamma_{\text{out}}} r \left( \frac{1}{\text{Th}_S} \frac{d\delta S}{dz} + \frac{1}{\text{Da}} \delta S \right) \delta S dr d\tau.
\end{aligned}
$$

However, applying the boundary conditions for $\delta S$ in system (3.5) one has

$$ 
\begin{aligned}
& \int_0^t \int_\Omega r \frac{d\delta S}{d\tau} \delta S dr dz d\tau \\
&= -\underbrace{\frac{1}{\text{Da}} \int_0^t \int_\Omega r \frac{d\delta S}{dz} \delta S dr dz d\tau}_{(I)} - \frac{1}{\text{Da}} \int_0^t \int_{\Gamma_{\text{out}}} r \delta S^2 dr dz d\tau \\
&\quad - \frac{\sigma^2}{\text{Th}_S} \int_0^t \int_\Omega r \left( \frac{d\delta S}{dr} \right)^2 dr dz d\tau - \frac{1}{\text{Th}_S} \int_0^t \int_\Omega r \left( \frac{d\delta S}{dz} \right)^2 dr dz d\tau \\
&\quad - \mu(1) \int_0^t \int_\Omega r \delta S \delta B dr dz d\tau.
\end{aligned}
\quad (3.16) $$

The integral denoted by (I) in equation (3.16) can be rewritten as

$$ (I) = -\frac{1}{2\text{Da}} \int_0^t \int_{\Gamma_{\text{in}}} r \delta S^2 dr d\tau + \frac{1}{2\text{Da}} \int_0^t \int_{\Gamma_{\text{out}}} r \delta S^2 dr d\tau. $$
Thus, equation (3.16) leads to

$$ 
\begin{aligned}
& \frac{1}{2} \int_0^t \int_{\Omega} r \frac{d(\delta S^2)}{d\tau} dr dz d\tau + \frac{\sigma^2}{\text{Th}_S} \int_0^t \int_{\Omega} r \left(\frac{d\delta S}{dr}\right)^2 dr dz d\tau \\
& \qquad + \frac{1}{\text{Th}_S} \int_0^t \int_{\Omega} r \left(\frac{d\delta S}{dz}\right)^2 dr dz d\tau + \frac{1}{2\text{Da}} \int_0^t \int_{\Gamma_{\text{in}}} r \delta S^2 dr d\tau \\
& \frac{1}{2\text{Da}} \int_0^t \int_{\Gamma_{\text{out}}} r \delta S^2 dr d\tau \\
& = -\mu(1) \int_0^t \int_{\Omega} r \delta S \delta B dr dz d\tau.
\end{aligned} 
\quad (3.17) $$

By multiplying (3.17) by $2\pi$ and applying Young's inequality (with $\epsilon > 0$ to be chosen afterward), we obtain

$$ 
\begin{aligned}
& \frac{1}{2} \int_0^t \frac{d}{d\tau} (\|\delta S(\tau)\|_{L^2(\Omega^*)}^2) d\tau + \frac{\min(1, \sigma^2)}{\text{Th}_S} \int_0^t \|\nabla \delta S(\tau)\|_{L^2(\Omega^*)}^2 d\tau \\
& \qquad + \frac{1}{2\text{Da}} \int_0^t \|\delta S(\tau)\|_{L^2(\Gamma_{\text{out}}^*)}^2 d\tau \\
& \le \epsilon\mu(1) \int_0^t \|\delta S(\tau)\|_{L^2(\Omega^*)}^2 d\tau + \frac{\mu(1)}{4\epsilon} \int_0^t \|\delta B(\tau)\|_{L^2(\Omega^*)}^2 d\tau.
\end{aligned} 
$$

Considering $A = \min\{\frac{1}{\text{Th}_S}, \frac{\sigma^2}{\text{Th}_S}, \frac{1}{2\text{Da}}\}$, it follows that

$$ 
\begin{aligned}
& \frac{1}{2} \int_0^t \frac{d}{d\tau} (\|\delta S(\tau)\|_{L^2(\Omega^*)}^2) d\tau + A \int_0^t (\|\nabla \delta S(\tau)\|_{L^2(\Omega^*)}^2 + \|\delta S(\tau)\|_{L^2(\Gamma_{\text{out}}^*)}^2) d\tau \\
& \le \epsilon\mu(1) \int_0^t \|\delta S(\tau)\|_{L^2(\Omega^*)}^2 d\tau + \frac{\mu(1)}{4\epsilon} \|\delta B\|_{L^2(\Omega^* \times (0,t))}^2.
\end{aligned} 
\quad (3.18) $$

Now, applying Friedrich's inequality (see e.g., [25, Theorem 6.1]) to inequality (3.18) with $E = \Gamma_{\text{out}}^*$, there exists a constant $C$ depending on $\Omega^*$ and $\Gamma_{\text{out}}^*$ such that

$$ 
\begin{aligned}
& \frac{1}{2} \int_0^t \frac{d}{d\tau} (\|\delta S(\tau)\|_{L^2(\Omega^*)}^2) d\tau \le (\epsilon\mu(1) - \frac{A}{C}) \int_0^t \|\delta S(\tau)\|_{L^2(\Omega^*)}^2 d\tau \\
& \qquad + \frac{\mu(1)}{4\epsilon} \|\delta B\|_{L^2(\Omega^* \times (0,t))}^2.
\end{aligned} 
\quad (3.19) $$

Next, applying the Gronwall's inequality in its integral form, it follows that

$$ 
\|\delta S(t)\|_{L^2(\Omega^*)}^2 \le \underbrace{(\|\delta S_{\text{init}}\|_{L^2(\Omega^*)}^2 + \frac{\mu(1)}{2\epsilon} \|\delta B\|_{L^2(\Omega^* \times (0,t))}^2)}_{(:=m(t))} \exp \underbrace{(2(\epsilon\mu(1) - \frac{A}{C})t)}_{(:=\alpha)}.
$$

Since $\|\delta B(t)\|_{L^2(\Omega)}^2 \le K^2 \|\delta B_{\text{init}}\|_{L^\infty(\Omega)}^2$ for all $t > 0$ (see equation (3.15)), taking $\epsilon < \frac{A}{\mu(1)C}$ it follows that $\alpha < 0$. Thus,

$$ 
\|\delta S(t)\|_{L^2(\Omega^*)}^2 \le (\|\delta S_{\text{init}}\|_{L^2(\Omega^*)}^2 + \frac{\mu(1)}{2\epsilon} t K^2 \|\delta B_{\text{init}}\|_{L^\infty(\Omega^*)}^2) e^{\alpha t} \xrightarrow{t \to \infty} 0.
$$

□

**Remark 3.7.** To the best of our knowledge, if the inequality “<” is replaced by the equality “=” in condition (3.4), the stability analysis of the steady state (1,0)
requires different techniques from those used in the proof of Theorem 3.5. This case has not been tackled here.

Taking into account Theorem 3.5, we conjecture (supported by the numerical experiments presented in Section 4) that the following result holds:

**Conjecture 3.8.** If $\mu$ satisfies (A1) (respectively, (A2)), a sufficient condition for $(S_2^*, B_2^*)$ to be a linearly asymptotically stable steady state of system (2.4) is that

$$ \mu(1) > \frac{\mathrm{Th}_B}{(2\mathrm{Da})^2} + \frac{(\beta_1(\mathrm{Da}, \mathrm{Th}_B))^2}{\mathrm{Th}_B} \quad (3.20) $$

(respectively,

$$ 1 > \frac{\mathrm{Th}_B}{(2\mathrm{Da})^2} + \frac{(\beta_1(\mathrm{Da}, \mathrm{Th}_B))^2}{\mathrm{Th}_B}. \quad (3.21) $$

**Remark 3.9.** In terms of the variables with dimensions appearing in system (2.1), conditions (3.20) and (3.21) are reformulated, respectively, as

$$ \mu(S_e) > \frac{u^2}{4D_B} + \frac{D_B}{H^2}(\tilde{\beta}_1(H, u, D_B))^2, $$

and

$$ \bar{\mu} > \frac{u^2}{4D_B} + \frac{D_B}{H^2}(\tilde{\beta}_1(H, u, D_B))^2. $$

**Remark 3.10.** From Theorem 3.5 and Conjecture 3.8, it follows that if $\mu$ satisfies (A2) and $\mu(1) < 1$, there is bistability in system (2.4) when

$$ \mu(1) < \frac{\mathrm{Th}_B}{(2\mathrm{Da})^2} + \frac{(\beta_1(\mathrm{Da}, \mathrm{Th}_B))^2}{\mathrm{Th}_B} < 1. $$

3.2.1. *Bounds for the flow rate assuring asymptotic stability of the steady states.* Conditions (3.4) and (3.21) include in their analytical expression the model parameters $Da$, $Th_B$ and $\mu(1)$, among which the flow rate $Da$ can be seen as a bioreactor control parameter. In this section, we present bounds for the parameter $Da$ assuring the asymptotic stability of the steady states $(1, 0)$ and $(S_2^*, B_2^*)$. To do so, we first define the following function.

**Definition 3.11.** For a fixed value $Th_B$, we define the function $f_{Th_B} : [0, +\infty) \to [0, +\infty)$

$$ f_{Th_B}(\mathrm{Da}) = \frac{\mathrm{Th}_B}{(2\mathrm{Da})^2} + \frac{(\beta_1(\mathrm{Da}, \mathrm{Th}_B))^2}{\mathrm{Th}_B}. $$

In Figure 2 we plot the value of functions $\beta_1(\mathrm{Da}, \mathrm{Th}_B)$ and $f_{Th_B}(\mathrm{Da})$ for $Th_B \in \{\frac{1}{5}, 1, 5\}$ and $Da \in [0, 2]$. For a fixed value $Th_B$, function $\beta_1(\cdot, \mathrm{Th}_B)$ is decreasing, bounded by $\pi$ (see the proof of Theorem 3.5 for a detailed explanation of this feature) and $\beta_1(\mathrm{Da}, \mathrm{Th}_B) \xrightarrow{\mathrm{Da} \to +\infty} 0$. One can also conclude that, for a fixed value $Th_B$, function $f_{Th_B}$ is decreasing, $f_{Th_B}(\mathrm{Da}) \xrightarrow{\mathrm{Da} \to 0} +\infty$ and $f_{Th_B}(\mathrm{Da}) \xrightarrow{\mathrm{Da} \to +\infty} 0$. Taking into account these properties of $f_{Th_B}$, we define the following variables.

**Definition 3.12.** We define

* $\mathrm{Da}^{\mathrm{W}}_{(2.4)}(\mathrm{Th}_B, \mu(1)) := (f_{\mathrm{Th}_B})^{-1}(\mu(1)).$

* $\mathrm{Da}^{\mathrm{NW}}_{(A1),(2.4)}(\mathrm{Th}_B, \mu(1)) := \mathrm{Da}^{\mathrm{W}}_{(2.4)}(\mathrm{Th}_B, \mu(1)).$

* $\mathrm{Da}^{\mathrm{NW}}_{(A2),(2.4)}(\mathrm{Th}_B) := (f_{\mathrm{Th}_B})^{-1}(1).$
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Graphical plots of functions $\beta_1(\text{Da}, \text{Th}_B)$ and $f_{\text{Th}_B}(\text{Da})$ (described in Definitions 3.4 and 3.11, respectively) for $\text{Th}_B \in \{\frac{1}{5}, 1, 5\}$ and $\text{Da} \in [0, 2]$.</figcaption></figure>

**Remark 3.13.** Following Theorem 3.5, if Conjecture 3.8 is true, it follows that

* If $\text{Da} < \text{Da}^{\text{W}}_{(2.4)}(\text{Th}_B, \mu(1))$, then the equilibrium state (1, 0) of system (2.4) is linearly asymptotically stable.

* If $\mu$ satisfies (A1) and $\text{Da} > \text{Da}^{\text{NW}}_{(A1),(2.4)}(\text{Th}_B, \mu(1))$, then the equilibrium state $(S_2^*, B_2^*)$ of system (2.4) is linearly asymptotically stable.

* If $\mu$ satisfies (A2) and $\text{Da} > \text{Da}^{\text{NW}}_{(A2),(2.4)}(\text{Th}_B)$, then the equilibrium state $(S_2^*, B_2^*)$ of system (2.4) is linearly asymptotically stable.

## 4. NUMERICAL EXPERIMENTS

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 3. In Section 4.1, we study the sensitivity of variables $\text{Da}^{\text{W}}_{(2.4)}(\text{Th}_B, \mu(1))$ and $\text{Da}^{\text{NW}}_{(A2),(2.4)}(\text{Th}_B)$, defined in Section 3.2.1, regarding the model parameters. Then, In Section 4.2, we carry out the numerical implementation of system (2.4)-(2.5) in order to check the interest of these functions. Finally, in Section 4.3, we compare the results of the stability analysis of systems (1.2) and (2.4).

In this section, the value of functions $\text{Da}^{\text{W}}_{(2.4)}(\text{Th}_B, \mu(1))$ and $\text{Da}^{\text{NW}}_{(A2),(2.4)}(\text{Th}_B)$ is approximated numerically using a self-implemented *Dichotomy method* (see, e.g. [21]). Moreover, for each pair $(\text{Th}_B, \text{Da})$, the value of $\beta_1(\text{Th}_B, \text{Da})$ (see Definition 3.4) was computed by using the MATLAB function (see www.mathworks.com/help/symbolic/vpasolve.html).

### 4.1. Sensitivity to model parameters.
In this section, we perform the sensitivity analysis of $\text{Da}^{\text{W}}_{(2.4)}(\text{Th}_B, \mu(1))$ with respect to the nondimensional parameters $\text{Th}_B$ and $\mu(1)$ (the sensitivity analysis of $\text{Da}^{\text{NW}}_{(A2),(2.4)}(\text{Th}_B)$ can be obtained with a similar methodology).
4.1.1. Sensitivity with respect to $\mu(1)Lg$. Taking into account that

$$Da_{(2.4)}^{W}(Th_{B}, \mu(1)) = (f_{Th_{B}})^{-1}(\mu(1))$$

and $f_{Th_B}$ is decreasing, one concludes that, for any fixed value $Th_B$, the function $Da_{(2.4)}^W(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/Da) suitable to avoid washout also increases (see, e.g, [8, 15, 44]).

4.1.2. Sensitivity with respect to $Th_B Lg$. To easily analyze the of $Da_{(2.4)}^W(Th_B, \mu(1))$ with respect to $Th_B$, we aim to approximate $Da_{(2.4)}^W(Th_B, \mu(1))$ by using the following variables:

* $\overline{Da}_{(2.4)}^W(Th_B, \mu(1)) := \frac{1}{2}\sqrt{\frac{Th_B}{\mu(1)}}$. This should be a good approximation of $Da_{(2.4)}^W(Th_B, \mu(1))$ assuming that the second term of the right hand side of condition (3.4) is negligible.

* $\overline{Da}_{(2.4)}^W(Th_B, \mu(1)) := (g_{Th_B})^{-1}(\mu(1))$, where $g_{Th_B} : [0, +\infty) \to [0, +\infty)$,

$$g_{Th_B}(Da) = \frac{(\beta_1(Th_B, Da))^2}{Th_B}.$$

This should be a good approximation of $Da_{(2.4)}^W(Th_B, \mu(1))$ assuming that the first term of the right hand side of condition (3.4) is negligible. Since $\beta_1(Th_B, Da) < \pi$ (see the proof of Theorem 3.5 for a detailed explanation of this fact), if $Th_B\mu(1) > \pi^2$, then the function $\overline{Da}_{(2.4)}^W(Th_B, \mu(1))$ is not defined. We approximate numerically $\overline{Da}_{(2.4)}^W(Th_B, \mu(1))$ applying the same methodology that the one used to approximate numerically $Da_{(2.4)}^W(Th_B, \mu(1))$, described above.

Figure 3 illustrates the difference between $Da_{(2.4)}^W(Th_B, \mu(1))$, $\overline{Da}_{(2.4)}^W(Th_B, \mu(1))$ and $\overline{Da}_{(2.4)}^W(Th_B, \mu(1))$ when $\mu(1) = 0.5$ and $Th_B \in [5 \cdot 10^{-3}, 5 \cdot 10^3]$. We observe that $\overline{Da}_{(2.4)}^W(Th_B, 0.5)$ approximates $Da_{(2.4)}^W(Th_B, 0.5)$ for values smaller than $\log(Th_B) = -2 (Th_B \approx 0.1)$ while $\overline{Da}_{(2.4)}^W(Th_B, 0.5)$ approximates $Da_{(2.4)}^W(Th_B, 0.5)$ for values larger than $\log(Th_B) = 6 (Th_B \approx 400)$.

The comparison between the functions $Da_{(2.4)}^W(Th_B, \mu(1))$, $\overline{Da}_{(2.4)}^W(Th_B, \mu(1))$ and $\overline{Da}_{(1.2)}^W(Th_B, \mu(1))$, shown in Figure 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 $Th_B \ge 10^4$, the function $\overline{Da}_{(2.4)}^W(Th_B, \mu(1))$ can be used as an approximation of $Da_{(2.4)}^W(Th_B, \mu(1))$; and if $Th_B \le 0.1$, the function $\overline{Da}_{(2.4)}^W(\mu(1))$ can be used as an approximation of $Da_{(2.4)}^W(Th_B, \mu(1))$.

Taking into account the approximations of $Da_{(2.4)}^W(Th_B, \mu(1))$ presented above and Figure 3, the sensitivity of $Da_{(2.4)}^W(Th_B, \mu(1))$ with respect to $Th_B$ reads as follows:

* If $Th_B \le 0.1$, the variable $Da_{(2.4)}^W(Th_B, \mu(1))$ is not sensible to parameter $Th_B$. Indeed, small values of $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
<figure><img src="image_5.png" /><figcaption>FIGURE 3. Comparison between the functions $Da_{(2.4)}^{W}(\text{Th}_B, 0.5)$, $\overline{Da}_{(2.4)}^{W}(\text{Th}_B, 0.5)$ and $\widehat{Da}_{(2.4)}^{W}(\text{Th}_B, 0.5)$ (depicted with solid, dotted and dashed lines, respectively), when $\text{Th}_B \in [5 \cdot 10^{-3}, 5 \cdot 10^3]$.</figcaption></figure>

will see in Section 4.3, if $\text{Th}_B \le 0.1$, the dynamics of the bioreactor can be modeled with ordinary differential equations.

• If $\text{Th}_B > 0.1$, the variable $Da_{(2.4)}^W(\text{Th}_B, \mu(1))$ seems to increase with parameter $\text{Th}_B$. This outcome is physically reasonable, since as parameter $\text{Th}_B$ increases (equivalently, the diffusion coefficient decreases) the flow rate (1/Da) should be chosen smaller to favor the reaction between the substrate and the biomass (see [8]).

• If $\text{Th}_B \ge 10^4$, the variable $Da_{(2.4)}^W(\text{Th}_B, \mu(1))$ is quadratically proportional to $\text{Th}_B$.

## 4.2. Numerical validation of the results.
In this section, we check the properties given in Remark 3.13 for the threshold values $Da_{(2.4)}^W(\text{Th}_B, \mu(1))$ and $Da_{(A2),(2.4)}^{\text{NW}}(\text{Th}_B, \mu(1))$ by using the numerical solution of system (2.4)-(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 [36, 38]). 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 $Da_{(2.4)}^W(\text{Th}_B, \mu(1))$ and $Da_{(A2),(2.4)}^{\text{NW}}(\text{Th}_B, \mu(1))$, we define the following variables:

• $\widetilde{Da}_{(2.4)}^W(\text{Th}_B, \mu(1)) := \sup\{\text{Da: the numerical solution of system (2.4)-(2.5)}$ (with parameters $\text{Th}_B$, $\text{Da}$, $\text{Th}_S = \text{Th}_B$, $\sigma = 1$, $\mu$ the nondimensional Monod function with $K_S = \frac{1-\mu(1)}{\mu(1)}$, $S_{\text{init}} = 0.1$ and $B_{\text{init}} = 0.9$) approaches asymptotically the steady state $(1, 0)\}$.
• $\widetilde{\mathcal{D}}_{(A2),(2.4)}^{\text{NW}}(\text{Th}_B) := \inf\{\text{Da}: \text{the numerical solution of system (2.4)-(2.5) (with parameters Th}_B, \text{Da, Th}_S = \text{Th}_B, \sigma = 1, \mu \text{ the nondimensional Haldane function with } \mu^*/\bar{\mu} = 1.7071, K_S = 0.3536 \text{ and } K_I = 2.8284) \text{ approaches asymptotically a steady state different from } (1,0)\}$.

We approximate numerically $\widetilde{\mathcal{D}}_{(2.4)}^{\text{W}}(\text{Th}_B, \mu(1))$ and $\widetilde{\mathcal{D}}_{(A2),(2.4)}^{\text{NW}}(\text{Th}_B)$ by using again a self-implemented *Dichotomy method*. Figure 4(a) illustrates the difference between $\widetilde{\mathcal{D}}_{(2.4)}^{\text{W}}(\text{Th}_B, \mu(1))$ and $\widetilde{\mathcal{D}}_{(2.4)}^{\text{NW}}(\text{Th}_B, \mu(1))$ when $\text{Th}_B \in [5 \cdot 10^{-3}, 1.5 \cdot 10^2]$ and $\mu(1) = 0.5$.

Similarly, Figure 4-(b) shows the difference between $\mathcal{D}_{(A2),(2.4)}^{\text{NW}}(\text{Th}_B)$ and $\widetilde{\mathcal{D}}_{(A2),(2.4)}^{\text{NW}}(\text{Th}_B)$ when $\text{Th}_B \in [5 \cdot 10^{-3}, 1.5 \cdot 10^2]$. We point out that these comparisons were also performed with $\widetilde{\mathcal{D}}_{(2.4)}^{\text{W}}(\text{Th}_B, \mu(1))$ and $\widetilde{\mathcal{D}}_{(A2),(2.4)}^{\text{NW}}(\text{Th}_B)$ defined using other model parameters $\sigma$, $\text{Th}_S$ and $\mu$ and similar results were obtained.

<figure><img src="image_6.png" /><figcaption>(a) Comparing $\log(\mathcal{D}_{(2.4)}^{\text{W}}(\text{Th}_B, 0.5))$ (solid) and $\log(\widetilde{\mathcal{D}}_{(A2),(2.4)}^{\text{NW}}(\text{Th}_B))$ (dashed) line) and $\log(\widetilde{\mathcal{D}}_{(2.4)}^{\text{W}}(\text{Th}_B, 0.5))$ (dashed lines). (b) Comparing $\log(\mathcal{D}_{(A2),(2.4)}^{\text{NW}}(\text{Th}_B))$ (solid line) and $\log(\widetilde{\mathcal{D}}_{(A2),(2.4)}^{\text{NW}}(\text{Th}_B))$ (dashed line).</figcaption></figure>

<caption>FIGURE 4. Numerical validation of results.</caption>

In Figure 5, we plot the steady-state solution $(S_2^*, B_2^*)$ of system (2.4), computed numerically when $\text{Th}_B = \text{Th}_S = e^4$, $\text{Da} = e^2$, $\sigma = 1$, $S_{\text{init}} = 0.1$, $B_{\text{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(\text{Th}_B) = 4$ and $\log(\text{Da}) = 2$, the equilibrium solution $(S_2^*, B_2^*)$ is linearly asymptotically stable (see Figure 4-(a)). Notice that the same steady-state solution can be obtained with nonhomogeneous initial conditions (for instance, $S_{\text{init}}(r, z) = rz$ and $B_{\text{init}}(r, z) = r(1-z)$).

The bistability of system (2.4), stated in Remark 3.10, is perceivable when numerically solving system (2.4). For instance, if $\text{Th}_B = 0.01$, $\text{Da} = 1.5$ and $\mu(1) = 0.5$, we observe that the solution of system (2.4) (computed with parameters $\sigma = 1$, $\text{Th}_S = 0.01$ and $\mu$ the nondimensional Haldane function with $\mu^*/\bar{\mu} = 1.7071$, $K_S = 0.0529$ and $K_I = 0.4235$) approaches $(1,0)$ if we choose $S_{\text{init}} = 0.9$ and $B_{\text{init}} = 0.1$, while it approaches a different equilibrium (similar to the one represented in Figure 5) solution if we set $S_{\text{init}} = 0.1$ and $B_{\text{init}} = 0.9$.

### 4.3. Comparison with the stability analysis of system (1.2)Lg.
In this section, we compare the stability analysis conditions associated to the ODE and PDE
<figure><img src="image_7.png" /><figcaption>FIGURE 5. Representation of the steady-state solution $(S_2^*, B_2^*)$ of (2.4) computed numerically when $Th_B = Th_S = e^4$, $Da = e^2$, $\sigma = 1$, $S_{\text{init}} = 0.1$, $B_{\text{init}} = 0.9$ and $\mu$ being the nondimensional Monod Function with $K_S = 1$ (so that $\mu(1) = 0.5$).</figcaption></figure>

systems (1.2) and (2.4), respectively. As done in Section 3.2.1 for system (2.4), we define the following variables:

**Definition 4.1.**

* $Da_{(1,2)}^W(\mu(1)) := 1/\mu(1).$

* $Da_{(A1),(1,2)}^{NW}(\mu(1)) := Da_{(1,2)}^W(\mu(1)).$

* $Da_{(A2),(1,2)}^{NW} := 1.$

**Remark 4.2.** According to Remark 2.2 and Definition 4.1, the stability analysis of system (1.2) (shown in Section 1) can be rewritten as

* If $Da < Da_{(1,2)}^W(\mu(1))$, then the equilibrium solution $(1,0)$ of system (1.2) is asymptotically stable.

* If $\mu$ satisfies (A1) and $Da > Da_{(A1),(1,2)}^{NW}(\mu(1))$, then the equilibrium solution $(S_2^*, B_2^*)$ of system (1.2) is asymptotically stable.

* If $\mu$ satisfies (A2) and $Da > Da_{(A2),(1,2)}^{NW}(1)$, then the equilibrium solution $(S_2^*, B_2^*)$ of system (1.2) is asymptotically stable.

Figure 6 illustrates the difference between the variable $Da_{(A2),(2,4)}^{NW}(Th_B)$ and the constant $Da_{(A2),(1,2)}^{NW} = 1$ (and the difference, when $\mu(1) = 0.5$, between the variable $Da_{(2,4)}^W(Th_B, \mu(1))$ and the constant $Da_{(1,2)}^W(\mu(1)) = 2$). In both cases $Th_B \in [5 \cdot 10^{-3}, 1.5 \cdot 10^2]$. Notice that the area limited between the curves $Da_{(A2),(2,4)}^{NW}(Th_B)$ and $Da_{(2,4)}^W(Th_B, 0.5)$ is the region of bistability of system (2.4) (see Remark 3.10).

We observe that $\log(Da_{(A2),(2,4)}^{NW}(Th_B)) \approx 0$ for values smaller than $\log(Th_B) \approx -2$ ($Th_B \approx 0.1$). Similarly, for the particular case when $\mu(1) = 0.5$, we observe that $\log(Da_{(2,4)}^W(Th_B, 0.5)) \approx \log(2)$ also for values smaller than $\log(Th_B) \approx -2$ ($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 $Th_B < 0.1$, the
<figure><img src="image_8.png" /><figcaption>FIGURE 6. Comparisong between $\log(Da_{(2.4)}^W(Th_B, 0.5))$, $\log(Da_{(A2),(2.4)}^{NW}(Th_B))$ (solid lines) and constant values $\log(Da_{(1.2)}^W(0.5)) = \log(2)$, $\log(Da_{(A2),(1.2)}^{NW}) = 0$ (dashed line) when $Th_B \in [5 \times 10^{-3}, 1.5 \times 10^2]$.</figcaption></figure>

stability results obtained for the ODE and PDE systems (1.2) and (2.4) are similar. This result is consistent with the physics of the problem. Indeed, small values of $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 [8]).

## 5. 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 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 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 Da (Damkhöler Number,
relating reaction and advective rates), $Th_B$ (Thiéle 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éle Modulus ($Th_B$) lower than 0.1, but not valid for values of Thiéle Modulus above this value.

**Remark 5.1.** 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 [10] and published in the same issue of this article.

**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ía" 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 MIS-TEA, Montpellier, France.

## REFERENCES

[1] J. F. Andrews; A mathematical model for the continuous culture of microorganisms utilizing inhibitory substrates. *Biotechnology and Bioengineering*, 10(6):707-723, 1968.

[2] R. Aris; Phenomena of multiplicity, stability and symmetry. *Annals of the New York Academy of Sciences*, 231(1):86-98, 1974.

[3] R. Aris; *The Mathematical Theory of Diffusion and Reaction in Permeable Catalysts: Vol. 1: The theory of the steady state*. Oxford University Press, USA, 1975.

[4] N. H. Asmar; *Partial Differential Equations with Fourier Series and Boundary Value Problems*. Pearson Prentice Hall, 2005.

[5] J. E. Bailey and D. F. Ollis; *Biochemical Engineering Fundamentals*. Chemical engineering. McGraw-Hill, 1986.

[6] A. W. Bush, A. E. Cool; The effect of time delay and growth rate inhibition in the bacterial treatment of wastewater. *Journal of Theoretical Biology*, 63(2):385-395, 1976.

[7] M. Crespo, B. Ivorra, A.M. Ramos; Existence and uniqueness of solution of a continuous flow bioreactor model with two species. *Real Academia de Ciencias Exactas, Físicas y Naturales. Serie A: Matemáticas*, 110(2):357-377, 2016.

[8] M. Crespo, A. M. Ramos, B. Ivorra, A. Rapaport; Modeling and optimization of activated sludge bioreactors for wastewater treatment taking into account spatial inhomogeneities. *Journal of Process Control*, 54:118-128, 2017.

[9] P. V. Danckwerts; Continuous flow systems. *Chemical Engineering Science*, 2(1):1-13, 1953.

[10] J. I. Díaz; Remarks on a precedent paper by Crespo, Ivorra and Ramos on the stability of bioreactor processes. This issue of EJDE-2017.

[11] J. I. Díaz, A. C. Casal; On the principle of pseudo-linearized stability: Applications to some delayed nonlinear parabolic equations. *Nonlinear Analysis*, 63(5):e997-e1007, 2005.

[12] D. Dochain; *Automatic Control of Bioprocesses*. ISTE. Wiley, 2010.

[13] D. Dochain and P. Vanrolleghem; *Dynamical Modelling and Estimation in Wastewater Treatment Processes*. IWA Publishing, 2001.

[14] A. K. Dramé; A semilinear parabolic boundary-value problem in bioreactors theory. *Electronic Journal of Differential Equations*, 2004.
[15] A. K. Dramé, C. Lobry, J. Harmand, A. Rapaport, F. Mazenc; Multiple stable equilibrium profiles in tubular bioreactors. *Mathematical and Computer Modelling*, 48(1112):1840–1853, 2008.

[16] A. E. Fischer, J. E. Marsden; Linearization stability of nonlinear partial differential equations. *Proceedings of Symposia in Pure Mathematics*, 27:219–263, 1975.

[17] P. Gajardo, J. Harmand, H.C. Ramírez, A. Rapaport; Minimal time bioremediation of natural water resources. *Automatica*, 47(8):1764–1769, 2011.

[18] J. P. Grivet; Nonlinear population dynamics in the chemostat. *Computing in Science Engineering*, 3, 2001.

[19] J. K. Hale; *Asymptotic Behavior of Dissipative Systems*. Mathematical Surveys and Monographs. American Mathematical Society, 2010.

[20] D. Herbert, R. Elsworth, R. C. Telling; The continuous culture of bacteria; a theoretical and experimental study. *Journal of general microbiology*, 1956.

[21] J. A. Infante, J. M. Rey; *Métodos Numéricos. Teoría, problemas y prácticas con MATLAB*. Ediciones Pirámide. Grupo Anaya, 2015.

[22] S. Ishihara, M. Otsuji, A. Mochizuki; Transient and steady state of mass-conserved reaction-diffusion systems. *Phys. Rev. E*, 75:015203, 2007.

[23] A. W. Knapp; *Advanced Real Analysis*. Cornerstones. Birkhäuser Boston, 2005.

[24] H. J. Kuiper; Invariant sets for nonlinear elliptic and parabolic systems. *Journal on Mathematical Analysis*, 11(6):1075–1103, 1980.

[25] L. Lanzani, Z. W. Shen; On the Robin boundary condition for Laplace's equation in Lipschitz domains. *Communications in Partial Differential Equations*, 29(1 & 2): 91–109, 2004.

[26] W. Liu; *Elementary Feedback Stabilization of the Linear Reaction-Convection-Diffusion Equation and the Wave Equation*. Mathématiques et Applications. Springer Berlin Heidelberg, 2009.

[27] G. B. Marin; Advances in Chemical Engineering: Multiscale analysis. *Advances in Chemical Engineering*, 2005.

[28] R. H. Martin; Asymptotic behavior for semilinear differential equations in Banach spaces. *SIAM Journal on Mathematical Analysis*, 9(6):1105–1119, 1978.

[29] H. Matano; Convergence of solutions of one-dimensional semilinear parabolic equations. *J. Math. Kyoto Univ.*, 18(2):221–227, 1978.

[30] H. Matano; Asymptotic behavior and stability of solutions of semilinear diffusion equations. *Publ. Res. Inst. Math. Sci.*, 15:401–454, 1979.

[31] C. R. McGowin, D. D. Perlmutter; Tubular reactor steady state and stability characteristics. *I. Effect of axial mixing*. *AIChE Journal*, 17(4):831–837, 1971.

[32] J. Moreno; Optimal time control of bioreactors for the wastewater treatment. *Optimal Control Applications and Methods*, 20(3):145–164, 1999.

[33] Y. Morita, T. Ogawa; Stability and bifurcation of nonconstant solutions to a reaction diffusion system with conservation of mass. *Nonlinearity*, 23(6):1387, 2010.

[34] C. V. Pao; Asymptotic stability of reaction-diffusion systems in chemical reactor and combustion theory. *Journal of Mathematical Analysis and Applications*, 82(2):503–526, 1981.

[35] C. V. Pao; *Nonlinear Parabolic and Elliptic Equations*. Fems Symposium. Springer US, 1992.

[36] C. Parés, C. Vázquez, F. Coquel; *Advances in Numerical Simulation in Physics and Engineering: Lecture Notes of the XV 'Jacques-Louis Lions' Spanish-French School*. SEMA SIMAI Springer Series. Springer International Publishing, 2014.

[37] A. J. Perumpanani, J. A. Sherratt, P. K. Maini; Phase differences in reaction-diffusion-advection systems and applications to morphogenesis. *IMA Journal of Applied Mathematics*, 55: 19–33, 1995.

[38] A. M. Ramos; *Introducción al Análisis Matemático del Método de Elementos Finitos*. Editorial Complutense, 2012.

[39] V. S. H. Rao, P. Raja Sekhara Rao; Basic chemostat model revisited. *Differential Equations and Dynamical Systems*, 17(1): 3–16, 2009.

[40] A. Rapaport, I. Haidar, J. Harmand; Global dynamics of the buffered chemostat for a general class of response functions. *Journal of Mathematical Biology*, 71(1): 69–98, 2015.

[41] J. Siebert, S. Alonso, M. Bär, E. Schöll; Dynamics of reaction-diffusion patterns controlled by asymmetric nonlocal coupling as a limiting case of differential advection. *Phys. Rev. E*, 89: 052909, 2014.
[42] H. Smith, P. Waltman; *The theory of the Chemostat*. In *Cambridge studies in Mathematical biology*, volume 13. Cambridge: Cambridge University Press, 1995.

[43] N. A. S. Smith, S. L. Mitchell, A. M. Ramos; Analysis and simplification of a mathematical model for high-pressure food processes. *Applied Mathematics and Computation*, 226: 20–37, 2014.

[44] S. A. Vejtasa, R. A. Schmitz; An experimental study of steady state multiplicity and stability in an adiabatic stirred reactor. *AIChE Journal*, 16(3): 410–419, 1970.

[45] H. F. Weinberger; *A First Course in Partial Differential Equations with Complex Variables and Transform Methods*. Dover Publications, 1995.

[46] J. J. Winkin, D. Dochain, P. Ligarius; Dynamical analysis of distributed parameter tubular reactors. *Automatica*, 36(3): 349–361, 2000.

MARÍA CRESPO (CORRESPONDING AUTHOR)

UMR MISTEA - MATHÉMATIQUES, INFORMATIQUE ET STATISTIQUE POUR L'ÉNVIRONNEMENT ET L'AGRONOMIE, (INRA/SUPAGRO). 2, PLACE P.VIALA, 34060 MONTPELLIER, FRANCE

*E-mail address: maria.crespo-moya@umontpellier.fr*

BENJAMIN IVORRA

DEPARTAMENTO DE MATEMÁTICA APLICADA & INSTITUTO DE MATEMÁTICA INTERISCIPLINAR, UNIVERSIDAD COMPLUTENSE DE MADRID, PLAZA DE CIENCIAS, 3, 28040 MADRID, SPAIN

*E-mail address: ivorra@ucm.es*

ÁNGEL MANUEL RAMOS

DEPARTAMENTO DE MATEMÁTICA APLICADA & INSTITUTO DE MATEMÁTICA INTERISCIPLINAR, UNIVERSIDAD COMPLUTENSE DE MADRID, PLAZA DE CIENCIAS, 3, 28040 MADRID, SPAIN

*E-mail address: angel@mat.ucm.es*