## 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.

## 1 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 capacity 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 two-step 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.

## 2 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 by

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].

### 2.1 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 half-saturation constant associated with

*S*

_{1},

*K*

_{2}is the half-saturation constant associated with

*S*

_{2}, and K

*is the inhibition constant associated with*

_{I}*S*

_{2}. Finally,

*μ*

_{1 max}and

*μ*

_{2 max}denotes the maximum value for the acidogenic and methanogenic growth rates, respectively. The Eqs. (2) and (3) are reduced to the Bernard’s model at

*I*= 1 and

_{pH}*Θ*= 1, i.e.

*T*= 30

*C and pH=7.*

^{∘}#### 2.1.1 Temperature activity coefficient (*Θ*)

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*= 30

_{opt}*C,*

^{∘}*T*

_{min}= 5

*C, and*

^{∘}*T*

_{max}= 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 2 *et* 2*al*. [25]. The CTM temperature activity coefficient *Θ* is given by

where *T _{a}* = (

*T*−

_{opt}*T*)(

_{max}*T*+

_{opt}*T*− 2

_{min}*T*),

*T*= (

_{b}*T*−

_{opt}*T*)(

_{min}*T*−

*T*),

_{opt}*T*is the operation temperature,

*T*

_{min}and

*T*

_{max}are the lower and upper temperatures where growth does not occur, respectively, and

*T*is the temperature at which the maximum specific growth rate equals its optimal value.

_{opt}#### 2.1.2 pH inhibition factor (*I*_{pH})

_{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]

where, *pH _{UL}* and

*pH*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,

_{LL}*pH*=8.0 and

_{UL}*pH*=6.0 have been considered. The

_{LL}*I*factor behavior for the treatment system is shown in Fig. 2.

_{pH}## 3 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 four-dimensional state-space system is given by doing the left-hand 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.

## 4 Results

### 4.1 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 *ii*) decreasing function at *μ*_{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*= 30

_{opt}*C. At this temperature (*

^{∘}*T*) 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 (

_{opt}*I*) 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.

_{pH}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.

Parameter | Values of Ref. [1] | Unit |
---|---|---|

k_{1} | 42.14 | g COD/g X_{1} |

k_{2} | 116.5 | mmol VFA/g X_{1} |

k_{3} | 268 | mmol VFA/g X_{2} |

K_{1} | 7.1 | g/L |

K_{2} | 9.28 | mmol/L |

K_{I} | 256 | (mmol/L)^{2} |

15.6 | g/L | |

112.7 | mmol/L | |

0.0 | g/L | |

0.0 | g/L | |

α | 0.5 | – |

μ_{1 max} | 1.2 | d^{−1} |

μ_{2 max} | 0.74 | d^{−1} |

### 4.2 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.

#### 4.2.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. 5 shows the behavior of methanogenic biomass concentration *X*_{2} changing with *D*. In this figure, solid lines represent stable equilibrium and dotted lines represent unstable 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 meaning, 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.

#### 4.2.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. 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

*C, 16.57*

^{.25∘}*C and 16.13*

^{∘}*C, respectively. In Fig. 6(b), there are two LP denoted by LP*

^{∘}_{3}and LP

_{4}at

*T*= 38

*C and also there are three BP denoted by BP*

^{.44∘}_{4}, BP

_{5}and BP

_{6}at

*T*= 38

*C, 38.25*

^{.11∘}*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.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.2.3 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} atpH=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(b) there are two LP denoted by LP_{3} and LP_{4} at pH=8.41, as well as three BP denoted by BP_{4}, BP_{5} and BP_{6} at pH=8.31, 8.35 and 8.65, respectively. The behavior analysis of the system leads to identify the existence of two fold bifurcations LP_{1} and LP_{3} at pH=5.59 and pH=8.41, respectively. At these limit points, forward and backward continuation produces the transcritical bifurcations BP_{2} and BP_{4} at pH=5.69 and pH= 8.31.

