Next Article in Journal
Prediction of Cutting Material Durability by T = f(vc) Dependence for Turning Processes
Next Article in Special Issue
Not Just Numbers: Mathematical Modelling and Its Contribution to Anaerobic Digestion Processes
Previous Article in Journal
Development of Indicator of Data Sufficiency for Feature-based Early Time Series Classification with Applications of Bearing Fault Diagnosis
Previous Article in Special Issue
Effect of Delays on the Response of Microalgae When Exposed to Dynamic Environmental Conditions
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Mathematical Modeling and Stability Analysis of a Two-Phase Biosystem

1
Institute of Mathematics and Informatics, Bulgarian Academy of Sciences, Acad. G. Bonchev Str., Block 8, 1113 Sofia, Bulgaria
2
Institute of Microbiology, Bulgarian Academy of Sciences, Acad. G. Bonchev Str., Block 26, 1113 Sofia, Bulgaria
*
Author to whom correspondence should be addressed.
Submission received: 24 April 2020 / Revised: 30 June 2020 / Accepted: 2 July 2020 / Published: 6 July 2020
(This article belongs to the Special Issue Modelling and Optimal Design of Complex Biological Systems)

Abstract

:
We propose a new mathematical model describing a biotechnological process of simultaneous production of hydrogen and methane by anaerobic digestion. The process is carried out in two connected continuously stirred bioreactors. The proposed model is developed by adapting and reducing the well known Anaerobic Digester Model No 1 (ADM1). Mathematical analysis of the model is carried out, involving existence and uniqueness of positive and uniformly bounded solutions, computation of equilibrium points, investigation of their local stability with respect to practically important input parameters. Existence of maxima of the input–output static characteristics with respect to hydrogen and methane is established. Numerical simulations using a specially elaborated web-based software environment are presented to demonstrate the dynamic behavior of the model solutions.

1. Introduction

Anaerobic digestion (AD) is a multi-step biotechnological process widely applied for the treatment of wastes and wastewater, where usually hydrogen ( H 2 ) is not accumulated as an intermediary product [1,2,3]. Recently, the interest in H 2 production by modifying operating conditions of the AD has increased. This process has been denominated as dark fermentative H 2 production in contrast with light derived H 2 producing processes [4,5,6,7,8]. In traditional AD, H 2 is not detected as it is consumed immediately e.g., by hydrogenotrophic methanogens to produce methane C H 4 and carbon dioxide C O 2 [9]. Conversely, H 2 can be produced as main component in biogas by engineering the process conditions in a way that favors Clostridium and Enterobacter metabolism against the predominance of archaeal microflora [10,11]. However, the main limitation of dark fermentative H 2 production is the rather low energy recovery, cf. [7,8,12]. In order to completely utilize the organic acids produced during dark fermentation and improve the overall energy conversion efficiency, a two-stage AD (TSAD) concept consisting of a hydrogenic process followed by a methanogenic process has been suggested (see e.g., [12,13,14]).
A lot of models describing separately the fermentative hydrogen production (see, e.g., [15,16]) and the AD for methane production (cf. e.g., [17,18,19,20]) are known. However, only a few models of TSAD processes are available; see, for example, [21,22,23,24] and the references therein. In [25], both bio-hydrogen and bio-methane productions are optimized in two bioreactors regarding two operating parameters, organic loading rate, and dilution rate; however, the first bioreactor is running in semi-continuous and the second—in batch operation mode because the working volumes of both bioreactors are equal. In [26], steady-state analysis for the TSAD process is presented using a very simple model without hydrolysis and without taking into account the most important fact for this scheme—the energy recovery.
In this paper, we propose a mathematical model, describing the process of simultaneous H 2 and C H 4 production by AD of organic wastes in a cascade of two connected continuously stirred bioreactors. The biochemical processes in the first bioreactor include disintegration of organic wastes, hydrolysis, and acidogenesis with hydrogen production. The methane production from acetate (methanogenesis) is separated in the second bioreactor. The proposed model is developed by reducing the well known Anaerobic Digester Model No 1 (ADM1) basic structure elaborated by the International Water Association (IWA) [17]. Modeling of TSAD using ADM1 is also presented in [27], where, however, methane is obtained from both bioreactors.
Our mathematical model is presented for the first time in the conference paper [21]. The present article represents an extension of the latter, providing (i) more detailed assumptions used to construct the model, which are based on practical experiments (Section 3); (ii) more precise and deeper theoretical investigations of the model dynamics such as existence and uniqueness of positive and uniformly bounded solutions, more precise computation of practically important equilibrium points of the model and study of their local asymptotic stability (Section 4); (iii) more detailed investigation concerning the optimization of hydrogen and methane production on steady state operation (Section 5), and showing the existence of maxima, which is important for the applications; (iv) numerical simulations in Section 6, based on a newly developed, so called one-page application “(H2,CH4) Bioreactors ver. 1.2” in the SmoWeb platform.
The main goal of this article is to present rigorous mathematical analysis of the proposed model dynamics, to detect and provide some model-based predictions, which could be useful in designing and engineering real laboratory experiments. The theoretical studies are supported and complimented by a variety of computer simulations.

2. Process Description

In the TSAD system, relatively fast growing acidogens and H 2 producing microorganisms are developed in the first-stage hydrogenic bioreactor ( B R 1 with working volume V 1 , see Figure 1) and are involved in the production of volatile fatty acids (VFA–mainly acetate and butyrate, [8]) and H 2 . Some recent studies on the microbial community in the hydrogen producing B R 1 show that dominant are microorganisms of the genus Clostridium [28]. In B R 1 , p H is kept in the range 5.0–5.5, and acetate degrading methanogens are eliminated. On the other hand, the slow growing acetogens and methanogens are developed in the second-stage methanogenic bioreactor ( B R 2 with working volume V 2 , and p H in the range 6.5–8.5), in which the produced VFA are further converted to C H 4 and C O 2 by methane-forming bacteria from the genus Methanobacterium, Methanococcusect, see, for details, [9] and the references therein.
It is known that in the two-stage H 2 + C H 4 system the energy yields are up to 43% more in comparison to the traditional one-stage C H 4 production process (cf. [8,25]).
Assume that the volumes V 1 and V 2 of the bioreactors are constant. Let F 1 and F 2 be the inflows in the first and second bioreactor, respectively, and let F 1 = F 2 = F be valid. It is well known that the dilution rates D 1 and D 2 are defined as
D 1 = F V 1 , D 2 = F V 2 .
Then,
D 2 = V 1 V 2 D 1 : = γ D 1 within γ : = V 1 V 2 .
Later on in the paper (see Section 5), by analyzing points of maximum biogas production in both bioreactors, biohydrogen in B R 1 and biomethane in B R 2 , we find theoretical values for the hydraulic retention time (HRT) 1 / D 1 , max and 1 / D 2 , max , which show that the volume V 2 of the second bioreactor for methane production should be larger than the volume V 1 of the first bioreactor. Therefore, γ < 1 should be valid. We determine the constant γ using the proposed model equations.

3. Model Description

In this section, we present the mathematical model describing the process of simultaneous H 2 and C H 4 production by AD of organic wastes in a cascade of two continuously stirred bioreactors. The model is derived on the basis of the ADM1 basic structure as well as on our experience with the two-phase AD process with hydrogen and methane production, see [29]. The ADM1 is a complex model, describing all known (till its creation) microbiological and biochemical reactions occurring in AD of organic wastes. ADM1 is suitable mainly for process knowledge and simulation, but not appropriate for theoretical studies and model-based control design due to the plenty of input parameters which are difficult to obtain. The reduction of ADM1 is done here in order to simplify the latter and to provide an opportunity for studying the asymptotic stability of the dynamics with respect to two essential indicators: (i) the separation of the model equations into two groups, corresponding to the reactions taking place in each one of the two bioreactors, and (ii) the type of the organic wastes used.
For simplification of the model, the following assumptions have been accepted:
  • In B R 1 only the main VFA (acetate) is formed, whereas propionate and butyrate are omitted. In this way, the slow acetogenic phase has been omitted and consequently in B R 2 only methanogenesis occurs.
  • In both bioreactors, the existence of inhibition phenomena is not taken into account.
  • Balance equations of the hydrogen and methane in the liquid phases have been neglected because they are practically not dissolved in liquids.
  • Hydrogenotrophic bacteria, producing methane from hydrogen, are eliminated previously and do not exist in this process.
  • Equations describing balances of inorganic components and some biochemical equations have been neglected in view of simplifying the model.
  • p H is considered in the range between 5.0 to 5.5 in the first bioreactor B R 1 , and in the interval 6.5–8.5 in the second bioreactor B R 2 . Within these ranges, it is assumed that the negative influence of p H could be disregarded.
These ranges of p H were obtained during our own laboratory experiments with TSAD of waste from maize treatment (not published yet).
The proposed models of the two bioreactors are of balancing type for continuously stirred tank reactors (CSTR) like the original ADM1 model. For each bioreactor, the ADM1 equations are adapted to describe the dynamics of the most important variables for the AD of organic waste.
Following the above assumptions, we present now the mathematical model for hydrogen and methane production in the two bioreactors. The definitions of the phase variables and parameters in the model equations are given in Table 1 and Table 2 below.
The dynamics in the first bioreactor are described by the following set of 10 nonlinear ordinary differential equations:
d d t S s u ( t ) = D 1 ( S s u i n S s u ) + k h y d , c h X c h + f s u , l i k h y d , l i X l i μ s u a a , s u X s u a a
d d t S a a ( t ) = D 1 ( S a a i n S a a ) + k h y d , p r X p r μ s u a a , a a X s u a a
d d t S f a ( t ) = D 1 ( S f a i n S f a ) + f f a , l i k h y d , l i X l i μ f a X f a
d d t S a c ( t ) = D 1 S a c + ( 1 Y s u a a ) f a c , s u μ s u a a , s u + f a c , a a μ s u a a , a a X s u a a + 0.7 ( 1 Y f a ) μ f a X f a
d d t X c ( t ) = D 1 ( X c i n X c ) k d i s X c
d d t X c h ( t ) = ( D 1 + k h y d , c h ) X c h + f c h , x c k d i s X c
d d t X p r ( t ) = ( D 1 + k h y d , p r ) X p r + f p r , x c k d i s X c
d d t X l i ( t ) = ( D 1 + k h y d , l i ) X l i + f l i , x c k d i s X c
d d t X s u a a ( t ) = D 1 + Y s u a a ( μ s u a a , s u + μ s u a a , a a ) X s u a a
d d t X f a ( t ) = ( D 1 Y f a μ f a ) X f a ,
with gaseous output
Q h 2 = ( Y h 2 , s u μ s u a a , s u + Y h 2 , a a μ s u a a , a a ) X s u a a + Y h 2 , f a μ f a X f a ,
where Q h 2 is the hydrogen flow rate [ L / h ] . In (2)–(11), the biomass specific growth rates μ s u a a , s u , μ s u a a , a a and μ f a are of the form
μ s u a a , s u = μ s u a a , s u ( S s u , S a a ) = k m , s u a a S s u k s , s u a a + S s u · S s u S s u + S a a , μ s u a a , a a = μ s u a a , a a ( S s u , S a a ) = k m , s u a a S a a k s , s u a a + S a a · S a a S s u + S a a , μ f a = μ f a ( S f a ) = k m , f a S f a k s , f a + S f a .
We extend the model (2)–(11) by the following two equations describing the process in the second bioreactor (cf. [30]):
d d t S a c , c h 4 ( t ) = D 2 ( S a c S a c , c h 4 ) Y a c μ a c , c h 4 X a c
d d t X a c ( t ) = ( D 2 + μ a c , c h 4 ) X a c
with gaseous output
Q c h 4 = Y c h 4 , a c μ a c , c h 4 X a c
where Q c h 4 is the methane flow rate [ L / h ] , and
μ a c , c h 4 = μ a c , c h 4 ( S a c , c h 4 ) = k m , a c S a c , c h 4 k s , a c + S a c , c h 4 .
All model coefficients in the Equations (2)–(13) are assumed to be positive. The dilution rates D 1 and D 2 are the decision or control variables. The numerical values in the rightmost column of Table 2 are taken from [31]. In the model variables and parameters, the subscripts h 2 and c h 4 indicate hydrogen ( H 2 ) and methane ( C H 4 ), respectively.

