Risk index to monitor an anaerobic digester using a dynamic model based on dilution rate, temperature, and pH

Abstract In this work, a risk index to monitor an anaerobic digestion process is established. The index is derived from the dynamical analysis of the system and it is composed by two conditions ensuring an optimum operating conditions inside bioreactor. The analysis is performed using bifurcation theory in order to study the effects of dilution rate, temperature and pH parameters at behavior system. For these purpose an extended version of the mathematical model proposed by Bernard [1] involving the dependence of growth kinetics on temperature and pH is used. The analysis shows both microbial growth and the performance of the bioprocess, when is highly affected by these parameters, indeed the washout condition occurs by combining a fold bifurcation and a transcritical bifurcation. From bifurcation diagrams both safety and optimal operation regions of the bioreactor are defined. Consequently, based on two conditions of stability a risk index to monitor the bioprocess on-line is proposed. The index is evaluated via numerical simulations allowing to detect system destabilization. The results obtained in this work may provide a useful methodology to monitor an anaerobic process and to guarantee the optimal performance of the bioreactor given that the model is sufficiently accurate.


Introduction
Anaerobic digestion is a biological processes used for the stabilization of organic matter by the action of microbial populations. The anaerobic digestion processes have been used since 70's for the treatment of both municipal and industrial solid and liquid waste. From the 80's, industrial wastewater treated by anaerobic digestion began to grow and extend worldwide, such that, in 2007, the overall number of anaerobic reactors treating industrial wastewater reached more about two thousand references and kept on increasing since then [2].
Compared to the aerobic treatment, the main advantages of the anaerobic digestion treatment system are the low amounts of sludge generated and methane production with energy value [3]. Nevertheless, the anaerobic digestion process is an attractive waste treatment for organic materials degradation, its proper operation requires some expertise and care, due to the highly non-linear dynamics of this process. Its efficient operation is affected by several factors such as: load disturbances, system uncertainties, limited online measurement information, constraints on variables, and uncertain kinetics [4]. The last factor is strongly affected by temperature and pH, both being keys for biomass growth during the anaerobic digestion process. Hence, understanding of the complex nonlinear behavior of the anaerobic digestion system allows a qualitative description of the stability of the system. Dynamic analysis of anaerobic digestion processes has been an important and active research area to determine better operational conditions, on-line monitoring, and design control strategies [5]. In this way, simple models such as the model proposed by Bernard [1] have been widely taken due to their relative simplicity and high ca-pacity in order to reproduce the dynamical behavior of main operational parameters of the process [5]. Nowadays, many simplified models have been proposed in the literature in order to describe anaerobic digestion process, depending on the aim, e.g., process understanding, dynamic simulation, optimization, or control [6].
For instance, Shen et al. [7] have investigated the stability of stationary solutions of an anaerobic digestion model. The authors show a sistematically analyses of the system and they argued that this kind of stability analysis provides insight and guidance for anaerobic digestion reactor design, operation and control [7]. Then, Hess and Bernard [8,9] proposed a criterion to evaluate the operational risk from the analysis of the nonlinear system. Dimitrova and Krastanov [10], and Rincón et al. [11] studied an open-loop buckle dynamic analysis and stability of the system varying the dilution rate parameter with the purpose of development control strategies for the system. Sbarciog et al. [3] proposed a methodology to estimate the separatrix between the stable attraction basins of the equilibrium. Their analysis leads to achieve the proper operation of the system. Recently, Benyahia et al. [5] proposed a generic methodology for the stability analysis of a twostep bioprocess system.
In consequence, the latest research of anaerobic systems is focused on the analysis of the equilibrium and local stability of the system against variations in the dilution factor. However, the anaerobic digestion process is affected by other factors, e.g., ammonia concentration, acclimation, sulfate, heavy metals, pH, temperature, presence of other ions [12]. Although they have been well known biological community, the dynamical analysis of anaerobic digestion systems under variation of temperature and pH, they have not been considered yet. Based on the background described above and the knowledge gap identified, the purpose of present work is to establish a risk index to monitor an anaerobic digestion process. In this way, optimal operating conditions by analyzing an anaerobic digestion process under variations in temperature, pH and dilution rate were determined. Then, equilibrium and stability analysis using by bifurcation theory are performed at different operational conditions. A simplified mathematical model proposed by Bernard et al. [1] was considered, assuming that the process occurs in two main stages: acidogenesis and methanogenesis. Also, mathematical expressions of temperature and pH activity coefficients are included in the Monod and Haldane bacterial growth functions, the cardinal temperature model (CTM) [13] for temperature and Angelidaki expression [14] for pH inhibition. The influence of temperature and pH on the two main reaction kinetics (acidogenic and methanogenic) and the dynamic behavior of the bioprocess were analyzed. From bifurcation analysis both safety and optimal operation regions of the bioreactor are defined, and two conditions are established in order to guarantee their optimum performance. These conditions lead to propose a risk index that allows to identify in a preventive way a possible destabilization of the bioreactor. It may be used for on-line monitoring from the measurement of dilution factor, temperature, and pH parameters.
The remainder of this paper is organized as follows. Sec. 2 describes the anaerobic digestion process assuming two main stages, acidogenic and methanogenic. The mathematical model accounting for the effects of temperature and pH inhibition on growth rates is discussed. Sec. 3 details the bifurcation parameters, operational methodology and the software packages used. Then, results are presented and discussed in Sec. 4. The dynamic analysis of the system using bifurcation theory is discussed in Sec. 4.2. Sec. 4.3 showing the results and the global behavior of the system oriented to on-line monitoring, and a risk index is proposed and evaluated in Sec. 4.4. Finally, in Sec. 5 the main conclusions of this work are provided.

