Early warning indicators for mesophilic anaerobic digestion of corn stalk: a combined experimental and simulation approach

Background Monitoring and providing early warning are essential operations in the anaerobic digestion (AD) process. However, there are still several challenges for identifying the early warning indicators and their thresholds. One particular challenge is that proposed strategies are only valid under certain conditions. Another is the feasibility and universality of the detailed threshold values obtained from different AD systems. In this article, we report a novel strategy for identifying early warning indicators and defining threshold values via a combined experimental and simulation approach. Results The AD of corn stalk (CS) was conducted using mesophilic, completely stirred anaerobic reactors. Two overload modes (organic and hydraulic) and overload types (sudden and gradual) were applied in order to identify early warning indicators of the process and determine their threshold values. To verify the selection of experimental indicators, a combined experimental and simulation approach was adopted, using a modified anaerobic bioconversion mathematical model (BioModel). Results revealed that the model simulations agreed well with the experimental data. Furthermore, the ratio of intermediate alkalinity to bicarbonate alkalinity (IA/BA) and volatile fatty acids (VFAs) were selected as the most potent early warning indicators, with warning times of 7 days and 5–8 days, respectively. In addition, IA, BA, and VFA/BA were identified as potential auxiliary indicators for diagnosing imbalances in the AD system. The relative variations for indicators based on that of steady state were observed instead of the absolute threshold values, which make the early warning more feasible and universal. Conclusion The strategy of a combined approach presented that the model is promising tool for selecting and monitoring early warning indicators in various corn stalk AD scenarios. This study may offer insight into industrial application of early warning in AD system with mathematical model. Electronic supplementary material The online version of this article (10.1186/s13068-019-1442-7) contains supplementary material, which is available to authorized users.


Background
Anaerobic digestion (AD), as an efficient technology for organic waste treatment, has been widely adopted worldwide [1]. Among others, straw has a great potential to serve as a feedstock for anaerobic methane (CH 4 ) production, due to its abundance and suitable bioconversion characteristics [2]. In China, straws are produced at a high annual rate of approximately 1 billion metric tons [3], and about 30% of those were underutilized [3,4]. In addition, the methanogenic potential of some major straws is 2.86-3.78 × 10 5 Nm 3 CH 4 /kg VS [5]. However, straw is consisting of cellulose, hemicellulose, and lignin, which is a typical high carbon-to-nitrogen (C/N) substrate [2,6]. And excessive volatile fatty acids (VFAs), which are intermediates of AD, may be produced when feedstock overloading occurs, especially with high C/N (e.g., > 30) [2,5,7]. Besides, the complex structure of lignocellulose is difficult for microbial cellulolytic enzymes to access, limiting degradation [8]. Therefore, some previous studies reported that the moderate organic loading