4. Investigation of the Model Solutions

This section presents a mathematical analysis of the dynamics (2)–(13). First, we establish existence, positivity, uniform boundedness, and uniqueness of the model solutions. The results are presented in Theorem 1. The interested reader can find a detailed proof of Theorem 1 in Appendix A.1 of the Appendix A.
Theorem 1.
If 0 < Y s u a a < 1 and 0 < Y f a < 1 , then all solutions of the dynamics (2)–(13) are positive and uniformly bounded, and thus exist for all time t [ 0 , + ) .
Note that the assumptions 0 < Y s u a a < 1 and 0 < Y f a < 1 in Theorem 1 are not restrictive. Since Y s u a a and Y f a are yield coefficients, their values are always less than 1 in practice; see last column in Table 2.
In the next two subsections, we find the simplest solutions of the model, namely the equilibrium points, and investigate their local asymptotic stability.

4.1. Equilibrium Points of the Model

The equilibrium points of the model are solutions of the algebraic equations obtained from (2)–(13) by setting the right-hand sides equal to zero:
D 1 ( S s u i n S s u ) + k h y d , c h X c h + f s u , l i k h y d , l i X l i μ s u a a , s u X s u a a = 0
D 1 ( S a a i n S a a ) + k h y d , p r X p r μ s u a a , a a X s u a a = 0
D 1 ( S f a i n S f a ) + f f a , l i k h y d , l i X l i μ f a X f a = 0
D 1 S a c ( 1 Y s u a a ) f a c , s u μ s u a a , s u + f a c , a a μ s u a a , a a X s u a a 0.7 ( 1 Y f a ) μ f a X f a = 0
D 1 ( X c i n X c ) k d i s X c = 0
( D 1 + k h y d , c h ) X c h f c h , x c k d i s X c = 0
( D 1 + k h y d , p r ) X p r f p r , x c k d i s X c = 0
( D 1 + k h y d , l i ) X l i f l i , x c k d i s X c = 0
D 1 Y s u a a ( μ s u a a , s u + μ s u a a , a a ) X s u a a = 0
( D 1 Y f a μ f a ) X f a = 0 .
D 2 ( S a c S a c , c h 4 ) Y a c μ a c , c h 4 X a c = 0
( D 2 μ a c , c h 4 ) X a c = 0 .
We are looking for nonnegative solutions of (14)–(25) since only these equilibrium points correspond to the physical and biological meaning of the model variables. The equilibrium points are computed as functions of the control variables, the dilution rates D 1 and D 2 .
Using Equation (18), we find the equilibrium component X ¯ c with respect to the phase variable X c , which is obviously positive:
X ¯ c = D 1 X c i n D 1 + k d i s .
Furthermore, Equations (19)–(21) imply the equilibrium components X ¯ c h , X ¯ p r and X ¯ l i respectively:
X ¯ c h = f c h , x c k d i s X ¯ c D 1 + k h y d , c h , X ¯ p r = f p r , x c k d i s X ¯ c D 1 + k h y d , p r , X ¯ l i = f l i , x c k d i s X ¯ c D 1 + k h y d , l i .
Equations (16) and (23) are used to determine the steady state components S ¯ f a and X ¯ f a with respect to the phase variables S f a and X f a . Obviously, Equation (23) suggests that the trivial solution X f a 0 always exists, but is not of practical interest. Assuming that X f a 0 and using the expression of μ f a , we obtain from (23) the following steady state component S ¯ f a
S ¯ f a = D 1 k s , f a Y f a k m , f a D 1 .
It is clear that S ¯ f a exists and S ¯ f a > 0 if and only if
D 1 < Y f a k m , f a .
Equation (16) then delivers the equilibrium component for X f a ,
X ¯ f a = D 1 ( S f a i n S ¯ f a ) + f f a , l i k h y d , l i X ¯ l i μ f a ( S ¯ f a ) .
Obviously, X ¯ f a from (29) satisfies X ¯ f a > 0 if and only if
S ¯ f a < S f a i n + 1 D 1 f f a , l i k h y d , l i X ¯ l i .
Remark 1.
The quantity S f a i n + 1 D 1 f f a , l i k h y d , l i X ¯ l i on the right-hand side of (30) can be considered as a worst-case upper bound of the total concentration of fatty acids (LCFA) S f a , which accumulation in B R 1 might occur as a result of some imbalances in the degradation phase, leading to wash-out of the LCFA degraders X f a [32].
The inequality (30) is equivalent with
D 1 3 k s , f a + S f a i n + D 1 2 ( k s , f a + S f a i n ) ( k h y d , l i + k d i s ) Y f a k m , f a S f a i n + D 1 f f a , l i k h y d , l i f l i , x c k d i s X c i n Y f a k m , f a S f a i n ( k h y d , l i + k d i s ) + k s , f a + S f a i n Y f a k m , f a k h y d , l i k d i s ( S f a i n + f f a , l i f l i , x c X c i n ) < 0 .
The left hand-side cubic polynomial in (31) always possesses one positive real root D 1 ( 1 ) , so the inequality (30) is satisfied if D 1 ( 0 , D 1 ( 1 ) ) . Therefore, the equilibrium components S ¯ f a and X ¯ f a exist and satisfy S ¯ f a > 0 and X ¯ f a > 0 if D 1 0 , min { Y f a k m , f a , D 1 ( 1 ) } holds true.
If D 1 = D 1 ( 1 ) , then X ¯ f a becomes equal to zero; this means that D 1 = D 1 ( 1 ) is a (transcritical) bifurcation value of the parameter D 1 , leading to wash-out of the biomass X f a for D 1 > D 1 ( 1 ) . We shall pay special attention to the case X ¯ f a = 0 motivating our considerations later on in the paper (see Section 5). If X ¯ f a = 0 , then we obtain the following steady state component for S f a
S ¯ f a 0 = S f a i n + 1 D 1 f f a , l i k h y d , l i X ¯ l i ,
which always exists and is positive if D 1 > D 1 ( 1 ) .
Here, and in what follows, the superscript 0 indicates X ¯ f a = 0 in the steady state components.
Furthermore, Equations (14), (15) and (22) are used to compute the steady state components S ¯ s u , S ¯ a a , and X ¯ s u a a in case positive solutions do exist. It is clear from Equation (22) that the trivial solution X s u a a 0 always exists, but we exclude this case for biological evidence. Let X s u a a 0 . Then, Equation (22) reduces to
D 1 Y s u a a μ s u a a , s u ( S s u , S a a ) + μ s u a a , a a ( S s u , S a a ) = 0 .
Using Equation (15), we find
X s u a a = D 1 ( S a a i n S a a ) + k h y d , p r X ¯ p r μ s u a a , a a ( S s u , S a a ) .
The equilibrium component with respect to X s u a a should be positive and this will be fulfilled if
S a a < S a a i n + 1 D 1 k h y d , p r X ¯ p r .
A similar remark to Remark 1 can be made to the above inequality (35) concerning the concentrations of the amino acids S a a and of the degraders X s u a a .
Substituting X s u a a from (34) into Equation (14) and using Equation (33), we obtain the following nonlinear algebraic system with respect to S a a and S s u :
D 1 ( S s u i n S s u ) + k h y d , c h X ¯ c h + f s u , l i k h y d , l i X ¯ l i μ s u a a , a a ( S s u , S a a ) D 1 ( S a a i n S a a ) + k h y d , p r X ¯ p r μ s u a a , s u ( S s u , S a a ) = 0 D 1 Y s u a a μ s u a a , s u ( S s u , S a a ) + μ s u a a , a a ( S s u , S a a ) = 0 .
It is straightforward to see that system (36) is equivalent to the following one:
A 1 S a a ( S a a + S s u ) A 2 S s u 2 ( 1 + S a a ) A 3 S a a S s u ( S a a + S s u ) = 0 B 1 ( S a a 2 + S s u 2 ) + B 2 S a a S s u ( S a a + S s u ) B 3 S a a S s u B 4 ( S a a + S s u ) = 0 ,
where
A 1 = D 1 S s u i n + D 1 k d i s X c i n D 1 + k d i s k h y d , c h f c h , x c D 1 + k h y d , c h + f s u , l i k h y d , l i f l i , x c D 1 + k h y d , l i , A 2 = D 1 S a a i n + D 1 k h y d , p r f p r , x c k d i s X c i n ( D 1 + k d i s ) ( D 1 + k h y d , p r ) , A 3 = D 1 k s , s u a a ; B 1 = k s , s u a a ( Y s u a a k m , s u a a D 1 ) , B 2 = Y s u a a k m , s u a a D 1 , B 3 = 2 D 1 k s , s u a a , B 4 = D 1 k s , s u a a 2 .
All coefficients A i , i = 1 , 2 , 3 , and B j , j = 1 , 2 , 3 , 4 , depend on D 1 (the dilution rate in the first bioreactor B R 1 ). Assume that there exists a value D 1 ( 2 ) > 0 of D 1 such that system (37) possesses positive solutions S ¯ a a and S ¯ s u for all D 1 ( 0 , D 1 ( 2 ) ) . Denote by D 1 ( 3 ) a third possible upper bound for D 1 , defined by
D 1 ( 3 ) = min { D 1 > 0 : S ¯ a a < S a a i n + 1 D 1 k h y d , p r X ¯ p r } .
Then, the equilibrium component X ¯ s u a a will exist and will be positive for D 1 0 , min { D 1 ( 2 ) , D 1 ( 3 ) } :
X ¯ s u a a = D 1 ( S a a i n S ¯ a a ) + k h y d , p r X ¯ p r μ s u a a , a a ( S ¯ s u , S ¯ a a ) .
Note that an equality S ¯ a a = S a a i n + 1 D 1 k h y d , p r X ¯ p r for some value of D 1 will lead to X ¯ s u a a = 0 and may cause wash-out of the biomass X s u a a , which is practically undesirable, and we shall exclude this case.
Equation (17) is used to find the steady state component
S ¯ a c = 1 D 1 ( 1 Y s u a a ) f a c , s u μ s u a a , s u ( S ¯ s u , S ¯ a a ) + f a c , a a μ s u a a , a a ( S ¯ s u , S ¯ a a ) X ¯ s u a a + 0 . 7 ( 1 Y f a ) μ f a ( S ¯ f a ) X ¯ f a .
The latter is positive if Y s u a a < 1 and Y f a < 1 are fulfilled.
If X ¯ f a = 0 , then the equilibrium component with respect to S a c becomes
S ¯ a c 0 = 1 D 1 ( 1 Y s u a a ) f a c , s u μ s u a a , s u ( S ¯ s u , S ¯ a a ) + f a c , a a μ s u a a , a a ( S ¯ s u , S ¯ a a ) X ¯ s u a a ;
S ¯ a c 0 exists and is positive if Y s u a a < 1 and D 1 D 1 ( 1 ) , min { D 1 ( 2 ) , D 1 ( 3 ) } hold true.
The above calculations deliver the following equilibrium points of the model (2)–(11):
E 1 ( D 1 ) = ( S ¯ s u , S ¯ a a , S ¯ f a , S ¯ a c , X ¯ c , X ¯ c h , X ¯ p r , X ¯ l i , X ¯ s u a a , X ¯ f a ) if   0 < D 1 < min { Y f a k m , f a , D 1 ( 1 ) , D 1 ( 2 ) , D 1 ( 3 ) } ; E 1 0 ( D 1 ) = ( S ¯ s u , S ¯ a a , S ¯ f a 0 , S ¯ a c 0 , X ¯ c , X ¯ c h , X ¯ p r , X ¯ l i , X ¯ s u a a , X ¯ f a = 0 ) if   D 1 ( 1 ) < D 1 < min { D 1 ( 2 ) , D 1 ( 3 ) } .
Consider now Equations (12) and (13) to compute the equilibrium components with respect to X a c and S a c , c h 4 . Thereby, we assume that S ¯ a c and S ¯ a c 0 are known; see (40) and (41), respectively.
Equation (25) has the obvious solution X a c 0 , which is not of practical interest. Assuming that X a c 0 , the equality μ a c , c h 4 ( S a c , c h 4 ) = D 2 delivers the following steady state component for S a c , c h 4
S ¯ a c , c h 4 = k s , a c D 2 k m , a c D 2 .
Obviously, there exists S ¯ a c , c h 4 > 0 if D 2 < k m , a c holds true. Then, the corresponding equilibrium component X ¯ a c is
X ¯ a c = 1 Y a c ( S ¯ a c S ¯ a c , c h 4 ) .
The latter is positive if
S ¯ a c > S ¯ a c , c h 4 .
If we consider S ¯ a c 0 instead of S ¯ a c we obtain
X ¯ a c 0 = 1 Y a c ( S ¯ a c 0 S ¯ a c , c h 4 ) ,
which is positive if
S ¯ a c 0 > S ¯ a c , c h 4 .
The inequalities (44) and (46) give relationships between D 1 and D 2 because S ¯ a c and S ¯ a c 0 depend also on D 1 ; see (40) and (41), respectively.
Therefore, the extended model (2)–(13) can possess the following two equilibrium points, parameterized on D 1 and D 2 :
E ( D 1 , D 2 ) = E 1 ( D 1 ) ( S ¯ a c , c h 4 , X ¯ a c ) , E 0 ( D 1 , D 2 ) = E 1 0 ( D 1 ) ( S ¯ a c , c h 4 , X ¯ a c 0 ) .
We summarize the above calculations in the following theorem.
Theorem 2.
Let 0 < Y s u a a < 1 and 0 < Y f a < 1 be fulfilled in the model Equations (2)–(13). Assume that there exists a value D 1 ( 2 ) > 0 such that positive solutions S ¯ a a and S ¯ s u of system (37) do exist for D 1 ( 0 , D 1 ( 2 ) ) . Then,
(i) 
the equilibrium point E ( D 1 , D 2 ) exists if 0 < D 1 < min { Y f a k m , f a , D 1 ( 1 ) , D 1 ( 2 ) , D 1 ( 3 ) } and 0 < D 2 < min { k m , a c , min { D 2 > 0 : S ¯ a c > S ¯ a c , c h 4 } } ;
(ii) 
the equilibrium point E 0 ( D 1 , D 2 ) exists if D 1 ( 1 ) < D 1 < min { D 1 ( 2 ) , D 1 ( 3 ) } and 0 < D 2 < min { k m , a c , min { D 2 > 0 : S ¯ a c 0 > S ¯ a c , c h 4 } } .
Theorem 2 can be specified using the coefficient values in Table 2. Numerical calculations produce the following values:
Y m k m , f a = 0 . 009 , D 1 ( 1 ) 0 . 00843 , D 1 ( 2 ) 0 . 0942 .
The left plot in Figure 2 visualizes the equilibrium component S ¯ a a , which exists for D 1 ( 0 , D 1 ( 2 ) ) , as well as the upper bound S a a i n + 1 D 1 k h y d , p r X ¯ p r of S ¯ a a according to (35). Obviously, (35) is satisfied for D 1 ( 0 , D 1 ( 2 ) ) so that D 1 ( 2 ) D 1 ( 3 ) . The right plot in Figure 2 shows the equilibrium component X ¯ s u a a ( D 1 ) , which vanishes at D 1 = D 1 ( 2 ) .
Furthermore, we have k m , a c = 0 . 0167 and
S ¯ a c > S ¯ a c , c h 4 D 2 ( 0 , 0 . 0156 ) = : ( 0 , D 2 ( 2 ) ) .
This relation is demonstrated in the left plot of Figure 3. The steady state component S ¯ a c = S ¯ a c ( D 1 ) (solid line) exists for D 1 ( 0 , D 1 ( 1 ) ) . The dash line represents the equilibrium component S ¯ a c , c h 4 = S ¯ a c , c h 4 ( D 2 ) . The horizontal dot line passes through the points ( D 1 ( 1 ) , S ¯ a c ( D 1 ( 1 ) ) ) and ( D 2 ( 2 ) , S ¯ a c , c h 4 ( D 2 ( 2 ) ) . The vertical dash-dot lines pass through D 1 ( 1 ) and D 2 ( 2 ) . This implies the following relation:
E ( D 1 , D 2 ) exists D 1 ( 0 , D 1 ( 1 ) ) = ( 0 , 0 . 00843 ) , D 2 ( 0 , D 2 ( 2 ) ) = ( 0 , 0 . 0156 ) .
Similarly, the right plot in Figure 3 visualizes the relation
S ¯ a c 0 > S ¯ a c , c h 4 D 2 ( 0 , 0 . 0152 ) = : ( 0 , D 2 ( 1 ) )
because the horizontal dot line passes through the points ( D 1 ( 1 ) , S ¯ a c 0 ( D 1 ( 1 ) ) ) and ( D 2 ( 1 ) , S ¯ a c , c h 4 ( D 2 ( 1 ) ) . This leads to the following relation:
E 0 ( D 1 , D 2 ) exists D 1 ( D 1 ( 1 ) , D 1 ( 2 ) ) = ( 0 . 00843 , 0 . 0942 ) , D 2 ( 0 , D 2 ( 1 ) ) = ( 0 , 0 . 0152 ) .
The equilibrium E ( D 1 , D 2 ) corresponds to coexistence of all microorganisms in both bioreactors B R 1 and B R 2 , whereas E 0 ( D 1 , D 2 ) is related to wash-out of the LCFA degrader X f a in B R 1 for hydrogen production.

4.2. Local Asymptotic Stability of the Equilibrium Points

We shall first investigate the local asymptotic stability of the coexistence equilibrium E 1 ( D 1 ) for D 1 0 , D 1 ( 1 ) and of the wash-out steady state (within X ¯ f a = 0 ) E 1 0 ( D 1 ) for D 1 D 1 ( 1 ) , D 1 ( 2 ) , i.e., of the equilibrium points corresponding to the model (2)–(11) of the first bioreactor B R 1 for hydrogen production.
It is well known that an equilibrium point is locally asymptotically stable, if all eigenvalues of the Jacobi matrix evaluated at this equilibrium have negative real parts, cf. e.g., [33]. The eigenvalues of this Jacobi matrix coincide with the roots of the corresponding characteristic polynomial.
Detailed calculations of the characteristic polynomials P ( E 1 ; λ ) and P ( E 1 0 ; λ ) with respect to the equilibrium points E 1 = E 1 ( D 1 ) and E 1 0 = E 1 0 ( D 1 ) are given in Appendix A.2 of the Appendix A. We present them here for reader’s convenience:
P ( E 1 ; λ ) = ( D 1 λ ) ( D 1 k d i s λ ) ( D 1 k h y d , c h λ ) × ( D 1 k h y d , p r λ ) ( D 1 k h y d , l i λ ) × λ 2 + D 1 + d μ f a d S f a ( S ¯ f a ) X ¯ f a λ + D 1 d μ f a d S f a ( S ¯ f a ) X ¯ f a × R ( E 1 ; λ ) , D 1 ( 0 , D 1 ( 1 ) ) ; P ( E 1 0 ; λ ) = ( D 1 λ ) 2 ( D 1 k d i s λ ) ( D 1 k h y d , c h λ ) × ( D 1 k h y d , p r λ ) ( D 1 k h y d , l i λ ) ( D 1 + Y f a μ f a ( S ¯ f a 0 ) λ ) × R ( E 1 0 ; λ ) , D 1 ( D 1 ( 1 ) , D 1 ( 2 ) ) ,
where R ( E 1 ; λ ) is the characteristic polynomial of the sub-matrix Δ ( 3 ) evaluated at the steady state components ( S ¯ s u , S ¯ a a , X ¯ s u a a ) for D 1 ( 0 , D 1 ( 1 ) ) , and R ( E 1 0 ; λ ) is the characteristic polynomial of the sub-matrix Δ ( 3 ) evaluated at the steady state components ( S ¯ s u , S ¯ a a , X ¯ s u a a ) for D 1 ( D 1 ( 1 ) , D 1 ( 2 ) ) , see Appendix A.2 in the Appendix A.
The linear quantities in P ( E 1 ; λ ) = 0 and P ( E 1 0 ; λ ) = 0 produce negative real eigenvalues of the Jacobi matrices. In the first polynomial P ( E 1 ; λ ) , the quadratic equation
λ 2 + D 1 + d μ f a d S f a ( S ¯ f a ) X ¯ f a λ + D 1 d μ f a d S f a ( S ¯ f a ) X ¯ f a = 0
possesses two real negative roots. The multiplier D 1 + Y f a μ f a ( S f a 0 ) λ in the second polynomial P ( E 1 0 ; λ ) delivers a negative root λ = D 1 + Y f a μ f a ( S ¯ f a 0 ) because, in this case, D 1 > Y f a μ f a ( S ¯ f a 0 ) holds true.
It remains to show that R ( E 1 ; λ ) = 0 and R ( E 1 0 ; λ ) = 0 have roots with negative real parts. These two problems are solved numerically using the coefficient values in Table 2. The real parts of the eigenvalues are presented graphically in Figure 4 and Figure 5, and are obviously negative. Therefore, the steady states E 1 = E 1 ( D 1 ) , D 1 ( 0 , D 1 ( 1 ) ) and E 1 0 = E 1 0 ( D 1 ) , D 1 ( D 1 ( 1 ) , D 1 ( 2 ) ) , are locally asymptotically stable.
Consider again the polynomial P ( E 1 ; λ ) . We know that at D 1 = D 1 ( 1 ) the equilibrium component X ¯ f a becomes equal to zero, and thus P ( E 1 ; λ ) = 0 has a root λ = 0 . This means that a transcritical bifurcation at D 1 = D 1 ( 1 ) occurs, leading to wash-out of the LCFA degraders X f a . As a result, the steady state E 1 ( D 1 ) no longer exists for D 1 > D 1 ( 1 ) , and E 1 0 ( D 1 ) remains the locally asymptotically stable equilibrium for D 1 ( D 1 ( 1 ) , D 1 ( 2 ) ) .
The equilibrium points of the model (12)–(13) of the second bioreactor B R 2 for methane production are known to be locally asymptotically stable for any admissible value of D 2 (see, e.g., [20]). Indeed, the characteristic polynomials of the Jacobi matrices of (12)–(13) evaluated at the equilibrium points ( S ¯ a c , c h 4 , X ¯ a c ) and ( S ¯ a c , c h 4 , X ¯ a c 0 ) , have the form
λ 2 + A λ + B = 0 a n d λ 2 + A 0 λ + B 0 = 0 ,
where
A = D 2 + Y a c d μ a c , c h 4 d S a c , c h 4 ( S ¯ a c , c h 4 ) X ¯ a c > 0 , B = Y a c μ a c , c h 4 ( S ¯ a c , c h 4 ) d μ a c , c h 4 d S a c , c h 4 ( S ¯ a c , c h 4 ) X ¯ a c > 0 , A 0 = D 2 + Y a c d μ a c , c h 4 d S a c , c h 4 ( S ¯ a c , c h 4 ) X ¯ a c 0 > 0 , B 0 = Y a c μ a c , c h 4 ( S ¯ a c , c h 4 ) d μ a c , c h 4 d S a c , c h 4 ( S ¯ a c , c h 4 ) X ¯ a c 0 > 0 .
The positive signs of A, B and A 0 , B 0 imply that the above quadratics possess two roots with negative real parts. This proves the local asymptotic stability of the equilibrium points E ( D 1 , D 2 ) for D 1 ( 0 , D 1 ( 1 ) ) , D 2 ( 0 , D 2 ( 2 ) ) and E 0 ( D 1 , D 2 ) for D 1 ( D 1 ( 1 ) , D 1 ( 2 ) ) , D 2 ( 0 , D 2 ( 1 ) ) .

5. Optimization of Hydrogen and Methane Flow Rates

In this section, we shall find values of the dilution rates D 1 and D 2 , so that the hydrogen and methane flow rates Q h 2 and Q c h 4 evaluated at the model equilibrium points achieve maximal values.
We compute Q h 2 on the two sets of steady states E 1 ( D 1 ) , D 1 ( 0 , D 1 ( 1 ) ) and E 1 0 ( D 1 ) , D 1 ( D 1 ( 1 ) , D 1 ( 2 ) ) :
Q h 2 = Q h 2 ( D 1 ) = Y h 2 , s u μ s u a a , s u ( S ¯ s u , S ¯ a a ) X ¯ s u a a + Y h 2 , a a μ s u a a , a a ( S ¯ s u , S ¯ a a ) X ¯ s u a a + Y h 2 , f a μ f a ( S ¯ f a ) X ¯ f a , D 1 0 , D 1 ( 1 ) ; Q h 2 0 = Q h 2 0 ( D 1 ) = Y h 2 , s u μ s u a a , s u ( S ¯ s u , S ¯ a a ) X ¯ s u a a + Y h 2 , a a μ s u a a , a a ( S ¯ s u , S ¯ a a ) X ¯ s u a a , D 1 D 1 ( 1 ) , D 1 ( 2 )
Recall that the superscript 0 in the equilibrium E 1 0 ( D 1 ) indicates X ¯ f a = 0 .
The functions Q h 2 ( D 1 ) and Q h 2 0 ( D 1 ) depend on the decision (manipulated) input D 1 and are called input–output static characteristics with respect to the hydrogen production.
The important (from practical point of view) question is whether the functions Q h 2 ( D 1 ) and Q h 2 0 ( D 1 ) are unimodal with respect to D 1 in the admissible intervals for D 1 . For example, the function Q h 2 ( D 1 ) is called unimodal if there exists a unique (admissible) point D 1 , max , such that Q h 2 ( D 1 ) possesses maximum Q h 2 , max = Q h 2 ( D 1 , max ) , Q h 2 ( D 1 ) is strongly increasing if D 1 < D 1 , max and Q h 2 ( D 1 ) is strongly decreasing if D 1 > D 1 , max .
Figure 6 (left plot) presents the graphs of the input–output static characteristics Q h 2 ( D 1 ) for D 1 ( 0 , D 1 ( 1 ) ) (solid line) and Q h 2 0 ( D 1 ) for D 1 ( D 1 ( 1 ) , D 1 ( 2 ) ) (dash line). The function Q h 2 ( D 1 ) achieves its maximum at D 1 , max 0.0072 ( 0 , D 1 ( 1 ) ) and Q h 2 , max 0.115 [ L / h ] . The maximum of Q h 2 0 ( D 1 ) is achieved at D 1 , max 0 0.0439 ( D 1 ( 1 ) , D 1 ( 2 ) ) and Q h 2 , max 0 0.169 [ L / h ] . Therefore, the absolute maximum of Q h 2 ( D 1 ) Q h 2 0 ( D 1 ) , D 1 ( 0 , D 1 ( 2 ) ) coincides with the maximum of Q h 2 0 ( D 1 ) . From a practical point of view, this means that the maximal flow rate of H 2 could be achieved at a operation mode, related to wash-out of the LCFA degraders X f a in B R 1 .
A similar question of maximizing the biogas production at steady states has been discussed in [34] for the well known four-dimensional AM2 model. There, two maxima of the input–output static characteristics with respect to methane flow rate have also been determined, related to coexistence or wash-out of the acidogenic bacteria in the bioreactor.
Using the methane flow rate Q c h 4 , we compute further the input–output static characteristics Q c h 4 and Q c h 4 0 on the set of the steady states E ( D 1 , D 2 ) and E 0 ( D 1 , D 2 ) respectively, namely:
Q c h 4 = Q c h 4 ( D 1 , D 2 ) = Y c h 4 , a c μ a c , c h 4 ( S ¯ a c , c h 4 ) X ¯ a c , Q c h 4 0 = Q c h 4 ( D 1 , D 2 ) = Y c h 4 , a c μ a c , c h 4 ( S ¯ a c , c h 4 ) X ¯ a c 0 .
The latter are parameterized on the two inputs D 1 and D 2 , since X ¯ a c and X ¯ a c 0 depend on both D 1 and D 2 ; see (43) and (45), respectively.
With D 1 = D 1 , max 0 0.0439 , we compute the equilibrium components S ¯ a c , c h 4 and X ¯ a c 0 according to (42) and (45), as well as Q c h 4 0 , all of them as functions of D 2 .
Figure 6 (right plot) visualizes the graph of Q c h 4 0 ( D 1 , max 0 , D 2 ) for D 2 ( 0 , D 2 ( 1 ) ) ; the latter is a unimodal function, taking its maximum at D 2 , max 0 0.00991 [ 1 / h ] with Q c h 4 , max 0 = Q c h 4 ( D 1 , max 0 , D 2 , max 0 ) 0.0388 [ L / h ] .
The latter results are theoretical and predicted by the model studies. Experimental data that confirm these results by considering H 2 and C H 4 volumetric production under similar conditions for which our model was developed (i.e., two bioreactors in a cascade with a similar ratio of working volumes and operating in continuous mode) can be found in [35].
Using the presentation (1), we define and compute the constant γ and in this way the relationship between the volumes V 1 and V 2 of the two bioreactors:
γ = D 2 , max 0 D 1 , max 0 0.226 V 2 4.42 V 1 .
The last relation implies that the volume of the second bioreactor B R 2 for methane production should be at least 4.4 times greater than the volume of the first bioreactor B R 1 for hydrogen production.

6. Dynamic Behavior of the Model Solutions

The dynamic model has been implemented as one-page application “(H2,CH4) Bioreactors ver. 1.2” in a web-based simulation software SmoWeb, see http://platform.sysmoltd.com/. SmoWeb is an open computational platform in Python and a library of applications is built on top of this platform. It does not require installation of any software by the user, since all computations are performed in a web cloud. Python is used as the implementation language because it is a powerful modern general purpose object-oriented language. Nowadays, Python has become an inseparable part of the scientific computing due to the vast abundance of open source libraries and interfaces to major software tools. A number of one-page applications have already been and are continuously being developed on top of the functionality provided by the platform, including bioprocess Modeling as well.
All figures below are developed using the one-page application “(H2,CH4) Bioreactors ver. 1.2”. The plots illustrate the dynamic behavior of the model solutions when D 1 , D 2 and the input concentration of the composites X c i n take different values from Table 3 and Table 4, respectively.
Figure 7 presents the time evolution of the flow rates Q h 2 ( t ) and Q c h 4 ( t ) when D 1 and D 2 take values from Table 3. Figure 8 and Figure 9 display the model solutions by variable D 1 and D 2 given in Table 3. It is seen that, if D 1 ( 0 , D 1 ( 1 ) ) = ( 0 , 0 . 00843 ) and D 2 ( 0 , D 2 ( 2 ) ) = ( 0 , 0 . 0156 ) , then the model solutions tend to the corresponding components of the locally stable coexistence equilibrium E ( D 1 , D 2 ) . When D 1 ( D 1 ( 1 ) , D 1 ( 2 ) ) = ( 0 . 00843 , 0 . 0942 ) and D 2 ( 0 , D 2 ( 1 ) ) = ( 0 , 0 . 0152 ) , then the solutions approach the locally stable wash-out equilibrium (with X ¯ f a = 0 ) after sufficiently large time t [ h ] ; in particular, one can see in the left plot of Figure 9 that X f a ( t ) approaches 0 for sufficiently large time t [ h ] .
Figure 10, Figure 11 and Figure 12 demonstrate the dynamics of Q h 2 ( t ) , Q c h 4 ( t ) as well as of the model solutions, when D 1 and D 2 are fixed, D 1 = 0.005 [ 1 / h ] , D 2 = 0.00113 [ 1 / h ] , but the input concentration of composites X c i n varies taking values from Table 4. The model solutions tend to the corresponding components of the coexistence equilibrium E ( D 1 , D 2 ) , which is the locally stable steady state. These numerical simulations demonstrate the fact that uncertainties in the input X c i n do not affect the stability of the model dynamics.

7. Discussion

In this paper, we propose a mathematical model describing the process of simultaneous production of hydrogen and methane by anaerobic digestion of organic wastes in a cascade of two connected continuously stirred bioreactors with different volumes. The proposed model is developed by adapting and reducing the universal Anaerobic Digester Model No 1 (ADM1). The equations in our model are separated in two groups, corresponding to the processes in the two bioreactors, for hydrogen and for methane production, respectively. The simplifications made in the two models naturally lead to some limitations compared to the original ADM1. However, as already shown, the latter also does not cover all possible reactions, e.g., possible inclusion in the model of the activity of the so called syntrophicacetate oxidizers [36]. Some new, more complex models have also been developed [37], which however have the same main disadvantage as ADM1—great complexity and thus impossibility to carry out some mathematical analytical studies.
The investigation of the hydrogen and methane flow rates, Q h 2 and Q c h 4 respectively, on steady state operation, shows existence of maxima for some values of the dilution rates D 1 and D 2 in the two bioreactors B R 1 and B R 2 . The local maximum for D 1 in B R 1 reflects the maximal yield of hydrogen from lipids at very long hydraulic retention time HRT = 1/ D 1 , max 1 / 0.0072 139 [hours] 5.8 [days]. This maximum is much less than the maximum yield of hydrogen derived from monosaccharides and aminoacids (proteins) for the accepted combination of coefficients and initial conditions reflecting the biodegradation of organic wastes. In [35] under similar conditions, results close to our theoretical model-based predictions are obtained.
The model also allows us to find the optimal ratio between the volumes V 1 and V 2 of the two bioreactors subject to the same optimization goal. Many studies in the literature report on TSAD processes; however, most of them are operated manually and the ratio of working volumes of CSTRs is not discussed. In some references, a ratio of 10 is accepted without comments. In [38], a TSAD system consisting of two CSTRs operating at mesophilic conditions are used to investigate the effect of HRT on hydrogen and methane production; the ratio of working volumes of the bioreactors is equal to 8 (without comments). In [39], optimization of TSAD of separately collected municipality bio-waste is carried out in pilot scale CSTRs at thermophilic temperature, using recirculation of the digestate in the second bioreactor to maintain p H in the first bioreactor at optimal value. There, the ratio of working volumes of the bioreactors is equal to 3.8 (without comments). In [40], optimal loading rate is obtained providing maximum H 2 and C H 4 productions in TSAD of cassava wastewaters using specific thermophilic bioreactors and constant recycling ratio of 1:1 with automatic control of p H in the hydrogenic bioreactor. The ratio of working volumes of both bioreactors is equal to 6 (without explanations). The review article [41] is the first one that combines optimization approaches for three possible products from AD—methane, hydrogen and volatile fatty acids (VFA), taking into account different process parameters and types of bioreactors, including acidogenesis and methanogenesis separation in two different bioreactors. However, the ratio of the bioreactors working volumes is not given and discussed.
Using the adapted values of the coefficients in our model, we obtain an optimal ratio of the working volumes of both bioreactors equal to 4.42 using a criterion for maximizing the bioenergy production. This value falls into the range of the above reported values; however, it is theoretically derived.

8. Conclusions

In this paper, we consider a mathematical model describing a biotechnological process of hydrogen and methane production by anaerobic digestion. The process is carried out in two connected continuously stirred bioreactors, thus the model equations are separated in two groups, describing the processes in each bioreactor. We establish existence and uniqueness of positive and uniformly bounded solutions of the model dynamics. Two equilibrium points of the model are computed and their local asymptotic stability is studied with respect to the practically important input parameters, the dilution rates D 1 and D 2 in the two bioreactors. The first equilibrium point corresponds to coexistence of substrates and biomass, and the second one is related to wash-out of the LCFA degraders in the first bioreactor. We show that there exists a value of the dilution rate D 1 , at which the coexistence equilibrium disappears as a result of bifurcation of the steady state and the wash-out equilibrium remains the only locally asymptotically stable equilibrium. The investigation of the hydrogen and methane flow rates, Q h 2 and Q c h 4 , on steady state operation, shows existence of maxima at some values of the dilution rates D 1 and D 2 . It is worth noting that the maximal hydrogen flow rate is achieved at the steady state where the LCFA degrader X f a is washed out. This is established here by taking numerical values for the model parameters from the literature [31]. Nevertheless, this fact could help in the design of an AD process to achieve maximal production of either hydrogen and methane. The future work envisages extending the model of the second bioreactor by equations describing an acetogenesis phase, where other VFA (propionate, butyrate and valeriate) will be transformed into acetate, as well as including inhibitory kinetics of the microorganisms growth. The elaborated web-based software will allow us to check different model-based optimization and control strategies that will contribute to designing and engineering a real process.

Author Contributions

M.B.: software development, simulation and visualization, writing and editing; N.D.: theoretical investigation of the model—existence and uniqueness of solutions, stability analysis of equilibrium points, writing and editing; I.S.: model development and validation, discussion and evaluation of the results, editing. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

This research has been partially supported by the Bulgarian Science Fund under contract No. DFNI–E02/13. The work of the second author has been partially supported by Grant No. BG05M2OP001-1.001-0003, financed by the Science and Education for Smart Growth Operational Program (2014–2020) in Bulgaria and co-financed by the European Union through the European Structural and Investment Funds. The authors are grateful to the anonymous reviewers for the useful comments that improved the content of this article.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Appendix A.1. Proof of Theorem 1

Equation (6) is linear with respect to X c and can be solved explicitly:
X c ( t ) = X c ( 0 ) e ( D 1 + k d i s ) t + D 1 X c i n D 1 + k d i s .
The solution X c ( t ) obviously exists for all t 0 , is positive, and is uniformly bounded. In practice, it is reasonable to consider only positive initial conditions X c ( 0 ) > 0 .
Assume that X c h ( 0 ) = 0 . Then, Equation (7) implies X ˙ c h ( 0 ) = f c h , x c k d i s X c ( 0 ) > 0 and so there exists ε > 0 such that X c h ( t ) > 0 for each t ( 0 , ε ) . Let X c h ( 0 ) 0 and suppose that there exists t ¯ c h > 0 such that X c h ( t ) > 0 for each t ( 0 , t ¯ c h ) and X c h ( t ¯ c h ) = 0 . Then, X ˙ c h ( t ¯ c h ) 0 , but it follows from (7) that X ˙ c h ( t ¯ c h ) > 0 , a contradiction. Thus, X c h ( t ) > 0 for all t 0 .
Similar arguments as above can be applied to X p r and X l i from Equations (8) and (9), respectively, so that X p r ( t ) > 0 and X l i ( t ) > 0 for all t 0 .
Equation (10) implies that, if X s u a a ( 0 ) = 0 , then X s u a a ( t ) = 0 for all t [ 0 , ) due to the uniqueness of solutions of Cauchy problems. Assuming that X s u a a ( 0 ) > 0 , we obtain
X s u a a ( t ) = X s u a a ( 0 ) e 0 t D 1 + Y s u a a ( μ s u a a , s u ( S s u ( ξ ) , S a a ( ξ ) ) + μ s u a a , a a ( S s u ( ξ ) , S a a ( ξ ) ) ) d ξ > 0 f o r a l l t 0 .
Similarly, if X f a ( 0 ) = 0 , then it follows from Equation (11) that X f a ( t ) = 0 for all t [ 0 , ) . Thus, it is reasonable to consider only positive initial conditions X f a ( 0 ) > 0 . Then, we have
X f a ( t ) = X f a ( 0 ) e 0 t D 1 Y f a μ f a ( S f a ( ξ ) ) d ξ > 0 f o r a l l t 0 .
Assume now that X f a ( 0 ) > 0 , X s u a a ( 0 ) > 0 and X j ( 0 ) 0 for j { c , c h , p r , l i } .
If S s u ( 0 ) = 0 , then it follows from Equation (2) that
S ˙ s u ( 0 ) = D 1 S s u i n + k h y d , c h X c h ( 0 ) + f s u , l i k h y d , l i X l i ( 0 ) > 0 ,
thus there exists ε > 0 with S s u ( t ) > 0 for t ( 0 , ε ) . Let S s u ( 0 ) 0 . Assume that there exists t ¯ s u > 0 such that S s u ( t ) > 0 for t ( 0 , t ¯ s u ) and S s u ( t ¯ s u ) = 0 . Then, S ˙ s u ( t ¯ s u ) 0 , but, from (2), it follows that S ˙ s u ( t ¯ s u ) > 0 , a contradiction. Therefore, S s u ( t ) > 0 for all t > 0 . Using Equations (3) and (4), one can show in a similar way that S a a ( t ) > 0 and S f a ( t ) > 0 for all t 0 .
Denote
Σ 1 = X c h + f c h , x c X c , Σ 1 i n = f c h , x c X c i n .
Since X c h ( t ) > 0 for all t 0 , it follows from Equations (6) and (7) that
Σ ˙ 1 = X ˙ c h + f c h , x c X ˙ c = D 1 f c h , x c X c i n ( X c h + f c h , x c X c ) k h y d , c h X c h D 1 Σ 1 i n Σ 1 .
Multiplying both sides of the inequality by e D 1 t > 0 , we obtain e D 1 t Σ ˙ 1 + D 1 e D 1 t Σ 1 e D 1 t D 1 Σ 1 i n , which is equivalent to d d t e D 1 t Σ 1 e D 1 t D 1 Σ 1 i n , and further
0 t d d ξ e D 1 ξ Σ 1 ( ξ ) d ξ D 1 Σ 1 i n 0 t e D 1 ξ d ξ ,
thus finally
Σ 1 e D 1 t Σ 1 ( 0 ) + Σ 1 i n 1 e D 1 t .
The latter inequality means that lim sup t Σ 1 ( t ) Σ 1 i n . Since X c ( t ) is positive and bounded, and X c h ( t ) is positive, this implies that X c h ( t ) is uniformly bounded and exists for all t [ 0 , ) .
Denote Σ 2 = X p r + f p r , x c X c , Σ 2 i n = f p r , x c X c i n . Then, using the fact that X c h ( t ) 0 for all t 0 , we obtain from Equations (6) and (9)
Σ ˙ 2 = D 1 f p r , x c X c i n X p r f p r , x c X c k h y d , p r X p r D 1 Σ 2 i n Σ 2 ,
which means that lim sup t Σ 2 ( t ) Σ 2 i n . Since X c ( t ) is positive and bounded, the latter inequality implies that X p r ( t ) is also uniformly bounded, and exists for all t [ 0 , ) .
Define Σ 3 = X l i + f l i , x c X c and Σ 3 i n = f l i , x c X c i n . Then, using Equations (6) and (9) and the fact that X l i ( t ) 0 for all t 0 , we obtain
Σ ˙ 3 = D 1 f l i , x c X c i n X l i f l i , x c X c k h y d , l i X l i D 1 Σ 3 i n Σ 3 ,
which implies lim sup t Σ 3 ( t ) Σ 3 i n . Since X c ( t ) is positive and bounded, the latter inequality means that X l i ( t ) is also uniformly bounded, and exists for all t [ 0 , ) .
Consider now
Σ 4 = S f a + 1 Y f a X f a + f f a , l i X l i + f f a , l i f l i , x c X c , Σ 4 i n = S f a i n + f f a , l i f l i , x c X c i n .
Then, we obtain from Equations (6), (9), and (11)
Σ ˙ 4 = D 1 Σ 4 i n Σ 4 , thus Σ 4 ( t ) = Σ 4 ( 0 ) e D 1 t + Σ 4 i n 1 e D 1 t > 0 .
Since X c ( t ) and X l i ( t ) are positive and bounded, S f a ( t ) and X f a ( t ) are positive, the latter inequality implies that S f a ( t ) and X f a ( t ) are uniformly bounded and exists for all t [ 0 , ) .
Further consider
Σ 5 = S s u + S a a + 1 Y s u a a X s u a a + X c h + X p r + f s u , l i X l i + f c h , x c k d i s + f l i , x c f s u , l i + f p r , x c X c , Σ 5 i n = S s u i n + S a a , i n + f c h , x c k d i s + f l i , x c f s u , l i + f p r , x c X c i n .
Then, we obtain
Σ ˙ 5 = D 1 Σ 5 i n Σ 5 Σ 5 ( t ) = Σ 5 ( 0 ) e D 1 t + Σ 5 i n 1 e D 1 t .
Since X c ( t ) , X c h ( t ) , X l i ( t ) and X p r ( t ) are positive and bounded, the above presentation implies that S s u ( t ) , S a a ( t ) and X s u a a ( t ) are also uniformly bounded, and exist for all t [ 0 , ) .
Consider now Equation (5). Since S a c ( t ) serves as input into the second bioreactor, see Equation (12), it is reasonable to consider only positive initial values S a c ( 0 ) > 0 . The positivity of the phase variables included on the right-hand side of Equation (5) imply that S ˙ a c ( t ) D 1 S a c ( t ) and thus S a c ( t ) S a c ( 0 ) e D 1 t > 0 for all t [ 0 , ) .
Furthermore, the boundedness of the specific growth rate functions imply that there exists a constant Γ 1 > 0 such that S ˙ a c ( t ) D 1 S a c ( t ) + Γ 1 for each t 0 . Then, we have consecutively
d d t e D 1 t S a c ( t ) Γ 1 e D 1 t e D 1 t S a c ( t ) S a c ( 0 ) + Γ 1 0 t e D 1 ξ d ξ S a c ( t ) e D 1 t S a c ( 0 ) + Γ 1 D 1 1 e D 1 t ,
which means that lim sup t S a c ( t ) Γ 1 D 1 .
Therefore, all solutions of the model (2)–(11) are positive, uniformly bounded, and thus exist for all t 0 .
Consider now the models (12) and (13). Obviously, if X a c ( 0 ) = 0 , then X a c ( t ) = 0 for all t 0 . Therefore, it is reasonable to consider only positive initial conditions X a c ( 0 ) > 0 . Then,
X a c ( t ) = X a c ( 0 ) e 0 t D 2 + μ a c , c h 4 ( S a c , c h 4 ( ξ ) ) d ξ > 0 f o r a l l t 0 .
If S a c , c h 4 ( 0 ) = 0 , then S ˙ a c , c h 4 ( 0 ) = D 2 S a c ( 0 ) > 0 . Assume that there exists t ^ > 0 such that S a c , c h 4 ( t ) 0 for t ( 0 , t ^ ) but S a c , c h 4 ( t ^ ) = 0 . Then, there exists ε > 0 such that S ˙ a c , c h 4 ( t ) > 0 for each t ( t ^ ε , t ^ ) , and so
S a c , c h 4 ( t ^ ) = S a c , c h 4 ( t ) + t t ^ S ˙ a c , c h 4 ( ξ ) d ξ > 0 ,
which is a contradiction. Thus, S a c , c h 4 ( t ) > 0 for all t 0 .
Finally, denote Σ 6 = S a c , c h 4 + Y a c X a c . Since S a c ( t ) is positive and uniformly upper bounded, there exists a constant Γ 2 > 0 such that S a c ( t ) Γ 2 for all t 0 . Then, it is straightforward to see that
Σ ˙ 6 = D 2 ( S a c Σ 6 ) D 2 ( Γ 2 Σ 6 ) ,
which means that Σ 6 ( t ) Σ 6 ( 0 ) e D 2 t + Γ 2 ( 1 e D 2 t ) . Thus, lim sup t Σ 6 ( t ) Γ 2 . Since S a c , c h 4 and X a c are positive, it follows that they are bounded as well. This proves Theorem 1.

Appendix A.2. Characteristic Polynomials of the Equilibrium Points E1(D1) and E 1 0 (D1)

We introduce the following notations for simplicity and better readability:
g 1 = D 1 ( S s u i n S s u ) + k h y d , c h X c h + f s u , l i k h y d , l i X l i μ s u a a , s u X s u a a g 2 = D 1 ( S a a i n S a a ) + k h y d , p r X p r μ s u a a , a a X s u a a g 3 = D 1 ( S f a i n S f a ) + f f a , l i k h y d , l i X l i μ f a X f a g 4 = D 1 S a c + ( 1 Y s u a a ) f a c , s u μ s u a a , s u + f a c , a a μ s u a a , a a X s u a a + 0.7 ( 1 Y f a ) μ f a X f a g 5 = D 1 ( X c i n X c ) k d i s X c g 6 = ( D 1 + k h y d , c h ) X c h + f c h , x c k d i s X c g 7 = ( D 1 + k h y d , p r ) X p r + f p r , x c k d i s X c g 8 = ( D 1 + k h y d , l i ) X l i + f l i , x c k d i s X c g 9 = D 1 + Y s u a a ( μ s u a a , s u + μ s u a a , a a ) X s u a a g 10 = ( D 1 Y f a μ f a ) X f a .
The Jacobi matrix J related to g j , j = 1 , 2 , , 10 , with respect to ( S s u , S a a , S f a , S a c , X c , X c h , X p r , X l i , X s u a a , X f a ) has the form
J = g 1 S s u g 1 S a a g 1 S f a g 1 S a c g 1 X c g 1 X c h g 1 X p r g 1 X l i g 1 X s u a a g 1 X f a g 2 S s u g 2 S a a g 2 S f a g 2 S a c g 2 X c g 2 X c h g 2 X p r g 2 X l i g 2 X s u a a g 2 X f a g 10 S s u g 10 S a a g 10 S f a g 10 S a c g 10 X c g 10 X c h g 10 X p r g 10 X l i g 10 X s u a a g 10 X f a
Using the explicit dependance of g j , j = 1 , 2 , , 10 , on the above variables, we obtain
J = g 1 S s u g 1 S a a 0 0 0 g 1 X c h 0 g 1 X l i g 1 X s u a a 0 g 2 S s u g 2 S a a 0 0 0 0 g 2 X p r 0 g 2 X s u a a 0 0 0 g 3 S f a 0 0 0 0 g 3 X l i 0 g 3 X f a g 4 S s u g 4 S a a g 4 S f a g 4 S a c 0 0 0 0 g 4 X s u a a g 4 X f a 0 0 0 0 g 5 X c 0 0 0 0 0 0 0 0 0 g 6 X c g 6 X c h 0 0 0 0 0 0 0 0 g 7 X c 0 g 7 X p r 0 0 0 0 0 0 0 g 8 X c 0 0 g 8 X l i 0 0 g 9 S s u g 9 S a a 0 0 0 0 0 0 g 9 X s u a a 0 0 0 g 10 S f a 0 0 0 0 0 0 g 10 X f a
The characteristic polynomial corresponding to the Jacobi matrix J is defined by d e t ( J λ I 10 ) , where λ is any complex number and I 10 is the ( 10 × 10 ) –identity matrix.
We note that the Jacobian J = ( J i j ) , i , j = 1 , 2 , , 10 , is a sparse matrix, whose fourth column and fifth row contain only one nonzero element, i.e., J 44 = g 4 S a c = D 1 and J 55 = g 5 X c = D 1 k d i s , respectively. Then, we obtain
d e t ( J λ I 10 ) = ( D 1 λ ) ( D 1 k d i s λ ) × g 1 S s u λ g 1 S a a 0 k h y d , c h 0 f s u , l i k h y d , l i μ s u a a , s u 0 g 2 S s u g 2 S a a λ 0 0 k h y d , p r 0 μ s u a a , a a 0 0 0 g 3 S f a λ 0 0 f f a , l i k h y d , l i 0 μ f a 0 0 0 D 1 k h y d , c h λ 0 0 0 0 0 0 0 0 D 1 k h y d , p r λ 0 0 0 0 0 0 0 0 D 1 k h y d , l i λ 0 0 g 9 S s u g 9 S a a 0 0 0 0 g 9 X s u a a λ 0 0 0 g 10 S f a 0 0 0 0 g 10 X f a λ = ( D 1 λ ) ( D 1 k d i s λ ) ( D 1 k h y d , c h λ ) ( D 1 k h y d , p r λ ) ( D 1 k h y d , l i λ ) × g 1 S s u λ g 1 S a a 0 μ s u a a , s u 0 g 2 S s u g 2 S a a λ 0 μ s u a a , a a 0 0 0 g 3 S f a λ 0 μ f a g 9 S s u g 9 S a a 0 g 9 X s u a a λ 0 0 0 g 10 S f a 0 g 10 X f a λ = ( D 1 λ ) ( D 1 k d i s λ ) ( D 1 k h y d , c h λ ) ( D 1 k h y d , p r λ ) ( D 1 k h y d , l i λ ) × d e t ( Δ ( 5 ) λ I 5 ) ,
where I 5 is the ( 5 × 5 ) -identity matrix and Δ ( 5 ) denotes the matrix
Δ ( 5 ) = g 1 S s u g 1 S a a 0 μ s u a a , s u 0 g 2 S s u g 2 S a a 0 μ s u a a , a a 0 0 0 g 3 S f a 0 μ f a g 9 S s u g 9 S a a 0 g 9 X s u a a 0 0 0 g 10 S f a 0 g 10 X f a
Furthermore, we obtain
d e t ( Δ ( 5 ) λ I 5 ) = μ f a g 1 S s u λ g 1 S a a 0 μ s u a a , s u g 2 S s u g 2 S a a λ 0 μ s u a a , a a g 9 S s u g 9 S a a 0 g 9 X s u a a λ 0 0 0 g 10 S f a 0 + g 10 X f a λ g 1 S s u λ g 1 S a a 0 μ s u a a , s u g 2 S s u g 2 S a a λ 0 μ s u a a , a a 0 0 g 3 S f a λ 0 g 9 S s u g 9 S a a 0 g 9 X s u a a λ = μ f a g 10 S f a + g 10 X f a λ g 3 S f a λ × d e t ( Δ ( 3 ) λ I 3 ) ,
where I 3 is the ( 3 × 3 ) -identity matrix and Δ ( 3 ) denotes the matrix
Δ ( 3 ) = g 1 S s u g 1 S a a μ s u a a , s u g 2 S s u g 2 S a a μ s u a a , a a g 9 S s u g 9 S a a g 9 X s u a a .
Replacing the derivatives g 10 S f a , g 10 X f a and g 3 S f a by their corresponding explicit expressions, the characteristic polynomial d e t ( J λ I 10 ) is presented by
d e t ( J λ I 10 ) = ( D 1 λ ) ( D 1 k d i s λ ) ( D 1 k h y d , c h λ ) ( D 1 k h y d , p r λ )     × ( D 1 k h y d , l i λ ) μ f a Y f a X f a d μ f a d S f a + D 1 + Y f a μ f a λ D 1 X f a d μ f a d S f a λ     × d e t ( Δ ( 3 ) λ I 3 ) .
The corresponding partial derivatives in Δ ( 3 ) are as follows:
g 1 S s u = D 1 μ s u a a , s u S s u X ¯ s u a a g 1 S a a = μ s u a a , s u S a a X ¯ s u a a g 2 S s u = μ s u a a , a a S s u X ¯ s u a a g 2 S a a = D 1 μ s u a a , a a S a a X ¯ s u a a g 9 S s u = Y s u a a X ¯ s u a a μ s u a a , s u S s u + μ s u a a , a a S s u g 9 S a a = Y s u a a X ¯ s u a a μ s u a a , s u S a a + μ s u a a , a a S a a g 9 X s u a a = ( D 1 Y f a μ ( S ¯ f a ) .
(i). The characteristic polynomial of the coexistence equilibrium
E 1 = E 1 ( D 1 ) = ( S ¯ s u , S ¯ a a , S ¯ f a , S ¯ a c , X ¯ c , X ¯ c h , X ¯ p r , X ¯ l i , X ¯ s u a a , X ¯ f a ) > 0 f o r D 1 0 , D 1 ( 1 ) .
Denote by P ( E 1 ; λ ) the characteristic polynomial d e t ( J λ I 10 ) evaluated at the equilibrium E 1 , and by R ( E 1 ; λ ) the characteristic polynomial d e t ( Δ ( 3 ) λ I 3 ) evaluated at the equilibrium components ( S ¯ s u , S ¯ a a , X ¯ s u a a ) for D 1 0 , D 1 ( 1 ) . Taking into account that D 1 + Y f a μ f a ( S ¯ f a ) = 0 , we obtain for P ( E 1 ; λ ) the presentation
P ( E 1 ; λ ) = ( D 1 λ ) ( D 1 k d i s λ ) ( D 1 k h y d , c h λ ) × ( D 1 k h y d , p r λ ) ( D 1 k h y d , l i λ ) × λ 2 + D 1 + d μ f a d S f a ( S ¯ f a ) X ¯ f a λ + D 1 d μ f a d S f a ( S ¯ f a ) X ¯ f a × R ( E 1 ; λ ) .
(ii). The characteristic polynomial of the equilibrium
E 1 0 = E 1 0 ( D 1 ) = ( S ¯ s u , S ¯ a a , S ¯ f a 0 , S ¯ a c 0 , X ¯ c , X ¯ c h , X ¯ p r , X ¯ l i , X ¯ s u a a , X ¯ f a = 0 ) f o r D 1 D 1 ( 1 ) , D 1 ( 2 ) .
Denote by P ( E 1 0 ; λ ) the characteristic polynomial d e t ( J λ I 10 ) evaluated at the equilibrium E 1 0 , and by R ( E 1 0 ; λ ) the characteristic polynomial d e t ( Δ ( 3 ) λ I 3 ) evaluated at the equilibrium components ( S ¯ s u , S ¯ a a , X ¯ s u a a ) for D 1 D 1 ( 1 ) , D 1 ( 2 ) . Taking into account that X ¯ f a = 0 and that D 1 + Y f a μ f a ( S ¯ f a 0 ) 0 in this case, we obtain for P ( E 1 0 ; λ ) the presentation
P ( E 1 0 ; λ ) = ( D 1 λ ) 2 ( D 1 k d i s λ ) ( D 1 k h y d , l i λ ) ( D 1 k h y d , p r λ ) × ( D 1 k h y d , l i λ ) ( D 1 + Y f a μ f a ( S ¯ f a 0 ) λ ) × R ( E 1 0 ; λ ) .

References

  1. Ahring, B. (Ed.) Biomethanation, I and II; Springer: Berlin/Heidelberg, Germany, 2003. [Google Scholar]
  2. Deublein, D.; Steinhauser, A. Biogas from Waste and Renewable Resources; Wiley-VCH Verlag: Weinheim, Germany, 2008. [Google Scholar]
  3. Musa Egieya, J.; Čuček, L.; Zirngast, K.; Isafiade, A.J.; Pahor, B.; Kravanja, Z. Synthesis of biogas supply networks using various biomass and manure types. Comput. Chem. 2019, 122, 129–151. [Google Scholar] [CrossRef]
  4. Aceves-Lara, C.-A.; Latrille, E.; Steyer, J.-P. Optimal control of hydrogen production in a continuous anaerobic fermentation bioreactor. Int. J. Hydrog. Energy 2010, 35, I0710–I0718. [Google Scholar] [CrossRef]
  5. Ausiello, A.; Micoli, L.; Turco, M.; Toscano, G.; Florio, C.; Pirozzi, D. Biohydrogen production by dark fermentation of Arundo donax using a new methodology for selection of H2-producing bacteria. Int. J. Hydrog. Energy 2017, 42, 30599–30612. [Google Scholar] [CrossRef]
  6. Guo, X.M.; Trably, E.; Latrille, E.; Carrere, H.; Steyer, J.-P. Hydrogen production from agricultural waste by dark fermentation: A review. Int. J. Hydrog. Energy 2010, 36, 10660–10673. [Google Scholar] [CrossRef]
  7. Pakarinen, O.M.; Kaparaju, P.L.N.; Rintala, J.A. Hydrogen and methane yields of untreated, water-extracted and acid (HCl) treated maize in one- and two-stage batch assays. Int. J. Hydrog. Energy 2011, 36, 14401–144407. [Google Scholar] [CrossRef]
  8. Ruggeri, B.; Tommasi, T.; Sanfilippo, S. BioH2 & BioCH4 through Anaerobic Digestion (From Research to Full-Scale Applications); Springer: London, UK, 2015. [Google Scholar]
  9. Gerardi, M.H. The Microbiology of Anaerobic Digesters; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 2003. [Google Scholar]
  10. Cárdenas, E.L.M.; Zapata-Zapata, A.D.; Kim, D. Hydrogen Production from Coffee Mucilage in Dark Fermentation with Organic Wastes. Energies 2019, 12, 71. [Google Scholar] [CrossRef] [Green Version]
  11. Xu, L.; Wang, Y.; Shah, S.A.A.; Zameer, H.; Solangi, Y.A.; Walasai, G.D.; Siyal, Z.A. Economic viability and environmental efficiency analysis of hydrogen production processes for the decarbonization of energy systems. Processes 2019, 7, 494. [Google Scholar] [CrossRef] [Green Version]
  12. Prapinagsorn, W.; Sittijunda, S.; Reungsang, A. Co-Digestion of Napier Grass and Its Silage with Cow Dung for Bio-Hydrogen and Methane Production by Two-Stage Anaerobic Digestion Process. Energies 2018, 11, 47. [Google Scholar] [CrossRef] [Green Version]
  13. Al-Rubaye, H.; Karambelkar, S.; Shivashankaraiah, M.M.; Smith, J.D. Process simulation of two-stage anaerobic digestion for methane production. Biofuels 2019, 10, 181–191. [Google Scholar] [CrossRef]
  14. Farhat, A.; Miladi, B.; Hamdi, M.; Bouallagui, H. Fermentative hydrogen and methane co-production from anaerobic co-digestion of organic wastes at high loading rate coupling continuously and sequencing batch digesters. Environ. Sci. Pollut. Res. 2018, 25, 27945–27958. [Google Scholar] [CrossRef]
  15. Nasr, N.; Hafez, H.; Naggar, M.H.E.; Nakhla, G. Application of artificial neural networks for modeling of biohydrogen production. Int. J. Hydrog. Energy 2013, 38, 3189–3195. [Google Scholar] [CrossRef] [Green Version]
  16. Wang, J.; Wan, W. Kinetic models for fermentative hydrogen production: A review. Int. J. Hydrog. Energy 2009, 34, 3313–3323. [Google Scholar] [CrossRef]
  17. Batstone, D.J.; Keller, J.; Angelidaki, I.; Kalyuzhnyi, S.V.; Pavlostathis, S.G.; Rozzi, A.; Sanders, W.T.M.; Siegrist, H.; Vavilin, V.A. The IWA Anaerobic Digestion Model No. 1 (ADM1). Water Sci. Technol. 2002, 45, 65–73. [Google Scholar] [CrossRef] [PubMed]
  18. Dochain, D.; Vanrolleghem, P.A.V. Dynamical Modeling and Estimation in Wasterwater Treatement Process; IWA Publishing: London, UK, 2001. [Google Scholar]
  19. Neba, F.A.; Asiedu, N.Y.; Addo, A.; Morken, J.; Østerhus, S.W.; Seidu, R. Biodigester rapid analysis and design system (B–RADeS): A candidate attainable region-based simulator for the synthesis of biogas reactor structures. Comput. Chem. 2020, 132, 106607. [Google Scholar] [CrossRef]
  20. Simeonov, I. Dynamical modeling and estimation in wasterwater treatement process. In Contemporary Approaches to Modeling, Optimisation and Control of Biotechnological Processes; Tzonkov, S., Ed.; Prof. M. Drinov Acad. Publ. House: Sofia, Bulgaria, 2010. [Google Scholar]
  21. Borisov, M.; Dimitrova, N.; Simeonov, I. Mathematical Modeling of anaerobic digestion with hydrogen and methane production. In Proceedings of the 6th IFAC Conference on Foundations of Systems Biology in Engineering, the International Federation of Automatic Control, Magdeburg, Germany, 9–12 October 2016; IFAC–PapersOnLine (TuPP.2: 1–8). Volume 49, pp. 231–238. [Google Scholar]
  22. Chorukova, E.; Simeonov, I. Mathematical modeling of the anaerobic digestion in two-stage system with production of hydrogen and methane including three intermediate products. Int. J. Hydrog. Energy 2020, 45, 11550–11558. [Google Scholar] [CrossRef]
  23. Simeonov, I.; Chorukova, E. Mathematical Modeling of the Anaerobic Digestion with Production of Hydrogen and Methane. In Proceedings of the 4th International Conference on Water, Energy and Environment (ICWEE), Burgas, Bulgaria, 1–3 June 2016; pp. 32–38. [Google Scholar]
  24. Guellout, Z.; Clion, V.; Benguerba, Y.; Dumas, C.; Ernst, B. Study of the dark fermentative hydrogen production using modified ADM1 models. Biochem. Eng. J. 2018, 132, 9–19. [Google Scholar] [CrossRef]
  25. Schievano, A.; Tenca, A.; Lonati, S.; Manzini, E.; Adani, F. Can two-stage instead of one-stage anaerobic digestion really increase energy recovery from biomass? Appl. Energy 2014, 124, 335–342. [Google Scholar] [CrossRef]
  26. Donoso-Bravo, A.; Gajardo, P.; Sebbah, M.; Vicencio, D. Comparison of performance in an anaerobic digestion process: One-reactor vs two-reactor configurations. Math. Biosci. Eng. 2019, 16, 2447–2465. [Google Scholar] [CrossRef]
  27. Blumensaat, F.; Keller, J. Modeling of two-stage anaerobic digestion using the IWA Anaerobic Digestion Model No. 1 (ADM1). Water Res. 2005, 39, 171–183. [Google Scholar] [CrossRef]
  28. Najdenski, H.; Ilyin, V.; Angelov, P.; Hubenov, V.; Korshunov, D.; Kussovski, V.; Dimitrova, L.; Simeonov, I. Laboratory biodegradation of potential cellulose wastes generated during long-term manned space missions. Ecol. Eng. Env. Prot. 2019, 1, 71–78. [Google Scholar] [CrossRef]
  29. Denchev, D.; Hubenov, V.; Simeonov, I.; Kabaivanova, L. Biohydrogen production from lignocellulosic waste with anaerobic bacteria. In Proceedings of the 4th International Conference on Water, Energy and Environment (ICWEE), Burgas, Bulgaria, 1–3 June 2016; pp. 7–12. [Google Scholar]
  30. Simeonov, I. Mathematical modeling and parameters estimation of anaerobic fermentation process. Bioprocess Eng. 1999, 21, 377–381. [Google Scholar] [CrossRef]
  31. Rosen, C.; Jeppsson, U. Aspects on ADM1 Implementation within the BSM2 Framework; Preprint CODEN:LUTEDX/(TEIE-7224); Department of Industrial Electrical Engineering and Automation, Lund University: Lund, Sweden, 2006; pp. 1–35. [Google Scholar]
  32. Hess, J.; Bernard, O. Design and study of a risk management criterion for an unstable anaerobic wastewater treatment process. J. Process. Control 2008, 18, 71–79. [Google Scholar] [CrossRef]
  33. Wiggins, S. Introduction to Applied Nonlinear Dynamical Systems and Chaos; Springer: New York, NY, USA, 1990. [Google Scholar]
  34. Bayen, T.; Gajardo, P. On the steady state optimization of the biogas production in a two-stage anaerobic digestion model. J. Math. Biol. 2019, 78, 1067–1087. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  35. Li, Q.; Li, Y. Coproduction of hydrogen and methane in a CSTR-IC two-stage anaerobic digestion system from molasses wastewater. Water Sci. Technol. 2019, 79, 270–277. [Google Scholar] [CrossRef] [PubMed]
  36. Simeonov, I.; Karakashev, D. Mathematical Modeling of the anaerobic digestion including the syntrophic acetate oxidation. In Proceedings of the 7th Vienna International Conference on Mathematical Modeling (MATHMOD), Vienna, Austria, 14–17 February 2012; Volume 7, pp. 304–309. [Google Scholar]
  37. Borisov, M.; Denchev, D.; Simeonov, I. Mathematical Modeling of a two-stage anaerobic digestion process with hydrogen and methane production using ADM1. Ecol. Eng. Environ. Prot. 2020, 1, 18–29. [Google Scholar]
  38. Dareioti, M.A.; Kornaros, M. Effect of hydraulic retention time (HRT) on the anaerobic co-digestionof agro-industrial wastes in a two-stage CSTR system. Bioresour. Technol. 2014, 167, 407–415. [Google Scholar] [CrossRef]
  39. Cavinato, C.; Bolzonella, D.; Faton, F.; Cecchi, F.; Pavan, P. Optimization of two-phase thermophilic anaerobic digestion of biowaste for hydrogen and methane production through reject water recirculation. Bioresour. Technol. 2011, 102, 8605–8611. [Google Scholar] [CrossRef]
  40. Intanoo, P.; Rangsanvigit, P.; Malakul, P.; Chavadej, S. Optimization of separate hydrogen and methane production from cassava wastewater using two-stage upflow anaerobic sludge blanket reactor (UASB) system under thermophilic operation. Bioresour. Technol. 2014, 173, 256–265. [Google Scholar] [CrossRef]
  41. Khan, M.A.; Ngo, H.H.; Guo, W.S.; Liu, Y.; Nghiem, L.D.; Hai, F.I.; Deng, L.J.; Wang, J.; Wu, Y. Optimization of process parameters for production of volatile fatty acid, biohydrogen and methane from anaerobic digestion. Bioresour. Technol. 2016, 219, 738–748. [Google Scholar] [CrossRef]
Figure 1. Two-phases process of AD with production of hydrogen ( H 2 ) and methane ( C H 4 ).
Figure 1. Two-phases process of AD with production of hydrogen ( H 2 ) and methane ( C H 4 ).
Processes 08 00791 g001
Figure 2. Left plot: the equilibrium point S ¯ a a (solid line) and the upper bound of S ¯ a a (dash line) according to (35) as functions of D 1 . Right plot: the graph of the equilibrium point X ¯ s u a a . The vertical dash-dot lines pass through D 1 ( 2 ) .
Figure 2. Left plot: the equilibrium point S ¯ a a (solid line) and the upper bound of S ¯ a a (dash line) according to (35) as functions of D 1 . Right plot: the graph of the equilibrium point X ¯ s u a a . The vertical dash-dot lines pass through D 1 ( 2 ) .
Processes 08 00791 g002
Figure 3. Existence of the equilibrium component S ¯ a c , c h 4 for D 2 ( 0 , D 2 ( 2 ) ) (left plot) and for D 2 ( 0 , D 2 ( 1 ) ) (right plot).
Figure 3. Existence of the equilibrium component S ¯ a c , c h 4 for D 2 ( 0 , D 2 ( 2 ) ) (left plot) and for D 2 ( 0 , D 2 ( 1 ) ) (right plot).
Processes 08 00791 g003
Figure 4. Eigenvalues of Δ ( 3 ) evaluated at ( S ¯ s u , S ¯ a a , X ¯ s u a a ) for D 1 ( 0 , D 1 ( 1 ) ) . The solid box on the horizontal axis denotes D 1 ( 1 ) .
Figure 4. Eigenvalues of Δ ( 3 ) evaluated at ( S ¯ s u , S ¯ a a , X ¯ s u a a ) for D 1 ( 0 , D 1 ( 1 ) ) . The solid box on the horizontal axis denotes D 1 ( 1 ) .
Processes 08 00791 g004
Figure 5. Eigenvalues of Δ ( 3 ) evaluated at ( S ¯ s u , S ¯ a a , X ¯ s u a a ) for D 1 ( D 1 ( 1 ) , D 1 ( 2 ) ) . The solid box and solid circle on the horizontal axis denote D 1 ( 1 ) and D 1 ( 2 ) , respectively.
Figure 5. Eigenvalues of Δ ( 3 ) evaluated at ( S ¯ s u , S ¯ a a , X ¯ s u a a ) for D 1 ( D 1 ( 1 ) , D 1 ( 2 ) ) . The solid box and solid circle on the horizontal axis denote D 1 ( 1 ) and D 1 ( 2 ) , respectively.
Processes 08 00791 g005
Figure 6. The input–output static characteristics Q h 2 ( D 1 ) Q h 2 0 ( D 1 ) for D 1 ( 0 , D 1 ( 2 ) ) (left) and Q c h 4 0 ( D 1 , m a x 0 , D 2 ) for D 2 ( 0 , D 2 ( 1 ) ) (right). The vertical dash-dot lines in the left plot pass through D 1 ( 1 ) (solid box) and D 1 ( 2 ) (solid circle). In the right plot, D 2 ( 1 ) is marked by a sold box.
Figure 6. The input–output static characteristics Q h 2 ( D 1 ) Q h 2 0 ( D 1 ) for D 1 ( 0 , D 1 ( 2 ) ) (left) and Q c h 4 0 ( D 1 , m a x 0 , D 2 ) for D 2 ( 0 , D 2 ( 1 ) ) (right). The vertical dash-dot lines in the left plot pass through D 1 ( 1 ) (solid box) and D 1 ( 2 ) (solid circle). In the right plot, D 2 ( 1 ) is marked by a sold box.
Processes 08 00791 g006
Figure 7. Time evolution of Q h 2 and Q c h 4 with variable D 1 and D 2 from Table 3.
Figure 7. Time evolution of Q h 2 and Q c h 4 with variable D 1 and D 2 from Table 3.
Processes 08 00791 g007
Figure 8. Dynamic behavior of the model solutions with variable D 1 from Table 3.
Figure 8. Dynamic behavior of the model solutions with variable D 1 from Table 3.
Processes 08 00791 g008
Figure 9. Dynamic behavior of the model solutions with variable D 1 from Table 3.
Figure 9. Dynamic behavior of the model solutions with variable D 1 from Table 3.
Processes 08 00791 g009
Figure 10. Time evolution of Q h 2 and Q c h 4 with variable X c i n from Table 4.
Figure 10. Time evolution of Q h 2 and Q c h 4 with variable X c i n from Table 4.
Processes 08 00791 g010
Figure 11. Dynamic behavior of the model solutions with variable X c i n from Table 4.
Figure 11. Dynamic behavior of the model solutions with variable X c i n from Table 4.
Processes 08 00791 g011
Figure 12. Dynamic behavior of the model solutions with variable X c i n from Table 4.
Figure 12. Dynamic behavior of the model solutions with variable X c i n from Table 4.
Processes 08 00791 g012
Table 1. Model Variables
Table 1. Model Variables
Definition of the Model Variables
S s u     concentration of monosacharides [ g C O D / L ]
S a a     concentration of amino acids (AA) [ g C O D / L ]
S f a     concentration of fatty acids (LCFA) [ g C O D / L ]
S a c     concentration of total acetate in B R 1 [ g C O D / L ]
S a c , c h 4     concentration of total acetate in B R 2 [ g C O D / L ]
X c     concentration of composites [ g C O D / L ]
X c h     concentration of carbohydrates [ g C O D / L ]
X p r     concentration of proteins [ g C O D / L ]
X l i     concentration of lipids [ g C O D / L ]
X s u a a     concentration of sugar and AA degraders [ g / L ]
X f a     concentration of LCFA degraders [ g / L ]
X a c     concentration of acetate degraders [ g / L ]
D 1 dilution rate in the first bioreactor [ h 1 ]
D 2 dilution rate in the second bioreactor [ h 1 ]
Table 2. Model Parameters
Table 2. Model Parameters
Definition of the Model ParametersValues
S s u i n     input concentration of S s u [ g C O D / L ] 0.01
S a a i n     input concentration of S a a [ g C O D / L ] 0.001
S f a i n     input concentration of S f a [ g C O D / L ] 0.001
X c i n     input concentration of X c [ g C O D / L ] 50
f c h , x c stoichiometric parameter [–]0.2
f p r , x c stoichiometric parameter [–]0.2
f l i , x c stoichiometric parameter [–]0.3
f s u , l i stoichiometric parameter [–]0.05
f f a , l i stoichiometric parameter [–]0.95
f a c , s u stoichiometric parameter [–]0.41
f a c , a a stoichiometric parameter [–]0.4
Y s u a a stoichiometric parameter [–]0.1
Y a c stoichiometric parameter [–]27.3
Y f a stoichiometric parameter [–]0.06
Y h 2 , s u physicochemical parameter [ L 2 / g ] 0.7
Y h 2 , a a physicochemical parameter [ L 2 / g ] 0.7
Y h 2 , f a physicochemical parameter [ L 2 / g ] 0.7
Y c h 4 , a c physicochemical parameter [ L 2 / g ] 75
k d i s biochemical parameter [ h 1 ] 0.0208
k h y d , c h biochemical parameter [ h 1 ] 0.417
k h y d , p r biochemical parameter [ h 1 ] 0.417
k h y d , l i biochemical parameter [ h 1 ] 0.417
k m , s u a a biochemical parameter [ h 1 ] 1.25
k s , s u a a biochemical parameter [ g / L ] 0.5
k m , f a biochemical parameter [ h 1 ] 0.15
k m , a c biochemical parameter [ h 1 ] 0.0167
k s , f a biochemical parameter [ g / L ] 0.67
k s , a c biochemical parameter [ g / L ] 0.4
Table 3. Values of D 1 and D 2 .
Table 3. Values of D 1 and D 2 .
Time [ h ] 0–50005000–10,00010,000–15,00015,000–20,000
D 1 [ 1 / h ] 0.005 0.01 0.015 0.02
D 2 [ 1 / h ] 0.00113 0.00226 0.00339 0.00452
Table 4. Values of X c i n .
Table 4. Values of X c i n .
Time [ h ] 0–50005000–10,00010,000–15,00015,000–20,000
X c i n [ g C O D / L ] 40507560

Share and Cite

MDPI and ACS Style

Borisov, M.; Dimitrova, N.; Simeonov, I. Mathematical Modeling and Stability Analysis of a Two-Phase Biosystem. Processes 2020, 8, 791. https://0-doi-org.brum.beds.ac.uk/10.3390/pr8070791

AMA Style

Borisov M, Dimitrova N, Simeonov I. Mathematical Modeling and Stability Analysis of a Two-Phase Biosystem. Processes. 2020; 8(7):791. https://0-doi-org.brum.beds.ac.uk/10.3390/pr8070791

Chicago/Turabian Style

Borisov, Milen, Neli Dimitrova, and Ivan Simeonov. 2020. "Mathematical Modeling and Stability Analysis of a Two-Phase Biosystem" Processes 8, no. 7: 791. https://0-doi-org.brum.beds.ac.uk/10.3390/pr8070791

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop