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

#### 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 (CH_{4}—biogas) and carbonic gas (CO_{2}). Methane can be used as an energetic component because it offers good calorific power, and CO_{2} 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 (Figure 1). 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].

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 CO_{2} and CH_{4}, 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. Section 2 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 where The discontinuity boundary separates the two regions and and is given by where is a smooth scalar function with a nonzero gradient over . The boundary is a closed set, and we must have that over .

##### 2.1. Sliding Solutions

Following [11], for , we define where denotes the standard scalar product, and defines the crossing or sliding region. The crossing set is defined by which corresponds to the set of all points , where the two vectors have nontrivial normal components of identical sign.

We also have the sliding set , which complements in : 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 combination of the two fields for each nonsingular sliding point , where is the so-called the Filippov vector field [11]: where .

#### 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] , (mg/L), d^{−1}, , (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,
where 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 .

Accounting for the above, the proposed model is
where is the SLG separator efficiency (and thus ), 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). is the substrate concentration in the input flow (in mg/L COD); is the substrate concentration in the reactor (in mg/L COD); is the substrate yield coefficient (in mg COD/mg VSS); is the bacterial growth model; is the biomass concentration in the reactor (in mg/L VSS); 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 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 (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].

Figure 2 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.

Thus we will consider as parameters, defining the nonsmooth discontinuous approximation to the Monod curve.

We have and thus the final proposed model is the following: where corresponds to the SLG separator deficiency, and is the SLG separator efficiency, , , , , , and .

Each region is ruled by a system of differential equations, which can be discontinuous in the border of each region. When , , and , 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 slightly vary from the corresponding nominal values , then the differential equations have discontinuous right-hand side, and Filippov methods can be applied.

However, only analyses of regions corresponding to mode one () and mode two () 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, , and when this is zero, the sliding region will be included in .

The gradient for is given by , which is constant and different from zero for each .

When applying the Filippov method, following [11], we can distinguish some critical points.

* Singular points* are such that
*Pseudoequilibriums* are such that
*Boundary equilibriums* are such that
*Tangent points* are such that
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: where and .

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 Then we have where , .

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, , is varied to be , , and (in mg/L VSS). Note that since the source was a leachate, the value of is always different from zero, but must be considered because bioreactor wash out can occur. The axis variables in the figures are (the concentration of the substrate in the reactor) (in mg/L COD) and (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 to 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 to when or when changes from to for and when also changes from to , for . 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

Figures 3(a)–3(d) show the analytical results obtained from the algebraic computations. An input flow biomass concentration of (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.

(a) Analysis with |

(b) Analysis with |

(c) Analysis with |

(d) Analysis with |

Figures 4(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 (in mg/L VSS). Figure 4(b) shows that when had a value of , 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 (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.

(a) Numerical computation with |

(b) Numerical computation with |

(c) Numerical computation with |

(d) Numerical analysis with |

(e) Numerical computation with |

Therefore, when increasing the value of , 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 Parameter

Bioreactor analysis when varying parameter was also performed for input flow biomass concentrations with , , and . For and (in mg/L VSS), it was observed that no interesting dynamics were present. However, richer dynamics were obtained for (in mg/L VSS) (Figures 5(a) to 5(e)), where only one stable node exists in the lower region. When increasing parameter from to , a pseudo-equilibrium is created in the sliding segment, and when is incresed further, another node appears in the upper region. When approached , 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.

(a) Parameter |

(b) Parameter |

(c) Parameter |

(d) Parameter |

(e) Parameter |

Basins of attraction can also be computed. For example, for and (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.

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

#### Acknowledgments

This research was funded through project *Dynamic analysis of a continuous stirred tank reactor (CSTR) with biofilm formation for the treatment of residual waste*, Support Program for postgraduate theses, DIMA 2012 of Universidad Nacional de Colombia, and Project no. 20201007075 (VRI Universidad Nacional de Colombia).