# A Dynamic Analysis for an Anaerobic Digester: Stability and Bifurcation Branches

Authors: Alejandro Rincón; Julieta Villa; Fabiola Angulo; Gerard Olivar
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

This work presents a dynamic analysis for an anaerobic digester, supported on the analytical application of the indirect Lyapunov method. The mass-balance model considered is based on two biological reaction pathways and involves both Monod and Haldane representations of the specific biomass growth rates. The dilution rate, the influent concentration of chemical oxygen demand (COD), and the influent concentration of volatile fatty acids (VFA) are considered as stability parameters. Several characteristics are determined analytically for the normal operation equilibrium point: (i) equilibrium coordinates, (ii) parameter conditions that lead to positive values of the equilibrium state variables, (iii) parameter conditions for locally stable nature of the equilibrium, (iv) coordinates for the local bifurcation points—fold and transcritical—, and (v) coordinates of the crossing between bifurcation points. These factors are computed analytically and explicitly as expressions of the dilution rate and the influent concentrations of COD and VFA.


---

## Body

## 1. Introduction

The anaerobic digestion processes have several advantages with respect to aerobic digestion ones, for instance, lower sludge production and higher pollution degradation [1–4]. The main drawback of anaerobic digestion processes is the easy destabilization. Indeed, large variation of dilution rate or the influent organic load may lead to the so-called biomass washout. Washout phenomenon involves the inactivation of biomass and the accumulation of volatile fatty acids [1, 3–6]. In practice, the dilution rate range is constrained in order to avoid it (cf. [4, pp. 1615]).The dependence of desired or undesired operation of biological systems on model parameters is usually analyzed by means of the indirect Lyapunov method (see [7–21]). This method is based on the linearization of the vector field and involves the definition of the following factors: (i) equilibrium points, (ii) parameter constraints that lead to expected intervals of the equilibrium state variables, (iii) Jacobian matrix, eigenvalues, and stability of each equilibrium, (iv) coordinates of local bifurcation points, and (v) parameter regions that lead to either desired or undesired system operation. The analytical development of the dynamics has many advantages with respect to the graphical analysis, as can be noticed from [4, 8, 20, 22, 23]: (i) it allows a fast computation of equilibrium points, bifurcation points and other properties, (ii) it allows an in-depth, complete, and efficient understanding of the influence of model parameters on the equilibria, bifurcation points and on other characteristics, specially in the case of three or more bifurcation parameters. Indeed, when there are three or more bifurcation parameters, the graphical anlaysis is cumbersome and the analytical results can be more practical.In the case of anaerobic digesters, the application and results of the method depend on the mass-balance considered. A mass balance model based on two biological reaction pathways involves two biomass species and two substrate concentrations. For instance, the model in [1] is based on two pathways: acidogenesis and methanization. The acidogenesis pathway involves the concentration of acidogenic biomass and the concentration of COD, whereas the methanization pathway involves the concentration of methanogenic biomass and the concentration of VFA. Each of these four concentrations is associated with a differential equation. The dilution rate, the influent COD concentration, and the influent VFA concentration are considered as stability parameters (see [2, 4, 22–24]—model IV—, [25]). The nonlinear nature of the specific biomass growth rate leads to multiple equilibria, but only some branch sections correspond to realistic process operation, either normal process operation or biomass washout condition. One equilibrium branch section corresponds to normal process operation, and it involves positive concentrations. The other equilibrium branch section corresponds to biomass washout, as one of the two biomass concentrations is zero. These realistic equilibrium branch sections have the following features: (i) they involve nonnegative concentrations and (ii) they attract the state trajectories if their initial state values are realistic (see [4, 22]). The dynamic analysis results show the effect of the dilution rate and the influent substrate concentrations on the stability and values of the mentioned equilibrium branch sections. Moreover, the dynamic analysis can also indicate the constraints on the dilution rate that allow us to avoid the biomass washout and keep the digester under normal operation [4]. As it can be noticed from [22–24], the interval of the dilution rate that leads to positive and stable features of the equilibrium biomass concentration depends on the influent substrate concentrations.Some works on dynamic analysis for two biological reaction pathways are discussed in the following. In [4], an analytical dynamic analysis is carried out, with determination of several characteristics for the normal operation equilibrium point: (i) the equilibrium state variables were established in analytical way, considering the dilution rate constrained to an interval such that biomass washout is avoided, (ii) the positive values of the equilibrium state variables is proven, and (iii) the local stability of the equilibrium is also proven. Nevertheless, the parameter constraints that allow to avoid the biomass washout, and the coordinates of the bifurcation points were not established. In [22], the equilibrium stability is analyzed in analytical way, taking into account the fact that the dynamics of the acidogenesis pathway is independent of the dynamics of the methanization pathway. This fact implies that the stability of each pathway can be analyzed separately on the basis of its differential equations. They consider the case that the normal operation equilibrium of the first biological pathway is globally asymptotically stable, and other cases of this pathway are not considered. The following characteristics are determined, for the dynamics of the second biological pahway: (i) the equilibrium points, (ii) the stability of the equilibrium points corresponding to washout and nonwashout of methanogenic bacteria, in terms of the magnitude of the dilution rate. They prove the following: (i) the normal operation equilibrium point of the acidogenesis step is stable for certain interval of the dilution rate, (ii) the normal operation equilibrium point of the methanization step is globally stable or locally stable for certain intervals of the dilution rate, and (iii) the washout equilibrium point of the methanization step is locally stable or globally stable for certain intervals of the dilution rate. Nevertheless, (i) the results are only valid for the case that the nonwashout equilibrium of the first biological pathway is globally stable, (ii) the case that the non-washout equilibrium of the first biological pathway is locally stable, is not considered, and (iii) the bifurcation points are not analyzed. In [25], one- and two-parameter bifurcation diagrams are shown. It is assumed that the acidogenesis step is faster than the methanization step, which allows considering only the dynamics of the methanization step and simplifies the stability analysis. The drawback is that the concentration of VFA at the influent is not taken into account. In addition, the bifurcation coordinates are established numerically instead of analytically. In [2], a linear state transformation is used in order to obtain a canonical state space representation. This allows to reduce the order of the dynamic system, thus simplifying the dynamic analysis, specially the definition of the eigenvalues. The following characteristics are determined: (i) the equilibrium points and its eigenvalues and (ii) the parameter conditions for the stability or unstability of each equilibrium points. Nevertheless, the conditions of the dilution rate that lead to stable nature of the normal operation equilibrium point are not explicitly expressed in terms of the influent concentrations of COD and VFA. In addition, the bifurcation points and bifurcation coexistence points are not analyzed. In [24], a dynamic analysis is developed for several bioreactor models, based on two biological reaction pathways, being the anaerobic digestion model corresponding to the so-called model IV. A linear state transformation is carried out to obtain a canonical state space representation. This allows the reduction of the dynamic order, thus simplifying the stability analysis. The bifurcation parameter coordinates are established and the corresponding bifurcation diagrams are developed, indicating the intervals on the dilution rate and the influent COD concentration that lead to locally stable nature of the normal operation equilibrium point. Nevertheless, the concentration of intermediate substrate in the influent is defined as zero, which implies a zero influent VFA concentration. In addition, the state coordinates of the bifurcation points and bifurcation coexistence points are not established. In [26], an analytical dynamic analysis is developed, leading to the following results: (i) the six possible equilibrium points including that for normal operation are determined, (ii) the hyperbolic or nonhyperbolic nature and the local stability of each equilibrium point are determined for different parameter conditions, (iii) the bifurcation points are determined, and (iv) one-parameter bifurcation diagrams are developed for the six equilibrium points. In order to compute the stability, the cascade structure of the model was used; that is, the equations for the first biological reaction pathway do not depend on the equations for the second biological pathway. The conditions for the stability of the equilibrium points are defined through of lumped parameters that are functions of the influent concentrations and the dilution rate. Nevertheless, (i) the parameter constraints that allow avoiding the biomass washout are not established in an explicit way as expressions on the dilution rate and influent concentrations, (ii) the coordinates of the bifurcation points are not established in an explicit way as functions of the dilution rate and influent concentrations, and (iii) the coordinates of the crossing between bifurcations are not established analytically.In the aforementioned works, the major drawback is the lack of an explicit and analytical definition of parameter conditions for the realistic nature of the normal operation equilibrium point (NOP). In the present work, a dynamic analysis is developed for an anaerobic digester, using a dynamical model with the assumption of two biological reaction pathways. The effect of the influent concentration of volattile fatty acids (VFA) and chemical oxygen demand (COD) is taken into account. Several characteristics of the normal operation equilibrium operation point are determined analytically: (i) the conditions for the realistic nature of the equilibrium, including parameter conditions that lead to positive values of the equilibrium state variables and locally stable nature of the equilibrium point and (ii) the coordinates of bifurcation values and bifurcation coexistence points. Those results are explicitly expressed through the dilution rate, the influent COD concentration, and influent VFA concentration. This is the main contribution with respect to closely related previous works.The rest of the paper is organized as follows. In Section2, the plant model, the equilibrium point, and the contidions to guarantee the existence of the normal operation point are found. In Section 3 the conditions that lead to locally stable nature of the NOP are presented. Section 4 shows the computation of the coordinates of local bifurcation points and the crossing between bifurcations, as well as numerical simulations for different scenarios. Finally, in Section 5 conclusions are presented.