Mathematical model
The mathematical model proposed by Bernard et al. [1] for the treatment system of raw industrial wine distillery vinasses is included in present work. This model has been widely taken [3,5,8,11]. The model assumes that the anaerobic digestion process is based on two main reactions: acidogenesis and methanogenesis [1]. Acidogenesis is the stage in which organic substrates S 1 are catabolized by acidogenic bacteria X 1 , generating intermediates compounds with lower molecular weight than usually, such as volatile fatty acids (VFA). The VFA concentration (S 2 ) is degraded by methanogenic bacteria X 2 into methane and carbon dioxide. Based on hydraulic tests, the authors assumed that the reactor behaves like a perfectly mixed tank, and the biomass is uniformly distributed within the reactor [1]. Also, it is assumed likely activated sludge (active biomass in the reactor bed) is without solid liquid separation, then both, the hydraulic and solid retention time are the same.
The mathematical model by component is given bẏ where D is the dilution rate, α is the fraction of bacteria in the liquid phase, k 1 is the yield for substrate degradation, k 2 is the yield for VFA production, k 3 is the yield for VFA consumption, and the superscript ( 0 ) refers to the input concentrations. Besides, µ 1 and µ 2 are the Monod [15] and Haldane's [16] kinetic models growth rate , respectively. Haldane inhibition model is used for describing the inhibition of acidogenesis caused by the VFA (acetic, propionic, butyric, isoburyric, pentanoic and isopentanoic) [1,17].

Growth rate kinetics
The Monod and Haldane's kinetic models growth rate are modified in function of temperature and pH. Both temperature and pH parameters affect directly the survival of bacterial population. Then, mathematical expressions of kinetic models including the effect of temperature and pH are given by where µ 1 is Monod modified, µ 2 Haldane modified, Θ is the cardinal temperature model (CTM) temperature coefficient, I pH is the inhibition factor by pH, K 1 is the halfsaturation constant associated with S 1 , K 2 is the halfsaturation constant associated with S 2 , and K I is the inhibition constant associated with S 2 . Finally, µ 1 max and µ 2 max denotes the maximum value for the acidogenic and methanogenic growth rates, respectively. The Eqs.

