Comparison of performance in an anaerobic digestion process: one-reactor vs two-reactor configurations

1.   Introduction

Anaerobic digestion has become a widely used bioprocess to obtain bioenergy from second generation biomass, making it one of the most sustainable biotechnological processes. It comprises a series of biochemical reactions carried out by a variety of microorganisms. As it is usually done for a well-established and consolidated technology, such as anaerobic digestion, ongoing and future research studies focus on maximizing the reactor's performance, for instance, the organic matter removal capacity or the biogas production, by evaluating different reactor configurations or operating conditions. In the case of anaerobic digestion, the one-reaction mesophilic treatment is the most used system. However, the separation of the treatment in two reactions, i.e., hydrolytic/acidogenic followed by a methanogenic reactor, has shown interesting results ^[1,2].

The use of mathematical modeling for bioprocesses appeared to be a useful tool because it allows for the application of mathematical optimization theory. Nevertheless, as modeling necessarily implies an approximation of the studied process, there is difficulty in choosing the level of accuracy of the model, and it is quite admitted that the more precise the model is, the more difficult it is to study. In the last thirty years, complex models of anaerobic digestion have been developed and fitted with experimental data to predict, for example, the organic matter removal and biogas production rates. This is true in the case of the Anaerobic Digestion Model 1 (ADM1) ^[3], a model that consists of 29 dynamic state variables, which makes it very difficult to carry out qualitative analysis. Simpler models focusing on different subprocesses of ADM1 have been proposed, as in ^[4] (two reactions considered: acidogenesis and methanogenesis), ^[5] (three reactions considered: acidogenesis, acetogenesis, and methanogenesis), and ^[6] (four reactions considered: hydrolysis, acidogenesis, acetogenesis, and methanogenesis). Mathematical analyses have been carried out for these models in ^{[7,8,9,10,5,6]}.

In the case of wastewater, acidogenesis and methanogenesis are considered the key steps of the entire process, and one of the most used models is introduced in ^[4]. Concerning this model, in ^[10], the authors present a complete analysis developed in a generic manner related to the existence of steady states and their stability properties. The problem of optimizing biogas production for the model proposed in ^[4] is studied in ^{[12,13,7,8,9]}. This issue is also analyzed in ^[5] and ^[6], where models with more reactions for the anaerobic digestion process are considered. In these works, one observes that steady states can be characterized as follows: (ⅰ) the extinction (washout) of one or more of the involved microorganisms and (ⅱ) the coexistence of microorganisms. The coexistence steady states are referred to as nominal operating points, and they are typically viewed as the desired operating conditions.

Due to the fact that the acidogenic and methanogenic microorganisms possess different growth rates, pH optima has led to the development of a two-phase (or two-stage) anaerobic digestion configuration, wherein each reaction takes place separately in a different reactor. The product of the first reaction in the first reactor is then transferred to the second reactor (^[15,14]). To our knowledge, there have been no studies addressing the mathematical analysis of a single reactor compared to the mentioned sequential configuration, wherein two reactors are sequentially connected, and the reactions of the anaerobic digestion process are separated. The objective of this work is to use the two-reaction anaerobic digestion model introduced in ^[4] and to compare at steady state the two mentioned configurations with regard to the following criteria: stability issues and soluble organic matter removal performance.

The paper is organized as follows. In the next section, we introduce the two-stage anaerobic digestion model established in ^[4] and specifically the one-reactor and two-reactor models to be studied. In Section 3, we present the stability analysis for the steady states of the one-reactor model and two-reactor model, comparing the sets of parameters for these two models ensuring stability properties of the desired stable coexistence steady state. In Section 4, we compare the soluble organic matter removal performance of the two models analyzed.

2.   Two-stage anaerobic digestion models

A representation of the two-reaction anaerobic digestion model is based on the coupling of two main reactions called acidogenesis and methanization. These two reactions can be described by the so-called Anaerobic Digestion Model AM2 (see ^[4]) represented by the following dynamical system:

which is based on the chemostat model (see ^[16,11] or also ^[17]). Here, $x_A$ and $x_M$ denote the microorganism concentrations associated with the acidogenesis and methanogenesis processes, respectively; $s_A$ is the concentration of the organic substrate characterized by its COD (chemical oxygen demand) and measured in [g/l]; and $s_M$ is the total concentration of VFA (volatile fatty acids) measured in [mmol/l]. The dilution rate of the continuously operated bioreactor is denoted by $D$ (i.e., $D = Q/V$, where $Q$ is the input and output flow rate of water and $V$ is the constant volume of water present in the bioreactor). The parameter $s^{in}_{A}$ represents the input acidogenesis substrate concentration^*. Coefficients $k_j$, $j = 1, 2, 3$ are positive parameters called pseudostochiometric coefficients associated with the bioreactions, and the functions $\mu_A(\cdot)$ and $\mu_M(\cdot)$ are the so-called growth rate functions.

^* The input concentration of methanogenic substrate in the influent is usually considered negligible; thus, for simplicity, it is not considered in this work.

The dilution rate $D$ is the decision variable due to the fact that one can control the input and output flow rates, both being equal in a continuously operated bioreactor.

