Monitoring management criterion of an anaerobic digester using an explicit model based on temperature and pH

TL;DR

This paper presents a stability analysis of an anaerobic digester using an explicit model based on temperature and pH, identifying bifurcation points and proposing an online risk monitoring criterion for safe and optimal operation.

Contribution

It introduces a simplified explicit model incorporating temperature and pH effects, along with a bifurcation analysis and a novel online risk criterion for bioreactor management.

Findings

Washout occurs via fold and transcritical bifurcations.

The proposed risk criterion enables early detection of destabilization.

Operational safety and optimality regions are clearly defined.

Abstract

In this paper, the stability analysis of an anaerobic digestion process is presented. The analysis is performed using a simplified mathematical model, which includes explicit temperature and pH dependence on kinetic growth rates. A detailed nonlinear and bifurcation analyses are performed in order to study the effects of temperature and pH parameters on the process behaviour. In addition, both safety and optimal operation regions of the bioreactor are established. Based on bifurcation diagrams, it is observed that the washout condition occurs by combining a fold bifurcation and a transcritical bifurcation. Consequently, to prevent the washout and guarantee optimal operation conditions of the bioreactor, a risk criterion oriented to monitor the bioprocess on-line is proposed, allowing to detect the system destabilisation.

Tables

Table 1. Table 1: Inlet concentrations and nominal kinetic parameter values

Parameter	Values of Ref. [12]	Unit
$k_1$	42.14	g COD$/$g X1
$k_2$	116.5	mmol VFA$/$g X1
$k_3$	268	mmol VFA$/$g X2
$K_1$	7.1	g$/$L
$K_2$	9.28	mmol$/$L
$K_I$	256	(mmol$/$L)2
$S_1^0$	15.6	g$/$L
$S_2^0$	112.7	mmol$/$L
$X_1^0$	0.0	g$/$L
$X_2^0$	0.0	g$/$L
$\alpha$	0.5	–
${\mu }_{1\max }$	1.2	d-1
${\mu }_{2\max }$	0.74	d-1

Table 2. Table 2: Stability of equilibrium points of the system

Region	Condition	Equilibria	Lyapunov Stability	Hurwitz Stability Criterion
I	or	$(0,S_1^0,0,S_2^0)$	Stable	Ensured
II	or	$(0,S_1^0,0,S_2^0)$	Unstable	Not ensured
		$(X_1^{\ast },S_1^{\ast },0,{\hat{S}}_2)$	Stable	Ensured
III	or	$(0,S_1^0,0,S_2^0)$	Unstable	Not ensured
		$(X_1^{\ast },S_1^{\ast },0,{\hat{S}}_2)$	Stable	Not ensured
		$(X_1^{\ast },S_1^{\ast },X_2^{\ast },S_2^{\ast })$	Stable	Ensured
		$(X_1^{\ast },S_1^{\ast },X_2^{\ast },S_2^{\ast })$	Unstable	Not ensured
		$(0,S_1^0,X_2^{\ast },S_2^{\ast })$	Unstable	Not ensured
		$(0,S_1^0,X_2^{\ast },S_2^{\ast })$	Unstable	Not ensured
IV	or	$(0,S_1^0,0,S_2^0)$	Unstable	Not ensured
		$(X_1^{\ast },S_1^{\ast },0,{\hat{S}}_2)$	Stable	Ensured
		$(X_1^{\ast },S_1^{\ast },X_2^{\ast },S_2^{\ast })$	Stable	Ensured
		$(X_1^{\ast },S_1^{\ast },X_2^{\ast },S_2^{\ast })$	Unstable	Not ensured
		$(0,S_1^0,X_2^{\ast },S_2^{\ast })$	Unstable	Not ensured
V	or	$(0,S_1^0,0,S_2^0)$	Unstable	Not ensured
		$(X_1^{\ast },S_1^{\ast },0,{\hat{S}}_2)$	Unstable	Not ensured
		$(X_1^{\ast },S_1^{\ast },X_2^{\ast },S_2^{\ast })$	Stable	Ensured
		$(0,S_1^0,X_2^{\ast },S_2^{\ast })$	Unstable	Not ensured

Equations

$$\dot{X_1}$$

$$\dot{S_1}$$

$$\dot{X_2}$$

$$\dot{S_2}$$

$${\bar{\mu }}_1$$

$${\bar{\mu }}_2$$

$$\Theta =\frac{(T-T_{max})(T-T_{\text{min}}{)}^2}{(T_{opt}-T_{min})(T_b-T_a)},$$

$$I_{pH}=\frac{1+2\times {10}^{0.5(p{H}_{LL}-p{H}_{UL})}}{1+{10}^{(pH-p{H}_{UL})}+{10}^{(p{H}_{LL}-pH)}},$$