In 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*should be between 0.37 and 1.

_{pH}Region | Condition | Equilibria | Lyapunov Stability | Hurwitz Stability Criterion |
---|---|---|---|---|

I | Θ < 0.22 or | Stable | Ensured 1200 | |

I < 0.22_{pH} | ||||

II | 0.22 < Θ < 0.34 or | Unstable | Not ensured | |

0.22 < I < 0.34_{pH} | Stable | Ensured | ||

III | 0.34 < Θ < 0.37 or | Unstable | Not ensured | |

0.34 < I < 0.37_{pH} | Stable | Not ensured | ||

Stable | Ensured | |||

Unstable | Not ensured | |||

Unstable | Not ensured | |||

Unstable | Not ensured | |||

IV | 0.37 < Θ < 0.40 or | Unstable | Not ensured | |

0.37 < I < 0.40_{pH} | Stable | Ensured | ||

Stable | Ensured | |||

Unstable | Not ensured | |||

Unstable | Not ensured | |||

V | 0.40 < Θ ≤ 1.00 or | Unstable | Not ensured | |

0.40 < I ≤ 1.00_{pH} | Unstable | Not ensured | ||

Stable | Ensured | |||

Unstable | Not ensured |

### 4.3 Overloading tolerance of the bioprocess

#### 4.3.1 Acidogenic stage

Two cases were studied solving Eqs. (6a)-(6b) with *ΘI _{pH}* <0.37, the acidogenic stage exhibits a unique solution with physical meaning, a stable equilibrium point

*i.e*.

*ΘI*≥0.37 corresponding to a nontrivial solution (see B) of the set equations from the acidogenic stage,

_{pH}*i.e*.

holds.

#### 4.3.2 Methanogenic stage

In this section the methanogenic stage is studied. The simulation was carried out for

From the solution of these differential equations, three particular solution cases are obtained: the first case *ΘI _{pH}* <0.37, the methanogenic stage shows a stable equilibrium point

*ΘI*<0.64, the methanogenic stage shows three equilibrium points with physical meaning

_{pH}and *ΘI _{pH}* ≥ 0.64, the methanogenic stage shows two equilibrium points

### 4.4 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.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:

the acidogenic biomass (X

_{1}), methanogenic biomass (X_{2}), organic substrate (S_{1}) and VFA (S_{2}) concentrations are unknown,the kinetic parameters are known,

the values of dilution rate (

*D*), temperature (*T*) and pH are known, andthe 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* effects 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.

#### 4.4.1 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 affluent 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.

## 5 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 useful 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.

## 6 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 well-operated 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.

## Acknowledgement

A. M. A. 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 (CSIC-UPC) from the Universitat Politécnica de Catalunya.

**Conflict of Interest**Conflict of Interests: The authors declare that there is no conflict of interests regarding the publication of this paper.

## References

[1] O. Bernard, Z. Hadj-Sadok, D. Dochain, A. Genovesi, and S.P. Steyer. Dynamical model development and parameter identification for an anaerobic wastewater treatment process, biotechnology and bioengineering. *Biotechnol. Bioeng*. 75:143–438, 2001.10.1002/bit.10036Search in Google Scholar
PubMed

[2] J. Jimenez, E. Latrille, J. Harmand, A. Robles, J. Ferrer, D. Gaida, C. Wolf, F. Mairet, O. Bernard, V. A.-Gonzalez, H. M.-Acosta, D. Zitomer, D. Totzke, H. Spanjers, F. Jacobi, A. Guwy, R. Dinsdale, G. Premier, S. Mazhegrane, G. R.-Filippi, A. Seco, T. Ribeiro, A. Páuss, and J.-P. Steye. Instrumentation and control of anaerobic digestion processes: a review and some research challenges. *Rev. Environ. Sci. Biotechnol* 14(4):615–648, 2015.10.1007/s11157-015-9382-6Search in Google Scholar