## 2. Plant Model and Normal Operation Point

We consider the upflow anaerobic fixed bed reactor of [1]. The anaerobic digestion can be explained as follows: the biomass degrades the organic substrate to produce biogas (a mixture of CO2 and CH4) and for growth. The mass balance model is (cf. [27])
(1)
X
˙
1
=
X
1
μ
1
-
α
D
,
X
˙
2
=
X
2
μ
2
-
α
D
,
S
˙
1
=
S
1
in
-
S
1
D
-
k
1
μ
1
X
1
,
S
˙
2
=
S
2
in
-
S
2
D
+
k
2
μ
1
X
1
-
k
3
μ
2
X
2
,
(2)
μ
1
=
μ
1
max
⁡
S
1
S
1
+
K
S
1
,
μ
2
=
μ
2
max
⁡
S
2
S
2
+
K
S
2
+
S
2
/
K
I
2
,
where X
1 represents the concentration of acidogenic biomass, X
2 represents the concentration of methanogenic biomass, S
1 represents the concentration of organic substrate characterized by its chemical oxygen demand, and S
2 represents the concentration of volatile fatty acids (VFA). In addition, D represents the dilution rate and α the proportion of biomass not attached to the reactor. On current operation of the bioreactor the concentrations X
1, X
2, S
1, and S
2 are all positive, and the equilibrium point called NOP tends to (
X
1
eq
,
X
2
eq
,
S
1
eq
,
S
2
eq
). When the dilution rate D experiences an excessive increase, the washout phenomenon occurs. In this case X
1
=
0 or X
2
=
0 or both X
1
=
0 and X
2
=
0. Just the value that induces the changes of the behaviour between the NOP to washout operation point is the value for which the bifurcation occurs, and this phenomenon is called a bifurcation. In this undesired state the active biomass flows outside the reactor, such that biological degradation of the influent organic matter does not occur anymore. Thus, it is necessary to inoculate the biomass again, what may take several months. Nomenclature section shows the parameters used in this work. We refer the reader to [1, 3, 27] for further details on the process.The goal of the dynamic analysis is to establish the following characteristics of the NOP: (i) the parameter conditions that lead to positive values of the equilibrium state variables, (ii) the parameter conditions that lead to locally stable nature of the equilibrium point, and (iii) the coordinates of bifurcation points and bifurcation coexistence points. Those results are explicitly expressed in terms of the dilution rate (D) and influent concentrations of volatile fatty acids (VFA) and chemical oxygen demand (COD).The dynamical system (1) has six equilibrium points, but only one of them, the so-called normal operation point (NOP), describes the normal operation of the digester (see [23, 26]). The other ones play an important role in the system dynamics and their interactions explain the different behaviors in the bioreactor and its transitions from stable operation to washout phenomena. The equilibrium condition for the NOP is found by imposing a null the vector field. Considering X
1
≠
0 and X
2
≠
0, we have
(3)
μ
1
eq
-
α
D
=
0
,
(4)
μ
2
eq
-
α
D
=
0
,
(5)
(
S
1
in
-
S
1
eq
)
D
-
k
1
μ
1
eq
X
1
eq
=
0
,
(6)
(
S
2
in
-
S
2
eq
)
D
+
k
2
μ
1
eq
X
1
eq
-
k
3
μ
2
eq
X
2
eq
=
0
,
where
(7)
μ
1
eq
=
μ
1
max
⁡
S
1
eq
K
S
1
+
S
1
eq
,
μ
2
eq
=
μ
2
max
⁡
S
2
eq
K
S
2
+
S
2
eq
+
S
2
eq
/
K
I
2
.
The NOP has coordinates:
(8)
X
1
eq
=
S
1
in
-
S
1
eq
k
1
α
,
(9)
X
2
eq
=
S
2
in
-
S
2
eq
+
(
k
2
/
k
1
)
(
S
1
in
-
S
1
eq
)
α
k
3
,
(10)
S
1
eq
=
α
D
K
S
1
μ
1
max
⁡
-
α
D
,
(11)
S
2
eq
=
-
b
-
b
2
-
4
a
c
2
a
,
where
(12)
a
=
α
D
,
b
=
K
I
2
(
α
D
-
μ
2
max
⁡
)
,
c
=
K
I
2
K
S
2
α
D
.
Normal operation implies that the reactor concentrations are positive and D is positive, so that the combination of X
1
eq
>
0, X
2
eq
>
0, S
1
eq
>
0, S
2
eq
>
0, D
>
0, and (8)–(11) yields the restrictions on D and S
1
in. Since S
2
eq is a positive real number, it follows that -
b
>
0 and, from (11), that the condition given by inequation (13) has to be fulfilled:
(13)
b
2
-
4
a
c
≥
0
.The following proposition presents some conditions must be satisfied byD in order the NOP exists.Proposition 1.
Consider the definitions given in (12), conditions established in (13) and the equilibrium point provided by (8)–(11). The dilution rate D and the equilibrium concentration S
2
e
q satisfy
(14)
i
D
(
S
2
e
q
)
=
μ
2
max
⁡
α
S
2
e
q
K
S
2
+
S
2
e
q
+
S
2
e
q
/
K
I
2
(15)
(
ii
)
D
(
S
2
e
q
)
e
x
h
i
b
i
t
a
m
a
x
i
m
u
m
v
a
l
u
e
w
h
e
n
S
2
e
q
=
K
I
K
S
2
i
n
c
o
n
s
e
q
u
e
n
c
e
D
≤
μ
2
max
⁡
α
K
I
K
I
+
2
K
S
2
(16)
(
iii
)
F
o
r
N
O
P
,
S
2
e
q
:
=
-
b
-
b
2
-
4
a
c
2
a
∈
[
0
,
K
I
K
S
2
]
.Proof.
Replacing (7) in (4) it follows that
(17)
D
=
μ
2
max
⁡
α
S
2
eq
K
S
2
+
S
2
eq
+
S
2
eq
/
K
I
2
.
The expression above indicates that D has a maximum value for a certain value of S
2
eq. Differentiating (17) yields
(18)
∂
D
∂
S
2
eq
=
μ
2
max
⁡
α
K
S
2
-
S
2
eq
/
K
I
2
S
2
eq
+
K
S
2
+
S
2
eq
/
K
I
2
2
and making this expression zero yields
(19)
S
2
eq
=
K
I
K
S
2
.
Substituting it into (2) yields
(20)
D
max
⁡
=
μ
2
max
⁡
α
K
I
K
I
+
2
K
S
2
.
In addition inequation (13) is only fulfilled if D
<
D
max
⁡; then, S
2
eq has a maximum when b
2
-
4
a
c
=
0 which is given by (21)(21)
S
2
max
⁡
eq
=
-
b
2
a
.
Replacing (20) in (21), we obtain that the maximum value of S
2
eq for the NOP fulfills
(22)
S
2
max
⁡
eq
=
K
I
K
S
2
.Other conditions can be computed from (8)–(11). From (8) it is easy to see that the condition for the existence of NOP is S
1
eq
<
S
1
in (on contrary case, X
1
eq
≤
0); this expression combined with (10) implies
(23)
D
<
μ
max
⁡
1
α
S
1
in
K
S
1
+
S
1
in
.
Therefore, if D satisfies (23) then X
1
eq
>
0 and S
1
eq
<
S
1
in. In a similar way, with a bit more complicated algebraic steps it is possible to find conditions on S
2
eq and D such that the existence of NOP can be proven. Equation (9) yields S
2
eq
<
S
2
in. By combining this expression with (9), (10), and (11), we find conditions on D such that X
2
eq
>
0. These conditions are
(24)
D
<
min
⁡
μ
2
max
⁡
α
K
I
K
I
+
2
K
S
2
,
D
~
,
where D
~ satisfies (25),
(25)
K
S
1
α
D
~
μ
1
max
⁡
-
α
D
~
+
k
1
k
2
-
b
~
-
b
~
2
-
4
a
~
c
~
2
a
~
=
S
1
in
+
k
1
k
2
S
2
in
with a
~
=
α
D
~, b
~
=
K
I
2
(
α
D
~
-
μ
2
max
⁡
), c
~
=
K
I
2
K
S
2
α
D
~.According to the results found above, the final condition to guarantee the existence of NOP is(26)
D
<
min
⁡
μ
1
max
⁡
α
S
1
in
K
S
1
+
S
1
in
,
μ
2
max
⁡
α
K
I
K
I
+
2
K
S
2
,
D
~
,
where D
~ satisfies (25).