$$X_1\ge 0,X_2\ge 0,$$

$$S_1\le S_1^0,S_2\le S_2^0+\frac{k_2}{k_1}{S}_1^0.$$

$$S_2^{\ast }=\frac{-b\pm \sqrt{b^2-4ac}}{2a},$$

$$X_2^{\ast }=\frac{\alpha D(S_2^{\ast }{K}_I+S_2^{\ast }{S}_0^2+K_2{K}_I)-S_2^{\ast }{\hat{\mu }}_2{K}_I}{S_2^{\ast }D{\alpha }^2{k}_3},$$

$$S_1^{\ast }$$

$$X_1^{\ast }$$

$$S_2^{\ast }$$

$$X_2^{\ast }$$

$$S_1^{\ast }$$

$$X_1^{\ast }$$

$$S_2^{\ast }$$

$$X_2^{\ast }$$

$$\text{Condition 1:}\frac{{\hat{\mu }}_1}{D}>\alpha .$$

$$\text{Condition 2:}{\left(\frac{\alpha D-{\hat{\mu }}_2}{D}\right)}^2\ge \frac{4K_2{\alpha }^2}{K_I}.$$

$$SRI=\{\begin{array}{rl}& 0\text{for}\frac{{\hat{\mu }}_1}{D}>\alpha \&{\left(\frac{\alpha D-{\hat{\mu }}_2}{D}\right)}^2\ge \frac{4K_2{\alpha }^2}{K_I},\\ & 1\text{for}\frac{{\hat{\mu }}_1}{D}\le \alpha \&{\left(\frac{\alpha D-{\hat{\mu }}_2}{D}\right)}^2<\frac{4K_2{\alpha }^2}{K_I},\end{array}$$

$$\dot{x}=f(x),x\in {\mathbb{R}}^n,$$

$$A=\left[{\begin{array}[]{*{20}{c}}A_{11}&0&A_{13}&0\\ 0&A_{22}&0&A_{24}\\ A_{31}&0&A_{33}&0\\ A_{41}&A_{42}&A_{43}&A_{44}\end{array}}\right],$$

$${\mu }_1^{\prime }=\frac{{\mu }_{1\max }\Theta I_{pH}{K}_1}{{\left(K_1+S_1^{\ast }\right)}^2},$$

$${\mu }_2^{\prime }=\frac{{\mu }_{2\max }\Theta I_{pH}{K}_I\left(K_2{K}_I-{S_2^{\ast }}^2\right)}{{\left(K_2{K}_I+S_2^{\ast }{K}_I+{S_2^{\ast }}^2\right)}^2}.$$

$${\lambda }_{1,2}=\frac{1}{2}(A_{44}+A_{22})\pm \frac{1}{2}\sqrt{{(A_{44}-A_{22})}^2+4A_{42}{A}_{24}},$$

$${\lambda }_{3,4}=\frac{1}{2}(A_{33}+A_{11})\pm \frac{1}{2}\sqrt{{(A_{33}-A_{11})}^2+4A_{31}{A}_{13}}.$$

$${\lambda }^4+a_1{\lambda }^3+a_2{\lambda }^2+a_3\lambda +a_4=0,$$

Full text

††journal: Biochemical Engineering Journal

Nomenclature

[TABLE]

1 INTRODUCTION

The growing interest for the protection of the environment leads to the development of more detailed studies dealing with the analysis of systems used in the wastewater treatment. Anaerobic digestion is one of the most common biological process in the industrial applications and municipal wastewater treatment [1]. Since the eighties, 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 that two thousand references and kept on increasing since then [2].

The anaerobic digestion process is an attractive waste treatment that involves degradation of organic materials by the action of microbial populations. This process is characterized by the existence of highly non-linear dynamics and its efficient operation is affected by several factors such as: load disturbances, system uncertainties, limited online measurement information, constraints on variables, and uncertain kinetics [3]. The last factor is strongly affected by temperature and pH, both being key for biomass growth during the anaerobic digestion process. Hence, the understanding of the complex nonlinear behaviour of the anaerobic digestion system allows to give a qualitative description of the stability of the system.

In the last decade, dynamic analysis of anaerobic digestion process has been studied. Shen et al. [4] have investigated the stability of stationary solutions of an anaerobic digestion model. The authors argued that the dynamics of the system depend on kinetic model and key parameters such as dilution rate.

Also, the authors provided a guidance for anaerobic digestion reactor design, operation and control in terms of robustness and productivity. However, the effects of biomass concentrations and products on the dynamics were not reported. In 2008, Hess and Bernard [5] proposed a simplified criterion to evaluate the operational risk from the analysis of the nonlinear system. This criterion showed a suitable performance in a real plant. When applied, it was able to diagnose an operation strategy by detecting an early destabilization caused by accumulation of volatile fatty acids. Later, Hess [6] presented an extension of the work in [5]. The authors demonstrated with real data the efficiency of the criterion and also proposed a methodology to monitor in real-time the trajectory of the system.

The results were obtained from a general approach, thus the methodology could be applied even without knowing the model parameters. Dimitrova and Krastanov [7] studied the equilibrium points stability and performed a bifurcation analysis of the model varying the dilution rate parameter in order to study the system in open loop, and proposed a feedback control law for asymptotic stabilization of the closed-loop system. Rincón et al. [8] studied an open-loop dynamic analysis and stability of the system varying the dilution rate parameter. The results show that the washout condition occurs after a combination of both fold and transcritical bifurcations. They also used their obtained results to develop an adaptive control strategy. Sbarciog et al. [1] proposed a methodology to estimate the separatrix between the stable attraction basins of the equilibria. Their analysis was based on growth kinetic models with and without inhibition and the results obtained from the inputs and initial states led to achieve the proper operation of the system. Recently, Benyahia et al. [9] proposed a generic methodology for the stability analysis of a two-step bioprocess system.

Dynamic analysis of anaerobic systems based on simplified mathematical models have been an important and active research area to determine better operational conditions, on-line monitoring, and design control strategies [9]. The simplified mathematical models assume a rate-limiting step, and these have been widely used due to their “relative simplicity and high capacity in order to reproduce the dynamical behaviour of main operational parameters of the process”[9] and computational easiness. Nowadays, simplified models have been positioned as a suitable alternative to describe biotechnological processes [10]. Initially, simplified models of anaerobic systems have been studied using specific kinetics with limited number of functioning conditions. Recently, research of anaerobic systems is focused on the analysis of the equilibria and local stability of the system against variations in the dilution factor. Nevertheless, the anaerobic digestion process is affected by other factors, e.g., ammonia concentration, acclimation, sulfate, heavy metals, pH, temperature, presence of other ions [11]. However, the dynamical analysis of anaerobic digestion systems including factors such as inhibitory effects of temperature and pH on the kinetics rates have not been considered yet.

Based on the background described above and the knowledge gap identified, the purpose of the this work is to determine optimal operating conditions by analysing an anaerobic digestion process under variations in dilution rate, temperature and pH. In this way, equilibria and stability analysis using bifurcation theory are performed at different operational conditions. A simplified mathematical model proposed by Bernard et al. [12] 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 CTM model [13] for temperature and Angelidaki expression [14] for pH inhibition. The influence of temperature and pH on the kinetic reactions and the stability of the equilibrium points were determined. To take this analysis into account, an operational indicator of the processes is proposed, which could be used for on-line monitoring from the measurement of dilution factor, temperature, and pH parameters. These criterion allows to identify in a preventive way a possible destabilization of the bioreactor.

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 analytical equilibrium points are presented in Sec. 4.2. The dynamic analysis using bifurcation theory, and the qualitative analysis of the system are discussed in Sec. 4.3. Sec. 4.4 presents the results and the global behaviour of the system oriented to on-line monitoring, and a risk monitoring management is proposed and evaluated in Sec. 4.5. Finally, Sec. 5 provides the main conclusions of this work.

2 MATHEMATICAL MODEL

The mathematical model proposed by Bernard et al. [12] is considered in this paper. The model assumes that the anaerobic digestion process is based on two main reactions: acidogenic and methanogenic. Acidogenesis is the stage in which organic substrates $S_{1}$ are catabolized by acidogenic bacteria $X_{1}$, generating intermediates compounds with lower molecular weight, such as Volatile Fatty Acids (VFA) concentration $S_{2}$. Then, $S_{2}$ is degraded by methanogenic bacteria $X_{2}$ into methane and carbon dioxide. The anaerobic digestion scheme of the process is shown in Fig. 1. The mathematical model by component in the liquid phase is given by

[TABLE]

where $D$ is the dilution rate, $\alpha$ is the fraction of bacterias 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 initial input concentration. Besides, $\mu_{1}$ and $\mu_{2}$ are the Monod and Haldane growth rate kinetic model, respectively [15, 16].

Both temperature and pH parameters affect directly the survival of bacterial population. Low values of pH are an usual indicator of the process destabilization by VFA accumulation. Also, a gas solubility and microbial growth rates are directly affected by the temperature. A decreased temperature process reduces the metabolic activity of methanogenic bacteria. On the other hand, an increased temperature process takes the advantage of stimulating the rate of the microorganisms but also leads higher concentration of ammonia and gases. In this paper, the growth rate kinetic models were modified as a function of temperature and pH as follows.

2.1 Growth rate kinetics

Mathematical expressions of kinetic models including the effect of temperature and pH are given by

[TABLE]

where $\bar{\mu}_{1}$ is Monod modified, $\bar{\mu}_{2}$ Haldane modified, $\Theta$ is the CTM temperature coefficient, $I_{pH}$ is the inhibition factor by pH, $K_{1}$ is the half-saturation constant associated with $S_{1}$, $K_{2}$ is the half-saturation constant associated with $S_{2}$, and KI is the inhibition constant associated with $S_{2}$. Finally, $\mu_{1\max}$ and $\mu_{2\max}$ denotes the maximum value for the acidogenic and methanogenic growth rates, respectively.

2.1.1 Temperature activity coefficient ($\Theta$)

Several authors have found that the effect of temperature on kinetic parameters can be easily modelled by means of theArrhenius’s equation [17, 18, 19]. However, some disadvantages are observed at each step of the microbial population [20]. In this way, 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 (CTM). The applicability of CTM in anaerobic processes has been shown by Donoso-Bravo et al. [21]. The CTM temperature activity coefficient $\Theta$ is given by

[TABLE]

where $T_{a}=(T_{opt}-T_{max})(T_{opt}+T_{min}-2T)$, $T_{b}=(T_{opt}-T_{min})(T-T_{opt})$, ${T}$ is the operation temperature, $T_{\text{min}}$ and $T_{\text{max}}$ are the lower and upper temperatures where the growth rate does not occurs, respectively, and $T_{opt}$ is the temperature at which the maximum specific growth rate equals its optimal value. The optimal temperature range for the bacterial microorganism activity has been established in the range $[25,35]^{\circ}$C according to Tchobanoglous [22]. 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^{\circ}$C).