[3] M. Sbarciog, M. Loccufier, and E. Noldus. Determination of appropriate operating strategies for anaerobic digestion systems. *Biochem. Eng. J*. 51(3):180–188, 2010.10.1016/j.bej.2010.06.016Search in Google Scholar

[4] H.O. Méndez-Acosta, B. Palacios-Ruiz, V. Alcaraz-González, V. González-Álvarez, and J.P. García-Sandoval. A robust control scheme to improve the stability of anaerobic digestion processes. *J. Process Contr*. 20(4):375–383, 2010.10.1016/j.jprocont.2010.01.006Search in Google Scholar

[5] B. Benyahia, T. Sari, B. Cherki, and J. Harmand. Bifurcation and stability analysis of a two step model for monitoring anaerobic digestion processes. *J. Process. Contr*. 22(6):1008–1019, 2012.10.1016/j.jprocont.2012.04.012Search in Google Scholar

[6] A. Donoso-Bravo, J. Mailier, C. Martin, J. Rodríguez, C. Aceves, and A. Vande. Model selection, identification and validation in anaerobic digestion: a review. *Water Res*. 45(17):5347–5364, 2011.10.1016/j.watres.2011.08.059Search in Google Scholar
PubMed

[7] S. Shen, G.C. Premier, A. Guwy, and R. Dinsdale. Bifurcation and stability analysis of an anaerobic digestion model. *Nonlinear Dyn*. 48(4):391–408, 2007.10.1007/s11071-006-9093-1Search in Google Scholar

[8] J. Hess and O. Bernard. Design and study of a risk management criterion for an unstable anaerobic wastewater treatment process. *J. Process Contr*. 18(1):71–79, 2008.10.1016/j.jprocont.2007.05.005Search in Google Scholar

[9] J. Hess and O. Bernard. Advanced dynamical risk analysis for monitoring anaerobic digestion process. *Biotechnol. Progr*. 25(3):643–653, 2009.10.1002/btpr.120Search in Google Scholar
PubMed

[10] N. Dimitrova and M. Krastanov. Nonlinear stabilizing control of an uncertain bioprocess model. *Int. J. Appl. Math. Comput. Sci*. 19(3):441–454, 2009.10.2478/v10006-009-0036-0Search in Google Scholar

[11] A. Rincón, F. Angulo, and G. Olivar. Control of an anaerobic digester through normal form of fold bifurcation. *J. Process Contr*. 19(8):1355–1367, 2009.10.1016/j.jprocont.2009.04.006Search in Google Scholar

[12] Y. Chen, J.J. Cheng, and K.S. Creamer. Inhibition of anaerobic digestion process: a review. *Bioresource Technol*. 99(10):4044–4064, 2008.10.1016/j.biortech.2007.01.057Search in Google Scholar

[13] L. Rosso, J.R. Lobry, and J.P. Flandrois. An unexpected correlation between cardinal temperatures of microbial growth highlighted by a new model. *J. Theor. Biol*. 162(4):447–463, 1993.10.1006/jtbi.1993.1099Search in Google Scholar

[14] I. Angelidaki, L. Ellegaard, and B.K. Ahring. A mathematical model for dynamic simulation of anaerobic digestion of complex substrates: Focusing on ammonia inhibition. *Biotechnol. Bioeng*. 42:159–166, 1993.10.1002/bit.260420203Search in Google Scholar

[15] J. Monod. La technique de culture continué. théorie et applications. *Ann. Inst. Pasteu*. 79(4):390–410, Oct 1950.10.1016/B978-0-12-460482-7.50023-3Search in Google Scholar

[16] L.Y. Lokshina, V.A. Vavilin, R.H. Kettunen, J.A. Rintala, C. Holliger, and A.N. Nozhevnikova. Evaluation of kinetic co-efficients using integrated monod and haldane models for low-temperature acetoclastic methanogenesis. *Water Res*. 35(12):2913–2922, 2001.10.1016/S0043-1354(00)00595-9Search in Google Scholar

