# Dynamics of a Bioreactor with a Bacteria Piecewise-Linear Growth Model in a Methane-Producing Process

Authors: Luz A. Melo Varela; Simeón  Casanova Trujillo; Gerard Olivar Tost
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

This paper shows a study, both analytical and numerical, of a continuous-time dynamical
system associated with a simple model of a wastewater biorreactor. Nonsmooth phenomena
and border-collision local bifurcations appear when the main parameters (dilution and biomass concentration at the inflow) are varied. We apply the Filippov methods following Kuznetsov’s work.


---

## Body

## 1. Introduction

Currently, anaerobic methods are applied to reduce water contamination problems. With these methods, we can reduce the percentages of chemical oxygen demand (COD) and biological oxygen demand (BOD), which are measurements of water quality. Depuration through anaerobic treatments converts organic matter in wastewater into methane (CH4—biogas) and carbonic gas (CO2). Methane can be used as an energetic component because it offers good calorific power, and CO2 can be recirculated to the bioreactor to improve the percentages of biogas yield, thus decreasing organic loads. Organic loads contain high concentrations of organic matter that originates from water circulation in garbage, which dissolves the elements present in it when running through the waste. The result is an environmentally damaging liquid that contaminates the soil and superficial and subterranean waters in their path. For this reason, leachates, among others, are one of the most significant contaminating agents in a landfill as has been extensively discussed in the literature (see, e.g., [1–3]).In this paper, we are mainly interested in leachates since our experimental secondary data comes from this sort of wastewaters. The most used systems for leachate treatment are the so-called high-rate systems, such as the UASB reactor (Figure1). This bioreactor separates different phases: biological (sludge bed), liquid (sludge blanket), and gas (upper section). The wastewater enters the reactor through its lower section and exits through the upper section. The reactor has no filling to support biological growth. The sludge created in the reactor can be divided into two regions: region one, the sludge bed, and region two, the sludge blanket, which is composed of granules or particles in addition to the wastewater, as discussed in [4–6].Figure 1
UASB reactor.The upper section of the reactor contains the solid-liquid-gas (SLG) separator, which prevents the discharge of solids from the reactor and separates them from the produced gas and effluent liquid. This section acts as a sludge sedimenter and gas collector because the gases produced under anaerobic conditions cause internal recirculation, which helps in the creation and sustainment of bacteria. The upper piece (so-called screen) generates a low-turbulence region, where 99% of the sludge in the suspension settles and returns to the reactor. The screen also serves to recover the gas that exits through the center region, as discussed by Kjeldsen et al. [5]. Therefore the SLG separator is fundamental in order to maintain settled sludge, a clarified effluent (gas-free), and properly separated gases.Anaerobic leachates depuration uses a combination of several processes. One of the most important is the so-called anaerobic digestion, which is a fermentation of organic matter. This is performed by bacteria when oxygen is not present. Fermentation subproducts are a mixture of gases (mainly CO2 and CH4, biogas) and also some biomass, which is kept in the process. Anaerobic digestion can be applied to leachates from cattle, forest and agroindustry, and disposal from transformation industries, either one by one, or together (codigestion). The method reported in the paper could be applied not only to leachates but also to a methanisation process that treats glucose or acetate residues.Many mathematical models for bioreactors have been obtained and there is a hugh literature on this topic (see e.g., [7–9]). Also, control of an anaerobic digester through normal form of smooth fold bifurcation has been implemented. This method has faster convergence rate and lower error than traditional methods. The idea is designing a nonlinear controller taking advantage of the knowledge of the bifurcation scenario [10]. Thus bifurcation diagrams and analytical results are important for a good design and control.In this paper we take real data from the leachates in “Esmeralda” landfill, which is located in Manizales (Caldas), a medium-size city in Colombia. Numerical simulations were carried out with standard procedures in MATLAB.The rest of the paper is organized as follows. Section2 is devoted to an overview of Filippov and Kuznetsov theory on nonsmooth systems and nonsmooth bifurcations. This theory will be applied to our model. In Section 3 we describe the model mathematically and we perform some basic algebraic computations following nonsmooth theory. Results are shown in Section 4, and they are compared with analytical computations. Also, nonsmooth bifurcations are reported. Finally, conclusions are stated in Section 5.