Temperature activity coeflcient (Θ)
Gas solubility and microbial growth rates are directly affected by the temperature. A decreased temperature process reduces the metabolic activity of methanogenic bacteria and the reactor's efficiency. On the other hand, an increased-temperature process stimulates the bacterial growth rates but it also leads to higher concentration of ammonia and gases. An anaerobic digestion process is carried out in three ranges of temperature (ambient, mesophilic, and thermophilic). This process under either mesophilic or thermophilic conditions results of higher metabolic rates. Then, the high temperatures could lead high risk for inhibition caused by ammonia toxicity, as well as higher energy requirements than usually [18]. In our case, the anaerobic digestion reactor is operated in the mesophilic range. According with Tchobanoglous [19] the optimal temperature for the bacterial microorganism activity is given in the range [25,35] ∘ C. Also, the metabolic activity of methanogenic bacteria is reduced significantly at the psychrophilic temperature range. In some cases, the growth rate activity can be inactivated (T < 15 ∘ C). Then, according to the system and the characteristics of the wastewater, the following parameters are assumed: T opt = 30 ∘ C, T min = 5 ∘ C, and Tmax = 40 ∘ C. The CTM coefficient behavior for the anaerobic treatment system of vinasse is shown in Fig. 1.
The effect of the temperature on the growth rate activity on anaerobic digestion process have been widely studied by some authors [20,21]. Several authors have found that the effect of temperature on kinetic parameters could be modelled by means of the Arrhenius's equation [22][23][24]. Rosso et al. [13] proposed a suitable temperature model in order to describe the temperature influence in an anaerobic process. This model has been referred in the literature as a cardinal temperature model. The applicability of CTM in anaerobic processes has been shown by Donoso-Bravo et al. [25]. The CTM temperature activity coefficient Θ is given by where T is the operation temperature, T min and Tmax are the lower and upper temperatures where growth does not occur, respectively, and T opt is the temperature at which the maximum specific growth rate equals its optimal value.

pH inhibition factor (I pH )
In anaerobic digestion processes, pH effects on the growth of microorganisms as well as the bioreactor performance, since the presence of ammonia concentration, acclimation, sulphate, and VFA leads to acidity or alkalinity and determines the kind of microorganisms in the anaerobic digester. In effect, low values of pH indicate process destabilization by VFA accumulation. The optimal pH range for the microorganisms activity according to Sanchez [26] is between 6.8 and 7.4 interval. At this pH range, the tolerance of anaerobic microorganisms in anaerobic process is better, while it favor's optimum operation condition of methanogenic bacteria growth. At lower pH values (i.e., pH ≤ 4), the VFA production and methanogenic bacteria inhibition occurs by acidifying and, at higher pH values (i.e., pH ≥ 8.2), the activity inhibition of the bacterial populations is presented in the anaerobic process. The pH inhibition factor I pH is given by [14,27] I pH = 1 + 2 × 10 0.5(pH LL −pH UL ) 1 + 10 (pH−pH UL ) + 10 (pH LL −pH) , (5) where, pH UL and pH LL are the upper and lower limits of pH, when the specific microbial growth rate is reduced 50% of its initial value without inhibition, respectively. In present work, pH UL =8.0 and pH LL =6.0 have been considered. The I pH factor behavior for the treatment system is shown in Fig. 2.

Materials and methods
Considering that the behavior of the Monod and Haldane kinetic models is well known in widespread applications as fermentative processes [28], in present paper the behavior of these models under the effects of temperature and pH was analyzed and discussed. Several authors [5,8,11] have reported the analysis of the existence and stability of the equilibrium points of the anaerobic digestion model of Bernard et al. [1] termed as AM2 model. However, the dynamic analysis of the system in function of the temperature and pH has not yet been reported. Here, analytical equilibrium points of the system are computed taking into account the effects of temperature and pH on the growth rates expressions. The equilibrium analysis of the fourdimensional state-space system is given by doing the lefthand side of Eq. (1) equal to zero: Then, based on the equilibrium analysis, a dynamical analysis of the system and a sensitive numerical analysis as a function of D, T and pH parameters are performed. The bifurcation analysis of the system is carried out using by bifurcation theory. The software package MATCONT [29] was used for the numerical simulations. MATCONT solves the equilibrium values of X 1 , S 1 , X 2 and S 2 from an starting point from the differential equations given by Eq. (6). First, all the parameters were kept fixed other than D, which was taken as a bifurcation parameter. Subsequently, the same procedure was performed taking T and pH as a bifurcation parameters. The bifurcation diagrams were analyzed by regions in order to study the number and nature of the equilibrium in the range of the operating parameters. The stability analysis around the equilibrium point was carried out by the Lyapunov's indirect method [30]. The result of its application leads to a necessary but not a sufficient conditions for the stability of the system. Therefore, the Hurwitz criteria [31] was applied in order to make it accuracy (see A).
Once the system has been described using the bifurcation theory with the considered bifurcation parameters (D, T and pH), the study was focused on the identification of critical values and operational optimal conditions. The global behavior of the system is studied by stages involving the effects of temperature and pH using the operational ranges defined by the bifurcation analysis. As mentioned by some authors, the acidogenic stage can be studied independently of the methanogenic stage [5,8,11]. Then, the behavior of each stage was studied separately based on dynamical analysis results, and the conditions that ensure the stability of the system were established.
Finally, the last step in the proposed analysis was the establishment of a risk index criteria, monitoring the bioreactor online. Numerical simulations using the experimental data reported by Bernard et al. [1] in an anaerobic up-flow reactor were performed, with the purpose of evaluating the risk-index applicability. All these results are presented in Section 4.

