ANALYSIS OF AN ANAEROBIC DIGESTION MODEL IN LANDFILL WITH MORTALITY TERM

. We study a mathematical model of anaerobic digestion with biomass recirculation, dedicated to landﬁll problems, and analyze its asymptotic be- havior. We show that the global attractor is composed of an inﬁnity of non-hyperbolic equilibria. For non-monotonic growth functions, this set is non connected, which impacts the performances of the bioprocess.

1. Introduction. The anaerobic digestion process is a natural biological process of decomposition of organic matter by microorganisms (bacteria) that are activated under anaerobic conditions, that is to say without oxygen. It is characterized by a succession of complex reactions both in parallel and in series. In the long term, the organic matter is transformed into biogas, a mixture mainly composed of methane and carbon dioxide. The main stages of this process are hydrolysis, acidogenesis, acetogenesis and methanogenesis. Models such as the "ADM1" allow a detailed description of this process, cf. [14]. However, such complex models are not well suited for mathematical analysis. It is why a number of simpler models have been investigated in the literature over theses last years [8,11,23,12,21,22]. When dealing with the digestion of wastewater, it is recognized that the limiting step is the methanoganesis. In such a case, modelling include one-, two-or three-step models. Of particular interest is the model by Bernard et al., 2001 ([6]) which proposes to model the anaerobic digestion process as a two-step process involving both the acidogenesis (using a Monod kinetics) and the methanogenesis processes (using a Haldane kinetics). This model, for which the mathematical analysis has been conducted by Benyahia et al., 2012, cf. [5], is very popular, notably for control purposes, since it remains of moderate complexity while being quite easy to calibrate to predict process behaviour with satisfying performances, [9,2] or still [15]. It has also been the basis for proposing a systematic way to link simple models to ADM1 predictions, cf. [14]. In this work, we consider the model proposed by M. Rouez [23]. He is considering the anaerobic process in two steps that are Hydrolysis and Methanogenesis. Such description is well suited as long as the acidogenesis is not the limiting step, which is the case for the digestion of solid waste, as considered in this paper. In landfills, there is always a part of the raw material that has no access to oxygen. In addition, we consider here mortality. When a biological process is running in continuous mode, it is the rule rather than the exception not to consider mortality terms, cf. for instance [16] or [17], and such term can have a great influence on the process behaviour, cf. for instance [13]. However, when working in a closed environment -as it is the case in landfills -the mortality can no longer be neglected: part of the mortality of microorganism then returns to the slowly biodegradable material, itself being further hydrolyzed in rapidly biodegradable material. The role of the mortality term and the growth function on the performances of the overall process, notably in terms of biogas production, is of primer importance in landfill applications. The main objective of the present paper is to give deeper insights and predictions of this role with the help of a mathematical model. The paper is organized as follows. In the next section, model and assumptions are presented and discussed. Then, the dynamics is mathematically analyzed in the following section. In particular, we show that there exists an infinity of equilibria to which solutions converge, depending on initial conditions. Finally, numerical simulations with different value of the death rate and various initial conditions are presented in a dedicated section, before discussions and conclusions are drawn.