[17] H.H. Gavala, I. Angelidaki, and B.K. Ahring. *Kinetics and Modeling of Anaerobic Digestion Process* pages 57–93. Springer Berlin Heidelberg, Berlin, Heidelberg, 2003.10.1007/3-540-45839-5_3Search in Google Scholar
PubMed

[18] L. Deng, H. Yang, G. Liu, D. Zheng, Z. Chen, Y. Liu, X. Pu, L. Song, Z. Wang, and Y. Lei. Kinetics of temperature effects and its significance to the heating strategy for anaerobic digestion of swine wastewater. *Applied Energy* 134:349 – 355, 2014.10.1016/j.apenergy.2014.08.027Search in Google Scholar

[19] G. Tchobanoglous, F. Burton, and H.D. Stensel. *Wastewater engineering: treatment and reuse* Metcalf & Eddy, Inc., New York, 2003.Search in Google Scholar

[20] R. H. Kettunen and J. A. Rintala. The effect of low temperature (5-29 * ^{∘}*c) and adaptation on the methanogenic activity of biomass.

*Appl Microbiol Biotechnol*48:570–576, 1997.10.1007/s002530051098Search in Google Scholar PubMed

[21] S.G. Pavlostathis and E. Giraldo-Gomez. Kinetics of anaerobic treatment: a critical review. *Crit. Rev. Env. Sci. Tec*. 21(5-6):411–490, 1991.10.1080/10643389109388424Search in Google Scholar

[22] R. Crites and G. Tchobanoglous. *Small and decentralized wastewater management systems* McGraw-Hill Companies Inc., New York, 1998.Search in Google Scholar

[23] G.C. Banik, T. Viraraghavan, and R.R. Dague. Low temperature effects on anaerobic microbial kinetic parameters. *Environ. Technol*. 19(5):503–512, 1998.10.1080/09593331908616706Search in Google Scholar

[24] Q. Yuan, R. Sparling, and J.A. Oleszkiewicz. VFA generation from waste activated sludge: effect of temperature and mixing. *Chemosphere* 82(4):603–607, 2011.10.1016/j.chemosphere.2010.10.084Search in Google Scholar
PubMed

[25] A. Donoso-Bravo, W.M.K.R.T.W. Bandara, H. Satoh, and G. Ruiz-Filippi. Explicit temperature-based model for anaerobic digestion: Application in domestic wastewater treatment in a UASB reactor. *Bioresource Technol*. 133:437–442, 2013.10.1016/j.biortech.2013.01.174Search in Google Scholar
PubMed

[26] E. Sánchez, R. Borja, P. Weiland, L. Travieso, and A. Martin. Effect of temperature and pH on the kinetics of methane production, organic nitrogen and phosphorus removal in the batch anaerobic digestion process of cattle manure. *Bioprocess Eng*. 22(3):247–252, 2000.10.1007/s004490050727Search in Google Scholar

[27] D.J. Batstone, J. Keller, I. Angelidaki, S.V. Kalyuzhnyi, S.G. Pavlostathis, A. Rozzi, W.T.M. Sanders, H. Siegrist, and V.A. Vavilin. The IWA anaerobic digestion model no 1 (ADM1). *Water Sci. Technol*. 45:65–73, 2002.10.2166/wst.2002.0292Search in Google Scholar

[28] V.M. Trejos, J. Fontalvo, and M.A. Gómez. Mathematical description and stability analysis of fermentative processes. *Dyna* 76(158):111–121, 2009.Search in Google Scholar

[29] A. Dhooge, W. Govaerts, Yu.A. Kuznetsov, W. Mestrom, A.M. Riet, and B. Sautois. *MATCONT and CL MATCONT: Continuation Toolboxes in MATLAB* Gent University and Utrech University, 2006.10.4249/scholarpedia.1375Search in Google Scholar