Open Access
Biotechnology for Biofuels *Correspondence: duanna@cau.edu.cn 1 College of Water Resources and Civil Engineering, China Agricultural University, Beijing 100083, China Full list of author information is available at the end of the article rates (OLRs) were very vital to avoid system acidification. The thermophilic AD process of straw can only be stably operated under relatively low organic loading rates (OLRs) and below 2 kg VS/(m 3 day) [9]. Meanwhile, Li et al. [10] suggested that the mesophilic anaerobic codigestion of rice stalk with cow manure should be operated at an OLR of 3-6 kg VS/(m 3 day). Ward et al. [11] found that biogas projects using straw as substrate also tend to be controlled at suboptimal OLRs to prevent process inhibition.
Besides optimal operational parameters like OLR, reliable early warning and regulation systems are also favorable for AD process. Previously, many studies have been carried out to explore feasible warning indicators in different AD systems. Some examples for several indicators proposed include VFA concentrations, alkalinity, biogas composition, specific intermediate metabolite (like glycerol, aromatic compound, etc.) concentrations, microbial community composition, and enzyme activity [12][13][14][15][16][17]. In addition, some coupled indicators, such as the ratio between intermediate and partial alkalinities or VFA concentrations and bicarbonate alkalinity (BA), showed better performance than individual indicators [18]. However, when comparing the results of previous studies, it must be pointed out that the proposed strategies are only valid under certain conditions, as some parameters may have different sensitivities to environmental fluctuations in different AD systems. For instance, Castellano et al. [19] suggested that hydrogen (H 2 ) concentrations have a high discriminatory ability for process state identification. On the contrary, Kleyböcker et al. [20] showed that H 2 partial pressure was not an ideal indicator, because of its unstable responses under organic overload conditions in an AD system treating rapeseed oil. In addition, the indicator threshold values might vary in AD systems with different substrates and operating conditions. In a specific case, Pullammanappallil et al. [21] suggested that the critical value of propionic acid was 2750 mg/L, while Holm-Nielsen et al. [22] found it to be 1500 mg/L. Conversely, propionic acid did not show early warning in some trials [20]. Consequently, finding effective indicators and rational threshold values for early warning and inhibition diagnosis is a challenge when working with AD systems.
Compared with traditional early warning methods, which only monitor process indicators by chemical pathway, modeling the AD process can provide a flexible and rapid solution for comparing and evaluating large numbers of such indicators, and it provides the possibility of automated warning for industrial applications. Mathematical models have long been used for the simulation of various AD scenarios, with several computer-aided implementations in existence [23,24]. Unlike simple models that are mainly used for calculating theoretical biogas and CH 4 yields, complex bioconversion models can be used to generate insights about process kinetics, microbial growth inhibition, substrate conversion, and product generation rates, to mention a few of their functionalities. In addition, these models can also handle extensive amounts of numerical data and provide qualitative or quantitative comparisons between the measured and simulated datasets. Hence, using these tools for the evaluation of early warning indicators appears to be a promising approach.
The aims of the present study were therefore to (1) evaluate the reaction of various process parameters during AD process of corn stalk fed to continuously stirred tank reactors under different overload modes (organic or hydraulic retention time (HRT) overload; and gradual or sudden overload); (2) compare simulation using a proven bioconversion model and experimental results at the same operational conditions; (3) identify and evaluate response parameters that are sufficiently sensitive to environmental disturbances; and (4) define threshold values for the sensitive indicators identified in (3), for the different overload conditions tested.

Gaseous parameters
The CH 4 yield and content of the two reactors were seen to have stabilized in the full-load phase, at 0.20 L CH 4 /g VS (62.20% of CH 4 ) on average with OLR of 1.50-2.24 g VS/(L day) in R1 and at 0.20 L CH 4 /g VS (58.73% of CH 4 ) with OLR between 1.87 and 2.24 g VS/(L day) in R2. During the gradual organic overload phase in R1 (day 101-113), both the CH 4 yield and CH 4 content decreased significantly to 0.02 L CH 4 /g VS and 51.76% on day 113. Regarding the sudden overload phase of R1, the CH 4 yield showed a stepwise rise with elevated OLR, up to 0.37 L CH 4 /g VS on day 155 with OLR at 3.37 g VS/(L day), and then sharply decreased to 0.12 L CH 4 /g VS (day 161). With respect to the gradual HRT overload phase of R2, the CH 4 yield showed an acidification response under OLR of 2.81-3.74 g VS/(L day). Subsequently, CH 4 production increased slowly and then decreased sharply after day 161. For each overload condition, the CH 4 content decreased slightly in the beginning and then returned to steady-state level (Fig. 1a1, a2). This result was consistent with a previous study, where the CH 4 content did not decrease significantly, as long as the pH was higher than 5.5 [1]; moreover, it reported that raising OLR did not result in a significant change in CH 4 content. However, a clear shift in populations of archaea from acetotrophic to hydrogenotrophic methanogenesis may have occurred [25].
On the other hand, solids evidently accumulated in the reactor, probably due to straw floating and inadequate stirring, which is approximately 5-6% in the last phase of the experiment (Additional file 1: Fig. S1). Solids accumulated could lead to gas-liquid phase transfer delay and oversaturation of CH 4 in the liquid phase [26]. Thereby, above gaseous parameters were relatively insensitive. Meanwhile, the ratio of methane and carbon dioxide (CH 4 /CO 2 ) showed more intense response to perturbations (Figs. 1b and 2b), where sudden drops were Fig. 1 A comparison of experimental and simulated CH 4 yields (a1, a2) and CH 4 /CO 2 ratios (b1, b2) related to the laboratory experiments. Data presented in subplots a1 and b1 represent R1, while those in subplots a2 and b2 represent R2 Fig. 2 A comparison of experimental and simulated acetic acid (a1, a2), propionic acid (b1, b2), butyric acid (c1, c2), and total VFA (d1, d2) concentrations related to the laboratory experiments. Data presented in subplots a1, b1, c1, and d1 represent R1, and those in subplots a2, b2, c2, and d2 represent R2 observed on days 112 and 161 in R1, and on day 161 in R2. Further, it realized a remarkable early warning potential in comparison with CH 4 yield (8 days earlier) under sudden overload, but barely provided any early warning under gradual overload. Previously, H 2 concentration was also suggested as a useful variable for detection of disturbances in both carbohydrate-and protein-based wastewaters [27]. However, in the current study, H 2 was only detected in the overload stage on days 112 and 155 in R1, and on day 159 in R2. In addition, these fluctuations were faint and short (1-2 days), with the maximum H 2 content of 0.17% being reached in the sudden overload phase of R1. It is speculated that the pH and total solid (TS) concentration of the effluent can affect the discrimination of H 2 .