2.1.2 pH inhibition factor ($I_{pH}$)

In anaerobic digestion processes, the pH affects 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. The pH inhibition factor $I_{pH}$ is given by [23, 14]

[TABLE]

where, $p{H_{UL}}$ and $p{H_{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. The optimal pH range for the microorganisms activity according to Sanchez [24] is given between 6.8 and 7.4. At this pH range, the tolerance of anaerobic microorganisms in anaerobic process is better, while it favours optimum operation condition of methanogenic bacteria growth. At lower pH values (i.e., pH $\leq$ 4), the VFA production and methanogenic bacteria inhibition occurs by acidifying and, at higher pH values (i.e., pH $\geq$ 8.2), the activity inhibition of the bacterial populations is presented in the anaerobic process.

3 MATERIALS AND METHODS

Although the behaviour of the Monod and Haldane kinetic models is well known in differents aplications as fermentative processes [25], in this paper the behaviour of these models under the effects of temperature and pH was analysed and discussed. In addition, several authors [5, 8, 9] have reported the analysis of the existence and stability of the equilibrium points of Bernard’s model [12]. 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 equilibria analysis of the four-dimensional state-space system is given by doing the left-hand side of (1) equal to zero ($\dot{X_{1}}=\dot{X_{2}}=\dot{S_{1}}=\dot{S_{2}}=0$). Then, 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 bifurcation theory. The software package MATCONT [26] is used for numerical simulations.MATCONT solves the equilibrium values of $X_{1}$, $S_{1}$, $X_{2}$ and $S_{2}$ from differential equations in (1) by continuation from an starting point satisfying the condition $\dot{X_{1}}=\dot{S_{1}}=\dot{X_{2}}=\dot{S_{2}}=0$. First, all the parameters were keep 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 analysed by regions in order to study the number and nature of the equilibria in the range of the operating parameters. The stability analysis around the equilibrium point was carried out by the Lyapunov’s indirect method [27]. The result of its application leads to a necessary but not a sufficient conditions for the stability of the system. Therefore, the Hurwitz criterion [28] was applied in order to ensure it (see A).

Once the system has been described using the bifurcation theory with the considered bifurcation parameters ($D$, $T$ and pH), the study of was focused on the identification of critical values and operational optimal conditions. The global behaviour 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 [8, 9, 5]. Then, the behaviour 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 criterion to monitor the bioreactor online. A risk index was proposed in order to monitor online the process. Numerical simulations using the experimental data reported by Bernard et al. [12] in an anaerobic upflow were performed, with the purpose of evaluating the risk-index applicability. All the results are presented in Section 4.

4 RESULTS

4.1 Effects of pH and T on growth rate kinetics

For the sake of simplicity, identical notation in all figures is used. In Fig. 2, the influence of temperature and pH in the microbial Monod growth rate ($\bar{\mu}_{1}$) as a function of the organic substrate $S_{1}$ is shown. As can be observed in this figure, the growth rate is an increasing function when ${S_{1}}\geq 0$, and it exhibits a maximum value at ${\bar{\mu}_{1}}\left({+\infty}\right)={\mu_{1\max}}$. In Fig. 2, the behaviour of the growth rate model is obtained in a range of temperature between 10∘C and 30∘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. 2. 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. 3 shows the behaviour of the Haldane growth rate as a function of VFA concentration, temperature and pH. This growth rate model is used in the methanogenic stage and shows the following behaviour: $i)$ increasing function in the range $0\geq S_{2}\geq S^{\text{max}}_{2}$ with a maximum value at $\bar{\mu}_{2}(S^{\text{max}}_{2})$ and $ii)$ decreasing function at $S_{2}>S^{\text{max}}_{2}$. Fig. 3 shows the Haldane growth rate ($\bar{\mu}_{2}$), changing with the temperature factor ($\Theta$) at 10, 15, 20, 25, 30∘C temperatures. Here, a maximum curve of values is observed at 30∘C corresponding to the optimal temperature. Behind 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. 3. Here, the maximum curve of biomass growth rate is obtained at the optimal value of pH (i.e., pH=7.0) and the behaviour 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. [12] (see Table 1).