## 3. Stability Analysis of the Normal Operation Point (NOP)

This section shows the analytical conditions that lead to locally stable nature of the NOP (8)–(11). The conditions are explicitly and analytically expressed as functions of the dilution rate and influent concentrations. The stability of the normal equilibrium point is determined by indirect Lyapunov method [28–30]. The Jacobian of the system (1)-(2) is given by
(27)
J
=
J
11
J
12
J
13
J
14
J
21
J
22
J
23
J
24
J
31
J
32
J
33
J
34
J
41
J
42
J
43
J
44
.
When the Jacobian matrix is evaluated at the NOP (i.e., (8)–(11)) the nonzero J
i
j are
(28)
J
11
=
μ
1
eq
-
α
D
,
(29)
J
13
=
μ
1
′
X
1
eq
,
(30)
J
22
=
μ
2
eq
-
α
D
,
(31)
J
24
=
μ
2
′
X
2
eq
,
(32)
J
31
=
-
k
1
μ
1
eq
,
(33)
J
33
=
-
D
-
k
1
μ
1
′
X
1
eq
,
(34)
J
41
=
k
2
μ
1
eq
,
(35)
J
42
=
-
k
3
μ
2
eq
,
(36)
J
43
=
k
2
μ
1
′
X
1
eq
,
(37)
J
44
=
-
D
-
k
3
μ
2
′
X
2
eq
being
(38)
μ
1
′
=
μ
1
max
⁡
K
S
1
K
S
1
+
S
1
eq
2
,
(39)
μ
2
′
=
μ
2
max
⁡
K
S
2
-
S
2
eq
/
K
I
2
K
S
2
+
S
2
eq
+
S
2
eq
/
K
I
2
2
.
Note that μ
1
′
>
0 and μ
2
′
>
0 for S
2
eq
∈
[
0
,
K
I
K
S
2
] (see (16)). Similarly, J
11
=
0 by (3) and J
22
=
0 by (4). The eigenvalues of the linearized system are (see [4, 23])
(40)
λ
1
=
t
1
2
+
t
1
2
2
-
d
1
,
(41)
λ
2
=
t
1
2
-
t
1
2
2
-
d
1
,
(42)
λ
3
=
t
2
2
+
t
2
2
2
-
d
2
,
(43)
λ
4
=
t
2
2
-
t
2
2
2
-
d
2
,
where
(44)
t
1
=
J
11
+
J
33
=
J
33
,
(45)
t
2
=
J
22
+
J
44
=
J
44
,
(46)
d
1
=
J
11
J
33
-
J
13
J
31
=
-
J
13
J
31
,
(47)
d
2
=
J
22
J
44
-
J
24
J
42
=
-
J
24
J
42
.
To establish the conditions that lead to stable nature of the NOP, the eigenvalues (40)–(43) and and equilibrium conditions (3)–(6) are used. Recall that an equilibrium is locally stable if the real part of all its eigenvalues are negative [28, pp. 55]. From (40) and (41) it follows that if t
1
<
0 and d
1
>
0, then Re
(
λ
1
)
<
0 and Re
(
λ
2
)
<
0, whereas (42) and (43) indicate that if t
2
<
0 and d
2
>
0, then Re
(
λ
3
)
<
0 and Re
(
λ
4
)
<
0. We will establish the range for D that assures these eigenvalue conditions, considering each term independently.Proposition 2.
Consider the definitions given in (12), the equilibrium point provided by (8)–(11), eigenvalues (40) and (41), and the term μ
1
′ defined in (38). If D fulfills
(48)
0
<
D
<
μ
1
max
⁡
α
S
1
i
n
K
S
1
+
S
1
i
n
,
then t
1
<
0 and d
1
>
0 and consequently Re
(
λ
1
)
<
0 and Re
(
λ
2
)
<
0.Proof.
Combining (44) and (33), we obtain
(49)
t
1
=
-
D
-
k
1
μ
1
′
X
1
eq
.
Now, combining (46), (29), and (32), we obtain
(50)
d
1
=
μ
1
′
X
1
eq
k
1
μ
1
eq
.
As D
>
0, μ
1
′
>
0, and X
1
eq
>
0 if D
<
(
μ
1
max
⁡
/
α
)
(
S
1
in
/
(
K
S
1
+
S
1
in
)
), then it follows from (40) and (41) that Re
(
λ
1
)
<
0 and Re
(
λ
2
)
<
0.Proposition 3.
Consider the definitions given in (12), the equilibrium point provided by (8)–(11), and the term μ
2
′ defined in (39). If
(51)
0
<
D
<
min
⁡
μ
2
max
⁡
α
K
I
K
I
+
2
K
S
2
,
D
~
,
where D
~ satisfies (25), then t
2
<
0 and d
2
>
0 and consequently Re
(
λ
3
)
<
0 and Re
(
λ
4
)
<
0.Proof.
Combining (45) and (37), we obtain
(52)
t
2
=
-
D
-
k
3
μ
2
′
X
2
eq
.
Now, combining (47), (31), and (35), we obtain
(53)
d
2
=
μ
2
′
X
2
eq
k
3
μ
2
eq
.
As 0
<
S
2
eq
<
K
I
/
K
S
2 according to (16), μ
2
′
>
0 if 0
<
S
2
eq
<
K
I
/
K
S
2, D
<
(
μ
2
max
⁡
/
α
)
(
K
I
/
(
K
I
+
2
K
S
2
)
) according to (15), and X
2
eq
>
0 if conditions given by (24) and (25) are satisfied; then, it follows from (42) and (43) that Re
(
λ
3
)
<
0 and Re
(
λ
4
)
<
0.According to the results obtained above, the final condition to NOP is stable is given by(54)
D
<
min
⁡
μ
1
max
⁡
α
S
1
in
K
S
1
+
S
1
in
,
μ
2
max
⁡
α
K
I
K
I
+
2
K
S
2
,
D
~
,
where D
~ satisfies (25).Remark 4.
The final condition to guarantee the NOP stability is the same as the final condition to assure the existence of NOP.Remark 5.
The final condition to guarantee the NOP stability indicates that (i) ifS
1
in is overly small, then the dilution rate D is constrained to small values even if S
2
in is large and (ii) if S
2
in is overly small, then the dilution rate D is constrained to small values even if S
1
in is large. If S
1
in and S
2
in satisfy
(55)
μ
1
max
⁡
α
S
1
in
S
1
in
+
K
S
1
>
μ
2
max
⁡
α
K
I
K
I
+
2
K
S
2
then the dilution rate D is only limited by the constant value (
μ
2
max
⁡
/
α
)
(
K
I
/
(
K
I
+
2
K
S
2
)
). Therefore, it is desirable to have large values of S
1
in, and S
2
in, as it would allow to have larger values of the dilution rate D.Remark 6.
The final condition to guarantee the NOP stability indicates that there are three ways whereby the stable nature of the NOP can be lost: (i) with small values ofS
1
in, (ii) with small values of S
2
in, (iii) with large values of D.


