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

1 Facultad de Ingenieŕıa y Arquitectura, Universidad Católica de Manizales, Carrera 23 No. 60-63, Manizales 170002, Colombia 2Departamento de Fı́sica y Matemáticas, Universidad Autónoma de Manizales, Antigua Estación del Ferrocarril, Manizales 170001, Colombia 3 Departamento de Ingenieŕıa Eléctrica, Electrónica y Computación-Percepción y Control Inteligente-Bloque Q, Facultad de Ingenieŕıa y Arquitectura, Universidad Nacional de Colombia-Sede Manizales, Campus La Nubia, Manizales 170003, Colombia


Introduction
The anaerobic digestion processes have several advantages with respect to aerobic digestion ones, for instance, lower sludge production and higher pollution degradation [1][2][3][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][4][5][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][8][9][10][11][12][13][14][15][16][17][18][19][20][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.

Mathematical Problems in Engineering
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][23][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][23][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 socalled 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 Section 2, 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.

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 CO 2 and CH 4 ) and for growth. The mass balance model is (cf. [27])̇1 where 1 represents the concentration of acidogenic biomass, 2 represents the concentration of methanogenic biomass, 1 represents the concentration of organic substrate characterized by its chemical oxygen demand, and 2 represents the concentration of volatile fatty acids (VFA). In addition, represents the dilution rate and the proportion of biomass not attached to the reactor. On current operation of the bioreactor the concentrations 1 , 2 , 1 , and 2 are all positive, and the equilibrium point called NOP tends to ( eq 1 , eq 2 , eq 1 , eq 2 ). When the dilution rate experiences an excessive increase, the washout phenomenon occurs. In this case 1 = 0 or 2 = 0 or both 1 = 0 and 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 ( ) 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 1 ̸ = 0 and 2 ̸ = 0, we have where eq 1 = 1 max The NOP has coordinates: Normal operation implies that the reactor concentrations are positive and is positive, so that the combination of eq 1 > 0, eq 2 > 0, eq 1 > 0, eq 2 > 0, > 0, and (8)-(11) yields the restrictions on and in 1 . Since eq 2 is a positive real number, it follows that − > 0 and, from (11), that the condition given by inequation (13) has to be fulfilled: The following proposition presents some conditions must be satisfied by in order the NOP exists. (12), conditions established in (13) and the equilibrium point provided by (8)- (11). The dilution rate and the equilibrium concentration