Effects of pH and T on growth rate kinetics
In this section, the behavior of the kinetic growth rates as a function of temperature and pH is described. The analysis is based in the assumption that both acidogenic and methanogenic microorganisms are affected by the same mechanism either from deviations in temperature and pH instead. For the sake of simplicity, identical notation in all figures is noted. In Fig. 3, the influence of temperature and pH in the microbial Monod growth rate (µ 1 ) as a function of the organic substrate S 1 is shown. The growth rate is an increasing function when S 1 ≥ 0, and it exhibits a maximum value at µ 1 (+∞) = µ 1 max . In Fig. 3(a), the growth rate as a function of the temperature is plotted between 10 ∘ C and 40 ∘ C. As can be observed in this figure, at low temperatures the growth rate is reduced dramatically and at high temperatures as 30 ∘ C, the maximum curve of the growth rate is obtained. The growth rate as a function of the pH is shown in Fig. 3(b). The numerical simulations are performed in the range of pH 5.0-9.0 in order to compare the bioreactor performance at different pH values. In the case of pH, as values move away from neutral pH (pH=7), a remarkable decrease of the growth rate is observed through the organic substrate values S 1 .  Fig. 4 shows the behavior of the Haldane growth rate as a function of VFA concentration (S 2 ), temperature and pH. This growth rate model is used in the methanogenic stage and shows the following behavior: i) increasing function in the range 0 ≥ S 2 ≥ S max 2 with a maximum value at µ 2 (S max 2 ) and ii) decreasing function at S 2 > S max 2 . Fig.  4(a) shows the Haldane growth rate (µ 2 ), changing with the temperature factor (Θ) at 10,15,20,25,30,35, 38 ∘ C temperatures. Here, a maximum curve of values is observed at T opt = 30 ∘ C. At this temperature (T opt ) the system exhibits an optimum conditions for the bacterial microorganisms activity. Below of these temperatures, the tendency of the curves is decreasing. The effect of the pH factor (I pH ) over the Haldane growth rate model is depicted in Fig. 4(b). Here, the maximum curve of biomass growth rate is obtained at the optimal value of pH (i.e., pH=7.0) and the behavior of the growth rate is decreasing at values away from the neutrality.
The simulations were performed using the kinetic parameters values reported by Bernard et al. [1] (see Table 1). The results show that the qualitatively behavior of the kinetics growth rates followed a similar trend, however, their quantitatively response and dependence to effects of temperature and pH are different.

Bifurcation analysis
As can be seen in B, the analytical solutions of the equilibrium points are function of D, T and pH parameters. In this section, the results of the bifurcation analysis assuming D, T and pH as bifurcation parameters of the system are shown. In all cases, simulations were performed using the inlet concentrations and parameters values reported in Table 1.