VFA parameters
As far as the full-load phase is concerned, the total VFA concentrations of the two reactors were found stable, with average values of 0.49 g/L (R1) and 0.60 g/L (R2), respectively. However, the total VFA concentration increased rapidly during the overload phase. Acetic acid was the most abundant acid, accounting for approximately 40% of the total VFA, hence it dominates total VFA concentration changes. Its sharp increase occurred on day 109 (from 0.21 to 0.39 g/L) and day 147 (from 0.23 to 0.59 g/L) in R1, and on day 155 (from 0.34 to 0.52 g/L) in R2. The sudden rise of total VFA was also occurred on day 109 (from 0.37 to 0.80 g/L) and day 147 (from 0.56 to 0.95 g/L) in R1, and on day 155 (from 0.82 to 2.15 g/L) in R2. The quick accumulation of acetic acid indicated an imbalance between the acid-forming phase and methaneforming phase of the digestion process, which was in line with the changes seen in CH 4 yield and CH 4 /CO 2 (Fig. 1). Propionic acid concentration was fairly stable during the full-load phase, especially in R2 (Fig. 2c1, c2). Its sharp increase was observed on day 109 (from 0.06 to 0.30 g/L) and day 153 (from 0.32 to 0.50 g/L) in R1, and on day 155 (from 0.18 to 0.98 g/L) in R2. It is noticeable that in R2, the sharp increases of VFAs concentration were always followed by the gradual returning to its steady-state concentration, and propionic acid showed the slowest recovery ( Fig. 2), which is in agreement with a previous study [22], but different to the finding of Boe et al. [28] stating that acetate revealed recovery while propionate was persistent. In the case of methanogenic populations, possible adaptation to the higher VFA concentrations may have occurred, which was probably the reason for reduced acidification [29]. What's more, a latest research reported that there is a kind of methanogen could uptake acetate but also propionate directly, which might be one of the reasons as well [30]. The findings about slow propionic acid reduction were also supported by the findings from Ahring et al. [31], who showed that propionate acid degraders are the slowest growing and most sensitive VFA-degrading microorganisms in the AD process.
High total alkalinity (TA) (about 6000 g CaCO 3 /L) was detected in both reactors at the beginning, due to the high alkalinity contained in the inoculum originating from the full-scale biogas plant. During the experiment, TA was maintained at an average of 1814.3 mg CaCO 3 /L (R1) and 1871.9 mg CaCO 3 /L (R2) (Fig. 3a1, a2). The result was also in accordance with a previous study, in that TA could remain stable until the pH fell below 4.3 as a consequence of high VFA concentration [32]. BA suddenly dropped on days 109 and 158 in R1, and day 154 in R2 ( Fig. 3c1, c2). Conversely, intermediate alkalinity (IA) suddenly rose on days 105 and 157 in R1, and day 153 in R2 ( Fig. 3b1, b2). Although both indicators showed fluctuation under overload conditions, IA proved to be more sensitive, since it got out of balance earlier.
From another perspective, coupled indexes appeared to be more sensitive than single indicators, thus their distinct fluctuations made the identification of acidification response easier. BA/TA decreased sharply on days 105 and 158 in R1, and on day 154 in R2 ( Fig. 3d1, d2). Conversely, IA/BA showed a sharp increase on days 105 (from 0.28 to 0.53) and 157 (from 0.44 to 0.54) in R1, and day 153 (from 0.48 to 0.68) in R2 (Fig. 3e1, e2). Unfortunately, all alkalinity indicators showed obvious signs of response delay under sudden overload. It is similar to the report of Li et al. [33] where the warning times of IA/ PA (where PA is the partial alkalinity and is analogous to BA) and BA/TA were shortened by 6 days and 2 days, respectively. Ahring et al. [12] also reported that most indicators were suitable for detecting gradual overloads, but were too slow to respond under sudden overload. Significant increase of VFA/BA was found on days 113 and 151 in R1, and day 155 in R2 ( Fig. 3f1, f2), and was confirmed by the work of Li et al. [33], who also found VFA/BA increase due to higher OLR in mesophilic AD of vegetable waste.

