Uncertainty and sensitivity analysis for reducing greenhouse gas emissions from wastewater treatment plants

The interest towards greenhouse gas (GHG) emissions from wastewater treatment plants (WWTPs) has considerably increased during the last decade (Kampschreur et al. 2009; Corominas et al. 2012; Mannina et al. 2016a). WWTPs can be a source of GHGs emissions as: direct (due to the biological processes), indirect internal (due to electric or thermal energy consumption) and indirect external (due to sludge disposal or chemicals production) emission (IPCC 2013). The waste and wastewater sector accounts for about 3% of the global GHG emissions (IPCC 2013). In order to reduce GHG emissions from WWTPs, researchers and professionals have three main aims: (i) mitigate pollutants emissions (both liquid and gaseous) (including GHGs); (ii) maintain the required quality of the treated wastewater; (iii) limit as much as possible the operational cost (Flores-Alsina et al. 2014; Mannina et al. 2016a).

In view of achieving these aims, several efforts have been provided in literature for establishing mathematical tools able to predict WWTP behaviour (in terms of GHG and liquid pollutants) (Kyung et al. 2015; Mannina et al. 2016a, 2016b, 2016c). Several modelling methods have been proposed to include GHG assessment (e.g. empirical models, process-based models often related to mass balance, dynamic mechanistic models often associated with life cycle assessment at plant wide-scale) (among others, Flores-Alsina et al. 2014; Bisinella de Faria et al. 2016; Mannina et al. 2017, 2018a). For each of the aforementioned modelling methods, plant-wide modelling approach can offer a straightway and effective solution for assisting in developing strategies aimed at reducing GHG emissions and improving environmental protection. Plant-wide mathematical models adopting dynamic mechanistic approaches are characterized by providing accurate predictions; however, they are more time demanding than simple comprehensive process-based models (Corominas et al. 2012). Therefore, despite the advantages of dynamic mechanistic models in terms of accurateness of the predictions, for a rapid GHG estimation, the simple comprehensive process-based models are suggested (Mannina et al. 2016a). However, it has to be mentioned that simple comprehensive process-based models are often based on a great number of assumptions and their response can be highly uncertain. Therefore, uncertainty analysis may help in obtaining model confidence and improve the model predictions. However, sensitivity and uncertainty analysis have rarely been performed in the GHG WWTP modelling studies with the aim to identify the key source of uncertainty (Behera et al. 2020). The existing studies mainly consider only the water line treatment (such as Mannina et al. 2017). Sweetapple et al. (2013) have performed an uncertainty analysis of a GHG plant-wide model, however they adopted a dynamic mechanistic modelling approach. Further, Mannina et al. (2016b, 2016c) have performed an uncertainty analysis of a GHG plant-wide comprehensive process-based model but without considering all the possible N₂O emissions pathways.

Bearing in mind the above introduction, in this study, sensitivity and uncertainty analysis of a new simple processes-based model have been performed in view of identifying the key sources of uncertainty.

The model adopted here is based on chemical oxygen demand (COD), total suspended solids (TSS) mass-balance and nitrogen kinetic/mass-balance. The model consists of 29 model factors (divided into kinetic, stoichiometric, influent fractionation and operational factors) and 21 state variables. Tables 1 and 2 summarize the detailed description of model factors and state variables, respectively. The model allows to assess the total equivalent CO₂ (CO_2eq) emissions (kgCO_2,eq/day or kgCO_2,eq/treated volume) from a WWTP as the sum between direct (DIR,CO_2eq) and indirect (IND,CO_2eq) emissions. The model considers uniform and constant influent features over time. Further, non-biodegradable compounds are considered conservative.

The model adopted here introduces the following main innovative aspects related to up-to-date available literature plant-wide models: (i) kinetic/mass-balance based model regarding nitrogen; (ii) includes the nitrification as a two-step process; (iii) includes the N₂O formation during nitrification and denitrification both in dissolved and off-gas forms. More specifically, the autotrophic biomass is divided into autotrophic ammonia oxidizing bacteria (AOB) and nitrite oxidizing biomass (NOB) in view of modelling the N₂O formation processes during nitrification. The secondary effluent ammonia and nitrite concentration are calculated according to the mass balance analysis of a well-mixed activated sludge reactor (Metcalf & Eddy 2003). Thus, the effluent ammonia concentration depends on the AOB kinetics parameters and on the sludge retention time. The dissolved N₂O concentration inside the aerobic reactor is modelled according to the relationship proposed by Wu et al. (2014), while the dissolved and gaseous N₂O concentration inside and from the anoxic reactor is modelled according to the relationships proposed by and Yan et al. (2017). According to Yan et al. (2017) the gaseous N₂O emitted from the anoxic reactor depends on the carbon to nitrogen ratio. The model allows to assess DIR,CO_2eq and IND,CO_2eq emissions at a plant-wide scale, considering the contribution due to the water line (headworks (HD), primary settler (PS), activated sludge process (ASP), secondary settler (SS), treated effluent discharge (EFF)) and to the sludge line (anaerobic digestion (AD), dewatering (D), biosolids disposal (TB) and energy recovery (ER) due to the biogas combustion) (Figure 1).