## 4. Determination of Bifurcation Points and Their Crossing

As we stated before, there are six equilibria in the system and the evolution of these equilibria determines the bifurcations and the behavior of the system when one or more parameters are varied. In Table1 the values of the equilibrium points when bifurcations occur are shown. The interaction of these equilibria raise three different bifurcations which are called TB(X
1), FB(X
2), and TB(X
2). The bifurcation TB(X
1) is of transcritical type and involves the intersection between the equilibrium corresponding to NOP and the equilibrium corresponding to inactivated acidogenic biomass (X
1
=
0). The bifurcation FB(X
2) is of fold type and involves collision between the equilibrium corresponding to NOP and one equilibrium that does not have physical meaning. The bifurcation TB(X
2) is of transcritical type and involves the intersection between the equilibrium corresponding to NOP and the one corresponding to inactivated methanogenic biomass (X
2
=
0). The mentioned bifurcations have common points where the parameters and state variables have the same values (see Table 2).Table 1
Equilibrium points.
Bifurcation
Equilibrium points:a, b, and c were defined in (12).
TB(X
1)
X
1
eq = 0
X
2
eq = 1
α
k
3
(
S
2
in
-
-
b
-
b
2
-
4
a
c
2
a
)
S
1
eq = S
1
in
S
2
eq = -
b
-
b
2
-
4
a
c
2
a
Where
S
1
in = α
D
K
S
1
μ
1
max
⁡
-
α
D, and D
=
μ
1
max
⁡
S
1
in
α
(
S
1
in
+
K
S
1
)
FB(X
2)
X
1
eq = 1
α
k
1
(
S
1
in
-
μ
2
max
⁡
K
I
K
S
1
μ
1
max
⁡
(
2
K
S
2
+
K
I
)
-
μ
2
max
⁡
K
I
)
X
2
eq = 1
α
k
3
(
S
2
in
+
k
2
k
1
S
1
in
-
K
I
K
S
2
-
k
2
k
1
μ
2
max
⁡
K
I
K
S
1
μ
1
max
⁡
(
2
K
S
2
+
K
I
)
-
μ
2
max
⁡
K
I
)
S
1
eq = μ
2
max
⁡
K
I
K
S
1
μ
1
max
⁡
(
2
K
S
2
+
K
I
)
-
μ
2
max
⁡
K
I
S
2
eq = K
I
K
S
2
Where
D = μ
2
max
⁡
K
I
α
(
2
K
S
2
+
K
I
)
TB(X
2)
X
1
eq = S
1
in
μ
1
max
⁡
-
α
D
(
S
1
in
+
K
S
1
)
α
k
1
(
μ
1
max
⁡
-
α
D
)
X
2
eq = 0
S
1
eq = α
D
K
S
1
μ
1
max
⁡
-
α
D
S
2
eq = -
b
-
b
2
-
4
a
c
2
a
Where
S
1
in = α
D
K
S
1
μ
1
max
⁡
-
α
D
+
k
1
k
2
(
-
b
-
b
2
-
4
a
c
2
a
-
S
2
in
)Table 2
Intersections of the bifurcation branches.
Bifurcations
Common values of the equilibrium points
TB(X
1) and FB(X
2)
X
1
eq = 0
X
2
eq = S
2
in
-
K
I
K
S
2
α
k
3
S
1
eq = μ
2
max
⁡
K
I
K
S
1
K
I
(
μ
1
max
⁡
-
μ
2
max
⁡
)
+
2
μ
1
max
⁡
K
S
2
S
2
eq = K
I
K
S
2
Where
S
1
in = μ
2
max
⁡
K
I
K
S
1
K
I
(
μ
1
max
⁡
-
μ
2
max
⁡
)
+
2
μ
1
max
⁡
K
S
2
FB(X
2) and TB(X
2)
X
1
eq = 1
k
2
α
(
K
I
K
S
2
-
S
2
in
)
X
2
eq = 0
S
1
eq = μ
2
max
⁡
K
I
K
S
1
K
I
(
μ
1
max
⁡
-
μ
2
max
⁡
)
+
2
μ
1
max
⁡
K
S
2
S
2
eq = K
I
K
S
2
Where
D = μ
2
max
⁡
K
I
α
(
2
K
S
2
+
K
I
) and
S
1
in = μ
2
max
⁡
K
I
K
S
1
K
I
(
μ
1
max
⁡
-
μ
2
max
⁡
)
+
2
μ
1
max
⁡
K
S
2
+
k
1
k
2
(
K
I
K
S
2
-
S
2
in
)In this section, the determination of the bifurcation points as well as crossing points between bifurcations TB(X
1) and FB(X
2) and bifurcations FB(X
2) and TB(X
2) is presented.


### 4.1. Determination of Bifurcation Points