[30] J.J. Slotine and W. Li. *Applied Nonlinear Control* Prentice-Hall Englewood Cliffs, New Jersey, 1991.Search in Google Scholar

[31] L. J. S. Allen. *An Introduction to Mathematical Biology* Pearson Education, New Jersey, 2007.Search in Google Scholar

## A Stability analysis of the system

## Theorem 1

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

*where x are the state variables and**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 λ*_{1}*, λ*_{2},· · · *, λ _{n} of A satisfy Re*(

*λ*)< 0.

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

*μ*_{1} with respect to *μ*_{2} with respect to

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}

*λ*+

*a*= 1, · · · ,

_{n}, where the coefficients a_{i}are real constants, i*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 positiveThus, the characteristic equation of the Jacobian matrix related to the system is given by

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*(*λ _{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

## 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

### Equilibrium 1: (washout condition)

The washout condition is given by the absence of biomass in both acidogenic and methanogenic stages, i.e.,

### Equilibrium 2: (washout by acidogenic biomass)

The washout condition is given by biomass absence in the acidogenic stage, i.e.,

and

where *a* = *αD*, *c* = *αDK*_{2}*K _{I}*, and

*μ*

_{2 max}

*ΘI*. Moreover,

_{pH}*S*

_{2}has a real solution if

*a*≠ 0 and the discriminant is positive or zero,

*Δ*≥ 0,

*i.e*.,

*b*

^{2}−4

*ac*≥ 0.

### Equilibrium 3: (washout by methanogenic biomass)

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

where

### Equilibrium 4: (nontrivial solution)

The analytical nontrivial solution is given by

where, *αDK*_{2}*K*_{I}k_{1}, and *ϕ* =

D | dilution rate, day^{−1} |

I_{pH} | inhibition factor by pH, dimensionless |

k_{1} | yield for substrate degradation, g COD/g X_{1} |

k_{2} | yield for VFA production, mmol VFA/g X_{1} |

k_{3} | yield for VFA consumption, mmol VFA/g X_{2} |

K_{1} | half-saturation constant associated with S_{1}, g/L |

K_{2} | half-saturation constant associated with S_{2}, mmol/L |

K_{I} | inhibition constant associated with S_{2}, mmol^{2} /L ^{2} |

pH | hydrogen potential, dimensionless |

pH_{UL} | upper limit of pH, dimensionless |

pH_{LL} | lower limit of pH, dimensionless |

T | operational temperature, C^{∘} |

T_{min} | minimum temperature, C^{∘} |

T_{max} | maximum temperature, C^{∘} |

T_{opt} | optimal temperature, C^{∘} |

S_{1} | outlet organic substrate concentration, g L^{−1} |

inlet organic substrate concentration, g L^{−1} | |

S_{2} | outlet volatile fatty acid concentration,mmol L^{−1} |

inlet volatile fatty acid concentration, mmol L^{−1} | |

inlet acidogenic biomass concentration, g L^{−1} | |

X_{1} | acidogenic biomass concentration, g L^{−1} |

inlet methanogenic biomass concentration, g L^{−1} | |

X_{2} | methanogenic biomass concentration, g L^{−1} |

α | proportion of dilution rate for bacteria, dimensionless |

μ_{1} | acidogenic growth rate kinetic modified, day^{−1} |

μ_{2} | methanogenic growth rate kinetic modified, day^{−1} |

μ_{1 max} | maximum value for acidogenic growth rate kinetic, day^{−1} |

μ_{2 max} | maximum value for methanogenic growth rate kinetic, day^{−1} |

Θ | temperature activity coefficient, dimensionless |

**Received:**2017-03-20

**Revised:**2018-10-04

**Accepted:**2019-01-11

**Published Online:**2019-08-31

© 2020 A.M. Alzate-Ibañez et al., published by De Gruyter

This work is licensed under the Creative Commons Attribution 4.0 International License.