Effect of dilution rate (D)
Bifurcation analysis for substrate and biomass concentration has been carried out assuming D as a bifurcation parameter of the system. The continuation in D reaches two limit points (LP) and three branch points (BP).  Table 1. stable equilibrium. This system exhibits two LP (LP 1 and LP 2 ) at D = 1.0719. Besides, LP 2 is associated to washout of methanogenic bacteria. Similarly, BP 1 at D = 0.9100 (washout of acidogenic biomass), BP 2 at D = 0.9720 (washout) and BP 3 at D = 1.6493 (washout condition). These bifurcation points are qualitatively equivalent to the results reported by Benyahia et al. [5]. The bifurcation diagrams from the dynamical analysis of the system exhibit a phenomenon denoted as backward bifurcation, which involves the existence of a fold bifurcation (LP 1 ) at D = 1.0719 and a transcritical bifurcation (BP 1 ) at D = 0.9100. On the other hand, LP 2 shows a second saddle-node bifurcation, which occurs when the critical equilibrium has a null eigenvalue. At the limit point, backward continuation produces the branch point (BP 2 ) at D = 0.9720. Hence, the bacterial population becomes extinct. Moreover, the bifurcation diagrams show that, depending on the bifurcation parameter value, some equilibrium points may exist or not. As can be seen in Fig. 5, the system is divided and analyzed into four regions as a function of D values.
In region I (D < 0.9100), the system shows three equilibrium points with physical meaning, one of them is stable and the others are unstable. In region II (0.9100 < D < 0.9720), the system shows four equilibrium points with physical meaning, two of them is stable and the others are unstable. In region III (0.9720 < D < 1.0719), the system shows five equilibrium points with physical mean-ing, two of them are stable and the remaining points are locally asymptotically unstable. At region III, although a stable operational node is presented, it is possible that disturbances affecting the system lead to the other stable equilibrium, which corresponds to the washout condition. Then, in region IV (D > 1.0719) there are two equilibrium points with physical meaning, one of them is stable and the other one is unstable. Although optimal operating ranges and stability of the system against variations in D have been extensively studied [5,8,11], a suitable operation range of the bioreactor in Region I at D < 0.91d −1 is defined.