4.2 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:

[TABLE]

[TABLE]

Thus, the following analytical expressions for the equilibrium are obtained, all of them assuming inlet biomass concentration equal zero ($X_{1}^{0}=X_{2}^{0}=0$).

Equilibrium 1: (washout condition)

The washout condition is given by the absence of biomass in both acidogenic and methanogenic stages, i.e., ${S_{1}^{*}}=S_{1}^{0}$, ${X_{1}^{*}}=0$, ${S_{2}^{*}}=S_{2}^{0}$, ${X_{2}^{*}}=0$.

Equilibrium 2: (washout by acidogenic biomass)

The washout condition is given by biomass absence in the acidogenic stage, i.e., ${S_{1}^{*}}=S_{1}^{0}$ and ${X_{1}^{*}}=0$. Thus, the equilibrium points associated to the methanogenic stage are

[TABLE]

and

[TABLE]

where $a=\alpha D$, $b=(\alpha D-\hat{\mu}_{2})K_{I}$, $c=\alpha DK_{2}K_{I}$, and $\hat{\mu}_{2}=\mu_{2\max}\Theta{I_{pH}}$. Moreover, $S_{2}$ has a real solution if $a\neq 0$ and the discriminant is positive or zero, $\Delta\geq 0$, i.e., $b^{2}-4ac\geq 0$.

