This paper examines the importance of influent fractionation, kinetic, stoichiometric and mass transfer parameter uncertainties when modeling biogas production in wastewater treatment plants. The anaerobic digestion model no. 1 implemented in the plant-wide context provided by the benchmark simulation model no. 2 is used to quantify the generation of CH_{4}, H_{2} and CO_{2}. A comprehensive global sensitivity analysis based on (i) standardized regression coefficients (SRC) and (ii) Morris' screening's (MS's) elementary effects reveals the set of parameters that influence the biogas production uncertainty the most. This analysis is repeated for (i) different temperature regimes and (ii) different solids retention times (SRTs) in the anaerobic digester. Results show that both SRC and MS are good measures of sensitivity unless the anaerobic digester is operating at low SRT and mesophilic conditions. In the latter situation, and due to the intrinsic nonlinearities of the system, SRC fails in decomposing the variance of the model predictions (*R*^{2} < 0.7) making MS a more reliable method. At high SRT, influent fractionations are the most influential parameters for predictions of CH_{4} and CO_{2} emissions. Nevertheless, when the anaerobic digester volume is decreased (for the same load), the role of acetate degraders gains more importance under mesophilic conditions, while lipids and fatty acid metabolism is more influential under thermophilic conditions. The paper ends with a critical discussion of the results and their implications during model calibration and validation exercises.

## INTRODUCTION

The benchmark simulation models (BSMs) have been designed to develop, test and evaluate local and/or plant-wide wastewater treatment plant (WWTP) control and monitoring strategies using dynamic models to realize good effluent quality at low operational cost (Gernaey *et al.* 2014). The success of these tools is best illustrated by the more than 400 international papers, presentations and theses on work related to the benchmark systems which have been published to date (Gernaey *et al.* 2014). Aside from control/monitoring strategies, the BSM has also been proven to be an effective tool in teaching and studying WWTP processes. However, despite the broad spectrum of control/monitoring strategies that the BSM tools allow, most of the studies have been focused on the water line (activated sludge section) putting sludge treatment somewhat aside.

Anaerobic digestion (AD) is a key process for sludge treatment and its operation has significant impact on the overall (good) plant performance. Since biogas is the final product of the AD process, biogas production could be used to show how well the digester is performing. High biogas production along with a high methane content generally indicates that organic matter from the wastewater is broken down adequately and therefore the amount of sludge to be disposed of is reduced. Although AD has been around for a long time, poor performance and sometimes digester failure are still frequent in plants worldwide (Chen *et al.* 2008). If biogas production unexpectedly drops, it usually indicates inhibition of one or several of the anaerobic degradation phases. The hypothesis behind this paper is that performing sensitivity analysis on the AD model when simulating scenarios of normal operation and collapse may assist in identifying factors contributing to process inhibitions and even failure.

Thus, the objective of the study is to apply and compare two different global sensitivity analysis (GSA) methods to estimate how the predicted biogas production (CH_{4}, H_{2} and CO_{2}) in the anaerobic digester is affected by the influent waste composition and model parameters of the AD model No. 1 (ADM1) (Batstone *et al.* 2002). The results of this study will: (i) improve understanding of the model and its behavior; (ii) help gain additional knowledge about the model parameter(s) which influence the simulation results of the biogas production the most at different operating conditions; (iii) provide process engineers with useful information about the parameters that should be focused on during calibration exercises; and (iv) help in the design of model-based monitoring or control strategies during anaerobic digester operation in order to avoid process failure.

## METHODS

### Wastewater treatment plant and simulation program

The WWTP under study has the same layout as the International Water Association (IWA) BSM2 proposed by Jeppsson *et al.* (2007) and Nopens *et al.* (2010) (see Figure 1). In this layout, the influent wastewater passes through the primary clarifier for primary settling, the overflow of which is biologically treated in the activated sludge tanks with a modified Ludzack–Ettinger configuration consisting of two anoxic tanks followed by three aerobic tanks. The biologically treated water is then finally passed through the secondary clarifier, whose overflow is the WWTP effluent. The solids from the primary clarifier together with the thickened sludge from the secondary clarifier are anaerobically digested to produce gas and stabilized sludge. The liquids collected from the sludge thickener and the sludge dewatering unit are recycled to the primary clarifier.

A dynamic influent file is generated based on Gernaey *et al*. (2011) to describe the influent wastewater characteristics. Different potential control handles can also be included such as valves and flow controllers for aeration, carbon source, internal recirculation, return sludge, etc. However, the open-loop version (without control) of the model was used for the simulations in this study. Further information about the models used to simulate BSM2 can be found in Jeppsson *et al.* (2007).

The model focused on in this study is the ADM1. This model categorizes the considered biochemical processes into five steps of disintegration, hydrolysis, acidogenesis, acetogenesis and methanogenesis. In this study, the effects of the stoichiometric, kinetic, mass transfer and influent fractionation-related parameters of the ADM1 on the simulated production of digester biogas are determined by performing a global sensitivity analysis.

#### Scenarios simulated

The AD process is simulated within the BSM2 plant-wide model in four different scenarios with different solids retention times (SRTs) and temperature regimes detailed in Table 1.

Scenario | Liquid volume (m^{3}) | Gas volume (m^{3}) | Solids retention time (d) | Temperature (°C) | |
---|---|---|---|---|---|

A | 3,400 | 300 | 20 | Mesophilic | 35 |

B | 900 | 75 | 5 | 35 | |

C | 3,400 | 300 | 20 | Thermophilic | 55 |

D | 900 | 75 | 5 | 55 |

Scenario | Liquid volume (m^{3}) | Gas volume (m^{3}) | Solids retention time (d) | Temperature (°C) | |
---|---|---|---|---|---|

A | 3,400 | 300 | 20 | Mesophilic | 35 |

B | 900 | 75 | 5 | 35 | |

C | 3,400 | 300 | 20 | Thermophilic | 55 |

D | 900 | 75 | 5 | 55 |

#### Model output/evaluation criteria

The effects of varying the ADM1 parameters leading to variations in production of methane (CH_{4}), hydrogen gas (H_{2}) and carbon dioxide (CO_{2}) are evaluated. These gaseous components were chosen as model output criteria since these are the final outputs of the AD process and can provide useful information on how to analyze its performance for different scenarios.

### Uncertainty framing and global sensitivity analysis methods

A total of 57 parameters are considered in this study, consisting of kinetic, stoichiometric, mass transfer and influent fraction-related parameters. Default parameter values at 35 and 55 °C are listed in Table A1 in the Appendix (available online at http://www.iwaponline.com/wst/071/029.pdf). In the list of kinetic parameters, the ratio of half saturation constants to Monod specific uptake rates is considered in order to account for parameter correlations (Benedetti *et al.* 2010).

#### Uncertainty framing

To structure the different parameter groups, two uncertainty classes are distinguished: low and medium uncertainty. This parameter uncertainty quantification is performed through literature and expert review. Suggested parameter values for both mesophilic and thermophilic conditions of the ADM1, activated sludge model (ASM)/ADM/ASM interfaces and BSM2 parameters are taken from Batstone *et al.* (2002), Nopens *et al*. (2009) and Jeppsson *et al.* (2007). Expert knowledge and process understanding is used to determine the uncertainty ranges (Omlin *et al.* 2001). The first uncertainty class (C1), corresponding to low uncertainty of 10%, includes yield (kinetic) and mass transfer parameters (number of parameters = 11). The second uncertainty class (C2), corresponding to medium uncertainty of 25%, includes stoichiometric, kinetic and influent fractionation-related parameters (number of parameters = 46).

Uncertainties are assigned to the model parameters according to a uniform probability distribution. Uniform probability distributions are assumed for all the parameters due to lack of relevant information about the real distribution (Freni & Mannina 2010; Beven & Binley 1992).

The stoichiometric parameters are chosen accordingly, warranting that each set of yield coefficients on a specific substrate ([*f_h2_su*, *f_bu_su*, *f_pro_su*, *f_ac_su*], [*f_h2_aa*, *f_va_aa*, *f_bu_aa*, *f_pro_aa*, *f_ac_aa*], [*f_sI_xc*, *f_xI_xc*, *f_ch_xc*, *f_pr_xc*, *f_li_xc*]) (see Table A1) should have a sum of 1. This is done for each set by expressing the last parameter as the difference between 1 and the sum of the other parameters.

#### Propagation of uncertainty

Parameter uncertainty is covered by using the Latin hypercube sampling (LHS) and Morris sampling methods.

The LHS method (Iman & Conover 1982; Helton & Davis 2003) ensures that the random numbers are representative of the real parameters' variability. Each parameter is sampled, assuming a uniform distribution with the default value used as the mean and a corresponding lower and upper limit, depending on the uncertainty class assigned. Morris' sampling (Morris 1991), on the other hand, uses a one-at-a-time approach, which means that in each successive simulation of the model only one input parameter is given a new value until all input values are changed. It facilitates a global sensitivity analysis by repeating this method *r* times (a typical value for *r* is between 10 and 20), each time using a different set of input values.

In this case study, 1,000 samples are generated for each parameter creating a LHS matrix for the standardized regression coefficients (SRCs) method while a suggested set of samples (total of 870) for Morris' screening (MS) obtained.

#### Monte Carlo simulations

In the SRC method, Monte Carlo (MC) simulations are carried out, where each simulation has a set of parameters randomly selected by LHS. The uniform parameter distributions used for the Monte Carlo simulations were generated in a similar manner for each of the considered four scenarios (see Table 1). On the other hand, in the MS method, each simulation takes a set of parameters from the sampling matrix wherein only one factor is changed from one parameter set to the next.

The samples obtained from both LHS and Morris' samplings are introduced as parameters in the BSM2 open-loop version implemented in MATLAB/Simulink. The steady-state BSM2 was run for 200 days using the ode15s (Gear type) solver. Dynamic simulations were performed for 28 days using the ode45 (Runge–Kutta type) solver. This is mainly because the stiff solver ode15s does not handle well dynamic influents and step changes. Simulation results produced an evaluation matrix where the AD products (CO_{2}, H_{2}, CH_{4}) are considered.

#### Sensitivity analysis

Model inputs and parameters that significantly affect model outputs are determined through sensitivity analysis. There are a number of methods available for local/global sensitivity analysis (Saltelli *et al.* 2005). In this study, two (global) methods: (i) SRCs and (ii) MS were selected based on (i) implementation complexity and (ii) reasonable computational demand. Detailed descriptions of these two GSA methods can be found in the Appendix (available online at http://www.iwaponline.com/wst/071/029.pdf).

## RESULTS

Figures 2 and 3 show the results of the GSA for (i) the produced gaseous components CH_{4}, H_{2} and CO_{2}; (ii) each evaluated scenario (A, B, C & D). Specifically, Figure 2 shows the squared sum of the SRCs, highlighting (only) the significant parameters (SRC^{2} > 0.1). All other parameters having a squared sum of SRCs less than 0.1 are lumped together. Figure 3 depicts the mean, *μ* (abscissa), and the standard deviation, *σ* (ordinate), of the elementary effects. The *μ* is used to detect input factors with an important overall influence on the output while *σ* is used to detect factors involved in interactions with other factors or whose effects are nonlinear. The screening criterion is based on selecting those parameters that remain outside the wedge formed by lines corresponding to twice the standard error of the mean (SEM) (i.e. influential parameters).

In both Scenarios A (mesophilic, SRT = 20d) and C (thermophilic, SRT = 20d), *frxs_adm* (anaerobically degradable fraction of activated sludge biomass) is the dominating factor for the uncertainty in CH_{4} production. For CO_{2} production, besides *frxs_adm*, the ratio of lipids to carbohydrates in the remaining influent biodegradable particulate matter is essential and strongly contributes to model output uncertainty. For CO_{2} production, besides *frxs_adm*, the ratio of lipids to carbohydrates in the remaining influent biodegradable particulate matter, after the amount of protein has been calculated in relation to nitrogen availability, is essential (defined as fraction of lipids = *frlixs* and fraction of carbohydrates = 1 − *frlixs*), and strongly contributes to model output uncertainty. As this ratio is defined on a chemical oxygen demand (COD) basis the available biodegradable COD in the digester will not change and consequently the CH_{4} production is not affected. However, the carbon content of carbohydrates is assumed to be about 50% higher than that of lipids (expressed as kmol C·(kg COD)^{−1}). Therefore a low value of *frlixs* means that more carbon is available per COD unit, and the carbon that is not part of the produced CH_{4} must end up in the pool of inorganic carbon (*S*_{IC}) and consequently the quantity of stripped CO_{2} will increase. A higher value of *frlixs* will naturally have an adverse effect on CO_{2} production. In the case of H_{2}, the kinetic parameters and stoichiometric yield coefficients of hydrogen degraders, such as uptake rate (*k_m_h2*) and affinity constant (*K_s_h2*), are the most important parameters. This is mainly because these two parameters are the dominating ones for hydrogenotrophic methanogenesis.

In Scenario B (mesophilic, SRT = 5d), the kinetics of acetate degraders become the most important parameters influencing biogas production (CH_{4}, H_{2} and CO_{2}). This is because the uptake rate of acetate degraders (*k_m_ac*) mainly determines aceticlastic methanogenesis. Hence, when the digestion process does not proceed sufficiently fast, particularly for low SRTs, there is: (i) an accumulation of acetate; (ii) a decrease of the pH, which affects the uptake of acetate and hydrogen, as well as the carbonate equilibria; and (iii) finally a collapse of the whole digester process (see correlation among *k_m_ac* and CH_{4} in Figure 4B). This is in agreement with common knowledge that the value of the biodegradable fraction of the activated sludge biomass is significant for biogas production. A high biodegradable fraction generally results in higher biogas production while the opposite normally leads to process instability especially in scenarios where SRT is low and the digester operates at mesophilic conditions.

In Scenario D (thermophilic, SRT = 5d), the number of significant parameters on the predicted CH_{4} and CO_{2} uncertainty is increased. Influent fractionation (*frxs_adm*) is still the main factor responsible for CH_{4} and CO_{2} uncertainty. Nevertheless, the role of fatty acid metabolism (*k_m_fa*), the anaerobic digester liquid volume (*V_liq*), which is closely related to SRT, and the fraction of lipids in the total biodegradable COD (*frlixs*) also gain significance. Again, this is mainly because those parameters affect acetogenesis from long chain fatty acids (LCFAs) and the system's total pool of inorganic carbon (*S*_{IC}). No substantial variation can be observed in influential parameters with respect to H_{2} production, wherein uptake rate (*k_m_h2*) and affinity constant (*K_s_h2*) are still the most influential parameters and only the yield of biomass on hydrogen (*Y_h2*) is starting to gain importance.

## DISCUSSION

The phenomena observed in this study highlight the explicit nonlinearities of the digestion process. Using CH_{4} production as an example, Figure 4B shows that at low *k_m_ac* values the anaerobic digester collapses and there is no more biogas production. In addition, Scenario B generates *R*^{2} values <0.5. Under these circumstances, methods based on surrogate linear models fail and the variance decomposition cannot be accomplished for the three compounds. In that case, an alternative method is to be used, for instance one based on elementary effects. It should be highlighted though, that for this particular case both SRC and MS (Figures 2 and 3) identify the same parameters. However, the prediction of the multivariate model created using SRC would be unreliable (Sin 2012).

Another interesting point that Figure 4 shows is the effect of single/multiple parameters. For example, in Scenarios A and C, the effect of a single parameter is clear. On the other hand, the wider dispersion of B and D demonstrates the effect of multiple parameters. The high values of *R*^{2}, 0.92–0.99 for Scenarios A, C and D, means that the gas production per day can be assumed as linear functions of the model inputs under the defined conditions.

It is also of interest to note the studies from Cosenza *et al.* (2013) and Ruano *et al.* (2011) mentioning the non-stability of results of MS depending on the chosen repetition number *r*. To investigate this, MS simulations were performed with a larger sample size and it was verified that the results of the MS for this study converge, i.e. the same significant parameters were identified in both simulations using a different number of repetitions (870 and 1,740).

Another observation is that MS is able to determine more significant parameters than the SRC method, allowing identification of important factors using fewer samples than what is typically used in the SRC method. According to Cosenza *et al.* (2013), the Morris method tends to overestimate the number of important factors and the same conclusion was reached in this study when comparing the Morris method with the SRC method, wherein more parameters were identified as significant by MS. A drawback of the MS method is that although it is able to rank parameters it cannot quantify the rankings. Parameters involved in interactions are detected using MS; however, these interactions are not quantifiable. In addition, the method cannot distinguish between nonlinear effects and interaction effects. Other methods, such as the extended Fourier Amplitude Sensitivity Test (e-FAST) method (Saltelli *et al*. 1999), may be used for such purposes.

Another important conclusion is that the results only reflect the way the model is structured. In the way the model is coded, there is no clear distinction between primary and secondary sludge, and this explains why *frxs_adm* is the most significant parameter because this is the only parameter representing the biodegradability of the anaerobic digester influent (both primary and secondary sludge).

Valuable information gained from this study is useful in a practical sense: (1) when designing the operation of the AD unit; (2) when developing model-based control strategies of the wastewater treatment plant; and/or (3) for calibration purposes.

Having realized that the AD process is susceptible to failure when operating at short SRTs and under mesophilic conditions, either this combination of circumstances should be avoided or suitable control strategies should be implemented to stabilize anaerobic digester operation under such conditions. By analyzing the parameters that highly influence biogas production, which is closely related to the AD performance, control strategies can be developed in order to avoid digester failure. For instance, a digester could have a controller that regulates the addition of actively digesting sludge (Ward *et al*. 2008). In this way, acetate degraders, whose role is critical at low SRT and mesophilic conditions, would be given an advantage and acetate accumulation avoided. Other control strategies that can be implemented are feed-rate control (von Sachs *et al*. 2003) to circumvent reactor overloading and pH control (Strömberg *et al.* 2013; Steyer *et al.* 1999) to regulate biogas production. Control systems that could be implemented for the above purposes can be as simple as on-off controllers or could be much more sophisticated, such as artificial neural networks (Ward *et al*. 2008).

Another practical application of the results of this study is related to model calibration. Knowing which parameters influence an AD plant's biogas production the most, process engineers are guided in the proper direction, i.e. which parameters they should invest more time/effort/money in determining. Process modelers calibrating AD models of systems with high SRTs should focus on determining influent fractionation (i.e. biodegradability and composition). On the other hand, the kinetic characteristics of acetate degraders should be focused on when dealing with systems having low SRTs, and particularly those operating under mesophilic conditions.

In addition to the aforementioned, it should be pointed out that the results of the uncertainty analysis presented in this paper are only as good as the underlying assumptions. The selection of the simulation model (Batstone *et al*. 2002; Gernaey *et al*. 2014) is crucial and will determine the values of the biogas production. In the same way, the definition of the problem, the scope of the study and the definition of the input/output (uncertainty factors) strongly affect the results and conclusions of the analysis. The reader is reminded that the description of important factors such as the characteristics of input ranges, the assigned probability distribution functions (uniform) and the sampling methodology (correlated/non-correlated), are kept rather general in this case study and therefore imply a certain degree of subjectivity. A modification in any of these factors could lead to a somewhat different interpretation of the results. The results of the analysis have to be interpreted in the context in which the problem has been formulated (Flores-Alsina *et al*. 2008).

## CONCLUSIONS

Several parameters of ADM1, with reference to biogas production, vary at different operating conditions: different temperature regimes and SRT. To determine these significant parameters, global sensitivity analysis was performed on the model using two methods: (i) SRCs and (ii) MS. Results of both methods are in agreement with each other.

In AD units operating with long SRTs, influent fractionation is determined to be the most influential parameter with respect to CH

_{4}and CO_{2}production, while for shorter SRTs at mesophilic conditions the role of acetate degraders becomes important. For H_{2}production, the role of hydrogen degraders is most influential but shifts to acetate degraders for mesophilic conditions with short SRTs.Knowledge of these highly influential parameters is useful for understanding the AD process including situations of process failure and for model calibration and validation exercises.

## ACKNOWLEDGEMENTS

Ms Solon and Dr Flores-Alsina acknowledge the People Program (Marie Curie Actions) of the European Union's Seventh Framework Programme FP7/2007-2013 under REA agreement 289193 (SANITAS) and REA agreement 329349 (PROTEUS), respectively. This article reflects only the authors' views and the European Union is not liable for any use that may be made of the information contained therein.