Simulation results using the BioModel
The results of the experimental reactor operation simulations are presented in terms of CH 4 yield and CH 4 /CO 2 (Fig. 1), and individual and total VFA concentrations ( Fig. 2). From a qualitative point of view and by looking at the fits between experimental and simulated data curves, it appears that the model was mostly successful in capturing the overall gas and VFA production trends of the two experiments with high accuracy. At the same time, by evaluating the goodness of these fits using statistical measures, relative root mean squared error (rRMSE) and mean absolute percentage error (MAPE) values provide a more detailed overview of simulation accuracy (Table 1). By comparing the visual and statistical measures, biogas and CH 4 yield, along with butyric and acetic acid concentration simulations appear to be the most accurate, with the lowest rRMSE and MAPE only slightly above the feasible range (Table 1, values for "Day 0-165"). This deviance is due to a few sections of the dataset, where the simulation-measurement fit was not satisfactory. For example, the gas yield and VFA concentration levels during the reactor startup periods were not matched in the simulations. However, reactor startup periods inherently involve significant stochasticity and, depending on the substrates and reactor history, potential microbial growth lag [34], which are hard to simulate with fixed kinetic models. Deeming this period irrelevant from the perspective of early  (Table 1, values for "Day 30-165") have improved significantly. On the other hand, VFA simulations showed varying rates of change, with those of R1 being more positive than R2. A more important difference between the measured and simulated values was, on the other hand, seen in R1, around days 110-120. Here, the drop in CH 4 yield and the consecutive process failure were slower in the simulation results compared to the laboratory experience, by approximately 2 days. At the same time, the drop was preceded by a slight increase in the respective yields, potentially owing to an initial positive response of simulated acetoclastic methanogens to an increased OLR. This means that in the model, the immediate response of microbial groups to an increase in substrate availability is a proportional increase in their productivity, which is in most cases followed by their negative response to the gradual accumulation of certain inhibitory compounds. In the case of acetoclastic methanogens, model inhibition was assumed by free ammonia and the saturation of volatile fatty acids, of which the latter was seen both experimentally and during the simulations. However, as the levels of acid saturation were significantly lower in the model simulations than in the physical reactors, inhibition during the simulations appeared to take effect slower than seen experimentally. Nevertheless, the simulations for both R1 and R2 were in general found to be in good agreement with experimental data, therefore they could be used for comparing the experimentally defined early warning indicators with their simulated counterparts. In the short term, this provided an additional method for verifying the quality of indicator selection. In addition, the long-term benefits of such simulations lie in their rapid generality and continuous interpretability: both contributing to the reduction of necessary analytic measurements and process data density. Thus, by means of simulated early warning indicators, monitoring, and forecasting, the fate of the experimental processes could be improved significantly.

Early warning indicators and threshold values The procedure of screening early warning indicators
Experimental and simulated data were used to screen potentially optimal warning indicators, based on the proposition of three important and mostly qualitative criteria. Firstly, optimal indicators had to show high sensitivity to changes, in the sense that there was enough warning time between an indicator's response point and process failure. Stable acidification response was another vital aspect, as excessive indicator sensitivity to acids may have led to false assumptions. Finally, the low cost of monitoring is essential in practice [20,35], therefore indicators had to be measured economically.
Based on the above selection criteria, reference points for reactor failure are needed for measuring warning time of different parameters. In the current study, these points of failure were declared when at least one of two events happened. The first one concerned the reduction of reactor pH to less than 6.4, which is well below the optimal pH range commonly reported for methanogenic archaea [36]. Meanwhile, the other event was the significant reduction of CH 4 yield (relative standard deviation (RSD) > 20%). pH and CH 4 yield are the most intuitive and easily measured indicators, but relatively insensitive for process destabilization, therefore, their changes were also recommended as the symbol of process failure in previous study [37].
In order to quantify the changes of different indicators at the points of failure, their variation amplitudes were calculated. These amplitudes were expressed on a percentage basis and relative to their steady-state values. Specifically, by comparing the variation amplitude of abrupt changes throughout the different overload modes and types applied in the experimental reactors, maximum and minimum values of each indicator were defined.