One-parameter bifurcation diagrams can be obtained from Table1. This results in a curve called equilibrium branch. The crossing, collision, or splitting of these branches corresponds to a local bifurcation, either fold, transcritical, or pitchfork. The local bifurcations involve the vanishing of the real part of at least one eigenvalue, which allows determining the bifurcation coordinates [29].Now, we establish local bifurcations of the NOP, considering the dilution rateD, the feed concentration S
1
in, and the feed concentration S
2
in as bifurcation parameters. The real part of at least one eigenvalue becomes zero when a local bifurcation occurs, which allows us to determine the bifurcation coordinates. From (41) and (43) it follows that λ
2
=
0 if d
1
=
0 and λ
4
=
0 if d
2
=
0. These two conditions are analyzed in the following. We begin by combining d
1
=
0 with (46):
(56)
-
J
13
J
31
=
0
,
α
k
1
μ
1
′
X
1
eq
=
0
and d
2
=
0 with (47):
(57)
-
J
24
J
42
=
0
,
α
k
3
μ
2
′
X
2
eq
=
0
.
Consequently, we have three possibilities:
(58)
(
1
)
X
1
eq
=
0
,
(59)
(
2
)
μ
2
′
=
0
,
(60)
(
3
)
X
2
eq
=
0
.
Combining the above conditions with the conditions for the equilibrium (3)–(6), we obtain the variable states and parameter values for each bifurcation.The bifurcation type, either fold or transcritical, can be determined on the basis of the behavior of the equilibrium branches. The fold bifurcation involves the collision and disappearance of two equilibrium branches. In some cases, it is related to the presence of a square root in the equilibrium coordinates, and the bifurcation coordinates can be determined by equating the square root to zero. Contrarily, the transcritical bifurcation involves the crossing of two equilibrium branches. Its coordinates can be determined by analyzing the intersection of the equilibrium branches (cf. [31, pp. 525]).Notice that the equilibrium (8)–(11) involves a square root in (11), which indicates the presence of a fold bifurcation. Imposing a null value for the square root, it yields
(61)
D
=
μ
2
max
⁡
α
K
I
K
I
+
2
K
S
2
or
D
=
μ
2
max
⁡
α
K
I
K
I
-
2
K
S
2
but only the first possibility satisfies condition (15). The graphical analysis of the bifurcations, and the fact that the parameter value (61) leads to a fold bifurcation, will be used to determine the bifurcation type.Bifurcation Branch TB(
X
1
). Combining (58) with (3)–(6), yields the coordinates of bifurcation TB(X
1):
(62)
X
1
eq
=
0
,
(63)
X
2
eq
=
1
α
k
3
S
2
in
-
-
b
-
b
2
-
4
a
c
2
a
,
(64)
S
1
eq
=
S
1
in
,
(65)
S
2
eq
=
-
b
-
b
2
-
4
a
c
2
a
,
(66)
S
1
in
=
α
D
K
S
1
μ
1
max
⁡
-
α
D
,
(67)
D
=
μ
1
max
⁡
S
1
in
α
(
S
1
in
+
K
S
1
)
,
where a, b, c are defined in (12).Bifurcation Branch FB(
X
2
). From (59) and (39) it follows that
(68)
K
S
2
-
S
2
eq
K
I
2
=
0
(69)
⟹
S
2
eq
=
K
I
K
S
2
.
From (4) it follows that D
=
(
1
/
α
)
μ
2
eq. Using (2) and (69), we obtain
(70)
D
=
μ
2
max
⁡
α
K
I
2
K
S
2
+
K
I
.
Notice that this value is independent of S
1
in, S
2
in, and is the same as the upper bound appearing in (15). Combining (69), (70), (3)–(6), and (59) yields the coordinates of bifurcation FB(X
2):
(71)
X
1
eq
=
1
α
k
1
S
1
in
-
μ
2
max
⁡
K
I
K
S
1
μ
1
max
⁡
(
2
K
S
2
+
K
I
)
-
μ
2
max
⁡
K
I
,
(72)
X
2
eq
=
1
α
k
3
k
2
k
1
μ
2
max
⁡
K
I
K
S
1
μ
1
max
⁡
2
K
S
2
+
K
I
-
μ
2
max
⁡
K
I
S
2
in
+
k
2
k
1
S
1
in
-
K
I
K
S
2
h
h
h
h
h
-
k
2
k
1
μ
2
max
⁡
K
I
K
S
1
μ
1
max
⁡
2
K
S
2
+
K
I
-
μ
2
max
⁡
K
I
,
(73)
S
1
eq
=
μ
2
max
⁡
K
I
K
S
1
μ
1
max
⁡
(
2
K
S
2
+
K
I
)
-
μ
2
max
⁡
K
I
,
(74)
S
2
eq
=
K
I
K
S
2
,
(75)
D
=
μ
2
max
⁡
K
I
α
(
2
K
S
2
+
K
I
)
.Bifurcation Branch TB(
X
2
). Combining (3)–(6) with (60), we obtain
(76)
X
1
eq
=
S
1
in
μ
1
max
⁡
-
α
D
(
S
1
in
+
K
S
1
)
α
k
1
(
μ
1
max
⁡
-
α
D
)
,
(77)
X
2
eq
=
0
,
(78)
S
1
eq
=
α
D
K
S
1
μ
1
max
⁡
-
α
D
,
(79)
S
2
eq
=
-
b
-
b
2
-
4
a
c
2
a
,
(80)
S
1
in
=
α
D
K
S
1
μ
1
max
⁡
-
α
D
+
k
1
k
2
-
b
-
b
2
-
4
a
c
2
a
-
S
2
in
,
where a, b, and c are defined in (12).Remark 7.
Notice that the coordinates of the bifurcation parameters, provided by (67), (75), and (80), are the same as the upper bounds that guarantee stability of NOP, established in (54). The reason is that the bifurcation points bound the locally stable nature of the NOP, and the parameter coordinates for the bifurcation points are parameter bounds that allow to keep the system under normal operation.Bifurcation Diagrams. Figure 1 shows four diagrams: a 2-parameter bifurcation diagram (in the upper left position) and three 1-parameter bifurcation diagrams. For the 2-parameter bifurcation diagram D and S
1
in are the bifurcation parameters. The vertical line in this subfigure is the bifurcation branch FB(X
2); it is a fold bifurcation branch. The solid curved line corresponds to the bifurcation branch TB(X
1); it is a transcritical bifurcation branch. The dash-dot curved line (very close to FB(X
2)) corresponds to the bifurcation branch TB(X
2); it is a transcritical bifurcation branch. The region D
>
(
μ
2
max
⁡
/
α
)
K
I
/
(
K
I
+
2
K
S
2
) located at the right of the vertical solid line implies the non-existence of the NOP. The 1-parameter bifurcation diagrams are computed using D as bifurcation parameter and fixing S
1
in according to the value shown by the arrows in the 2-parameter diagram.Bifurcation diagrams for the NOP. The 2-parameter bifurcation diagram is in upper left position. The corresponding 1-parameter bifurcation diagrams for the specific values ofS
in according to arrows A
1, A
2, and A
3 are in the other positions.
(a)
(b)
(c)
(d)From Proposition1 and (15) and (75), it follows that the value D
=
(
μ
2
max
⁡
/
α
)
K
I
/
(
K
I
+
2
K
S
2
) is not only the parameter coordinate for the bifurcation branch FB(X
2), but also the limit for the existence of the NOP. The arrows displayed in Figure 1(a) and noted as A
1, A
2, and A
3 allow us to define if the bifurcation is either transcritical or fold: (i) arrow A
1 indicates that bifurcation TB(X
1) is transcritical, bifurcation TB(X
2) is transcritical and bifurcation FB(X
2) is fold, (ii) arrow A
3 confirms that bifurcation FB(X
2) is fold and bifurcation TB(X
1) is transcritical, and (iii) arrow A
2 shows the crossing of transcritical and fold bifurcations.Thus,(i)
TB(X
1)—described by (62) to (67) with X
1
eq
=
0—is a transcritical bifurcation that involves the intersection between the equilibrium corresponding to normal operation and the equilibrium corresponding to X
1
eq
=
0. For larger values of D the concentration X
1 would converge to zero, which is known as washout of acidogenic biomass.
(ii)
Bifurcation FB(X
2)—described by (71) to (75)—corresponds to a fold bifurcation where the equilibrium branch corresponding to NOP collides with the other equilibrium branch. The collision and disappearance can only be noticed in the branches of X
2
eq and S
2
eq. For larger values of D, the concentration X
2 would converge to zero, which is known as washout of methanogenic biomass.
(iii)
Bifurcation TB(X
2)—described by (76) to (80) with X
2
eq
=
0—corresponds to a transcritical bifurcation that involves the intersection between the equilibrium corresponding to NOP and the equilibrium corresponding to X
2
eq
=
0. For larger values of D the concentration X
2 would converge to zero, which is known as washout of methanogenic biomass.Figure2 is a zoomed version close to TB(X
2) of Figure 1. This figure has the following features: (i) the region for local stability of the NOP is denoted by NO, (ii) the region for local stability of an equilibrium corresponding to acidogenic biomass washout is denoted by W/X
1, (iii) the region for local stability of an equilibrium corresponding to methanogenic biomass washout is denoted by W/X
2, (iv) the region for local stability of an equilibrium corresponding to washout of both acidogenic and methanogenic biomass is denoted by W/(
X
1
&
X
2
), (v) the region for the fold bifurcation FB(X
2) corresponds to the vertical line, (vi) the regions for the transcritical bifurcations TB(X
1) and TB(X
2) (continuous and dash-doted lines, resp.) are curved lines, (vii) the two big points correspond to coexisting bifurcations: the upper one corresponds to the coexistence of FB(X
2) and TB(X
1) whereas the lower one corresponds to the coexistence of FB(X
2) and TB(X
2), and (viii) in the region D
>
(
μ
2
max
⁡
/
α
)
K
I
/
(
K
I
+
2
K
S
2
), located at the right of the vertical solid line, the NOP does not exist.Figure 2
Two-parameter bifurcation diagram.From Figure1 it can be observed that if S
2
in is fixed there are two ways whereby the local stability of the normal operation condition can be lost: (i) by increasing the value of D, (ii) by decreasing the value of S
1
in, which agrees with Remark 6. In the same figure there are three cases (A
3, A
2, and A
1) in which the occurrence of local bifurcation is analyzed. For some values of D
>
0 and S
1
in there is a locally stable equilibrium point with X
1
eq
>
0, X
2
eq
>
0, S
1
eq
>
0 and S
2
eq
>
0, which is associated with the NOP. For a given value of S
1
in, there is a value of the dilution rate, D
=
D
*, for which a local bifurcation occurs, related to a change of stability (see Figure 1). Further increase of D, satisfying D
>
D
*, leads to a locally stable equilibrium point that involves either X
1
eq
=
0 or X
2
eq
=
0 or both, while S
1
eq
>
0 and S
2
eq
>
0. These equilibria are associated with washout.Figures3, 4, and 5 show the verification of the two-parameter bifurcation diagram of Figure 1. In each case, NOP is chosen by selecting suitable values of S
1
in and D. In order to obtain the diagram, D is increased so that a bifurcation occurs. The difference between the three cases is the characteristics of the first bifurcation attained. For the case A
3 the bifurcation is related to the washout of the methanogenic biomass (X
2). For the case A
2 the bifurcation is related to the washout of both the acidogenic and the methenogenic biomasses (X
1 and X
2), and for the case A
1 the bifurcation is related to the washout of the acidogenic biomass (X
1).Bifurcation diagrams and transient responses forX
1 and X
2 states. Case A
3.
(a)
(b)
(c)
(d)Bifurcation diagrams and transient responses forX
1 and X
2 states. Case A
2.
(a)
(b)
(c)
(d)Bifurcation diagrams and transient responses forX
1 and X
2 states. Case A
1.
(a)
(b)
(c)
(d)In Figures3, 4, and 5, the upper left part shows the bifurcation diagram for the chosen value of S
1
in, with D as bifurcation parameter, and indicates the considered case, either case A
3, A
2, or A
1. The figures surrounding the bifurcation diagram show the transient evolution of the state variables X
1 and X
2. Figure 3 shows the results for the case A
3 with S
1
in
=
6.5
gCOD
/
L and D
∈
{
0.8,1.1,1.3
} days−1 which correspond to NOP, washout operation for X
2, and washout operation for both X
1 and X
2, respectively. Figure 4 shows the results for the case A
2, with S
1
in
=
5.73
gCOD
/
L, and D
∈
{
0.7,1.2
} days−1 which correspond to NOP, and washout operation for both X
1 and X
2, respectively. Figure 5 shows the results for the case A
1, with S
1
in
=
0.5
gCOD
/
L, and D
∈
{
0.1,0.8,1.07,1.2
} days−1 which correspond to NOP for the first value, washout operation for X
2, for the second value, and washout operation for both X
1 and X
2, for third and fourth values.