2.
Model and assumptions. The hydrolysis step of transformation of organic matter of concentration is an important phase in the biodegradation process because it is a substrate preparation step. It is modeled by the first order equation where X(t) is the concentration of organic matter at time t and K h is the hydrolysis constant.
On the other hand, the methanogenesis is the last step in the anaerobic digestion process and leads to the production of biogas. Based on the principle of mass conservation and the fundamental relations of biological kinetic (that are the growth rate of bacteria and the use rate of the substrate), this step is modeled as follow where S(t) and B(t) represent the concentrations of soluble organic matter and methanogenic biomass at time t. The symbols µ, Y , and K d are successively the specific growth rate, the rate of use of the substrate and the mortality rate. Parameters f 2 and (1 − f 2 ) are the stoichiometric coefficients, representing the parts of soluble organic matter transformed into carbon dioxide and methane during methanogenesis step, of cumulative concentrations at time t denoted by the variables [CO 2 ](t) and [CH 4 ](t).
In the present work, we consider that during the process the death of methanogenic bacteria constitute a substrate (given by a proportion αK d B of the death biomass, where α is a constant parameter) to the hydrolysis step, which brings a source term in the equation (1). This means that only a fraction of the biomass mortality is used again as a substrate in the methanogenesis step, as described in [6,4,10]). This couples the hydrolysis and methanogenesis steps as represented in the model below, which makes the mathematical analysis of the model not straightforward. In addition, we introduce the stoichiometric coefficients f 1 and (1 − f 1 ) to represent the parts of biomass X transformed into organic matter S and carbon dioxide CO 2 during the "hydrolysis / acidogenesis" process. The overall process is depicted on Figure 1. It is modeled by the following dynamic system (that was already proposed by Rouez in [23] but for a specific function µ and coefficient α = 1).
The total biogas G(t) produced by the process at time t is the sum of carbon dioxide [CO 2 ](t) and methane [CH 4 ](t). Notice that in absence of initial substrate S or initial biomass B, the model (2) does not produce any biogas. We shall see that this is not the case for model (3), provided that there is initial matter X.
The specific growth rate function satisfy the following properties.
Many functions that satisfy this hypothesis are met in the literature. The most popular ones are (see Figure 2) : 1. the Monod law [18], which is related to a growth saturation or limitation : where µ max is the maximum growth rate and K S the half-saturation constant. 2. the Haldane law [1], which is also characterized by an inhibition phenomenon for large values of the substrate concentration: where K I is the inhibition constant. The Haldane expression is often considered to be more appropriate to the anaerobic process. Notice however that the Haldane function can be seen as a generalization of the Monod one on a fixed interval for large values of K I .
, f 1 and f 2 fulfill the following conditions. 1. The proportion of nutrient recycling α cannot exceed 1 2. The mortality rate K d is a positive parameter which is below the maximum growth rate 3. The rate of use of the substrate is a strictly positive parameter such that 4. The stoichiometric coefficients parameters f 1 and f 2 are strictly positive and satisfy 0 < f 1 < 1 and 0 < f 2 < 1 (9) Condition (7) means that the choice of bacteria and operating conditions, which impact the growth and death rates (such as temperature, pH...), are such that the bacterial growth is possible. For convenience, we define the set which has non-empty interior, under Hypotheses 1 and 2. For the Monod expression, one has and for the Haldane expression, the set E has two connected components: which is positive under Hypothesis 2.
3. Study of the asymptotic behavior. The dynamics (3) has a cascade structure. We thus begin by the study of the asymptotic behavior of the sub-system (3a). Proof. Let us first show that for any initial condition (X 0 , S 0 , B 0 ) in R 3 + , the forward solution (X(t), S(t), B(t)) remains in the non-negative orthant R 3 + . • If B 0 = 0, then the solution verifies B(t) = 0 for any t > 0, whatever are X 0 , S 0 . Therefore, by uniqueness of solutions of the Cauchy problem, a solution B(·) cannot cross the plane B = 0 neither reaches 0 in finite time. Then one has for any t ≥ 0, from which one deduces the inequality X(t) ≥ X 0 e −K h t for any t ≥ 0. This proves that the solution X(t) is non-negative for any positive t, and cannot reach 0 in finite time, whatever is ( for any t ≥ 0. By comparison of solutions of scalar ordinary differential equations (see e.g. [24]), one has S(t) ≥ S(t) for any t > 0, where S is solution of the differential equation dS dt = − 1 Y µ(S)B(t) with S(0) = S 0 . As µ(0) = 0, S(t) = 0 for any t > 0 is solution when S 0 = 0. By uniqueness of the solutions (the function µ being Lipschitz continuous), a solution S(·) cannot cross the plane S = 0 neither reaches 0 in finite time.
We deduce that S(t) is also non-negative for any t > 0, and cannot reaches 0 in finite time, whatever is (X 0 , S 0 , B 0 ) in R 3 + . We show now the system (3a) is dissipative. Consider the "storage" function From equations (3a), one has By Hypothesis 2, one has αY f 1 < 1 and as X(·) is non-negative, we deduce that V is non-increasing. Being bounded from below by 0, V (t) converges to a limit V ∞ ≥ 0 when t tends to +∞. This implies that the variables X(t), S(t) and B(t) are bounded. Moreover, one has from which we deduce that d 2 V /dt 2 is bounded and thus dV /dt is uniformly continuous on R + . By Barbalat's Lemma [3], we obtain that dV (t)/dt converges to 0 when t tends to +∞. Therefore, we obtain lim t→+∞ X(t) = 0.
Finally, as V (t) converges to a limit when t tends to +∞, we obtain Let us show that S cannot be equal to 0 when X 0 or S 0 is non null. If not, S(t) has to tends to zero and therefore there exists T > 0 such that µ(S(t)) ≤ f 1 Y αK d for any t > T . Then one has When f 1 X(0) + S(0) > 0, the variable X(·) or S(·) cannot reaches 0 in finite time, as shown previously, and one has then f 1 X(t) + S(t) ≥ f 1 X(T ) + S(T ) > 0 for any t > T , which is in contradiction with the fact that f 1 X(t) + S(t) tends to 0 when t tends to +∞. Finally, from the fact that V is non increasing, we obtain Let us now characterize the set of equilibria that can be reached by positive solutions.
Proof. If there exists a solution of (3a) such that S(·) converges asymptotically to S > 0 that does not belong to E, there exists T > 0 such that one has but then the solution B(·) > 0 verifies With standard comparison theorem for differential inequalities (see e.g. [24]), we conclude that B(·) cannot converge asymptotically to 0, and thus a contradiction with the results of Proposition 1.
The system (3a) admits a continuum of equilibria which are not hyperbolic. One cannot conclude about their stability by studying the single linearization of the dynamics. The Center Manifold Theorem [7] could be used but it turns out to be not enough informative for our problem. However, we obtain the following result.
Proposition 3. Assume that Hypotheses 1 and 2 are fulfilled. For each steady state E = (0, S , 0) with S ∈ int E, there exists an invariant two-dimensional manifold M in R 3 + such that any solution of (3a) with initial condition in M converges asymptotically to E.
Proof. Let us fix S > 0 such that µ(S ) < K d . For solutions with B(t) > 0, t ∈ [0 + ∞), consider the variable Then, a direct computation gives We can then write the system (3a) in R 2 + × R + equivalently on the domain with g(Z, X, B) := µ(S − f 1 X + (Z − β)B).
One can check that the domain D is (positively) invariant by the dynamics (13). Moreover, the dynamics (13) is well defined for B = 0 and regular on the set D (which is also invariant). Any trajectory (Z(·), X(·), B(·)) in D matches a trajectory (X(·), S(·), B(·)) of (3a) where However, trajectories of (13) with B(·) = 0 does not necessarily match a trajectory of (3a). As µ(S ) = K d , one can check that 0 is the only equilibrium of (13) in D.
The Jacobian matrix at 0 is given by Under the condition µ(S ) < K d , 0 is thus an hyperbolic equilibrium with a two-dimensional stable manifold S and a one-dimensional unstable manifold U (see e.g. [19]). Clearly, one has U = R × {0} × {0}. Moreover, any solution of (13) with B(0) = 0 verifies B(t) = 0 for any t > 0 and lim t→+∞ X(t) = 0, which implies that the function converges to µ(S ) when t → +∞, whatever is Z(0). Therefore, for any Z(0), X(0), there exists T > 0 such that and then one obtains from (13) that Z(t) tends to infinity when t → +∞. We deduce that any point (Z, X, B) ∈ S \ {0} has to be such that B = 0. Finally, we conclude that any trajectory of (13) in the two-dimensional invariant manifold M := S ∩ D matches a trajectory of (3a) and converges asymptotically to 0, that is (X(t), S(t), S(t)) converges asymptotically to (0, S , 0). Propositions 2 and 3 together allow to state the following result. Finally, let us consider the sub-system (3b).
Proposition 4. Under Hypotheses 1, 2, for any non-negative vector (X 0 , S 0 , B 0 ), the solutions of (3) verify where S is the asymptotic value of S(·) and the coefficients a, b, c, d are positive numbers given by the following expressions.
Proof. Consider again the function V defined in (11). Equation (12) gives From the first equation of (3a), we obtain also and with (16), we obtain Then, from the last equation of (3a), one has We conclude, from equations (3b), that [CO 2 ](t) and [CH 4 ](t) converges asymptotically to finite values. The integration of equations (3a) between t = 0 and t = +∞ leads to the following system of equations whose solution is unique, given by the expressions because 1−αY f 1 > 0 under Hypothesis 2, which provide the formulae (14), (15).

Remark 1.
One can check that when B 0 = 0, there is no production of methane (as expected) because the trajectory stays on the half line and one has then cX 0 + d(S 0 − S ) = 0. The production of dioxide is then equal to (a + αb)X 0 whatever is S 0 .
When the set E is not connected, as an union of disjoint intervals E = I 1 ∪I 2 ∪· · · , the state space can be split into a family {B i } of attraction basins of the subsets {0} × I i × {0}. These basins conduct the system to different levels of performances. For instance, with the Haldane law, there are two basins B − , B + leading to equilibria with S in [0, λ − ] or in [λ + , +∞). The separating surface S of B − and B + is numerically investigated in the next section. 4. Application. Proposition 4 shows that the production of the total biogas G is impacted by the final value S of the remaining soluble matter, which is itself related to the death rate K d (as S belongs to the set E).
For the Haldane growth function, the process leads either to a relatively large production of biogas (when S < λ − ) or to a relatively small production of it (when S > λ + ), depending on the initial condition. The difference between these two situations get larger when the death rate K d is small, differently to the Monod case which is more robust.
Because the Haldane function is often more realistic (and the Monod function can be seen as a particular case for large values of the parameter K I ), we present here simulations only for the Haldane one, with µ m = 0.3, K S = 160 and K I = 10 and other model parameters reported in the following Table.   Table 1.
With these parameters, the corresponding Haldane function is plotted in Figure 3.  Let us fix the death rate K d = 0.02. For such value, we obtain λ − = 12.5 and λ + = 127.4 (see equations (10)). Consider S 0 = 0 and B 0 = 2. Depending on the initial condition of X, S converges either towards a value smaller than λ − or larger than λ + . Such a behaviour can be seen with five initial conditions of X in Figure  4 (X 0 = 340 to X 0 = 360 with an increment of 5). As expected, there is a sudden change on the asymptotic value of S when X 0 passes a threshold: for X 0 = 340, 345 and 350, S converges towards values smaller than λ − while it converges towards values larger than λ + for X 0 = 355 and 360.
As K d increases, λ − increases while λ + decreases and for very small values of K d , λ + may be very large. In Figure 5, we plotted the value of the threshold on X 0 for which the asymptotic value S is under λ − or over λ + for a range of values  Table 1 The biogas production can be seen as an increasing function of X 0 as long as S converges under λ − (since asymptotically there remains a low amount of substrate) and an increasing one as soon as the threshold on X 0 is crossed over (when S converges over λ + ). In Figure 6, we plotted the total production of biogas production for K d = 0.02 as a function of X 0 from 300 to 400 (recall that in this case the threshold value for X 0 is about 357), which shows the discontinuity.  Figure 6. Biogas production as a function of X 0 for K d = 0.02 Finally, we have plotted in Figure 7 the values λ − , λ + as functions of K d , which show the distance between the two attractors. In addition, the biogas productions obtained on both sides of the switching value of X 0 , denoted by Biogas − and Biogas + , have been plotted. As expected, for small values of K d , the difference between the two extremes values of the biogas productions are high. This study reveals interesting insights for practitioners. 1. When S 0 is small (or null), considering that an inhibition cannot occur because it concerns only high values of S , and adopting a Monod function instead of an Haldane one (or a large value of K I ), could lead to wrong predictions. The biogas production could be poor under high initial concentration of organic matter X 0 , because it conducts the system to large asymptotic values of S (see Figures 4 and 6). 2. When the death rate K d is high (but not too much to fulfill condition (7)), the performances on biogas production are quite sensitive to the initial load of organic matter X 0 , because the threshold on X 0 is low (see Figure 5). 3. When the death rate K d is low, the system has good performances on the condition that the initial load X 0 does exceed the threshold, which is relatively large. However, if it is exceeded, the performances collapse dramatically (see Figure 7).

5.
Conclusion. With the help of Barbalat's Lemma and the Stable and Unstable Manifolds Theorem, we have shown that each trajectory of the system is bounded and converges to one of the non-hyperbolic equilibria. For non-monotonic growth functions, such as the Haldane function, the global attractor is non-connected. In this case, we have shown that the performances in terms of biogas production is discontinuous with respect to the initial condition. For the Haldane law, a too high initial load of organic matter could be penalizing because there could be a significant quantity of residual soluble matter, especially when the death rate is small. A strategy could be then to fractionate the load of organic matter over the time, i.e. have a smaller initial load and re-introduce the remaining quantity of matter to be treated later (in one or several times). This leads to a control problem for choosing optimally the proportion of splitting and the corresponding re-introduction time(s) to obtain the best performances. This could be the matter of a future work. Another possibility is to play with the recirculation rate (in a completely mixed system), as in [20].
Moreover, the introduction of spatial heterogeneity in the model, in terms of a system of reaction diffusion p.d.e. instead of a system of o.d.e. would be a reasonable extension of the present work, taking into consideration that spatial heterogeneity is usually observed in real landfills. One expects to find an infinite number of stationary solutions, as for the o.d.e. model. A relevant question of interest for practitioners would be then to study the effect of the diffusion on the performances of the biogas production, in particular for the case of the Haldane growth function.