In this study, the value of the specific CO₂ equivalent emission for energy consumption (SEnCO_2,eq) suggested by EIA (2009) has been adopted; the sum in brackets of Equation (3) represents the total energy demand of the plant (e_D).

The model has been applied to a conventional activated sludge (CAS) WWTP having the Ludzack-Ettinger (anoxic and aerobic biological reactors) configuration for nitrogen removal (Figure 1). The plant treats 60,000 m³ d⁻¹ of real wastewater and consists of a water line (headworks, primary settler, CAS units according to Ludzack-Ettinger configuration, secondary settler and disinfection unit) and a sludge line (anaerobic sludge digester with energy recovery, sludge dewatering unit) (Mannina et al. 2016b, 2016c) (Figure 1). A detailed description of the model can be found in Mannina et al. (2019, 2020).

The difference between S_Ti and S_i represents the interaction among the model factors. The Extended-FAST method requires nxNR simulations, where n is the number of factors and NR is the number of runs per model factor (NR = 500–1,000 according to Saltelli et al. 2004).

The uncertainty analysis has been performed by running Monte Carlo simulations varying the model factors selected as important through the GSA. More specifically, three uncertainty scenarios have been investigated by Monte Carlo simulations changing the following factors: scenario 1 – all model important factors; scenario 2 – model factors related to the influent features; scenario 3 – model factors related to the operational conditions.

For each scenario, results of the Monte Carlo simulations have been interpreted by analysing the cumulative distribution function (CDF). Specifically, for each model output, the CDF of the normalized value has been considered; the normalized value has been obtained by dividing each simulated value by the maximum one. Moreover, the comparison of the uncertainty analysis results among the model outputs taken into account for each scenario has been performed by comparing the value of the relative uncertainty bandwidth. This latter has been computed by dividing the width between the 5th and 95th percentiles for the 50th percentile.

TN_OUT includes non-oxidized ammonia, nitrite, nitrate and dissolved nitrous oxide while organic nitrogen is neglected. β_COD and β_TN have been considered equal to 1 and 20, respectively, according to literature (Maere et al. 2011).

For a detailed description of the symbol and variation range of each factor the reader is referred to Table 1. Due to the lack of knowledge about the distribution of the model factors, a uniform prior distribution was considered for each factor. The Extended-FAST method was applied using the sensitivity package developed by Pujol (2007) in the R environment (Development Core Team, 2007). For the Extended-FAST application, 1,000 NR have been considered, consequently 29,000 simulations have been run.

To classify important, non-influential and interacting factors, thresholds of the sensitivity measures were selected. Specifically, factors with S_i value greater than 0.01, at least for one model output, were classified as important. Interacting model factors were selected using the normalised index value (S_Ni), which corresponds to the ratio between the interaction of the i-th model factor related to one model output and the maximum value among the interactions for that model output. Factors with S_Ni greater than 0.5 for at least one model output were considered to be interacting. Model factors with S_Ni and S_i values lower than 0.5 and 0.01, respectively, were considered to be non-influential. The uncertainty analysis was performed by all the model factors classified as important or interacting. For each analysed scenario, 1,000 Monte Carlo simulations (by adopting a Latin Hypercube Sampling) have been performed.

In Figure 2, the Extended-FAST results for each analysed model output are reported. For sake of shortness, the values of S_i, S_Ti, S_Ti-S_i and S_Ni are summarized in Table 3.

The sum of S_i explains more than 94% of the total variance for all model outputs suggesting that the model is highly linear and additive. This statement is also confirmed by the value of the sum of S_Ti, which is always close to 1. This latter result suggests that a very low interaction among factors takes place. By applying the Extended-FAST method, 17 model factors resulted to be important (S_i > 0.01 and/or S_Ni > 0.5) for at least one model output.

By analyzing data reported in Figure 2, one can observe that six factors have significant impact on DIR,CO_2,eq. Specifically, npbCOD_IN, pCOD/VSS, nbsCOD_IN, pbCOD_IN, pbCOD_PI and npbCOD_PI (S_i equal to 0.78, 0.064, 0.051, 0.081, 0.011 and 0.011, respectively) have an S_i value higher than 0.01 for DIR,CO_2,eq. Among these factors, three (npbCOD_IN, nbsCOD_IN, and pbCOD_IN) are related to the influent wastewater fractionation, two (pbCOD_PI and npbCOD_PI) are related to the primary effluent fractionation and one (pCOD/VSS) refers to particulate COD and volatile suspended solids.

Influent fractionation factors are also strongly interacting; indeed, the S_Ni value for npbCOD_IN and pbCOD_IN is equal to 1 and 0.67, respectively (Figure 2, Table 3). The influent fractionation factors influence the bCOD availability to heterotrophic biomass growth and, consequently, the direct CO₂ produced during the biomass respiration value.

For example, the higher the nbsCOD_IN fraction and the lower the availability of substrate to be degraded during the biomass metabolism (both during aerobic and anoxic conditions). Hence, CO₂ produced during the biomass respiration is reduced as a result of the conservative nature of nbsCOD_IN. Further, even the physical processes occurring inside the primary settler (as demonstrated by the importance of pbCOD_PI and npbCOD_PI) can influence the COD availability for the biomass respiration and consequently the DIR,CO_2,eq.

The indirect emissions (IND,CO_2,eq) are mostly influenced by three main factors: the oxygen transfer efficiency (OTE) (S_i equal to 0.679), fractionation factors (nbsCOD_IN, pbCOD_IN and npbCOD_IN with S_i equal to 0.024, 0.132 and 0.167, respectively) and the yield coefficient for heterotrophic biomass growth (Y_H) (S_i equal to 0.015 (Figure 2, Table 3). Such results are consistent with literature which shows that energy consumption is mostly due to the activate sludge process (ASP): More specifically, ASP consumes are in the range of 60–70% (Mamais et al. 2015; Mannina et al. 2018b; Massara et al. 2018).

OTE is an important factor for the indirect emissions since regulates the air flow rate required to maintain the aerobic conditions inside the aerobic reactor. Similarly, nbsCOD_IN, pbCOD_IN and npbCOD_IN are important factors for indirect emissions since they regulate the availability of soluble COD required for the biological processes. For example, as the fraction of sCOD decreases the oxygen required for the aerobic processes decreases, thus influencing the power requirements of the aeration process and of the entire WWTP. Finally, the yield coefficient (Y_H) influences the amount of the air flow to be inserted in the aerobic reactor and consequently the power requirements.

The COD_OUT is mostly influenced by the maximum growth rate of heterotrophic biomass (μ), the half-saturation factor for heterotrophic biomass (k_s, S_i equal to 0.067) and the decay rate of heterotrophic biomass (k_d, S_i equal to 0.532). Among these factors μ and k_d also resulted to be interacting having an S_Ni value of 1 and 0.98, respectively. All model factors resulted to be important for COD_OUT have a key role on regulating the heterotrophic active biomass inside the aerobic reactor, thus influencing the COD consumption and consequently the COD concentration in the effluent wastewater.

Eleven model factors mostly related to the influent fractionation, heterotrophic, AOB and NOB resulted to be important for the total nitrogen at the effluent (TN_OUT) (Figure 2, Table 3). Among these factors, it is important to observe that Y_H (S_i equal to 0.289)_, which is directly related to the heterotrophic biomass, is the most important factor for TN_OUT (with an S_i value equal to 0.29). This result is related to the anoxic activity of heterotrophic biomass which strongly influences the amount of nitrogen removed from the system and consequently the TN_OUT. Conversely, the model interacting factor for TN_OUT is the sludge retention time of the conventional activated process (SRT_ASP) (S_i and S_Ni value equal to 0.01 and 1, respectively). Indeed, SRT_ASP is a key factor for different processes occurring inside the system. Among these processes, SRT_ASP also regulates the autotrophic biomass growth within the aerobic reactor (nitrification). The nitrification process is responsible for the nitrate availability to be denitrified and consequently removed from the system.

Finally, the EQI is influenced by 13 model factors. Among these factors, five are related to the influent COD fractionation (nbsCOD_IN, pbCOD_IN, npbCOD_IN, pbCOD_PI and npbCOD_PI), five to the biomass activity (k_d, Y_H, k_dN,AOB, Y_N,AOB and μ_N,NOB), one to the nitrogen content in the secondary sludge (iNVSS_SS) and two to the operational conditions (SRT_DIG, SRT_ASP) (Figure 2, Table 3). Among these factors, Y_H has the highest S_i value (0.28), indicating that the role of heterotrophic bacteria is relevant in terms of pollutants discharge (both nitrogen and COD) (Figure 2, Table 3), while the influence of SRT_ASP is only due to the interaction contribution. Indeed, despite the value of S_i for SRT_ASP (related to EQI) being lower than 0.01, this factor proved to be important due its high interaction. Indeed, the S_TI-S_i and S_Ni value resulted to be equal to 0.03 and 1, respectively. The high interaction of SRT_ASP is mainly due to the complexity of the model in terms of biological processes which are strongly regulated by SRT_ASP.

The uncertainty analysis was performed by considering the three scenarios described earlier. For each scenario, Monte Carlo simulations have been performed by varying important model factors grouped according to the scenarios.

In Figure 3, the cumulative distribution functions (CDFs) of DIR,CO_2,eq, IND,CO_2,eq, EQI, COD_OUT and TN_OUT of each scenario is reported.

In Figure 3, x-axes report the normalized value of each model output. The normalized value has been obtained by dividing each Monte Carlo output model value by the maximum one. From a visual inspection of Figure 3, one may observe that the width of the uncertainty bands, (calculated as difference between the 95th and 5th percentiles of the normalized values), changes with the model output and the analysed scenarios. This is mainly due to the fact that some of the model outputs entail a different level of complexity in terms of involved phenomena.

For scenario 1, where all important model factors have been varied, the uncertainty bands' widths of EQI (0.51), TN_OUT (0.51) and IND,CO_2,eq (0.38) is higher than that of COD_OUT (0.16) and DIR,CO_2,eq (0.32) (Figure 3(a)). This result is due to the fact that a greater number of processes influences EQI, TN and IND,CO_2,eq than the other model output. Therefore, the uncertainty of model factors varied during Monte Carlo simulations (operational, influent fractionation and kinetic/stoichimetric factors) strongly influence the EQI, TN and IND,CO_2,eq predictions. Conversely for scenario 2, the uncertainty band widths of TN_OUT (0.36), EQI (0.35), and COD_OUT (0.33) are higher than that of DIR,CO_2,eq (0.32) and IND,CO_2,eq (0.28) (Figure 3(b)). This result is mainly due to the fact that influent fractionation factors, as reported above, strongly influence TN_OUT, EQI and COD_OUT. Finally, for scenario 3, where only the influential factors related to the operational conditions were varied, the uncertainty bands' widths of IND,CO_2,eq (0.36) and DIR,CO_2,eq (0.24) are higher than that of EQI (0.18), TN_OUT (0.18), and COD_OUT (0.09) (Figure 3(c)). This result has relevant interest since underlines that the operational conditions may have an important role influencing both direct and indirect emissions predictions.

The results of Figure 3 show that the uncertainty band width decrease from scenario 1 to scenario 3, underlying that the reduction of the number of model factors varied from scenario 1 to scenario 3 reduce the uncertainty of model predictions. However, it is interesting to observe that from scenario 1 to scenario 2 the uncertainty bands width related to COD_OUT increases. This result is mainly due to the fact that the uncertainty of the influent fractionation factors, varied during in scenario 2, strongly influence the uncertainty of COD_OUT. Further, from scenario 2 to scenario 3, the uncertainty bands' width related to IND,CO_2,eq increases. Such a fact underlines the key role of operational conditions in influencing the IND,CO_2,eq predictions.

In order to provide a quantitative assessment of the model uncertainty and to make comparable the results among the model outputs, the relative uncertainty band width for each model output and scenario has been computed by dividing the width between the 5th and 95th percentiles to the 50th percentile. In Figure 4, the relative uncertainty band widths for each model output and scenario are reported.

By analysing Figure 4, one may observe that for scenarios 1 and 2, the highest uncertainty is related to TN due to the treated effluent (the relative uncertainty bandwidth is equal to 0.79 for scenario 1 and 0.48 for scenario 2), while for scenario 3 the highest uncertainty is related to IND,CO_2eq (0.47), thus, corroborating again that operational conditions mostly influence indirect emissions (scenario 3). The highest uncertainty of TN for scenarios 1 and 2 is mainly due to the uncertainty of kinetic, stoichiometric and fractionation factors which influence the biological processes and consequently the amount of nitrogen removed.

Data in Figure 4 also show a reduction of the total relative uncertainty from scenario 1 to scenario 3, therefore the uncertainty of model prediction decreases after fixing stoichimetric/kinetic model factors. This result suggest that an accurate estimation of stoichimetric/kinetic model factors has to be performed before applying the model.

The key findings of this study are summarised as follows:

• The sensitivity analysis reveals that model factors related to the influent wastewater and primary effluent COD fractionation exhibit a significant impact on direct, indirect and EQI model factors.

• The effluent concentration of COD and total nitrogen is mainly influenced by the heterotrophic stoichiometric/kinetic factors; revealing that for TN the role of denitrification (anoxic growth of heterotrophic bacteria) is relevant.

• The uncertainty analysis reveals that TN_OUT has the highest uncertainty in terms of relative uncertainty band for scenario 1 and scenario 2, thus revealing that uncertainty of all influential model factors and of influent fractionation factors has a relevant role on the total nitrogen prediction.

• Results of the uncertainty analysis show that the uncertainty of model prediction decreases after fixing stoichimetric/kinetic model factors.