### 4.2. Determination of Crossing Bifurcation Branches

Now, we analyze the crossing between bifurcations TB(X
1), FB(X
2), and TB(X
2). At the collision points, the bifurcation parameters and bifurcation variable states have the same values. There are two ways whereby the crossing points can be determined: (i) by combining the bifurcation coordinates (62),…,(67);  (71),…,(75);  (76),…,(80), (ii) by combining (58), (59), and (60) with equilibrium conditions (3)–(6).In order to determine the crossing between TB(X
1) and FB(X
2), the bifurcation coordinates (62),…,(67) and (71),…,(75) are combined, so that
(81)
X
1
eq
=
0
,
X
2
eq
=
S
2
in
-
K
I
K
S
2
α
k
3
,
S
1
eq
=
μ
2
max
⁡
K
I
K
S
1
K
I
(
μ
1
max
⁡
-
μ
2
max
⁡
)
+
2
μ
1
max
⁡
K
S
2
,
S
2
eq
=
K
I
K
S
2
,
D
=
μ
2
max
⁡
K
I
α
(
2
K
S
2
+
K
I
)
,
S
1
in
=
μ
2
max
⁡
K
I
K
S
1
K
I
(
μ
1
max
⁡
-
μ
2
max
⁡
)
+
2
μ
1
max
⁡
K
S
2
.Remark 8.
Notice that (i) the parameter coordinatesD and S
1
in are constant, (ii) the equilibrium values of the state variables X
1
eq, S
1
eq, and S
2
eq do not depend on D nor S
1
in nor S
2
in and (iii) the state variable X
2
eq depends on S
2
in but does not depend on D nor S
1
in.To determine the crossing between FB(X
2) and TB(X
2), the bifurcation coordinates (71),…,(75) and (76),…, (80) are combined, so that
(82)
X
1
eq
=
1
k
2
α
K
I
K
S
2
-
S
2
in
,
X
2
eq
=
0
,
S
1
eq
=
μ
2
max
⁡
K
I
K
S
1
K
I
(
μ
1
max
⁡
-
μ
2
max
⁡
)
+
2
μ
1
max
⁡
K
S
2
,
S
2
eq
=
K
I
K
S
2
,
D
=
μ
2
max
⁡
K
I
α
(
2
K
S
2
+
K
I
)
,
S
1
in
=
μ
2
max
⁡
K
I
K
S
1
K
I
μ
1
max
⁡
-
μ
2
max
⁡
+
2
μ
1
max
⁡
K
S
2
+
k
1
k
2
K
I
K
S
2
-
S
2
in
.Remark 9.
Notice that (i) the parameterD is constant, (ii) the parameter S
1
in depends on S
2
in but it does not depend on D, (iii) the state variables X
2
eq, S
1
eq, S
2
eq are constant, and (iv) the state variable X
1
eq depends on S
2
in but does not depend on D nor S
1
in.Figure2 shows the intersection between bifurcations by means of large black points.


## 4.1. Determination of Bifurcation Points