Screening early warning indicators
According to the experimental data, days 112 and 159 were identified as reference points for measuring the early warning times of different indicators under gradual and sudden overload conditions in R1, and this reference point was day 160 for R2. The abrupt change dates and values, and warning times of different indicators are shown in Tables 2  and 3. Almost all the parameters have performed certain response to the overload shock. Overall, gas phase indicators provided a delayed response in comparison with liquid phase indicators, only CH 4 /CO 2 have 6 days warning time under sudden overload. This could probably be due to the properties of corn stalk (CS), given the floating problem caused by low specific gravity can lead to stirring issues and eventually delays in mass transfer from the liquid to the gas phase [38,39].
It is obvious that parameters' sensitivity is influenced by the overload mode. For the gradual overload, the acidification responses of the various indicators did not show significant difference (3-7 days for individual VFA and 6-7 days for alkalinity) (Tables 2 and 3), regardless of the overload type. Compared with the gradual overload, indicators demonstrated varied acidification responses during sudden overload. More specifically, VFA-related indicators (6-12 days) had longer warning time, while alkalinity and its coupled indicators (1-2 days) had shorter warning time (Tables 2 and  3). This contradiction is probably related to the change of alkalinity. At the beginning of the AD process, the effluent taken from the full-scale plant for inoculation provided a high initial alkalinity in both reactors. Later, the gradual overload applied to R1 generated a large amount of VFA, most of which was consumed by alkalinity that delayed VFA accumulation. By contrast, as the buffer capacity of the AD system declined before the second (sudden) overload phase, VFA accumulated rapidly and showed a longer warning time.
Besides early warning time, stability and measurability of indicators are also important for monitoring biogas plants [19]. For example, IA and IA/BA, as well as BA and BA/TA, showed similar warning times in the two reactors (Tables 2 and 3). However, the average RSD values of these indicators during the full-load phase were 9.57 and 11.59 for IA, 8.23 and 10.31 for IA/BA, 6.26 and 6.94 for BA, and 2.83 and 3.45 for BA/TA, in R1 and R2, respectively. The RSD values of coupled indicators are relatively smaller, which indicates that the coupled indicators are more consistent before overload shock and implies less misjudgment as well. Furthermore, AD is an interrelated process, and the relationship between certain process variables might be the reason for coupled indicators showing better stability than individual ones. Although Li et al. [18] also proposed using coupled indicators to achieve early warning, present study showed that total VFA, acetic, and propionic acids can provide similar warning times in the reactors (Tables 2 and 3). Nevertheless, total VFA can be determined relatively easier than individual VFA in biogas plants [40], making total VFA a more widely acceptable variable for early warning indication. Fortunately, due to the development of detection technology, total VFA and alkalinity can both be monitored online by transducers, titration, or infrared spectroscopy, as well as via online sampling and gas chromatography [41][42][43][44].
In addition to the experiment-based identification of indicator variables with the highest early warning potential, results of the simulations were used to calculate the same variables numerically and evaluate their warning efficiency in comparison to the empirical values. This approach was unlike any previous work on the topic that the authors reviewed [28,[45][46][47][48][49]. Up until now, early warning indicators were selected mostly based on the offline or online monitoring of different values of VFA, pH, alkalinity, biogas fractions, trace elements, and their various combinations, although, in one instance, stable carbon isotopes of CH 4 were also named as potential indicators of process imbalance [50]. These indicator values could then be either evaluated comparatively, or used as inputs to sophisticated control systems. Despite their usually good early warning indicator potential, so far the majority of these solutions remained inapplicable by the industry, mainly due to their limited scope or significant costs involved [51]. Therefore, present simulation-based system could offer a competitive solution for the identification and monitoring of early warning indicators, through its validation by relevant experimental data, application flexibility, customizability at scale, and relatively low price point, when compared to laboratoryintensive processes. The experimental and simulated variables were compared separately for reactors R1 and R2 and in the gradual (Rg) and sudden overload (Rs) modes, thus creating three distinct indicator groups of interest, referred to as R1 g, R1 s, and R2 g from here onwards. In order to provide a reasonable basis for indicator comparison, reference values for identifying the early warning time in the simulations were chosen according to the experimental points defined earlier. Figure 4 shows a comparison of the experimental and simulated indicators for the three indicator groups described above. For each group and individual indicator, the orange bars enclosed with red (lower) and blue (higher) horizontal lines represent the ranges of percentagewise differences calculated between the reference (critical) indicator values and their highest values measured (exp) or simulated (sim) during the  (Tables 2 and 3). As an example, for the indicator group R2 g (Fig. 4c1, c2), the reference value of the experimental propionate indicator (Propionate_exp) was found to be 0.98 g/L on day 155; then during the five steady-state periods identified from the gradual mode operational data of R2, the highest experimental values measured were given as 0.20, 0.23, 0.22, 0.20, and 0.19 g/L, representing, respectively, 390, 326, 345, 390, and 416% changes between their values and the reference value of 0.98 g/L; finally, selecting the lowest (326%) and the highest (416%) changes, the minimum and maximum boundaries of the bar belonging to R2 g Propionate_exp (Fig. 4c1, c2) were given. The ranges for the simulated propionate indicator (Propionate_sim), along with all the other indicators were calculated in an identical manner. In cases where the minimum and maximum indicator values were equal, average indicator values were used for the calculation (see standalone red horizontal lines). Figure 4 shows that the experimental and simulated difference ranges agreed well for some indicators, while for others they were significantly different. Given that most of the indicators were partially or completely dependent on VFA concentrations or alkalinity (VFA and alkalinity being interrelated), changes in these two types of process variables could influence the results to a great extent. Regarding the generation of experimental and simulated VFA data points, a major influencing factor was how accurately the model considered the conversion of initial compounds to intermediate and terminal products, compared to reality. This implies that under experimental conditions, the conversion of complex organic substrates to VFA and then further to gases is a function of a series of stochastic and microbial community-driven events, while model simulations, which are relatively simplified descriptions of reality, are generated assuming structured kinetic equations and stoichiometric yield coefficients. Such fundamental differences could inherently lead to deviations between measured and simulated VFA concentrations, despite the good agreement of the experimental and simulated CH 4 yields (Fig. 1a1, a2) showing that the overall model mass balances are otherwise reliable. On the other hand, the simulation of the various alkalinity fractions posed considerable difficulties. Under laboratory conditions, the alkalinity of a sample is commonly measured through titrimetry, and it is the preferable method for the routine analysis of anaerobic digestion samples as well, due to its speed, simplicity and competitive price [52]. Nonetheless, it involves a significant level of uncertainty when used for the offline measurement of dissolved carbonate concentration in AD samples with inhomogeneous matrices [53]. Meanwhile, as the model cannot simulate titration, simplification was necessary.
It was decided that using the simulation results, the BA of the reactor would be calculated by using the bicarbonate and carbonate ion fractions calculated by the model, while IA would be represented by the sum of acetic, propionic and butyric acids, expressed in terms of acetic acid equivalents. However, this simplification meant that any inaccuracies in the simulation of these compounds, together with the uncertainties brought about by the titration method would potentially cause disagreements between experimental and simulated alkalinity results. For this reason, this factor was considered during the evaluation of the comparative results.
Considering the overall sensitivity, stability and measurability of indicators and based on the experimental and simulation results, IA/BA and VFA were selected as optimal early warning indicators. Further to that, IA, BA and VFA/BA were defined as auxiliary indicators for the diagnosis of the AD system treating CS. IA/BA was also suggested as warning indicator in other lab-scale research [32], and an industrial scale research [38], due to its sensitivity to pronounced changes under overloading conditions. VFA was usually recommended as warning indicator [12,18].

Threshold value
The abrupt change values and their change amplitude in the previously selected main (IA/BA and VFA) and auxiliary (IA, BA and VFA/BA) early warning indicators were investigated in current study. The variation amplitudes of abrupt changes in indicators monitored during experiments and simulation are shown in Fig. 4, and the accurate values were shown in Tables 2 and 3. The abrupt change value of IA/BA was below 0.7, which was 0.53 for organic gradual overload (R1 g), 0.54 for organic sudden overload (R1 s) and 0.68 for hydraulic gradual overload (R2 g), respectively. Meanwhile, Martín-González et al. [32] proposed the critical value of IA/PA to be 0.24 for municipal waste (37 °C), and in another study Ferrer et al. [38] found this number to be 0.72 for sewage sludge (55 °C). This divergence, however, may be the result of differences in feedstock composition and AD operating conditions. Consequently, compared with providing a determined threshold, observing the relative variation for indicators might be a more promising strategy, while evaluating the effectiveness of early warning indicators. This conclusion found further support in previous studies [12,18].
Accordingly, an acidification risk that requires attention would appear in present AD processes when the IA/ BA changed more than 10% in the experiment or 20% in the simulation, based on the steady-state data.
The abrupt change values of VFA were 0.80 g/L (R1 g), 0.95 g/L (R1 s) and 2.15 g/L (R2 g), respectively. Compared to the average values during the steady-state period, if the VFA increased more than 51% (experimental data) or 19% (simulation data), the biogas system would potentially be at risk from instability.
The abrupt change values of IA, BA and VFA/BA in this study were lower than 0.90, higher than 1.05 and lower than 0.80, respectively. When the IA value increased by nearly 20%, BA value decreased about 11% and VFA/BA value increased approximately 30%. This implied that the biogas system was imbalanced, thus attention should be paid and necessary actions might have to be taken to regulate the AD process.
While based on the assessment of several relevant publications, above discussion could potentially be extended by future analyses of the available literature on the early warning indicators in AD. These, together with experiments carried out and evaluated in a manner similar to the one hereby presented, could offer further model verification, deeper insights into the dynamic behavior of such interconnected processes and eventually provide optimized early warning indicators for biogas plants.

Conclusion
Monitoring and providing early warning are essential operations in the AD process. Using a mathematical model to simulate the selected experimental process variables provided good data fit and played a key role in the evaluation of the early warning indicators. Based on both experimental and simulated results, the optimal early warning indicators were identified to be IA/BA and VFA. Besides, IA, BA, and VFA/BA could be used as auxiliary indicators for diagnosing the AD system of CS. It is concluded that this modeling can be a promising tool for monitoring the change signals from early warning indicators and improving the standards of AD plant operation.

Feedstock and inoculum
Corn stalk was obtained from Weichang County, Hebei, China. The collected corn stalk was dried and smashed to approximately 3 mm. The inoculum, obtained from an anaerobic reactor in a wastewater treatment plant (Beijing, China), was acclimated at mesophilic temperature by feeding with pig manure for 2 weeks. The properties of corn stalk and inoculum were shown in Table 4.

Experimental setup
The experiment was carried out using two identical 20 L continuously stirred tank reactors (CSTR) having 17 L working volume. The reactors were maintained at mesophilic conditions (35 ± 1 °C) by a heating water bath (SY-200, Changfeng Instrument and Apparatus Company, Beijing, China) and were continuously mixed at a stirring speed of 60 rpm.
Reactor 1 (R1) operated at stepwise elevated OLR by increasing influent feedstock concentration. The experiment in R1 was divided into two phases, which were denoted as gradual overload (day 0-113 with OLR from 1.50 to 2.99 g VS/(L day)) and sudden overload (day 121-165 with OLR from 1.50 to 3.37 g VS/(L day)). Due to process inhibition, the AD system started recovering from day 114 to 120, during which period a re-inoculation was made, substrate feeding was stopped, and the reactor effluent was recycled as feed. R1 was operated at a fixed HRT of 25 days during the whole experiment. At the same time, Reactor 2 (R2) was operated at a gradually increasing OLR, by shortening HRT and keeping influent feedstock concentration at 6% TS. Both reactors were kept operating until the process completely failed. The operational parameters and periods of the experiment are presented in Table 5.
Produced biogas was collected in a gas bag, and the gas volume was measured by a gas flowmeter (LML-1, Changchun auto filter co., Ltd, Jilin, China). The effluent was drawn daily for the analysis of pH, VFA, TA, BA, IA, TS and VS.

Analytical methods
Total solid (TS) and volatile solid (VS) were measured according to the standard methods [54]. Crude protein was estimated by multiplying the total Kjeldahl nitrogen by 6.25 and the total Kjeldahl nitrogen was measured by a Kjeldahl apparatus (K1305A, Sonnen Automated Analysis Instrument Co., Ltd., Shanghai, China). Crude fiber was determined using a fiber analyzer (Model A220, ANKOM Technology Corporation, NY, USA). Organic elemental components of the corn stalk were determined using an elemental analyzer (Exeter Analytical, Inc. CE-440 Elemental Analyzer, Chicago, USA). Biogas composition (CH 4 , H 2 and carbon dioxide) was determined by a gas chromatograph (1490, Agilent Technologies, USA) equipped with a thermal conductivity detector as previously described [55]. Liquid samples were centrifuged at 4000 rpm for 10 min and then used for the chemical analyses. Before VFA analysis, samples were filtered through a 0.22 μm membrane. The VFA concentrations were measured by a high performance liquid chromatograph (LC-10A, Shimadzu Corporation, Kyoto, Japan), according to the method proposed by [55]. Alkalinity and pH were tested by an automatic potentiometric titrator (ZDJ-4B, Shanghai INESA Scientific Instrument Co., Ltd, China), with a glass and calomel electrode used as the indicator and reference electrode, respectively. For titration, 0.20 mol/L HCl was used as titrant, and the system was calibrated with anhydrous Na 2 CO 3 . IA, PA, and TA were determined using a three-point method [56], by recording the HCl consumption at the respective pH points of 5.75, 4.3, and 3.8, and converting those values to calcium carbonate (CaCO 3 ) equivalents (see Eq. (1)).
where 50.05 is a coefficient used to convert alkalinity units from mEq/L to mg CaCO 3 /L.BA is estimated by multiplying the PA by 1.25 according to Anderson and Yang [56]. To simplify the data analysis, only BA was analyzed in current study.

Statistical method
The state of the reactor was defined as steady state when the daily biogas production was within 10% variation, for at least 6 consecutive days [57]. The date of process failure was determined by the time point where a significant decrease (RSD > 20%) appeared in the CH 4 yield. Relative standard deviations (RSDs) were calculated according to Eq. 2, for a quantitative assessment of fluctuations in daily indicator values compared with the previous day.
In the above equation, S is the standard deviation of the measured indicator value compared to the previous day, and x is the average of the values at day i and i−1. The larger the RSD value, the greater the fluctuation. In the current study, RSD > 10% (yellow symbol in all figures) and RSD > 20% (red symbol in all figures) were identified as the signs of slightly and highly unstable process, respectively. Thereby, the date when RSD exceeds 20% was determined as the time point where sudden changes took place.
For the quantification of goodness of fit between simulations and experimental data, relative root meansquared error (rRMSE) and mean absolute percentage error (MAPE) were used, according to Eqs. 3 and 4. where y sim i and y exp i , respectively, represent a single simulated or measured (experimental) data value, while y exp is the average of all experimental data values and n is the total number of experimental data points available. RMSE is commonly used for AD model evaluation [58,59] and in certain cases, offers advantages over other measures, especially when sensitivity to large errors between experimental and simulated data points is required [60]. After dividing RMSE by the mean of the experimental data points involved, the resulting rRMSE can be compared through different variable datasets, without having to consider its specific units. While rRMSE is defined in the range of 0 (no error) to infinity (no fit) and can therefore take significantly different values for different variables, MAPE is expressed in percentages. By design, it takes values between 0 and 100% for simulated data points that are at most double as large as experimental ones, but is in general said to produce reliable simulations when MAPE is less than 50% [61]. Consequently, it offers a complementary metric for error assessment less distorted by extreme data points, and provides a holistic view of dynamic simulation regression [62].

Mathematical model simulations
For the simulation of the two experimental reactors and the different overload conditions, an advanced bioconversion model (BioModel) was used. The model was developed by Angelidaki et al. [63,64], and was later extended by Kovalovszki et al. [65] and Lovato et al. [66], considering various anaerobic co-digestion scenarios for validation. Compared to the extended model established earlier, however, two minor changes were made in the present model implementation. Of the two, the first involved the removal of ammonia as a microbial growth limiting substrate from the model's kinetic equations, while through the second, the acetic acid inhibition effect on acetolactic methanogens was replaced by a total VFA inhibition effect. Former change was considered reasonable, given that the single substrate in these experiments-corn stalk-is a negligible source of ammonia. Meanwhile, the argument for extending the range of VFA inhibition on methanogens lies in their sensitivity to undissociated organic acids [67], as well as the dilute reactor medium, which in general might limit the availability of physical shelter and granule-forming sites for methanogens [68,69]. The above described VFA inhibition effect was kinetically controlled similar to the original acetic acid inhibition effect, using inhibition constants with manually estimated values of 315 mg/L and (4) MAPE = 100 n n i=1 y exp i − y sim i y exp i , 365 mg/L for R1 and R2, respectively. By comparison, similar inhibition constants found in published literature take values on a wide range, from 10 s to 1000 s mg/L and considering either acetic and propionic acid as the main inhibitors, or all VFA collectively [64,[70][71][72][73]. The values reported depend largely on the specific reactor conditions and substrates, therefore the low values used in this work can be justified by the highly dilute reactor contents and the lack of solid buffers (potentially leading to faster acidification) and the general sensitivity of methanogenic archaea to the accumulation of acids. In addition, the slight difference in the magnitude of the values estimated appears negligible compared to the above literature sources, and can be attributed to reactor-specific conditions.