Equilibrium 3: (washout by methanogenic biomass)

The washout by absence of biomass in the methanogenic stage is given by

[TABLE]

where $\varUpsilon=k_{2}S_{1}^{0}+k_{2}K_{1}+S_{2}^{0}k_{1}$ and $\hat{\mu}_{1}=\mu_{1\max}\Theta I_{pH}$.

Equilibrium 4: (nontrivial solution)

The analytical nontrivial solution is given by

[TABLE]

where, $\gamma={k_{2}}(S_{1}^{0}+{K_{1}})+S_{2}^{0}{k_{1}}+{k_{1}}{K_{I}}$, $\xi=S_{2}^{0}{k_{1}}\hat{\mu}_{1}({K_{I}}+1)+{k_{2}}S_{1}^{0}\hat{\mu}_{1}$, $\varphi=S_{2}^{*}{k_{1}}{K_{I}}\hat{\mu}_{2}$ and $\psi=\alpha D{K_{2}}{K_{I}}{k_{1}}$, and $\phi=(\hat{\mu}_{1}-\alpha D)\varphi-\left({\hat{\mu}_{1}+\alpha D}\right)\psi$.Notice that the equilibrium points of the bioreactor are function of $D$, $T$ and pH parameters.

4.3 Bifurcation analysis

This section presents the results of bifurcation analysis assuming $D$, $T$ and pH as bifurcation parameters of the system. In all cases, the simulations were performed using the inlet concentrations and parameters values in Table 1. Notice that all parameters are real positive.

4.3.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). Fig. 4 shows the behaviour of methanogenic biomass concentration $X_{2}$ changing with $D$. In this figure, solid lines represent stable equilibria and dotted lines represent unstable equilibria. This system exhibits two LP (LP1 and LP2) at $D=1.0719$. Besides, LP2 is associated to washout of methanogenic bacteria. Similarly, BP1 at $D=0.9100$ (washout of acidogenic biomass), BP2 at $D=0.9720$ (washout) and BP1 at $D=1.6493$ (washout condition). These bifurcation points are qualitatively equivalent to the results reported by Benyahia et al. [9]. 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 (LP1) at $D=1.0719$ and a transcritical bifurcation (BP1) at $D=0.9099$. On the other hand, LP2 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 (BP2) 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. 4, the system is divided and analysed into four regions as a function of $D$ values. In region I ($D<0.9720$), the system shows three equilibrium points with physical meaning, one of them is stable and the others are unstable. In region II ($0.9720<D<1.0719$), the system shows five equilibrium points with physical meaning, two of them are stable and the remaining points are locally asymptotically unstable. In 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, 9], a suitable operation range of the bioreactor in Region I at $D<0.91d^{-1}$ is defined.