One-parameter bifurcation diagrams can be obtained from Table1. This results in a curve called equilibrium branch. The crossing, collision, or splitting of these branches corresponds to a local bifurcation, either fold, transcritical, or pitchfork. The local bifurcations involve the vanishing of the real part of at least one eigenvalue, which allows determining the bifurcation coordinates [29].Now, we establish local bifurcations of the NOP, considering the dilution rateD, the feed concentration S
1
in, and the feed concentration S
2
in as bifurcation parameters. The real part of at least one eigenvalue becomes zero when a local bifurcation occurs, which allows us to determine the bifurcation coordinates. From (41) and (43) it follows that λ
2
=
0 if d
1
=
0 and λ
4
=
0 if d
2
=
0. These two conditions are analyzed in the following. We begin by combining d
1
=
0 with (46):
(56)
-
J
13
J
31
=
0
,
α
k
1
μ
1
′
X
1
eq
=
0
and d
2
=
0 with (47):
(57)
-
J
24
J
42
=
0
,
α
k
3
μ
2
′
X
2
eq
=
0
.
Consequently, we have three possibilities:
(58)
(
1
)
X
1
eq
=
0
,
(59)
(
2
)
μ
2
′
=
0
,
(60)
(
3
)
X
2
eq
=
0
.
Combining the above conditions with the conditions for the equilibrium (3)–(6), we obtain the variable states and parameter values for each bifurcation.The bifurcation type, either fold or transcritical, can be determined on the basis of the behavior of the equilibrium branches. The fold bifurcation involves the collision and disappearance of two equilibrium branches. In some cases, it is related to the presence of a square root in the equilibrium coordinates, and the bifurcation coordinates can be determined by equating the square root to zero. Contrarily, the transcritical bifurcation involves the crossing of two equilibrium branches. Its coordinates can be determined by analyzing the intersection of the equilibrium branches (cf. [31, pp. 525]).Notice that the equilibrium (8)–(11) involves a square root in (11), which indicates the presence of a fold bifurcation. Imposing a null value for the square root, it yields
(61)
D
=
μ
2
max
⁡
α
K
I
K
I
+
2
K
S
2
or
D
=
μ
2
max
⁡
α
K
I
K
I
-
2
K
S
2
but only the first possibility satisfies condition (15). The graphical analysis of the bifurcations, and the fact that the parameter value (61) leads to a fold bifurcation, will be used to determine the bifurcation type.Bifurcation Branch TB(
X
1
). Combining (58) with (3)–(6), yields the coordinates of bifurcation TB(X
1):
(62)
X
1
eq
=
0
,
(63)
X
2
eq
=
1
α
k
3
S
2
in
-
-
b
-
b
2
-
4
a
c
2
a
,
(64)
S
1
eq
=
S
1
in
,
(65)
S
2
eq
=
-
b
-
b
2
-
4
a
c
2
a
,
(66)
S
1
in
=
α
D
K
S
1
μ
1
max
⁡
-
α
D
,
(67)
D
=
μ
1
max
⁡
S
1
in
α
(
S
1
in
+
K
S
1
)
,
where a, b, c are defined in (12).Bifurcation Branch FB(
X
2
). From (59) and (39) it follows that
(68)
K
S
2
-
S
2
eq
K
I
2
=
0
(69)
⟹
S
2
eq
=
K
I
K
S
2
.
From (4) it follows that D
=
(
1
/
α
)
μ
2
eq. Using (2) and (69), we obtain
(70)
D
=
μ
2
max
⁡
α
K
I
2
K
S
2
+
K
I
.
Notice that this value is independent of S
1
in, S
2
in, and is the same as the upper bound appearing in (15). Combining (69), (70), (3)–(6), and (59) yields the coordinates of bifurcation FB(X
2):
(71)
X
1
eq
=
1
α
k
1
S
1
in
-
μ
2
max
⁡
K
I
K
S
1
μ
1
max
⁡
(
2
K
S
2
+
K
I
)
-
μ
2
max
⁡
K
I
,
(72)
X
2
eq
=
1
α
k
3
k
2
k
1
μ
2
max
⁡
K
I
K
S
1
μ
1
max
⁡
2
K
S
2
+
K
I
-
μ
2
max
⁡
K
I
S
2
in
+
k
2
k
1
S
1
in
-
K
I
K
S
2
h
h
h
h
h
-
k
2
k
1
μ
2
max
⁡
K
I
K
S
1
μ
1
max
⁡
2
K
S
2
+
K
I
-
μ
2
max
⁡
K
I
,
(73)
S
1
eq
=
μ
2
max
⁡
K
I
K
S
1
μ
1
max
⁡
(
2
K
S
2
+
K
I
)
-
μ
2
max
⁡
K
I
,
(74)
S
2
eq
=
K
I
K
S
2
,
(75)
D
=
μ
2
max
⁡
K
I
α
(
2
K
S
2
+
K
I
)
.Bifurcation Branch TB(
X
2
). Combining (3)–(6) with (60), we obtain
(76)
X
1
eq
=
S
1
in
μ
1
max
⁡
-
α
D
(
S
1
in
+
K
S
1
)
α
k
1
(
μ
1
max
⁡
-
α
D
)
,
(77)
X
2
eq
=
0
,
(78)
S
1
eq
=
α
D
K
S
1
μ
1
max
⁡
-
α
D
,
(79)
S
2
eq
=
-
b
-
b
2
-
4
a
c
2
a
,
(80)
S
1
in
=
α
D
K
S
1
μ
1
max
⁡
-
α
D
+
k
1
k
2
-
b
-
b
2
-
4
a
c
2
a
-
S
2
in
,
where a, b, and c are defined in (12).Remark 7.
Notice that the coordinates of the bifurcation parameters, provided by (67), (75), and (80), are the same as the upper bounds that guarantee stability of NOP, established in (54). The reason is that the bifurcation points bound the locally stable nature of the NOP, and the parameter coordinates for the bifurcation points are parameter bounds that allow to keep the system under normal operation.Bifurcation Diagrams. Figure 1 shows four diagrams: a 2-parameter bifurcation diagram (in the upper left position) and three 1-parameter bifurcation diagrams. For the 2-parameter bifurcation diagram D and S
1
in are the bifurcation parameters. The vertical line in this subfigure is the bifurcation branch FB(X
2); it is a fold bifurcation branch. The solid curved line corresponds to the bifurcation branch TB(X
1); it is a transcritical bifurcation branch. The dash-dot curved line (very close to FB(X
2)) corresponds to the bifurcation branch TB(X
2); it is a transcritical bifurcation branch. The region D
>
(
μ
2
max
⁡
/
α
)
K
I
/
(
K
I
+
2
K
S
2
) located at the right of the vertical solid line implies the non-existence of the NOP. The 1-parameter bifurcation diagrams are computed using D as bifurcation parameter and fixing S
1
in according to the value shown by the arrows in the 2-parameter diagram.Bifurcation diagrams for the NOP. The 2-parameter bifurcation diagram is in upper left position. The corresponding 1-parameter bifurcation diagrams for the specific values ofS
in according to arrows A
1, A
2, and A
3 are in the other positions.
(a)
(b)
(c)
(d)From Proposition1 and (15) and (75), it follows that the value D
=
(
μ
2
max
⁡
/
α
)
K
I
/
(
K
I
+
2
K
S
2
) is not only the parameter coordinate for the bifurcation branch FB(X
2), but also the limit for the existence of the NOP. The arrows displayed in Figure 1(a) and noted as A
1, A
2, and A
3 allow us to define if the bifurcation is either transcritical or fold: (i) arrow A
1 indicates that bifurcation TB(X
1) is transcritical, bifurcation TB(X
2) is transcritical and bifurcation FB(X
2) is fold, (ii) arrow A
3 confirms that bifurcation FB(X
2) is fold and bifurcation TB(X
1) is transcritical, and (iii) arrow A
2 shows the crossing of transcritical and fold bifurcations.Thus,(i)
TB(X
1)—described by (62) to (67) with X
1
eq
=
0—is a transcritical bifurcation that involves the intersection between the equilibrium corresponding to normal operation and the equilibrium corresponding to X
1
eq
=
0. For larger values of D the concentration X
1 would converge to zero, which is known as washout of acidogenic biomass.
(ii)
Bifurcation FB(X
2)—described by (71) to (75)—corresponds to a fold bifurcation where the equilibrium branch corresponding to NOP collides with the other equilibrium branch. The collision and disappearance can only be noticed in the branches of X
2
eq and S
2
eq. For larger values of D, the concentration X
2 would converge to zero, which is known as washout of methanogenic biomass.
(iii)
Bifurcation TB(X
2)—described by (76) to (80) with X
2
eq
=
0—corresponds to a transcritical bifurcation that involves the intersection between the equilibrium corresponding to NOP and the equilibrium corresponding to X
2
eq
=
0. For larger values of D the concentration X
2 would converge to zero, which is known as washout of methanogenic biomass.Figure2 is a zoomed version close to TB(X
2) of Figure 1. This figure has the following features: (i) the region for local stability of the NOP is denoted by NO, (ii) the region for local stability of an equilibrium corresponding to acidogenic biomass washout is denoted by W/X
1, (iii) the region for local stability of an equilibrium corresponding to methanogenic biomass washout is denoted by W/X
2, (iv) the region for local stability of an equilibrium corresponding to washout of both acidogenic and methanogenic biomass is denoted by W/(
X
1
&
X
2
), (v) the region for the fold bifurcation FB(X
2) corresponds to the vertical line, (vi) the regions for the transcritical bifurcations TB(X
1) and TB(X
2) (continuous and dash-doted lines, resp.) are curved lines, (vii) the two big points correspond to coexisting bifurcations: the upper one corresponds to the coexistence of FB(X
2) and TB(X
1) whereas the lower one corresponds to the coexistence of FB(X
2) and TB(X
2), and (viii) in the region D
>
(
μ
2
max
⁡
/
α
)
K
I
/
(
K
I
+
2
K
S
2
), located at the right of the vertical solid line, the NOP does not exist.Figure 2
Two-parameter bifurcation diagram.From Figure1 it can be observed that if S
2
in is fixed there are two ways whereby the local stability of the normal operation condition can be lost: (i) by increasing the value of D, (ii) by decreasing the value of S
1
in, which agrees with Remark 6. In the same figure there are three cases (A
3, A
2, and A
1) in which the occurrence of local bifurcation is analyzed. For some values of D
>
0 and S
1
in there is a locally stable equilibrium point with X
1
eq
>
0, X
2
eq
>
0, S
1
eq
>
0 and S
2
eq
>
0, which is associated with the NOP. For a given value of S
1
in, there is a value of the dilution rate, D
=
D
*, for which a local bifurcation occurs, related to a change of stability (see Figure 1). Further increase of D, satisfying D
>
D
*, leads to a locally stable equilibrium point that involves either X
1
eq
=
0 or X
2
eq
=
0 or both, while S
1
eq
>
0 and S
2
eq
>
0. These equilibria are associated with washout.Figures3, 4, and 5 show the verification of the two-parameter bifurcation diagram of Figure 1. In each case, NOP is chosen by selecting suitable values of S
1
in and D. In order to obtain the diagram, D is increased so that a bifurcation occurs. The difference between the three cases is the characteristics of the first bifurcation attained. For the case A
3 the bifurcation is related to the washout of the methanogenic biomass (X
2). For the case A
2 the bifurcation is related to the washout of both the acidogenic and the methenogenic biomasses (X
1 and X
2), and for the case A
1 the bifurcation is related to the washout of the acidogenic biomass (X
1).Bifurcation diagrams and transient responses forX
1 and X
2 states. Case A
3.
(a)
(b)
(c)
(d)Bifurcation diagrams and transient responses forX
1 and X
2 states. Case A
2.
(a)
(b)
(c)
(d)Bifurcation diagrams and transient responses forX
1 and X
2 states. Case A
1.
(a)
(b)
(c)
(d)In Figures3, 4, and 5, the upper left part shows the bifurcation diagram for the chosen value of S
1
in, with D as bifurcation parameter, and indicates the considered case, either case A
3, A
2, or A
1. The figures surrounding the bifurcation diagram show the transient evolution of the state variables X
1 and X
2. Figure 3 shows the results for the case A
3 with S
1
in
=
6.5
gCOD
/
L and D
∈
{
0.8,1.1,1.3
} days−1 which correspond to NOP, washout operation for X
2, and washout operation for both X
1 and X
2, respectively. Figure 4 shows the results for the case A
2, with S
1
in
=
5.73
gCOD
/
L, and D
∈
{
0.7,1.2
} days−1 which correspond to NOP, and washout operation for both X
1 and X
2, respectively. Figure 5 shows the results for the case A
1, with S
1
in
=
0.5
gCOD
/
L, and D
∈
{
0.1,0.8,1.07,1.2
} days−1 which correspond to NOP for the first value, washout operation for X
2, for the second value, and washout operation for both X
1 and X
2, for third and fourth values.