## 2. An Overview of Filippov Systems

Nonsmooth systems (continuous piecewise-linear, continuous piecewise-smooth, discontinuous piecewise-smooth, and so on) have been studied in the literature [11–14].Through theory mainly developed by Filippov we can determine the solution of a system ruled by differential equations with discontinuous terms on the right-hand side. According to this method (Filippov method), the borders of all state-velocity vectors within the region of a point on a discontinuous surface must be complemented by a minimum convex set, and the state-velocity vector of the sliding motion must belong to this set, as discussed in [11].For a dynamical system in the state-space where Filippov method can be applied, and assuming only two regions separated by the discontinuity, we can write(1)z˙={F(1)(z),z∈R1,F(2)(z),z∈R2,
where
(2)R1={z∈ℝ2:H(z)<0},R2={z∈ℝ2:H(z)<0}.
The discontinuity boundary Σ separates the two regions R1 and R2 and is given by
(3)Σ={z∈ℝ2:H(z)=0},
where H(z) is a smooth scalar function with a nonzero gradient over Σ. The boundary Σ is a closed set, and we must have that F(1)≠F(2) over Σ.


### 2.1. Sliding Solutions

Following [11], for z∈Σ, we define
(4)σ(z)=〈Hz(z),F(1)(z)〉〈Hz(z),F(2)(z)〉,
where 〈,〉 denotes the standard scalar product, and σ(z) defines the crossing or sliding region. The crossing set Σc⊂Σ is defined by
(5)Σc={z∈Σ:σ(z)>0}
which corresponds to the set of all points z∈Σ, where the two vectors F(i)(z) have nontrivial normal components of identical sign.We also have the sliding setΣs, which complements Σc in Σ:
(6)Σs={z∈Σ:σ(z)≤0}.
The crossing set is open, whereas the sliding set is the union of the sliding closed segments and sliding isolated points.The Filippov method associates the following convex combinationg(z) of the two fields F(i)(z) for each nonsingular sliding point z∈Σ, where g(z) is the so-called the Filippov vector field [11]:
(7)g(z)=λF(1)(z)+(1-λ)F(2)(z),
where λ=〈Hz(z),F(2)(z)〉/〈Hz(z),F(2)(z)-F(1)(z)〉.


## 2.1. Sliding Solutions

Following [11], for z∈Σ, we define
(4)σ(z)=〈Hz(z),F(1)(z)〉〈Hz(z),F(2)(z)〉,
where 〈,〉 denotes the standard scalar product, and σ(z) defines the crossing or sliding region. The crossing set Σc⊂Σ is defined by
(5)Σc={z∈Σ:σ(z)>0}
which corresponds to the set of all points z∈Σ, where the two vectors F(i)(z) have nontrivial normal components of identical sign.We also have the sliding setΣs, which complements Σc in Σ:
(6)Σs={z∈Σ:σ(z)≤0}.
The crossing set is open, whereas the sliding set is the union of the sliding closed segments and sliding isolated points.The Filippov method associates the following convex combinationg(z) of the two fields F(i)(z) for each nonsingular sliding point z∈Σ, where g(z) is the so-called the Filippov vector field [11]:
(7)g(z)=λF(1)(z)+(1-λ)F(2)(z),
where λ=〈Hz(z),F(2)(z)〉/〈Hz(z),F(2)(z)-F(1)(z)〉.


## 3. Bioreactor Mathematical Model

We used the model proposed by Muñoz (2006) [15] for anaerobic digestion in a UASB reactor, but an approximation by straight lines to the bacterial growth model is applied (Monod kinetics). The system, originally smooth, is converted into a nonsmooth one because it is modeled by a system of differential equations with a discontinuous right-hand side. For this approximation method by straight lines, only substrate concentrations from 0 to 6000 (mg/L COD) have been considered. This is because, in addition to achieving a good approximation, the values only have physical sense in this region and are consistent with the design conditions of the UASB bioreactor.The assumptions in the model include that the operation temperature is 20∘C, the constants in the model are those observed by Muñoz (2006) [15] Ks=5522.3, (mg/L), μmax=1.32 d−1,  Y=3.35, (mg COD/mg VSS), and acidogenesis and metanogenesis are considered to be the only processes governed by Monod kinetics, which assumes that the bacterial growth follows Michaelis-Menten kinetics for processes catalyzed by enzymes. Therefore,
(8)μ(S)=μmaxSKs+S,
where Ks is the substrate semisaturation constant.According to this model, for the discussed biological process, the rate of microbial growth will asymptotically tend to the maximum valueμmax.Accounting for the above, the proposed model is(9)S˙=D(Sin-S)-Yμ(S)X,X˙=D(Xin-(1-η)X)+μ(S)X,
where η is the SLG separator efficiency (and thus α=1-η), D is the dilution factor (in d−1) and represents the influent volumetric flow per unit of reactor volume (the inverse of the hydraulic retention time). Sin is the substrate concentration in the input flow (in mg/L COD); S is the substrate concentration in the reactor (in mg/L COD); Y is the substrate yield coefficient (in mg COD/mg VSS); μ(S) is the bacterial growth model; X is the biomass concentration in the reactor (in mg/L VSS); Xin is the biomass concentration in the input flow (in mg/L VSS); and η is the sedimentation efficiency of the separator (SLG).We slightly modify this model by using an approximation toμ(S) by straight lines so that the originally smooth system is converted to a nonsmooth one. We take this approach after observing the experimental data in [15], which resembles much more to piecewise-linear than to a classical smooth Monod model. Zero is the minimum value for S (physically, negative concentrations cannot be observed) and 6000 [in mg/L COD] is the maximum (the maximum substrate concentration value at the input flow based on the operation conditions) [15].Figure2 corresponds to a continuous piecewise-linear approximation, but we will also consider discontinuous piecewise-linear approximations, taking into account the observed data in [15], as slight perturbations to the continuous one.Figure 2
Piecewise-linear continuous approximation of Monod growth model. Nonsmooth points correspond toS1=1800 and S2=4000.Thus we will consider0<p1<p2 as parameters, defining the nonsmooth discontinuous approximation μ-(S) to the Monod curve.We have(10)μ-(S)=μ1S1S,for0≤S<p1,μ-(S)=μ1+μ2-μ1S2-S1(S-S1),forp1≤S<p2,μ-(S)=μ2+μ3-μ2S3-S2(S-S2),forS≥p2,
and thus the final proposed model is the following:
(11)S˙=D(Sin-S)-Yμ-(S)X,X˙=D(Xin-(1-η)X)+μ-(S)X,
where α=1-η corresponds to the SLG separator deficiency, and η is the SLG separator efficiency, μ1=0.32, S1=1800, μ2=0.5544, S2=4000, μ3=0.7, and S3=6000.Each region is ruled by a system of differential equations, which can be discontinuous in the border of each region. Whenp1=p1*=1800, p2=p2*=4000, and p3=p3*=6000, corresponding to the continuous piecewise-linear approximation of the Monod curve, our system of differential equations is piecewise-continuous, and thus no sliding regions are possible. But when parameters pi slightly vary from the corresponding nominal values pi*, then the differential equations have discontinuous right-hand side, and Filippov methods can be applied.However, only analyses of regions corresponding to mode one (0≤S<p1) and mode two (p1<S<p2) were performed since these are the regions where a physically possible dynamic was observed when representing the positive equilibrium points. Thus, equations to be analyzed are (11) in modes one and two.In our case,H(z)=S-p1, and when this is zero, the sliding region will be included in S=p1.The gradient forH(z) is given by ∇Hz(z)=(0,1), which is constant and different from zero for each z∈Σ.When applying the Filippov method, following [11], we can distinguish some critical points.Singular points are z∈Σs such that
(12)〈Hz(z),F(1)(z)-F(2)(z)=0〉.Pseudoequilibriums are z∈Σs such that
(13)g(z)=0,F1,2(z)≠0.Boundary equilibriums are z∈Σs such that
(14)F(1)(z)=0,F(2)(z)=0.Tangent points are z∈Σs such that
(15)〈Hz(z),F(i)(z)〉=0,F1,2(z)≠0.
From now on, we will be interested in the case of discontinuous piecewise-linear approximations, leading to the Filippov solutions.


### 3.1. Algebraic Computations

Boundary-node nonsmooth bifurcations were observed in this model, which are obtained when a node approaches the switching surface and collides with a pseudoequilibrium. This is considered a local bifurcation.We consider again the two vector fields:(16)F(1)={D(Xin-αX)+ASXD(Sin-S)-YASXF(2)={D(Xin-αX)+[μ1+B(S-S1)]XD(Sin-S)-Y[μ1+B(S-S1)]X,
where A=μ1/S1 and B=(μ2-μ1)/(S2-S1).Some basic algebraic computations can be performed for this system and obtain the Filippov vector field. Then, for example, in order to obtain the pseudoequilibriums we have to impose(17)(D(Sin-S)-YX[μ1+B(S-S1)])×(D(Xin-αX)+ASX)-(D(Sin-S)-YASX)×(D(Xin-αX)+[μ1+B(S-S1)]X)=0.
Then we have
(18)X2(YmDα-YASDα)+X(DASK-YmDXin-DmK+YASDXin)=0,
where m=μ1+B(S-S1),  K=Sin-S.Solutions to (18) correspond to pseudoequilibriums.


## 3.1. Algebraic Computations

Boundary-node nonsmooth bifurcations were observed in this model, which are obtained when a node approaches the switching surface and collides with a pseudoequilibrium. This is considered a local bifurcation.We consider again the two vector fields:(16)F(1)={D(Xin-αX)+ASXD(Sin-S)-YASXF(2)={D(Xin-αX)+[μ1+B(S-S1)]XD(Sin-S)-Y[μ1+B(S-S1)]X,
where A=μ1/S1 and B=(μ2-μ1)/(S2-S1).Some basic algebraic computations can be performed for this system and obtain the Filippov vector field. Then, for example, in order to obtain the pseudoequilibriums we have to impose(17)(D(Sin-S)-YX[μ1+B(S-S1)])×(D(Xin-αX)+ASX)-(D(Sin-S)-YASX)×(D(Xin-αX)+[μ1+B(S-S1)]X)=0.
Then we have
(18)X2(YmDα-YASDα)+X(DASK-YmDXin-DmK+YASDXin)=0,
where m=μ1+B(S-S1),  K=Sin-S.Solutions to (18) correspond to pseudoequilibriums.


## 4. Results

We analyse the system given by (11), when the bacterial concentration in the input flow of the bioreactor, Xin,  is varied to be 0, 240, and 320 (in mg/L VSS). Note that since the source was a leachate, the value of Xin is always different from zero, but Xin=0 must be considered because bioreactor wash out can occur. The axis variables in the figures are S (the concentration of the substrate in the reactor) (in mg/L COD) and X (the concentration of the biomass in the reactor) (in mg/L VSS).In the following, we plot a series of figures that resulted from the analytical calculations for the critical points, including equilibrium points, pseudoequilibriums, tangent points, and singular sliding points, which served as references for the numerical analysis.Only figures corresponding toXin=240 were used (the average value of bacterial input in the bioreactor input flow, obtained by Muñoz (2006) [15], in its experimental part), since a similar behavior was observed for the other values.Parameterα seems to be very important since a boundary-node bifurcation occurred when α is varied from 0.15 to 0.18 when Xin=0 or when α changes from 0.19 to 0.24 for Xin=240 and when α also changes from 0.24 to 0.27, for Xin=320. This shows that the higher the biomass amount at the bioreactor input, the lower the sedimentation efficiency of the SLG separator. This lower efficiency blocks a good, previously stabilized sludge recirculation, which  affects the bacterial performance and presents sludge mixture with the treated effluent.When increasing the bacterial concentration at the bioreactor input, it is also expected that there will be different types of generated bacteria. This is due to the effluent coming from the landfill, other reactions, or the presence of toxic substances that avoid the proper functioning of the bioreactor.


### 4.1. Bifurcations with Parameterα

Figures3(a)–3(d) show the analytical results obtained from the algebraic computations. An input flow biomass concentration of 240 (in mg/L VSS) was chosen. The system evolution shows how the equilibrium in region one was moved towards the switching boundary, approaching a collision. This results in a boundary-node bifurcation. However, the pseudoequilibrium also moved within the sliding segment between the tangent points and finally disappears when it collides with the left tangent point.Algebraic computations when parameterα is varied. Equilibriums in *. Pseudoequilibriums in black square. Boundary points in red *. Tangent points in blue ∘. Singular points in green star.
(a)
Analysis withα=0.2
(b)
Analysis withα=0.22
(c)
Analysis withα=0.23
(d)
Analysis withα=0.24Figures4(a)–4(e) show the system evolution when applying the Filippov convex method through numerical simulations. A change in the phase portrait is presented when parameter α is varied, keeping Xin=240 (in mg/L VSS). Figure 4(b) shows that when α had a value of 0.203, the birth of a pseudo-equilibrium was observed in the sliding segment, and there was a node in the upper region. When the value of α was close to 0.23 (Figure 4(d)), the equilibrium in the lower region approached the switching surface and a collision of this equilibrium subsequently occurred. This process is similar to the one described before, generating a boundary-node bifurcation.Numerically computed phase portraits when parameterα is varied.
(a)
Numerical computation withα=0.198
(b)
Numerical computation withα=0.203
(c)
Numerical computation withα=0.22
(d)
Numerical analysis withα=0.23
(e)
Numerical computation withα=0.24Therefore, when increasing the value ofXin, the α value for which the boundary-node bifurcation occurred also increased. This result shows that the SLG separator efficiency is inversely related to the type of water treated, since, at higher bacterial concentrations in the effluent, the SLG sedimenter is less efficient. This fact allows a good organic matter conversion into methane because the suspended sludge does not return to the reactor when a deficient sedimentation occurs, which affects biogas production.


### 4.2. Bifurcations with ParameterD

Bioreactor analysis when varying parameterD was also performed for input flow biomass concentrations Xin with 0, 240, and 320. For 240 and 320 (in mg/L VSS), it was observed that no interesting dynamics were present. However, richer dynamics were obtained for Xin=0 (in mg/L VSS) (Figures 5(a) to 5(e)), where only one stable node exists in the lower region. When increasing parameter D from 2.0 to 2.45, a pseudo-equilibrium is created in the sliding segment, and when D is incresed further, another node appears in the upper region. When D approached 2.40, the equilibrium in the lower region approached the switching surface. Subsequently, a collision of the equilibrium point occurred within this boundary, which led to its catastrophic disappearance.Numerical computations varying parameterD.
(a)
ParameterD=2
(b)
ParameterD=2.05
(c)
ParameterD=2.3
(d)
ParameterD=2.4
(e)
ParameterD=2.45Basins of attraction can also be computed. For example, forD=2.1 and Xin=0 (shown in Figure 6), we observed that the pseudo-equilibrium only attracted the initial points that were in the switching surface, which was also recognized as an unstable sliding segment. The equilibriums within both regions are attractors.Figure 6
Numerically computed basins of attraction for parameterD=2.1.A comparison between numerically and algebraically computed pseudo-equilibriums and the corresponding bifurcations was performed. Error rates were less than 0.02%.


## 4.1. Bifurcations with Parameterα

Figures3(a)–3(d) show the analytical results obtained from the algebraic computations. An input flow biomass concentration of 240 (in mg/L VSS) was chosen. The system evolution shows how the equilibrium in region one was moved towards the switching boundary, approaching a collision. This results in a boundary-node bifurcation. However, the pseudoequilibrium also moved within the sliding segment between the tangent points and finally disappears when it collides with the left tangent point.Algebraic computations when parameterα is varied. Equilibriums in *. Pseudoequilibriums in black square. Boundary points in red *. Tangent points in blue ∘. Singular points in green star.
(a)
Analysis withα=0.2
(b)
Analysis withα=0.22
(c)
Analysis withα=0.23
(d)
Analysis withα=0.24Figures4(a)–4(e) show the system evolution when applying the Filippov convex method through numerical simulations. A change in the phase portrait is presented when parameter α is varied, keeping Xin=240 (in mg/L VSS). Figure 4(b) shows that when α had a value of 0.203, the birth of a pseudo-equilibrium was observed in the sliding segment, and there was a node in the upper region. When the value of α was close to 0.23 (Figure 4(d)), the equilibrium in the lower region approached the switching surface and a collision of this equilibrium subsequently occurred. This process is similar to the one described before, generating a boundary-node bifurcation.Numerically computed phase portraits when parameterα is varied.
(a)
Numerical computation withα=0.198
(b)
Numerical computation withα=0.203
(c)
Numerical computation withα=0.22
(d)
Numerical analysis withα=0.23
(e)
Numerical computation withα=0.24Therefore, when increasing the value ofXin, the α value for which the boundary-node bifurcation occurred also increased. This result shows that the SLG separator efficiency is inversely related to the type of water treated, since, at higher bacterial concentrations in the effluent, the SLG sedimenter is less efficient. This fact allows a good organic matter conversion into methane because the suspended sludge does not return to the reactor when a deficient sedimentation occurs, which affects biogas production.


## 4.2. Bifurcations with ParameterD

Bioreactor analysis when varying parameterD was also performed for input flow biomass concentrations Xin with 0, 240, and 320. For 240 and 320 (in mg/L VSS), it was observed that no interesting dynamics were present. However, richer dynamics were obtained for Xin=0 (in mg/L VSS) (Figures 5(a) to 5(e)), where only one stable node exists in the lower region. When increasing parameter D from 2.0 to 2.45, a pseudo-equilibrium is created in the sliding segment, and when D is incresed further, another node appears in the upper region. When D approached 2.40, the equilibrium in the lower region approached the switching surface. Subsequently, a collision of the equilibrium point occurred within this boundary, which led to its catastrophic disappearance.Numerical computations varying parameterD.
(a)
ParameterD=2
(b)
ParameterD=2.05
(c)
ParameterD=2.3
(d)
ParameterD=2.4
(e)
ParameterD=2.45Basins of attraction can also be computed. For example, forD=2.1 and Xin=0 (shown in Figure 6), we observed that the pseudo-equilibrium only attracted the initial points that were in the switching surface, which was also recognized as an unstable sliding segment. The equilibriums within both regions are attractors.Figure 6
Numerically computed basins of attraction for parameterD=2.1.A comparison between numerically and algebraically computed pseudo-equilibriums and the corresponding bifurcations was performed. Error rates were less than 0.02%.


## 5. Conclusions

When applying the Filippov convex method to a UASB, nonsmooth local boundary-node bifurcations with washout conditions in the reactor were shown. These bifurcations occur when parameterD or α is varied. When increasing the input biomass concentration, parameter α tends to increase the bioreactor inefficiency. This is due to the fact that it had more biomass in the input flow. Predator organisms, such as anaerobic ciliates or chemical products that generate biomass death in the reactor, can exist. In this case the leachates are not properly transformed into biogas. The comparison between the analytical section (algebraically computed from the Filippov vector field) and the numerical approximation yields an error close to 0.02%, which validates the performed calculations. The importance of parameter α was observed in the operation of the UASB bioreactor because the boundary-node bifurcation was present regardless of the biomass concentration at the bioreactor input. This shows that the equilibrium can be controlled with this parameter, either in region one or region two. The obtained results serve as a basis for bioreactor automatic control where a higher decontamination in the treated effluent and an improved conversion of the organic matter to biogas are expected.



---

*Source: 685452-2013-11-14.xml*