4.3.2 Effect of temperature parameter ($T$)

In this section, bifurcation analysis of the system is performed by using $T$ as a bifurcation parameter. In Fig. 5, 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. 5, the bifurcation diagram exhibits a symmetry in the range of temperatures analysed. In Fig. 5, there are two LP denoted by LP1 and LP2 at $T=15.52^{\circ}$C as well as three BP denoted by BP1, BP2 and BP3 at $T=13.25^{\circ}$C, 16.57*∘C and 16.13∘C, respectively. In Fig. 5, there are two LP denoted by LP3 and LP4 at $T=38.44^{\circ}$C and also there are three BP denoted by BP4, BP5 and BP6 at $T=38.11^{\circ}$C, 38.25∘C and 39.03∘*C, respectively. In Fig. 5, LP2 and LP4 are associated to washout by acidogenic biomass, and BP1, BP3, BP5 and BP6 points are associated to washout condition.

Fig. 5 shows two-fold bifurcations (LP1 and LP3), where the behaviour 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 BP2 and BP5. 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 equilibria 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.52∘C$<T<$16.13∘C and 38.25∘C$<T<$38.44∘C), shows six equilibrium points with physical meaning, two of them are stable and the others are unstable. In the case of Region IV (16.13∘C$>T>$16.57∘C and 38.11∘*C$>T>$38.25 *∘C), there are five equilibrium points with physical meaning, two of them are locally asymptotically stable and the others are unstable. Finally, in Region V (16.57∘C$>T>$38.11∘*C), 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.

4.3.3 Effect of potential hydrogen parameter (pH)

In this section, pH is used as a bifurcation parameter of the system behaviour. Here, four LP and six BP bifurcation points with physical meaning are observed. The methanogenic biomass concentration $X_{2}$ changing with the pH parameter is presented in Fig. 6. Fig. 6 shows two LP points denoted by LP1 and LP2 at pH=5.59 and also there are three BP denoted by BP1, BP2 and BP3 at pH=5.35, 5.69 and 5.63, respectively. In Fig. 6 there are two LP denoted by LP3 and LP4 at pH=8.41, as well as three BP denoted by BP4, BP5 and BP6 at pH=8.31, 8.35 and 8.65, respectively. The behavioural analysis of the system leads to identify the existence of two fold bifurcations LP1 and LP3 at pH$=$5.59 and pH$=$8.41, respectively. At these limit points, forward and backward continuation produces the transcritical bifurcations BP2 and BP4 at pH=5.69 and pH= 8.31.

In Fig. 6, 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 equilibria 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 $\Theta$ 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. Figures 5 and 6 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 (LP1 and LP3) and at transcritical bifurcations (BP2 and BP4). The transcritical bifurcations are associated to the washout condition. Notice that the values between LP1 and BP2, and LP3 and BP4, are quite close. Therefore, as limit operational values are considered the BP2 and BP4 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 $\Theta$ and $I_{pH}$ should be between 0.37 and 1.

4.4 Overloading tolerance of the bioprocess

4.4.1 Acidogenic stage

Equations ((1a))-((1b)) with $\dot{X_{1}}=\dot{S_{1}}=0$ are solved and two cases were studied. In the first case at $\Theta I_{pH}<$0.37, the acidogenic stage exhibits a unique solution with physical meaning, a stable equilibrium point (E${}_{1}^{0}$), i.e., S${}_{1}^{*}$=S${}_{1}^{0}$ and X${}_{1}^{*}=0$. According to Fig. 7, 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${}_{1}^{0}$). Moreover, E${}_{1}^{1}$ at $\Theta I_{pH}\geq$0.37 corresponding to a nontrivial solution (see Sec. 4.2) of the set equations from the acidogenic stage, i.e., S${}_{1}^{*}<S_{1}^{0}$ and X${}_{1}^{*}>$0 (see Fig. 7). Thus, from ((14)), it is defined that the stable operating point E${}_{1}^{1}$ in terms of the parameters is guaranteed as long as the condition

[TABLE]

holds.

4.4.2 Methanogenic stage