The growth functions or kinetics $\mu_A(\cdot)$ and $\mu_M(\cdot)$ are usually of a Monod and Haldane type (see Example 1 for the definitions), respectively, as in ^[4,7,8]. However, in this paper, we consider generic kinetics $\mu_A(\cdot)$ and $\mu_M(\cdot)$, satisfying the following qualitative properties:

Assumption 1. The function $\mu_A : [0, +\infty) \longrightarrow \mathbb{R}$ is concave, increasing, and continuously differentiable, with $\mu_A(0) = 0$ and $\lim \limits_{s_A \to +\infty} \mu_A(s_A) = \bar \mu_A$, for some $\bar \mu_A > 0$.

Assumption 2. The function $\mu_M : [0, +\infty) \longrightarrow \mathbb{R}$ is continuously differentiable, and for some $S^{max}_M > 0$, it is concave and increasing over $[0, S^{max}_M]$ and decreasing over $(S^{max}_M, +\infty)$, with $\mu_M(0) = 0$, $\mu_M(s_M) > 0$ for all $s_M > 0$, and also $\lim \limits_{s_M \to +\infty} \mu_M(s_M) = 0$.

Observe that under Assumptions 1 and 2, the set $\mathcal{D}': = [0, +\infty)^4$ is invariant under (2.1), which allows us to study only positive state variables.

Example 1. It is straightforward to check that if $\mu_A(\cdot)$ is a Monod function and $\mu_M(\cdot)$ is of Haldane type, then they satisfy Assumptions 1 and 2. Recall that the Monod function is defined (see ^[11]) by

and the Haldane function by

where $\bar \mu_A, k_A, \bar \mu_M, k_M$, and $k_I$ are positive parameters. In this case, the value $S^{max}_M$ in Assumption 2 is given by $S^{max}_M = \sqrt{k_M k_I}$.

We shall call the one-reactor model the system presented in (2.1), where the two reactions (acidogenesis and methanogenesis) of the anaerobic digestion process occur in one bioreactor (see Figure 1).

2.1.   Sequential two-reactor anaerobic digestion model

In this section, we introduce a variation of (2.1), considering the two reactions (acidogenesis and methanogenesis) of the anaerobic digestion process, separated in two bioreactors connected sequentially, with separation of the two biomass populations. Thus, the acidogenesis process takes place in the first reactor, where only the acidogenic microorganism is located. The output flow of the first bioreactor is filtered^† to the input flow to the second bioreactor, which consists only of acidogenic and methanogenic substrates. Finally, the methanogenic process takes place in the second bioreactor, where only the methanogenic microorganism is present as a biomass. The scheme of this sequential configuration is depicted in Figure 2.

^† We neglect the effect of the filter device on the flow rate.

The described configuration is represented by the following dynamical system:

The notations used in (2.2) are basically the same as those in (2.1). The subscripts $1$ and $2$, for the substrate concentrations, refer to which bioreactor these concentrations are in. Recall that in the first reactor, there is only acidogenic microorganisms, and in the second bioreactor, there is only methanogenic microorganisms. The notation $D$ stands for the global dilution rate of the system, which is given by the input flow divided by the sum of the volumes of each bioreactor, the volumes being denoted by $V_1$ and $V_2$, respectively. The parameter $r$ represents the ratio of the total volume $V = V_1+V_2$ attributed to the first tank, i.e., $r = \frac{V_1}{V_1+V_2} = \frac{V_1}{V} \in (0, 1)$.

The growth functions $\mu_A(\cdot)$ and $\mu_M(\cdot)$ in (2.2) are the same as those used in (2.1). Thus, as with the one-reactor model (2.1), one can prove straightforwardly that the set $\mathcal{D}: = [0, +\infty)^6$ is invariant under (2.2), reducing our study to only positive state variables.

3.   Stability analysis

In this section, we recall the stability analysis carried out in ^[10] for the one-reactor model (2.1), and we study the stability properties of equilibria corresponding to the two-reactor model (2.2).

3.1.   Stability analysis of the one-reactor model

In ^[10], the authors study the steady states of the one-reactor model (2.1), analyzing their local and global stability behaviors. We will first recall these results, which will be needed in the next sections. Let us first introduce some notations. It is well known that in a chemostat model, the existence and value of a steady state depends on the value of the dilution rate, the substrate input concentration and the growth functions considered. If $D < \bar \mu_A$, then we define

Notice that $s_A^{eq}$ is well defined because, from Assumption 1, the acidogenic growth function $\mu_A(\cdot)$ is increasing.

From Assumption 2, $D < \mu_M(S^{max}_M)$ ensures the existence of two solutions of the equation $\mu_M(s) = D$, which will be denoted by $s_{Ms}^{eq}$ and $s_{Mu}^{eq}$, such that $s_{Ms}^{eq} < s_{Mu}^{eq}$ (the subscript $s$ stands for stable, while $u$ stands for unstable). Then, assuming $D < \bar \mu_A$, we can define

It can be proven (see ^[10]) that the following are the only possible steady states $(x_A^*, s_A^*, x_M^*, s_M^*)$ of the one-reactor model (2.1):