## 4.2. Determination of Crossing Bifurcation Branches

Now, we analyze the crossing between bifurcations TB(X
1), FB(X
2), and TB(X
2). At the collision points, the bifurcation parameters and bifurcation variable states have the same values. There are two ways whereby the crossing points can be determined: (i) by combining the bifurcation coordinates (62),…,(67);  (71),…,(75);  (76),…,(80), (ii) by combining (58), (59), and (60) with equilibrium conditions (3)–(6).In order to determine the crossing between TB(X
1) and FB(X
2), the bifurcation coordinates (62),…,(67) and (71),…,(75) are combined, so that
(81)
X
1
eq
=
0
,
X
2
eq
=
S
2
in
-
K
I
K
S
2
α
k
3
,
S
1
eq
=
μ
2
max
⁡
K
I
K
S
1
K
I
(
μ
1
max
⁡
-
μ
2
max
⁡
)
+
2
μ
1
max
⁡
K
S
2
,
S
2
eq
=
K
I
K
S
2
,
D
=
μ
2
max
⁡
K
I
α
(
2
K
S
2
+
K
I
)
,
S
1
in
=
μ
2
max
⁡
K
I
K
S
1
K
I
(
μ
1
max
⁡
-
μ
2
max
⁡
)
+
2
μ
1
max
⁡
K
S
2
.Remark 8.
Notice that (i) the parameter coordinatesD and S
1
in are constant, (ii) the equilibrium values of the state variables X
1
eq, S
1
eq, and S
2
eq do not depend on D nor S
1
in nor S
2
in and (iii) the state variable X
2
eq depends on S
2
in but does not depend on D nor S
1
in.To determine the crossing between FB(X
2) and TB(X
2), the bifurcation coordinates (71),…,(75) and (76),…, (80) are combined, so that
(82)
X
1
eq
=
1
k
2
α
K
I
K
S
2
-
S
2
in
,
X
2
eq
=
0
,
S
1
eq
=
μ
2
max
⁡
K
I
K
S
1
K
I
(
μ
1
max
⁡
-
μ
2
max
⁡
)
+
2
μ
1
max
⁡
K
S
2
,
S
2
eq
=
K
I
K
S
2
,
D
=
μ
2
max
⁡
K
I
α
(
2
K
S
2
+
K
I
)
,
S
1
in
=
μ
2
max
⁡
K
I
K
S
1
K
I
μ
1
max
⁡
-
μ
2
max
⁡
+
2
μ
1
max
⁡
K
S
2
+
k
1
k
2
K
I
K
S
2
-
S
2
in
.Remark 9.
Notice that (i) the parameterD is constant, (ii) the parameter S
1
in depends on S
2
in but it does not depend on D, (iii) the state variables X
2
eq, S
1
eq, S
2
eq are constant, and (iv) the state variable X
1
eq depends on S
2
in but does not depend on D nor S
1
in.Figure2 shows the intersection between bifurcations by means of large black points.


## 5. Discussion and Conclusions

Although the system has six equilibria, we only consider the equilibrium associated to NOP and two other ones that interact with it. Only one of the six equilibrium points corresponds to the normal process operation and there are some conditions on the dilution rateD and the influent concentrations S
1
in and S
2
in under which such equilibrium has physical meaning: positive values of the equilibrium concentrations and stable nature.Analysis are shown that there are mainly three ways through which the stable nature of the NOP can be lost: with a large value of the dilution rate, with a low value of the influent COD concentration, and with a low value of the influent VFA concentration. Even more, the loss of the stable nature of the NOP occurs through either a fold or a transcritical bifurcation, where an equilibrium point with washout of some biomass population becomes stable. At the bifurcation points, the value ofD depends on S
1
in and S
2
in as follows: (i) it depends on S
1
in in the bifurcation TB(X
1), (ii) it does not depend on S
1
in nor S
2
in in bifurcation FB(X
2), and (iii) it depends on both S
1
in and S
2
in in bifurcation TB(X
2).For a givenS
1
in, the value of D corresponds to the bifurcation; it is a limit value that should not be exceeded, so that washout is avoided. The dilution rate D is defined as D
=
Q
/
V, where Q is the influent flow rate and V is the active reactor volume. The volume V remains constant, such that the variation of D is due to the variation of the flow rate Q, and limitation of D can be achieved by properly limiting the flow rate Q.For some parameter values, two bifurcations cross, and in such points the bifurcation coordinates have the same values.



---

*Source: 514797-2014-11-13.xml*