In this section the methanogenic stage is studied. The simulation was carried out for S${}_{1}^{*}<S_{1}^{0}$ and X${}_{1}^{*}>$0, which correspond to an stable operational point different to the washout condition. The curves in Fig. 8 show the phase-diagram results for the system of differential equations given by ((1c))-((1d)) with $\dot{X_{2}}=\dot{S_{2}}=0$. From the solution of these differential equations, three particular solution cases are obtained: the first case $\Theta I_{pH}<$0.37, the methanogenic stage shows a stable equilibrium point (E${}_{2}^{0}$), when S${}_{2}^{*}$=S${}_{2}^{0}$ and X${}_{2}^{*}=$0. This equilibrium point corresponds to the washout condition by absence of methanogenic biomass. The phase diagram is depicted in Fig. 8. In the second case 0.37$\leq\Theta I_{pH}<$0.64, the methanogenic stage shows three equilibrium points with physical meaning (E${}_{2}^{0}$, E${}_{2}^{1}$ and E${}_{2}^{2}$). The equilibrium points E${}_{2}^{1}$ and E${}_{2}^{2}$ corresponds to a nontrivial solution at S${}_{1}^{*}<S_{1}^{0}$ and S${}_{2}^{*}\leq S_{2}^{0}+\frac{k_{2}}{k_{1}}S_{1}^{0}$. As can be observed in Fig. 8, the equilibrium points E${}_{2}^{0}$ and E${}_{2}^{1}$ are asymptotically stable. Finally, in the third case $\Theta I_{pH}\geq 0.64$, the methanogenic stage shows two equilibrium points (E${}_{2}^{0}$ and E${}_{2}^{1}$). The results of this case are shown in Fig. 8. The stable equilibrium point E${}_{2}^{1}$ corresponds with a nontrivial solution, which guarantees a stable and optimal operation condition of the process. Then, the stability and operability of the system are guaranteed when the following system condition is taken into account:

[TABLE]

4.5 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 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.3 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 ((18)) and ((19)), respectively. 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$ effects is given by

[TABLE]

where both $\hat{\mu}_{1}$ and $\hat{\mu}_{2}$ are function of $T$ and pH.

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 (20) indicates that the behaviour of the reactor is in undesirable operational area with washout risk (SRI$=1$), and it is necessary to make operational decisions.

4.5.1 Numerical simulation using a risk index

The risk index is given by ((20)) while the simulations were computed by using $D$ and pH data from Bernard et al. [12]. The numerical results are shown in Fig. 9. 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 (18) and (19). Solid line corresponds to left side of (18)-(19). For the sake of simplicity, normalized values are used in Fig. 9.

The risk index computations are shown in Fig. 9. As can be seen in this figure, 0 is the value used to define stable and normal operational, and 1 is used for unstable condition. Here, the results are shown as a function of time in order to compare with the behaviour of experimental data reported by Bernard et al. [12]. The non-compliance with the stability conditions means that the system is conducted to destabilization, and possibly led to the washout. As one can notice, the on-line measurement of state variables ($S_{1}$, $X_{1}$, $S_{2}$ and $X_{2}$) is difficult, but in this case the risk index only depends on kinetics and the easiness easily of measuring on-line the parameters of the process.

A comparative analysis between the behaviour of the substrate and VFA concentrations reported by Bernard et al. [12], and the estimation of the risk index (see Fig. 9) 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 [5], high values of the risk index (close to one) are representative of regimes of acid accumulation.

In this case, a possible acidification process from the dilution rate and pH measurements were detected through the values of the risk index. Then, this index allows to prevent the destabilization of the reactor without direct measurement of the state variables and changes in the measured parameters. The proposed risk index could be applied in real processes using the adequate on-line instrumentation to measure the required parameters. Therefore, the measurements allowing to make a sensibility analysis of the risk index and evaluate its applicability on-line.

5 Conclusions

A mathematical model for an anaerobic digester reactor with two main reaction stages -methanogenic and acidogenic- was presented in this paper. The model takes into account the influence of some important variables such as temperature and pH. A dynamic analysis of non-trivial and biologically feasible steady state using bifurcation theory was performed. The systematic dynamical analysis was developed taking into account variations of temperature and pH on the growth rate models. Regions with multiple steady states, unstable operational points, washout conditions and transcritical bifurcations were found.

Furthermore, a detailed study of the bifurcation branches was performed. The qualitative and quantitative study of the equilibrium points allowed to identify normal and optimal operational conditions in which the operational condition was becoming stable and the washout is avoided.