Finally, to present the stability results, we define the sets $\mathcal{D'}: = [0, +\infty)^4$, $\mathcal{D'_*}: = (0, +\infty)\times[0, +\infty)^3$ and $\mathcal{D'}_{**}: = (0, +\infty)\times[0, +\infty)\times (0, +\infty)\times[0, +\infty)$.

Theorem 3.1 (^[10]). The stability behavior of steady states corresponding to the one-reactor model (2.1), depending on the values of $D$, $s^{in}_{A}$ and growth functions $\mu_A(\cdot)$ and $\mu_M(\cdot)$, is the following:

where $GAS$ stands for Globally Asymptotically Stable and $LAS$ stands for Locally Asymptotically Stable.

Usually, the stable coexistence steady states $e_{2s}$ (which exist for the range of dilution rates $D$ indicated in the above table) are referred to as nominal operating points, and they are typically viewed as the desired operating conditions. One of our objectives is to compare the conditions ensuring the stability (at least local) of this steady state between the one-reactor model and the two-reactor model. For this purpose, in the next section, we develop the same analysis carried out in ^[10] but for the two-reactor model (2.2).

3.2.   Stability analysis of the two-reactor model

To provide an equivalent result of Theorem 3.1 in the context of the two-reactor model, it is worth noting that the model consists of two usual chemostat models connected in series. Consequently, the existence of the steady states depends on the value of the dilution rate and on the growth functions (see, for instance, ^[18]). Therefore, we introduce the following notations. If $\frac{D}{r} < \bar \mu_A$, then we define

As emphasized before, $s^{eq}_{A1}$ is well defined from Assumption 1.

Assuming that $\frac{D}{1-r}\leq \mu_M(S^{max}_M)$, Assumption 2 implies the existence of two solutions of the equation $\mu_M(s) = \frac{D}{1-r}$, which will be denoted by $s^{eq}_{M2s}$ and $s^{eq}_{M2u}$ such that $s^{eq}_{M2s} < s^{eq}_{M2u}$. Then, if one has $\frac{D}{r} < \bar \mu_A$, one can define

The next result characterizes the existence of steady states under hypotheses on the dilution rate and the growth functions.

Proposition 1. The sequential model (2.2) admits at most four steady states, denoted in the form $(x_A^*, s_{A1}^*, s_{M1}^*, x_M^*, s_{A2}^*, s_{M2}^*)$:

Moreover,

1. If $\frac{D}{r}\geq\mu_A(s^{in}_{A})$, then the system only admits $E_0$.

2. If $\frac{D}{r} < \mu_A(s^{in}_{A})$ with $\frac{D}{1-r} > \mu_M(S^{max}_M)$ or with $\frac{D}{1-r}\leq\mu_M(S^{max}_M)$ and $s^{eq}_{M1} < s^{eq}_{M2s}$, then the system only admits $E_0$ and $E_1$.

3. If $\frac{D}{r} < \mu_A(s^{in}_{A})$ with $\frac{D}{1-r}\leq\mu_M(S^{max}_M)$ and $s^{eq}_{M2s}\leq s^{eq}_{M1}\leq s^{eq}_{M2u}$, then the system admits $E_0$, $E_1$ and $E_{2s}$.

4. If $\frac{D}{r} < \mu_A(s^{in}_{A})$ with $\frac{D}{1-r}\leq\mu_M(S^{max}_M)$ and $s^{eq}_{M1}\geq s^{eq}_{M2u}$, then the system admits the four steady states $E_0$, $E_1$, $E_{2s}$ and $E_{2u}$.

Proof. See the proof in the Appendix.

In the previous result, we use coexistence to refer to those steady states where the concentration of both microorganisms is positive, even though these populations are not in the same bioreactor.

Remark 1. It is worth noting that even if Proposition 1 is given for general growth functions satisfying Assumptions 1 and 2, explicit expressions of those functions are usually needed to know which of the four cases of the same proposition is fulfilled. Minding this fact, it can be observed that assumption $s^{eq}_{M2s}\leq s^{eq}_{M1}\leq s^{eq}_{M2u}$ in the third case in Proposition 1 can be substituted by $\mu_M(s^{eq}_{M1})\leq \frac{D}{1-r}$.

To establish the stability behaviors of the steady states in the two-reactor model, similar to the analysis in the previous section, we define the sets $\mathcal{D} = [0, +\infty)^6$, $\mathcal{D_*}: = (0, +\infty)\times[0, +\infty)^5$ and $\mathcal{D}_{**}: = (0, +\infty)\times[0, +\infty)^2\times (0, +\infty)\times[0, +\infty)^2$.

Theorem 3.2. The stability properties of steady states corresponding to the two-reactor model (2.2), depending on the values of $D$, $r$, $s^{in}_{A}$ and growth functions $\mu_A(\cdot)$ and $\mu_M(\cdot)$, are presented in the following table:

Proof. See the proof in the Appendix.

Regarding theorems 3.1 and 3.2, we prove in the next result that the bioprocess in one reactor (given by (2.1)) is more stable than the sequential process in two reactors (given by (2.2)) in the sense that, whenever stable coexistence exists in the two-reactor model and is locally or globally stable, then this steady state also exists in the one-reactor model, being locally or globally stable too. To prove this result, consider the following sets of parameters:

where $r \in (0, 1)$ is the parameter that characterizes the two-reactor model (2.2).

Proposition 2. For the sets defined by (3.1)- (3.4), one has

Proof. See the proof in the Appendix.

The previous result indicates that for parameters $(D, s_{A}^{in})$ for which the steady state $E_{2s}$ of model (2.2) is stable (at least locally), the corresponding steady state $e_{2s}$ of (2.1) (associated with $(D, s_{A}^{in})$) is also stable, independently of the distribution of volumes in the two-reactor model, with a distribution parameterized by $r \in (0, 1)$. In other words, there are more parameters $(D, s_{A}^{in})$ ensuring the stability of the coexistence steady state for the one-reactor model than for the two-reactor model. In Figure 3, we depict the sets $P_r \cup \hat{P}_r$ and $P \cup \hat{P}$ in the plane of parameters $s_{A}^{in}-D$ for different values of $r$. These sets were obtained considering the growth functions introduced in Example 1 and parameter values found in ^[4], also presented in Table 1.

We can obtain a similar result to that in Proposition 2 regarding the global stability of the stable coexistence in both models. For this purpose, let us introduce the following set:

where $\mu_M^{-1}(D)_+$ represents the greatest solution $s$ of the equation $\mu_M(s) = D$.

Proposition 3. For the sets defined in (3.2), (3.4), and (3.5), one has the following inclusion:

Proof. See the proof in the Appendix.

In the previous result, the intersection with the set $C$ is necessary because the inclusion of $\hat P_r \subseteq \hat P$ is not true. That is, there could exist some parameters $(D, s^{in}_{A})$ for which the stable coexistence is globally stable in the two-reactor model, but for these parameters, the stable coexistence steady state of the one-reactor model can be just locally stable.

In Figure 4, we depict sets $\hat P_r \cap C$ and $\hat P$ for different values of $r$ using the parameter values of Table 1. One can observe that these sets do not differ so much with sets presented in Figure 3.

4.   Comparison of the organic matter removal performance

Finally, in this section, we compare substrate concentrations at steady state in both models. First, regarding the possible steady states in the one-reactor and two-reactor models, the desired equilibria correspond to stable coexistence because, if one or two species are washed out, there is no removal of soluble organic matter. Therefore, to compare the soluble organic matter removal performance of the two configurations under study, it is sufficient to compare the substrate concentrations of the stable coexistence equilibria.

Proposition 4. For $r \in (0, 1)$, consider $(D, s^{in}_{A})$ in the set $P_r$ defined by (3.3), $e_{2s} = (x_A^{eq}, s_A^{eq}, X^{eq}_{Ms}, s_{Ms}^{eq})$ the stable coexistence steady state for the one-reactor model (2.1) and $E_{2s} = (X^{eq}_{A}, s^{eq}_{A1}, s^{eq}_{M1}, X^{eq}_{Ms}, s^{eq}_{A1}, s^{eq}_{M2s})$ the stable coexistence steady state for the two-reactor model (2.2), with steady states associated with parameters $(D, s^{in}_{A})$; then,

Proof. See the proof in the Appendix.

The above results simply establish that the organic matter removal capacity is better in the one-reactor model. Nevertheless, if separating the two microorganism populations in two reactors, allowing for improvement of the performances of both populations in term of their growth functions (e.g., adapting in each reactor the best pH and temperature for each microorganism), the organic matter removal capacity could be better in the two-reactor configuration. Indeed, if $\mu_A(\cdot)$ and $\mu_M(\cdot)$ are the growth functions for the microorganisms in the one-reactor configuration and $\hat \mu_A(\cdot)$ and $\hat \mu_M(\cdot)$ are the growth functions in the two-reactor configuration, with $\mu_A(\cdot) < \hat \mu_A(\cdot)$ and $\mu_M(\cdot) < \hat \mu_M(\cdot)$, one could find $r \in (0, 1)$ (distribution of the volume $V$ into two volumes $rV$ and $(1-r)V$) to improve the organic matter removal capacity, as is established in the next result.

Proposition 5. Consider $\varepsilon_A > 0$ and $\varepsilon_M > 0$ such that

and $r \in (0, 1)$ is such that

then, the organic matter removal performance is better for the two-reactor configuration (associated with the distribution of volume $r$ satisfying (4.2)).

Proof. See the proof in the Appendix.

In the above proposition, notice that the improvement of the growth functions in separating both populations, measured with $\varepsilon_A > 0$ and $\varepsilon_M > 0$, must not be too low. In fact, from (4.2), one deduces that $\varepsilon_A \varepsilon_M \ge 1$. For instance, if the improvements are equal, $\varepsilon_A = \varepsilon_M = \varepsilon$, then $\varepsilon \ge 1$, meaning that (from (4.1)) the growth functions of the microorganisms separated have to be at least double the growth functions when the populations are in the same reactor. In this case, defining $I(\varepsilon): = [1/(1+ \varepsilon), \varepsilon/(1+ \varepsilon)]$, the interval established in (4.2) (when $\varepsilon_A = \varepsilon_M = \varepsilon$), where $r$ has to belong, one has that $I(\varepsilon) \neq \emptyset$ if and only if $\varepsilon \ge 1$, also $I(1) = \{1/2\}$, and $I(\varepsilon) \to (0, 1)$ when $\varepsilon \to + \infty$, as is depicted in Figure 5.

5.   Conclusions

In this work, using the two-reaction anaerobic digestion model proposed in ^[4], we have studied two configurations: (ⅰ) a single bioreactor, where the acidogenesis and methanogenesis reactions occur inside, and (ⅱ) two bioreactors connected sequentially, where each reaction occurs separately in each reactor. We have proven that the set of operating conditions consisting of dilution rates $D$ and input acidogenic substrate concentrations $s^{in}_{A}$, which ensure that the coexistence steady state is stable (at least locally), is larger for the one-reactor model than the two-reactor model, implying that larger dilution rates can be applied in the one-reactor model. Furthermore, the soluble organic matter removal performance of the one-reactor model is better than that of the two-reactor model, leading to the conclusion that dividing the reactions in two reactors does not improve the stability of the anaerobic digestion process or the soluble organic matter removal capacity unless the improvement of the growth functions of both microorganisms populations is considerably important (established in a precise manner in Proposition 5) at the moment of separating them into two reactors. This may partly explain why this type of configuration has not reached the same level of development as that reached by the conventional one-phase reactor. In fact, two-stage systems are in decline, at least for municipal solid waste ^[19].

Acknowledgments

This research benefited from the financial support of FONDECYT grant (Chile) N 1160567 and Basal Project CMM Universidad de Chile. The fourth author was also funded by the CONICYT PFCHA/DOCTORADO BECAS CHILE/2017–21171813.

Conflict of interest

All authors declare no conflicts of interest in this paper.

Appendix

5.1.   Proof of Proposition 1

We are interested in finding the positive solutions of the algebraic system

that is, with $(x_A, s_{A1}, s_{M1}, x_M, s_{A2}, s_{M2}) \in \mathcal{D}$.

One can observe that the two first equations of (5.1) describe the steady state of an usual one specific chemostat with an increasing bounded growth function. Therefore, it is well known (see ^[18] for instance) that if $\frac{D}{r}\geq \mu_A(s_{A}^{in})$, then the only positive steady state for this submodel is $(0, s_{A}^{in})$, whereas if $\frac{D}{r} < \mu_A(s_{A}^{in})$, it admits $(0, s_{A}^{in})$ and $(X_A^{eq}, s_{A1}^{eq})$ as positive steady states.

Next, we can observe that if we consider $(0, s_{A}^{in})$ for the solution of the first two equations of (5.1), then one necessarily has that $(0, s_{A}^{in}, 0, 0, s_{A}^{in}, 0)$. Indeed, if $x_A = 0$, then it is straightforward that $s_{M1} = 0$ and $s_{A2} = s_{A}^{in}$. Moreover, multiplying the fourth equation of (5.1) by $k_3$ and then adding it to the sixth one yields $s_{M2} + k_3x_M = 0$, which implies that $s_{M2} = x_M = 0$.

Let us suppose now that $\frac{D}{r} < \mu_A(s_{A}^{in})$ and consider $(X_A^{eq}, s_{A1}^{eq})$ as the solutions of the first two equations of (5.1). Then, it is obvious that the third equation of (1) vanishes for $s_{M1} = s_{M1}^{eq}$, which is positive because $s_{A}^{in} \geq s_{A1}^{eq}$ by $\frac{D}{r} < \mu_A(s_{A}^{in})$ and the fact that $\mu_A(\cdot)$ is increasing.

Moreover, we observe that the fourth and sixth equations of (5.1) describe the steady state of an usual chemostat of one species under a Haldane-type growth function with substrate input $s_{M1}^{eq}$ (see ^[18]). Thus, if $\frac{D}{1-r} > \mu_M(S_{^M}^{max{\rm{ }}})$, one can define $s_{M2s}^{eq}, s_{M2u}^{eq}, X_{Ms}^{eq}, X_{Mu}^{eq}$ such that $s_{M2s}^{eq} < s_{M2u}^{eq}$, and we have the three following cases:

● if $X_{Ms}^{eq}\leq 0$, that is, $s_{M1}^{eq} \leq s_{M2s}^{eq}$, then the only steady state of the submodel is $(0, s_{M1}^{eq})$;

● if $X_{Mu}^{eq}\leq0$, that is, $s_{M2s}^{eq} \leq s_{M1}^{eq} \leq s_{M2u}^{eq}$, then $(0, s_{M1}^{eq})$ and $(X_{Ms}^{eq}, s_{M2s}^{eq})$ are the two steady states of the submodel;

● if $X_{Mu}^{eq} \geq 0$, that is $s_{M1}^{eq} \geq s_{M2u}^{eq}$ then $(0, s_{M1}^{eq})$, $(X_{Ms}^{eq}, s_{M2s}^{eq})$, and $(X_{Mu}^{eq}, s_{M2u}^{eq})$ are the three steady states of the model.

Finally, we observe that the fifth equation of (5.1) always leads to $s_{A2}^{eq} = s_{A1}^{eq}$, which concludes our determination of the steady states.

5.2.   Proof of Theorem 3.2

First, we compute the Jacobian matrix of the model 2.2:

and we evaluate this Jacobian in each of the four steady states, $E_0, E_1, E_{2s}, E_{2u}$ defined in Proposition 1, deducing the conditions for the stability of each one.

1. For $\lambda \in \mathbb{R}$, one has that

hence, the eigenvalues of $J(E_0)$ are $-\frac{D}{1-r}$, $-\frac{D}{r}$ and $\mu_A(s_{A}^{in}) - \frac{D}{r}$. This implies that if $\mu_A(s_{A}^{in}) < \frac{D}{r}$, then $E_0$ is LAS; if $\mu_A(s_{A}^{in}) > \frac{D}{r}$, then $E_0$ is unstable; and if $\mu_A(s_{A}^{in}) = \frac{D}{r}$, it is nonhyperbolic.

2. Suppose now that $\frac{D}{r} < \mu_A(s_A^{in})$ so that $E_1$ is well defined in $\mathcal{D}$. Then, one has that

where:

Therefore, the eigenvalues of $J(E_1)$ are $-\frac{D}{1-r}$, $-\frac{D}{r}$, $\mu_M(s_{M1}^{eq}) < \frac{D}{1-r}$ and also the eigenvalues of $A$. However, it can be observed that $A$ is the Jacobian matrix of the usual chemostat model with Monod-type growth function described by Equations 1 and 2 of (2.2) evaluated at $(X_A^{eq}, s_{A1}^{eq})$. Thus, it is well known that the eigenvalues of $A$ are with negative real part whenever $\frac{D}{r} < \mu_A(s_A^{in})$ (see ^[18]).

Consequently, we deduce that $E_1$ is LAS if $\mu_M(s_{M1}^{eq}) < \frac{D}{1-r}$, unstable if $\mu_M(s_{M1}^{eq}) > 0$ and nonhyperbolic if $\mu_M(s_{M1}^{eq}) = 0$.

Therefore, if $\mu_M(S_{^M}^{max{\rm{ }}}) < \frac{D}{1-r}$, then $\mu_M(s_{M1}^{eq}) < \frac{D}{1-r}$, and $E_1$ is LAS.

Suppose from now on that $\frac{D}{1-r} \leq \mu_M(S_{^M}^{max{\rm{ }}})$ so that $s_{M2s}^{eq}$ and $s_{M2u}^{eq}$ are well defined. If $s_{M1}^{eq} < s_{M2s}^{eq}$ or $s_{M1}^{eq} > s_{M2u}^{eq}$, then $\mu_M(s_{M1}^{eq}) < \frac{D}{1-r}$ and $E_1$ is LAS. If $s_{M2s}^{eq} < s_{M1}^{eq} < s_{M2u}^{eq}$, then $\mu_M(s_{M1}^{eq}) > \frac{D}{1-r}$ and $E_1$ is unstable. If $s_{M1}^{eq} = s_{M2s}^{eq}$ or $s_{M1}^{eq} = s_{M2u}^{eq}$, then $E_1$ is nonhyperbolic.

3. Suppose now that $s_{M1}^{eq} > s_{M2s}^{eq}$. Therefore, $E_{2s} \in \mathcal{D}$. Then, one has

where

First, we observe that if $s_{M1}^{eq} = s_{M2s}^{eq}$ or $s_{M2s}^{eq} = S_{^M}^{max{\rm{ }}}$, then one has, respectively, that $X_{Ms}^{eq} = 0$ or $\mu_M'(s_{M2s}^{eq}) = 0$, which implies that $det(B_s - \lambda I) = \lambda (\lambda + \frac{D}{1-r})$, leading to $E_{2s}$ being nonhyperbolic.

Therefore, let us suppose that $s_{M1}^{eq} > s_{M2s}^{eq}$ and $s_{M2s}^{eq} < S_{^M}^{max{\rm{ }}}$. Then, $\mu_M'(s_{M2s}^{eq}) > 0$, implying that $tr(B_s) < 0$ and $det(B_s) > 0$ and leading to $E_{2s}$ being LAS.

4. Finally, we suppose that $s_{M1}^{eq} \geq s_{M2u}^{eq}$ so that $E_{2u}\in\mathcal{D}$. Then, one has

where

As done in the previous case, it can be observed that $s_{M1}^{eq} = s_{M2u}^{eq}$ or $s_{M2u}^{eq} = S_{^M}^{max{\rm{ }}}$ implies that $E_{2u}$ is nonhyperbolic.

Therefore, we assume that $s_{M1}^{eq} > s_{M2u}^{eq}$ and $s_{M2u}^{eq} < S_{^M}^{max{\rm{ }}}$. This implies that $det(B_u) < 0$. Consequently we have that the discriminant of the characteristic polynomial of $B_u$ is positive, yielding that the greatest eigenvalue of $B_u$ is defined by $\lambda = \frac{-tr(B_u) + \sqrt{tr^2(B_u) - det (B_u)}}{2}$, which is positive. Therefore, $E_{2u}$ is unstable.

Now, we establish under what conditions steady states $E_0$, $E_1$, and $E_{2s}$ are GAS. For this purpose, the idea is to use the cascade structure of (2.2) to separately study the global asymptotic behavior of the state variables. Indeed, one can observe that (2.2) admits two independent subsystems: one composed by Equations 1 and 2 only, and the other, by Equations 1, 2, 3 and 5. This implies that once the initial condition $(x_A^0, s_{A1}^0)$ of the two first state variables are fixed, the state variables $s_{M1}$ and $s_{A2}$ satisfy the nonautonomous system

where $(x_A(t), s_{A1}(t))$ is the solution of the first independent subsystem associated with the initial condition $(x_A^0, s_{A1}^0)$.

Then, one can fix an initial condition $(s_{M1}^0, s_{A2}^0)$ of the third and fifth state variables and obtain that $x_M$ and $s_{M2}$ satisfy the nonautonomous system

where $(x_A(\cdot), s_{A1}(\cdot), s_{M1}(\cdot), s_{A2}(\cdot))$ is the solution of the second independent subsystem associated with the initial condition $(x_A^0, s_{A1}^0, s_{M1}^0, s_{A2}^0)$.

Let $(x_A(\cdot), s_{A1}(\cdot), s_{M1}(\cdot), x_M(\cdot), s_{A2}(\cdot), s_{M2}(\cdot))$ be the solution of (2.2) associated with an initial condition $(x_A^0, s_{A1}^0, s_{M1}^0, x_M^0, s_{A2}^0, s_{M2}^0) \in \mathcal{D} = [0, +\infty)^6$.

Suppose for the moment that $\frac{D}{r} > \mu_A(s_{A}^{in})$. We claim that $E_0$ is GAS on $\mathcal{D}$. Indeed, because $\frac{D}{r} > \mu_A(s_{A}^{in})$, it is well known that $(x_A(\cdot), s_{A1}(\cdot))$ converges to $(0, s_{A}^{in})$. Therefore, (5.2) is asymptotically autonomous (see ^[18]) with the limit system

which only admits $(0, s_{A}^{in})$ as steady state, which is GAS on $[0, +\infty)^2$. Therefore, one can apply the Thieme theorem (see Theorem A.9 of ^[18]) to obtain that $(s_{M1}(\cdot), s_{A2}(\cdot))$ converges to $(0, s_A^{in})$.

This in turn implies that (5.3) is also asymptotically autonomous with the limit system

which only admits $(0, 0)$ as steady state, which is GAS on $[0, +\infty)^2$. Therefore, one can apply the Thieme theorem once again to obtain that $(x_M(\cdot), s_{M2}(\cdot))$ converges to $(0, 0)$, which proves the claim.

Suppose now that $\frac{D}{r} < \mu_A(s_A^{in})$ with $x_A^0\not = 0$. Then, it is well known that $(x_A(\cdot), s_{A1}(\cdot))$ converges to $(x_A^{eq}, s_{A1}^{eq})$. Therefore, (5.2) is asymptotically autonomous with the limit system

which only admits $(s_{M1}^{eq}, s_{A1}^{eq})$ as steady state, which is GAS on $[0, +\infty)^2$. Therefore, by the Thieme theorem, one has that $(s_{M1}(\cdot), s_{A2}(\cdot))$ converges to $(s_{M1}^{eq}, s_{A1}^{eq})$.

Then, one obtains that (5.3) is once again asymptotically autonomous, but this time with the limit system

Then, assuming that $\frac{D}{1-r} > \mu_M(S_{^M}^{max{\rm{ }}})$, it is well known that $(0, s_{M1}^{eq})$ is the only steady state of (5.4), which is GAS on $[0, +\infty)^2$, implying by Thieme theorem that $E_1$ is GAS on $\mathcal{D}_*$ as wanted.

Assuming that $\frac{D}{1-r} < \mu_M(S_{^M}^{max{\rm{ }}})$ with $s_{M1}^{eq} < s_{M2s}^{eq}$ also implies that $(0, s_{M1})$ is the only steady state of (5.4), which is GAS on $[0, +\infty)^2$, also yields that $e_1$ is GAS on $\mathcal{D}_*$.

Finally, if $\frac{D}{1-r} < \mu_M(S_{^M}^{max{\rm{ }}})$ with $s_{M2s}^{eq} < s_{M1}^{eq} < s_{M2u}^{eq}$, it well known that supposing $x_M^0 \not = 0$ gives that $(X_{Ms}^{eq}, s_{M2s}^{eq})$ is the only steady state of (4), which is GAS on $(0, +\infty)\times[0, +\infty)$. Therefore, applying the Thieme theorem once again implies that $E_{2s}$ is GAS on $\mathcal{D}_{**}$ and concludes the proof.

5.3.   Proof of Proposition 2

For $r \in (0, 1)$, let us consider some values of $D, s^{in}_A$ such that $(D, s^{in}_A) \in P_r \cup \hat{P}_r$. This implies that the following conditions hold, according to Theorem 3.2:

(A1) $\frac{D}{r} < \mu_A(s^{in}_{A}),$

(A2) $\frac{D}{1-r} < \mu_M(S_{^M}^{max{\rm{ }}}),$

(A3) $s_{M2s}^{eq} < s_{M1}^{eq}.$

Additionally, let us recall that by Theorem 3.1, the following conditions must hold for $D, s^{in}_A$ such that $(D, s^{in}_A) \in P \cup \hat{P}$

(B1) $D < \mu_A(s^{in}_{A}),$

(B2) $D < \mu_M(S_{^M}^{max{\rm{ }}}),$

(B3) $s_{Ms}^{eq} < s_{Mw}^{eq}.$

Given that $D < \frac{D}{r}$ and $D < \frac{D}{1-r}$ if $r\in(0, 1)$, it is immediate that $(A1) \implies (B1)$ and $(A2) \implies (B2)$.

For the last implication, let us recall that:

The result of Proposition 4 indicates to us that $s_{A1}^{eq} > s_A^{eq}$ when $r\in (0, 1)$, and this in turn implies that $s_{Mw}^{eq} > s_{M1}^{eq}$. Alternately, Proposition 4 also implies that $s_{Ms}^{eq} < s_{M2s}^{eq}$. Taking into account these inequalities and assuming condition (A3), one has:

thus, $(A3) \implies (B3)$, completing the proof.

5.4.   Proof of Proposition 3

Let us assume that $(D, s_{A}^{in}) \in \hat{P}_r \cap C$. This implies that the following conditions hold:

(C1) $\frac{D}{r} < \mu_A(s^{in}_{A})$,

(C2) $\frac{D}{1-r} < \mu_M(S_{^M}^{max{\rm{ }}}),$

(C3) $s_{M2s}^{eq} < s_{M1}^{eq} < s_{M2u}^{eq},$

(C4) $s_{Mw}^{eq} < s_{Mu}^{eq}.$

Let us recall that for $(D, s_{A}^{in}) \in \hat{P}$, the following conditions must hold:

(D1) $D < \mu_A(s^{in}_{A})$

(D2) $D < \mu_M(S_{^M}^{max{\rm{ }}})$

(D3) $s_{Ms}^{eq} < s_{Mw}^{eq} < s_{Mu}^{eq}$

The result of Proposition 2 implies that $(C1)\implies (D1)$ and $(C2) \implies (D2)$. Additionally, Proposition $2$ implies that:

All this, together with condition (C4), implies (D3), thus proving the inclusion.

5.5.   Proof of Proposition 4

First, let us consider the steady states corresponding to acidogenic substrates in both models. We recall that the formula for these equilibria are, for the one-reactor and two-reactor models, respectively:

From Assumption 1, the function $\mu_A(\cdot)$ is increasing, and thus, given that, for $r \in (0, 1)$, $\frac{D}{r} > D$, then $s_A^{eq} < s_{A1}^{eq}$ must hold.

Now, let us consider the steady states corresponding to methanogenic substrates in both models. Considering the steady states in each case, for the one-reactor and two-reactor models, the formula for the equilibria are, respectively:

where $\mu_M^{-1}(d)$ represents the lowest solution of the equation $\mu_M(s) = d$.

Consider Assumption 2 and the fact that for the one-reactor model, we have that $s_{Ms}^{eq} \geq s_{Mu}^{eq}$, $\mu(s_{Ms}^{eq})\leq \mu(S_{^M}^{max{\rm{ }}})$ and $\mu(s_{Mu}^{eq}) \leq \mu(S_{^M}^{max{\rm{ }}})$; thus, $s_{Ms}^{eq}$ is in the increasing part of $\mu_M(\cdot)$.

For the two-reactor model, we have that $s_{M2s}^{eq} \ge _{M2u}^{eq}$, $\mu(s_{M2s}^{eq})\leq \mu(S_{^M}^{max{\rm{ }}})$ and $\mu(s_{M2u}^{eq}) \leq \mu(S_{^M}^{max{\rm{ }}})$. This implies that $s_{M2s}^{eq}$ is located in the increasing section of $\mu_M(\cdot)$.

Given that both steady states are in the increasing section of $\mu_M(\cdot)$ and given that $\frac{D}{1-r} > D$ when $r \in (0, 1)$, it is implied that $s_{M2s}^{eq} > s_{Ms}^{eq}$, thus concluding the desired result.

5.6.   Proof of Proposition 5

Using the same notation as in Proposition 4, we have that $s_A^{eq}$ and $s_{Ms}^{eq}$ are the acidogenic and methanogenic substrates (at steady states) for the one-reactor configuration and that $s_{A1}^{eq}$ and $s_{M2s}^{eq}$ are the corresponding concentrations in the two-reactor configuration. Then,

Alternately, from inequalities (4.2), one has

Therefore,

where the last inequality is obtained from the condition (4.1). Since $\mu_A(\cdot)$ is increasing, we deduce that $s_A^{eq} \ge s_{A1}^{eq}$.

Similarly,

As $s_{Ms}^{eq}$ and $s_{M2s}^{eq}$ are in the increasing part of the functions $\mu_M(\cdot)$ and $\hat \mu_M(\cdot)$, respectively, we obtain $s_{Ms}^{eq} \ge s_{M2s}^{eq}$, thus concluding that the organic matter removal performance is better for the two-reactor configuration because $s_A^{eq} \ge s_{A1}^{eq}$ and $s_{Ms}^{eq} \ge s_{M2s}^{eq}$.