Effect of temperature parameter (T)
In this section, bifurcation analysis of the system is performed by using T as a bifurcation parameter. In Fig. 6, methanogenic biomass concentration X 2 as a function of the bifurcation parameter T is shown. Here, bifurcation diagram of the system exhibits four LP and six BP with physical meaning. As can be seen in Fig. 6, the bifurcation diagram exhibits a symmetry in the range of temperatures analyzed. In Fig. 6(a), there are two LP denoted by LP 1 and LP 2 at T = 15.52 ∘ C as well as three BP denoted by BP 1 , BP 2 and BP 3 at T = 13.25 ∘ C, 16.57 ∘ C and 16.13 ∘ C, respectively. In Fig. 6(b), there are two LP denoted by LP 3 and LP 4 at T = 38.44 ∘ C and also there are three BP denoted by BP 4 , BP 5 and BP 6 at T = 38.11 ∘ C, 38.25 ∘ C and 39.03 ∘ C, respectively. In Fig. 6, LP 2 and LP 4 are associated to washout by acidogenic biomass, and BP 1 , BP 3 , BP 5 and BP 6 points are associated to washout condition. The "washout condition"means that there are not biomass in the bioreactor by absence of both, methanogenic and acidogenic bacterias. Fig. 6 shows two-fold bifurcations (LP 1 and LP 3 ), where the behavior of the equilibrium points switches from stable to unstable. At these limit points, a forward continuation leads to the extinction of biomass population of the branch points BP 2 and BP 5 . At these points the system exhibits an exchange stability at T =16.57, 38.25 ∘ C by transcritical bifurcation.
The analysis of the number and nature of the equilibrium is studied by regions, then, the system is divided into five regions (from I to V) delimited by the bifurcation points obtained. In Region I at T <13.25 ∘ C and T >39.04 ∘ C, the system shows one equilibrium point with physical meaning, being stable. In Region II (13.25 ∘ C< T <15.52 ∘ C and 38.44 ∘ C < T <39.04 ∘ C) the system has two equilibrium points with physical meaning, one of them is stable and the other one is unstable. Region III (15. , there are four equilibrium points with physical meaning and only one of them is stable. At this range of temperature, the bacterial growth is favoured, and a stable equilibrium point is guaranteed. These conditions favour a single equilibrium point, which corresponds to normal operating conditions. Therefore, the washout condition by temperature effects is avoided.

Effect of potential hydrogen parameter (pH)
In this section, pH is used as a bifurcation parameter of the system. Here, four LP and six BP bifurcation points with physical meaning are observed. In Fig. 7 the methanogenic biomass concentration X 2 changing with the pH parameter is presented. Fig. 7(a) shows two LP points denoted by LP 1 and LP 2 at pH=5.59 and also there are three BP denoted by BP 1 , BP 2 and BP 3 at pH=5.35, 5.69 and 5.63, respectively. In Fig. 7 Fig. 7, the system is divided into five regions (from I to V) as follows. In Region I (pH< 5.35 and pH> 8.65), the system present one equilibrium point with physical meaning and stable. In Region II (5.35 <pH< 5.59 and 8.41 >pH> 8.65), the system shows two equilibrium points with physical meaning, one of them is locally asymptotically stable and the other one unstable. In region III (5.59 <pH< 5.65 and 8.35 >pH> 8.41), there are six equilibrium points with physical meaning, two of them are stable and the others locally asymptotically unstable. In Region IV (5.65 >pH> 5.69 and 8.31 >pH> 8.35), the system exhibits five equilibrium points with physical meaning, two equilibrium are stable and the others unstable. Finally, in Region V (5.69 >pH> 8.31), four equilibrium points with physical meaning are determined, only one of them locally asymptotically stable. At this range of pH, the bacterial growth is favoured, and the washout condition is avoided. Then, in Region V suitable operational conditions for anaerobic digestion process are given.
The equilibrium points in terms of Θ and I pH are summarized in Table 2. In this table, the stability conditions and feasible equilibrium points are computed according to the established regions in the bifurcation analysis. Notice that, at Regions III and IV, the system exhibits two stable stationary solution competing. The former, corresponds to a nontrivial solution of the system while the latter with a washout by methanogenic biomass. Figs. 6 and 7 show how the equilibrium curves Branch 1 and Branch 2 cross each other and change their direction at the bifurcation points denoted by LP. Indeed, the Branch 1 shows the qualitative change of the system properties under quantitative parameter variation at a named fold bifurcation (LP 1 and LP 3 ) and at transcritical bifurcations (BP 2 and BP 4 ). The transcritical bifurcations are associated to the washout condition. Notice that the values between LP 1 and BP 2 , and LP 3 and BP 4 , are quite close. Therefore, as limit operational values are considered the BP 2 and BP 4 bifurcation points. At these bifurcation point values the conditions that leads to the system collapse are avoided. Then, from bifurcation analyses a normal operation value of temperature and pH in terms of Θ and I pH should be between 0.37 and 1.

Overloading tolerance of the bioprocess 4.3.1 Acidogenic stage
Two cases were studied solving Eqs. (6a)-(6b) withẊ 1 = S 1 = 0. In the first case at ΘI pH <0.37, the acidogenic stage exhibits a unique solution with physical meaning, a stable equilibrium point (E 0 1 ), i.e., S * 1 =S 0 1 and X * 1 = 0. According to Fig. 8(a), this case corresponds to a trivial solution to the washout condition by absence of acidogenic biomass. In the second case, the acidogenic stage exhibits two equilibrium points: one of them is a stable equilibrium point (E 1 1 ) and the other one an unstable equilibrium point (E 0 1 ). Moreover, E 1 1 at ΘI pH ≥0.37 corresponding to a nontrivial solution (see B) of the set equations from the acidogenic stage, i.e., S * 1 < S 0 1 and X * 1 >0 (see Fig. 8(b)). Thus, from Eq. (23), it is defined that the stable operating point E 1 1 in terms of the parameters is guaranteed as long as the condition Condition 1: holds.

Methanogenic stage
In this section the methanogenic stage is studied. The simulation was carried out for S * 1 < S 0 1 and X * 1 >0, which correspond to a stable operational point different to the washout condition. The curves in Fig. 9 show the phasediagram results for the system of differential equations given by Eqs. (6c)-(1d) withẊ 2 =Ṡ 2 = 0.

General operational indicator control-oriented risk index
Based on the previous dynamical analyses, it is shown that the stability of the system depends on certain conditions  Fig. 9: Phase diagram where E 0 2 corresponds to the washout of X 2 , E 1 2 is the operating stable equilibrium point of methanogenic stage and E 2 2 is the unstable equilibrium point of acidogenic stage. (a) ΘI pH <0.37 (S * 2 =S 0 2 ), (b) 0.37≤ ΘI pH <0.64 (S * 2 ≤ S 0 2 + k 2 that could be measurable in terms of the values of D, T and pH. In fact, from previous sections, it was determined that, while conditions defined in Sec. 4.2 guarantee the operability of the system, some other conditions guarantee the stability of the process. Therefore, an effective on-line monitoring and control parameter during the anaerobic digestion process is the accomplishment of Conditions 1 and 2 given by Eq. (7) and Eq. (8), respectively. The following assumptions are considered: (i) the acidogenic biomass (X 1 ), methanogenic biomass (X 2 ), organic substrate (S 1 ) and VFA (S 2 ) concentrations are unknown, (ii) the kinetic parameters are known, (iii) the values of dilution rate (D), temperature (T) and pH are known, and (iv) the dilution rate is saturated in the ranges of the bioreactor operation known from dynamic analysis (between known constant bounds).
Then, a risk index is proposed in order to warn in a preventative way a system destabilization and hence, the washout condition. A monitoring strategy that minimizes the risk of reactor destabilization due to T, pH and D ef- fects is given by where SRI is the sensitive risk index. If Conditions 1 and 2 are satisfied, there is no of risk (SRI= 0), therefore, the reactor is in a normal and stable operational area. Otherwise, if those conditions are not satisfied, the risk-index Eq. (9) indicates that the behavior of the reactor is in undesirable operational area with washout risk (SRI= 1). At this value of the SRI, operational decisions in order to avoid the washout are needed. Then, temperature and pH conditions should be kept into an optimum range in order to obtain high performance of the bacteria activity. Also, both bioreactor retention time and nutrients must be guaranteed for the growth of biomass, as well as, avoid toxic compounds into the bioreactor. It takes time for the reactor to reach its equilibrium condition (Condition 4). In the particular case of the washout condition, an inoculation is necessary.

Numerical simulation using a risk index
The risk index is given by Eq. (9) while the simulations were computed by using D and pH data from Bernard et al. [1]. The numerical results are shown in Fig. 10. In this figure, both Conditions 1 and 2 are included. Dotted line is the limit condition that corresponds to the constant value calculated from the right side of the inequalities at Eq. (7) and Eq. (8). Solid line corresponds to left side of Eqs. (7)- (8). For the sake of simplicity, normalized values are used in Fig. 10(a).
The risk index computations are shown in Fig. 10(b). As can be seen in this figure, risk-index equal zero is the value used to define stable and normal operational, and 1 is used for unstable condition. The non-compliance with the stability conditions means that the system is conducted to destabilization, and possibly led to the washout.
A comparative analysis between the behavior of the substrate and VFA concentrations reported by Bernard et al. [1], and the estimation of the risk index (see Fig.  10(b)) allows to see alerts when there is an increase in the VFA concentration. These increments are one of the most common causes of destabilization in bioreactors. As it has mentioned before by other authors [8], high values of the risk index (close to one) are representative of regimes of acid accumulation. Nevertheless, not all VFA have the same impact on the stability of the system. For example, the increase of the acetic acid concentration can be the result of a respective increase of the feed-stock loading, which does not imply a subsequent instability. On the other hand, the accumulation of propionic and butyric acids strongly indicate a potential instability (coupled with the increase of hydrogen production). In any case, all VFA affect the bioreactor performance and its stability.
In this case, a possible acidification process from the dilution rate and pH measurements were detected through the values of the risk index, and the index allows to prevent the destabilization of the reactor without direct measurement of the state variables and changes in the measured parameters. At this values of SRI equal to 1, operational decisions such as changes in the hydraulic retention time or an alkaline pretreatment before entering the a uent into the bioreactor are needed.
The proposed risk index could be applied in real processes using the adequate on-line instrumentation to measure the required parameters. The on-line measurements of the state variables (S 1 , X 1 , S 2 and X 2 ) are not required to calculate the risk index. Therefore, the measurements requirements are limited to dilution rate, pH and temperature, it allows to make a sensibility analysis of the risk index and evaluate its applicability on-line easily.

Conclusions
An efficient risk index criteria from the dynamical analysis of the anaerobic digestion process has been established. For the study of the dynamic behavior of nonlinear system, an extension of the phenomenological model initially proposed by Bernard [1] was considered. The mathematical model involving an explicit temperature and pH functions into the bacterial growth rates in order to evaluate the effect of these parameters in the microbial growth and hence on the bioprocess behavior.
The bifurcation analysis allowed to characterize changes in the qualitative behavior of the anaerobic digestion process by varying of the dilution rate, temperature and pH parameters. The dynamic behavior of the system exhibits regions with multiple steady states, unstable operational points, bifurcation points and washout conditions. The bifurcation diagrams showed normal and optimal operating regions in which the operational behavior was becoming stable and the hazardous conditions as washout are presented. Then, from the dynamic analysis have established two conditions that guaranteed a safety operating region to normal working of the bioreactor. These conditions were considered to established a risk index criteria to diagnosis and predict a future destabilization of the bioreactor from the dilution rate, temperature and pH measures.
This criterion could be used to monitor the operational conditions of the bioprocess on-line and in real time. The bioreactor could be monitored by on-line sensing devices in situ. In fact, the criterion proposed here is a use-ful tool for designing a real-time simulator for supervisory control of the process. For its implementation is only required a commercial instrumentation available and a visualization programming tool. The obtained results confirm the importance of dynamical analysis oriented to monitoring, optimization and control of anaerobic digestion processes. The risk-index based criterion proposed in this paper should be elaborated to decide and take the best actions and strategies in order to anticipate destabilization problems by temperature and pH in the bioreactor.

Future work
The anaerobic digestion process is a multi-parametric system, and additional analytical studies of the dynamical process could provide more insight into the process and strengthen the results in order to design a satisfying control strategies of these processes. Therefore, future research topics are oriented to study more complex systems with multi-stage units, recycle and multiplicity of steady states. For example, an increase in the ratio of carbon dioxide and methane could suggests the transition from a welloperated steady-state to a potential instability, therefore, future research could involving for the dynamical analysis the gaseous phase of the model. Thus studies will be oriented to analyze the behavior of the bioreactor by the effects of pH in function of the VFA and carbon dioxide concentration. Finally, biogas recovery from anaerobic digestion process reduces green house gas emission. Then, the study and implementation of more a sophisticated systems as the anaerobic membrane bioreactors are needed. Work in this direction is currently in progress.
Thus, the linearisation around x * leads to the following Jacobian matrix where A 11 = µ * 1 − αD, A 13 = µ , 1 X * 1 , A 22 = µ * 2 − αD, A 24 = µ , 2 X * 2 , A 31 = −k 1 µ * 1 , A 33 = −k 1 µ , 1 X * 1 − D, A 41 = k 2 µ * 1 , A 42 = −k 3 µ 2 , A 43 = k 2 µ , 1 X * 1 , and A 44 = −k 3 µ , 2 X * 2 − D. Moreover, µ * 1 , µ * 2 are growth rate kinetic model evaluated at S * 1 and S * 2 , respectively. Likewise, µ ′ 1 is the derivative of µ 1 with respect to S * 1 and µ ′ 2 is the derivative of µ 2 with respect to S * 2 both computed as follow: After perform some algebraic operations, the Jacobian eigenvalues are given by The stability analysis around the equilibrium point using the Lyapunov's indirect method is restricted to infinitesimal neighbourhoods of the equilibrium point [30], i.e., this method is a necessary condition but not sufficient in order to ensure the stability of the system. Hence, the Hurwitz criterion [31] is used to establish some conditions with the range of values from the parameters ensuring system stability conditions.

Theorem 2. Routh-Hurwitz Criterion.
Given the polynomial P(λ) = λ n + a 1 λ n−1 + · · · + a n−1 λ + an, where the coefficients a i are real constants, i = 1, · · · , n, define the n Hurwitz matrices using the coefficients a i of the characteristic polynomial. All of the roots of the characteristic polynomial have negative real part if and only if the determinants of all Hurwitz matrices are positive.
Thus, the characteristic equation of the Jacobian matrix related to the system is given by λ 4 + a 1 λ 3 + a 2 λ 2 + a 3 λ + a 4  In order to ensure that Re(λ i ) < 0, for all i, it is considered each term independently and it is verified, according to the Routh-Hurwitz simplified criterion, that a 1 > 0, a 3 > 0, a 4 > 0, and a 1 a 2 a 3 > a 2 3 + a 2 1 a 4 [31].

B Analytical solutions of equilibrium points
The solutions of the system must be real positive and the feasible values of the state variables X 1 , X 2 , S 1 , S 2 satisfy the following qualitative properties: Thus, the following analytical expressions for the equilibrium are obtained, all of them assuming inlet biomass concentration equal zero (X 0 1 = X 0 2 = 0).