In this way, a risk management criterion was established as a function of dilution rate, pH, and temperature measurements. The criterion is based on two conditions established from the dynamic analysis of the system. The purpose of the criterion is to generate an alert when the process evolves into an unstable operational area and at washout risk. This criterion could be used to monitor the operational conditions of the bioprocess on-line and in real time, i.e., the bioreactor could be monitored by on-line sensing devices in situ. In fact, the criterion proposed here is a useful tool for designing a real-time simulator for supervisory control of the process. The implementation in a real process will be require using a commercial instrumentation available and a visualization programming tool. The monitoring system requires the use of more economical instrumentation compared with the actual strategies used nowadays for the state variables measurements.

The results obtained 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 could be used to decide and take the best actions and strategies in order to anticipate future problems in the bioreactor.

Future research topics are oriented to study more complex systems with multi-stage units, recycle and multiplicity of steady states. The ranges of temperature and pH established here through the bifurcation analysis are the goal of the control design strategies. Work in this direction is currently in progress.

Acknowledgments

A.M. Alzate wants to acknowledge the program “Programa de becas para estudiantes sobresalientes de posgrado”from Universidad Nacional de Colombia - Sede Manizales and the Institut de Robótica i Informática Industrial (IRI) from the Universitat Politécnica de Catalunya.

Appendix A Stability analysis of the system

Theorem 1

Lyapunov’s Indirect Method. Consider a dynamical system defined by

[TABLE]

where $x$ are the state variables and $f:\mathbb{R}^{n}\longmapsto\mathbb{R}^{n}$ is smooth function. Suppose that it has an equilibrium $x^{*}$, and $A$ denotes the Jacobian matrix of $f(x)$ evaluated at $x^{*}$. Then, $x^{*}$ is asymptotically stable if all eigenvalues $\lambda_{1}$, $\lambda_{2}$,$\cdots$, $\lambda_{n}$ of $A$ satisfy $Re({\lambda_{i}})<0$.

Thus, the linearisation around $x^{*}$ leads to the following Jacobian matrix

[TABLE]

where $A_{11}={\mu_{1}^{*}}-\alpha D$, $A_{13}=\mu_{1}^{,}{X_{1}^{*}}$, $A_{22}={\mu_{2}^{*}}-\alpha D$, $A_{24}=\mu_{2}^{,}{X_{2}^{*}}$, $A_{31}={-{k_{1}}{\mu_{1}^{*}}}$, $A_{33}=-{k_{1}}\mu_{1}^{,}{X_{1}^{*}}-D$, $A_{41}={{k_{2}}{\mu_{1}^{*}}}$, $A_{42}={-{k_{3}}{\mu_{2}}}$, $A_{43}={{k_{2}}\mu_{1}^{,}{X_{1}^{*}}}$, and $A_{44}=-{k_{3}}\mu_{2}^{,}{X_{2}^{*}}-D$. Moreover, ${\mu_{1}^{*}}$, ${\mu_{2}^{*}}$ are growth rate kinetic model evaluated at $S_{1}^{*}$ and $S_{2}^{*}$, respectively. Likewise, $\mu_{1}^{\prime}$ is the derivative of $\mu_{1}$ with respect to $S_{1}^{*}$ and $\mu_{2}^{\prime}$ is the derivative of $\mu_{2}$ with respect to $S_{2}^{*}$ both computed as follow:

[TABLE]

[TABLE]

After perform some algebraic operations, the Jacobian eigenvalues are given by

[TABLE]

[TABLE]

The stability analysis around the equilibrium point using the Lyapunov’s indirect method is restricted to infinitesimal neighbourhoods of the equilibrium point [27],* i.e.*, this method is a necessary condition but not sufficient in order to ensure the stability of the system. Hence, the Hurwitz criterion [28] 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 polinomial $P(\lambda)=\lambda^{n}+a_{1}\lambda^{n-1}+\cdots+a_{n-1}\lambda+a_{n}$, where the coefficients $a_{i}$ are real constants, $i=1,\cdots,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

[TABLE]

where the coefficients are $a_{1}=-A_{44}-A_{33}-A_{22}-A_{11}$, $a_{2}=-A_{42}A_{24}+A_{44}(A_{33}+A_{22}+A_{11})-A_{31}A_{13}+A_{33}(A_{22}+A_{11})+A_{22}A_{11}$, $a_{3}=A_{42}A_{24}(A_{33}+A_{11})+A_{44}(A_{31}A_{13}-A_{33}A_{22}-A_{33}A_{11}-A_{22}A_{11})+A_{22}(A_{31}A_{13}-A_{33}A_{11})$, and $a_{4}=A_{42}A_{24}(A_{31}A_{13}-A_{33}A_{11})-A_{44}A_{22}(A_{31}A_{13}+A_{33}A_{11}).$ In order to ensure that $Re({\lambda_{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_{3}^{2}+a_{1}^{2}a_{4}$ [28].