Proposition 1. Consider the definitions given in
Proof. Replacing (7) in (4) it follows that The expression above indicates that has a maximum value for a certain value of eq 2 . Differentiating (17) yields and making this expression zero yields Substituting it into (2) yields In addition inequation (13) is only fulfilled if < max ; then, eq 2 has a maximum when 2 − 4 = 0 which is given by (21) eq Replacing (20) in (21), we obtain that the maximum value of eq 2 for the NOP fulfills Other conditions can be computed from (8)- (11). From (8) it is easy to see that the condition for the existence of NOP is eq 1 < in 1 (on contrary case, eq 1 ≤ 0); this expression combined with (10) implies Therefore, if satisfies (23) then eq 1 > 0 and eq 1 < in 1 . In a similar way, with a bit more complicated algebraic steps it is possible to find conditions on eq 2 and such that the existence of NOP can be proven. Equation (9) yields eq 2 < in 2 . By combining this expression with (9), (10), and (11), we find conditions on such that eq 2 > 0. These conditions are wherẽsatisfies (25), According to the results found above, the final condition to guarantee the existence of NOP is wherẽsatisfies (25).

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][29][30]. The Jacobian of the system (1)-(2) is given by = ( ) .
According to the results obtained above, the final condition to NOP is stable is given by wherẽ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) if in 1 is overly small, then the dilution rate is constrained to small values even if in 2 is large and (ii) if in 2 is overly small, then the dilution rate is constrained to small values even if in 1 is large. If in 1 and in 2 satisfy then the dilution rate is only limited by the constant value ( 2 max / )( /( +2√ 2 )). Therefore, it is desirable to have large values of in 1 , and in 2 , as it would allow to have larger values of the dilution rate .

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 Table 1 the values of the equilibrium points when bifurcations occur are shown. The interaction of these equilibria raise three different bifurcations which are called TB( 1 ), FB( 2 ), and TB( 2 ). The bifurcation TB( 1 ) is of transcritical type and involves the intersection between the equilibrium corresponding to NOP and the equilibrium corresponding to inactivated acidogenic biomass ( 1 = 0). The bifurcation FB( 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( 2 ) is of transcritical type and involves the intersection between the equilibrium corresponding to NOP and the one corresponding to inactivated methanogenic biomass ( 2 = 0). The mentioned bifurcations have common points where the parameters and state variables have the same values (see Table 2). In this section, the determination of the bifurcation points as well as crossing points between bifurcations TB( 1 )

Determination of Bifurcation Points.
One-parameter bifurcation diagrams can be obtained from Table 1. 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 rate , the feed concentration in 1 , and the feed concentration in 2 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)   Consequently, we have three possibilities: 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 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.

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 2parameter bifurcation diagram (in the upper left position) and three 1-parameter bifurcation diagrams. For the 2parameter bifurcation diagram and in 1 are the bifurcation parameters. The vertical line in this subfigure is the bifurcation branch FB( 2 ); it is a fold bifurcation branch. The solid curved line corresponds to the bifurcation branch TB( 1 ); it is a transcritical bifurcation branch. The dash-dot curved line (very close to FB( 2 )) corresponds to the bifurcation branch TB( 2 ); it is a transcritical bifurcation branch. The region > ( 2 max / ) /( + 2√ 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 as bifurcation parameter and fixing in 1 according to the value shown by the arrows in the 2-parameter diagram.
From Proposition 1 and (15) and (75), it follows that the value = ( 2 max / ) /( + 2√ 2 ) is not only the parameter coordinate for the bifurcation branch FB( 2 ), but also the limit for the existence of the NOP. The arrows displayed in Figure 1(a) and noted as 1 , 2 , and 3 allow us to define if the bifurcation is either transcritical or fold: (i) arrow 1 indicates that bifurcation TB( 1 ) is transcritical, bifurcation TB( 2 ) is transcritical and bifurcation FB( 2 ) is fold, (ii) arrow 3 confirms that bifurcation FB( 2 ) is fold and bifurcation TB( 1 ) is transcritical, and (iii) arrow 2 shows the crossing of transcritical and fold bifurcations.
Thus, (i) TB( 1 )-described by (62) to (67) with eq 1 = 0is a transcritical bifurcation that involves the intersection between the equilibrium corresponding to normal operation and the equilibrium corresponding Mathematical Problems in Engineering to eq 1 = 0. For larger values of the concentration 1 would converge to zero, which is known as washout of acidogenic biomass.
(ii) Bifurcation FB( 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 eq 2 and eq 2 . For larger values of , the concentration 2 would converge to zero, which is known as washout of methanogenic biomass.
(iii) Bifurcation TB( 2 )-described by (76) to (80) with eq 2 = 0-corresponds to a transcritical bifurcation that involves the intersection between the equilibrium corresponding to NOP and the equilibrium corresponding to eq 2 = 0. For larger values of the concentration 2 would converge to zero, which is known as washout of methanogenic biomass. acidogenic biomass washout is denoted by W/ 1 , (iii) the region for local stability of an equilibrium corresponding to methanogenic biomass washout is denoted by W/ 2 , (iv) the region for local stability of an equilibrium corresponding to washout of both acidogenic and methanogenic biomass is denoted by W/( 1 & 2 ), (v) the region for the fold bifurcation FB( 2 ) corresponds to the vertical line, (vi) the regions for the transcritical bifurcations TB( 1 ) and TB( 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( 2 ) and TB( 1 ) whereas the lower one corresponds to the coexistence of FB( 2 ) and TB( 2 ), and (viii) in the region > ( 2 max / ) /( + 2√ 2 ), located at the right of the vertical solid line, the NOP does not exist.
From Figure 1 it can be observed that if in 2 is fixed there are two ways whereby the local stability of the normal operation condition can be lost: (i) by increasing the value of , (ii) by decreasing the value of in 1 , which agrees with Remark 6. In the same figure there are three cases ( 3 , 2 , and 1 ) in which the occurrence of local bifurcation is analyzed. For some values of > 0 and in 1 there is a locally stable equilibrium point with eq 1 > 0, eq 2 > 0, eq 1 > 0 and eq 2 > 0, which is associated with the NOP. For a given value of in 1 , there is a value of the dilution rate, = * , for which a local bifurcation occurs, related to a change of stability (see Figure 1). Further increase of , satisfying > * , leads to a locally stable equilibrium point that involves either eq 1 = 0 or eq 2 = 0 or both, while eq 1 > 0 and eq 2 > 0. These equilibria are associated with washout. Figures 3,4,and 5 show the verification of the twoparameter bifurcation diagram of Figure 1. In each case, NOP is chosen by selecting suitable values of in 1 and . In order to obtain the diagram, is increased so that a bifurcation occurs. The difference between the three cases is the characteristics of the first bifurcation attained. For the case 3 the bifurcation is related to the washout of the methanogenic biomass ( 2 ). For the case 2 the bifurcation is related to the washout of both the acidogenic and the methenogenic biomasses ( 1 and 2 ), and for the case 1 the bifurcation is related to the washout of the acidogenic biomass ( 1 ).
In Figures 3, 4, and 5, the upper left part shows the bifurcation diagram for the chosen value of in 1 , with as bifurcation parameter, and indicates the considered case, either case 3 , 2 , or 1 . The figures surrounding the bifurcation diagram show the transient evolution of the state variables 1 and 2 . Figure 3 shows the results for the case 3 with in 1 = 6.5 gCOD/L and ∈ {0.8, 1.1, 1.3} days −1 which correspond to NOP, washout operation for 2 , and washout operation for both 1 and 2 , respectively.  shows the results for the case 2 , with in 1 = 5.73 gCOD/L, and ∈ {0.7, 1.2} days −1 which correspond to NOP, and washout operation for both 1 and 2 , respectively. Figure 5 shows the results for the case 1 , with in 1 = 0.5 gCOD/L, and ∈ {0.1, 0.8, 1.07, 1.2} days −1 which correspond to NOP for the first value, washout operation for 2 , for the second value, and washout operation for both 1 and 2 , for third and fourth values.

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 rate and the influent concentrations in 1 and in 2 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 of depends on in 1 and in 2 as follows: (i) it depends on in 1 in the bifurcation TB( 1 ), (ii) it does not depend on in 1 nor in 2 in bifurcation FB( 2 ), and (iii) it depends on both in 1 and in 2 in bifurcation TB( 2 ). For a given in 1 , the value of corresponds to the bifurcation; it is a limit value that should not be exceeded, so that washout is avoided. The dilution rate is defined as = / , where is the influent flow rate and is the active reactor volume. The volume remains constant, such that the variation of is due to the variation of the flow rate , and limitation of can be achieved by properly limiting the flow rate .
For some parameter values, two bifurcations cross, and in such points the bifurcation coordinates have the same values.