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

1 Departamento de Matemáticas y Estadı́stica, Facultad de Ciencias Exactas y Naturales, Universidad Nacional de Colombia-Sede Manizales, Carrera 27 No. 64-60 Manizales, Caldas, Colombia 2Departamento de Ingenieŕıa Eléctrica, Electrónica y Computación, Percepción y Control Inteligente-Bloque Q, Facultad de Ingenieŕıa y Arquitectura, Universidad Nacional de Colombia-Sede Manizales, Bloque Q, Campus La Nubia, Manizales 170003, Colombia


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 4biogas) 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][2][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][5][6].
The upper section of the reactor contains the solid-liquidgas (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 (socalled 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 socalled 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][8][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.
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  1 and  2 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  (1) ̸ =  (2) over Σ.

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]   = 5522.3,(mg/L),  max = 1.32 d −1 ,  = 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, 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  max .
Accounting for the above, the proposed model is where  is the SLG separator efficiency (and thus  = 1 − ),  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). in 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);  in 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 0 <  1 <  2 as parameters, defining the nonsmooth discontinuous approximation () to the Monod curve.
We have and thus the final proposed model is the following: where  = 1 −  corresponds to the SLG separator deficiency, and  is the SLG separator efficiency,  1 = 0.32,  1 = 1800,  2 = 0.5544,  2 = 4000,  3 = 0.7, and  3 = 6000.Each region is ruled by a system of differential equations, which can be discontinuous in the border of each region.When  1 =  * 1 = 1800,  2 =  * 2 = 4000, and  3 =  * 3 = 6000, corresponding to the continuous piecewiselinear 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 (0 ≤  <  1 ) and mode two ( 1 <  <  2 ) 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, () =  −  1 , and when this is zero, the sliding region will be included in  =  1 .
The gradient for () is given by ∇  () = (0, 1), which is constant and different from zero for each  ∈ Σ.
When applying the Filippov method, following [11], we can distinguish some critical points.

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  =  1 / 1 and  = ( 2 −  1 )/( 2 −  1 ).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  =  1 + ( −  1 ),  =  in − .

Results
We analyse the system given by (11), when the bacterial concentration in the input flow of the bioreactor,  in , is varied to be 0, 240, and 320 (in mg/L VSS).Note that since the source was a leachate, the value of  in is always different from zero, but  in = 0 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  in = 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  in = 0 or when  changes from 0.19 to 0.24 for  in = 240 and when  also changes from 0.24 to 0.27, for  in = 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.

Bifurcations with Parameter 𝛼. Figures 3(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.
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 = 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 boundarynode bifurcation.Therefore, when increasing the value of  in , 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.

Bifurcations with Parameter 𝐷.
Bioreactor analysis when varying parameter  was also performed for input flow biomass concentrations  in 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  in = 0 (in mg/L VSS) (Figures 5(a) to 5(e)), where only one stable node exists in the lower region.When increasing parameter  from 2.0 to 2.45, a pseudoequilibrium is created in the sliding segment, and when  is incresed further, another node appears in the upper region.When  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.
Basins of attraction can also be computed.For example, for  = 2.1 and  in = 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.
A comparison between numerically and algebraically computed pseudo-equilibriums and the corresponding bifurcations was performed.Error rates were less than 0.02%.

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 boundarynode bifurcation was present regardless of the biomass
Analysis with  = 0.24

Figure 3 :
Figure 3: 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.