## Abstract

Modelling conversion processes in sewers can help minimize odour and pipe corrosion issues, but model uncertainties and errors must be understood. In this study, the Wastewater Aerobic/Anaerobic Transformation in Sewers (WATS) model is implemented in two different frameworks; 1-D (CSTR-in-series) and computational fluid dynamics (CFD) to study the uncertainties due to model parameters and its mathematical form. The 1-D model is used to conduct uncertainty/sensitivity analysis using Monte Carlo simulations. Time-averaged outputs were represented using a general linearized model to quantify the importance of specific parameters. The sulfide formation rate per unit area of the biofilm is the most influential parameter. Parameters controlling anaerobic hydrolysis and fermentation are also significant. Uncertainty due to model structure is studied using CFD to explore the influences of non-homogeneous surface reactions and solids settling. These showed that the 1-D model provides a reasonable characterisation of the process for simple flows in pressure mains.

## HIGHLIGHTS

Monte Carlo simulations reveal the most important parameters in sewer model.

Sulfide formation, anaerobic hydrolysis, fermentation parameters are most significant.

Fluid dynamics simulations show biofilm reaction homogenization is appropriate.

Solids settling can be neglected in pressure mains with continuous flow.

### Graphical Abstract

## INTRODUCTION

There are a number of potential problems that can arise from the biological and chemical reactions that occur in sewer systems, most commonly odour nuisance and corrosion (Carrera *et al.* 2016). The build-up of hydrogen sulfide and volatile organic compounds is the major cause of odour, while formation of hydrogen sulfide and its subsequent oxidation to sulfuric acid on the moist walls of the sewer pipes leads to pipe corrosion. The cost of replacement and maintenance due to corrosion is significant for municipalities (Rootsey & Yuan 2011). Production of methane by methanogens in the sewer can cause explosion risk in confined spaces. Furthermore, methane is a potent greenhouse gas that is responsible for almost 20% of forcing in climate change models (Minami & Takata 1997). Guisasola *et al.* (2009) reports that methane discharge from sewers contributes significantly to the overall greenhouse gas emissions that are associated with wastewater transport and treatment.

Developing a more complete understanding of the processes occurring in sewer networks has the potential to enable improvements in the performance of wastewater treatment plants (WWTPs). For example, removal of readily biodegradable chemical oxygen demand (COD) can be achieved in the sewer system through conversion to slowly biodegradable organic matter stored in biomass, which is more readily removed by mechanical treatment. In other words, the sewer is not only a collection system; it is also a pre-treatment stage that can be exploited in different ways to improve the performance of the WWTP (Guo *et al.* 2019). When considering such strategies, it is important to consider downstream processes where it may be necessary to preserve and produce readily biodegradable substrate, since this is beneficial for denitrification and biological phosphorus removal processes in advanced WWTP designs (Hvitved-Jacobsen *et al.* 2002). Furthermore, some microbial activities, such as methanogenesis, consume the available organic carbon which may be required for other downstream processes (Gutierrez *et al.* 2009). In light of these complex and coupled interactions, process modelling becomes increasingly important.

Moving from an empirical approach (Pomeroy & Parkhurst 1978) to a conceptual understanding of the biological and chemical mechanisms of the in-sewer transformations allows for greater insight into the processes that occur. Microbial transformation processes in biological systems are, to a large extent, identical to the in-sewer transformations (Bjerre *et al.* 1998). Therefore, the mathematical model of the microbial processes described in the IAWQ Activated Sludge Model No. 1 (ASM1) (Henze *et al.* 1987) was adopted in the work of Bjerre *et al.* (1998) to model the aerobic processes in the wastewater collection system with two additional types of hydrolysis processes. Hvitved-Jacobsen *et al.* (1998) reported that certain modifications were necessary to adapt the ASM1 model to gravity sewer systems. First, the biomass decay rate was found to be unrealistic, so the biomass decay process was replaced by the energy maintenance requirement rate as a reasonable alternative process. Furthermore, the biofilm concentration was assumed to be constant and any growth of the biofilm was assumed to be released in the bulk water phase. Hvitved-Jacobsen *et al.* (1998) also included anaerobic processes in their model, based on the IAWQ Activated Sludge Model No. 2 (Gujer *et al.* 1995) where, even if aerobic conditions are assumed (as in gravity sewers), anaerobic regions could certainly exist within the sewer system (Carrera *et al.* 2016). The sulfur cycle was also included as it was necessary to describe sulfide formation and consumption (Hvitved-Jacobsen *et al.* 1999). The work of Hvitved-Jacobsen *et al.* (1999) is considered to be the foundation of the Wastewater Aerobic/Anaerobic Transformation in Sewers (WATS) model, which describes the aerobic and anaerobic transformation of carbon in the sewer, in addition to sulfide formation by the sulfate-reduction process. Yongsiri *et al.* (2003) later improved upon the original WATS model by introducing the emission of the sulfide from the bulk water phase, which is considered as the first step toward incorporating physical processes alongside the biological and chemical transformations. The application and parameter estimation of the WATS model were reported by Tanaka & Hvitved-Jacobsen (2002). However, the uncertainty and sensitivity analysis of the WATS model is scarce in the literature.

Model-based design, operation, and development of control strategies for wastewater systems management is prone to risk of not meeting regulatory standards or operating a system inefficiently as a result of model uncertainties, therefore uncertainty and sensitivity analysis is essential (Sin *et al.* 2009). Uncertainty analysis is concerned with the propagation uncertainty from different sources onto the global model output, while sensitivity analysis is concerned with the weight of each of the model inputs on the model output. According to Saltelli *et al.* (2008), sensitivity analysis in wastewater applications can be categorized into three classes: local methods, global methods, and screening methods. This work is focused on global sensitivity analysis, which is capable of providing information regarding the effect of model parameters on the model output over the space of all possible parameter values. Many global sensitivity analysis techniques have been described in the literature, including regression-based methods, e.g. the Standard Regression Coefficient (SRC) method (Sin *et al.* 2011; Flores-Alsina *et al.* 2012), and variance-based methods, e.g. the Fourier Amplitude Sensitivity Testing (FAST) and Extended-FAST methods (Chen *et al.* 2012; Cosenza *et al.* 2014).

Uncertainty analysis framing defines the sources of the uncertainty that are to be considered in the study, including the model assumptions and approximations, as well as the identification of the uncertainty range assigned for each of the model inputs. Sin *et al.* (2009) conducted uncertainty analysis under different model framings, concluding that the uncertainty of the different framings is almost additive if combined together. The main sources of uncertainty on model outputs, as described in Sin *et al.* (2009), arise from the uncertainty of the model inputs, the form (structure) of the mathematical model, and stochastic events. In the literature, sewers are modelled by coupling biochemical models with a hydraulic model, where the sewer pipe is commonly approximated as a series of Continuous Stirred Tank Reactors (CSTRs) (Sharma *et al.* 2008). One downside of the CSTR model is the spatial distribution of the species and their associated reactions cannot be taken into account, representing a model structure uncertainty. This can be particularly problematic for biofilm reactions, e.g. sulfide production by sulfate-reduction bacteria (SRB) (Mohanakrishnan *et al.* 2009), which are highly localized at the biofilm accumulated on the sewer walls. Similarly, settling of particulate species can generate spatial non-uniformities that cannot easily be captured using existing modelling approaches. This study is intended to shed light on the process model uncertainties that may arise from two common simplifications that are assumed in such models: (i) the homogenization of surface reactions using the ratio of wall surface area to volume and (ii) the omission of particulate settling.

A computational fluid dynamics (CFD) approach for modelling sewer systems offers some potential improvements over a CSTR model, since it is able to resolve the spatial variations in species concentration and reaction rates along with the transport of the different species by advection and diffusion. CFD solvers are robustly capable of solving the hydrodynamics of the sewer system, while the biological and chemical reactions embodied by the various kinetic models discussed previously can be implemented as additional coupled advection–diffusion transport equations. Of course, CFD simulation of pressure mains is computationally expensive, but can be useful for examining the details of the spatial distribution and the heterogeneity of the reactions in certain cross-sections of the pipes. As such, the uncertainties that accompany the usage of the CSTR model structure can be examined using CFD and can ultimately be used to further improve the lumped-parameter-based models. For gravity sewers, CFD could also play an important role in determining the mass transfer mechanisms (Teuber 2020), thereby improving modelling of hydrogen sulfide emission.

The aim of this work was to determine the sources of uncertainty in the modelling of sewer systems using the WATS model. These sources are associated with the biochemical model parameters (input uncertainty) and with the mathematical form of the model (model structure). The input uncertainty related to the model parameters is examined using the SRC method. The structural uncertainties that will be considered are: (i) the representation of the sewer pipe as a series of tanks, (ii) the homogenization of biofilm surface reactions, and (iii) exclusion of the physical processes of solids settling. The structural uncertainties are assessed using a detailed CFD model which allows for each of these factors to be included and their effects examined.

## MATERIALS AND METHODS

### Measurement data

Measurement data provided in Guo *et al.* (2018) are used for the calibration of the models in the present study. A 24-hour on-site sampling and measurement campaign on a force main system in California, USA was reported. The current study is concerned with one of three force mains that transport sewage from three different catchments to a WWTP. The main under consideration is 9.09 km in length and 40.64 cm in diameter. The flow rates and the water characteristics of the sewage were measured at the inlet and the outlet of the pipe during the sampling period. Soluble COD, volatile fatty acid (VFA), and total sulfide measurements are the main characteristics that are used for the calibration and validation (see Supplementary information).

### Biochemical model

The WATS model is the basis of most of the recently developed models (e.g. Sharma *et al.* 2008) to model the anaerobic processes that include fermentation, hydrolysis of particulate COD, and sulfide production. Therefore, the WATS biochemical model reported by Yongsiri *et al.* (2003) and Nielsen *et al.* (2005) is adapted in this work to model the biological and chemical conversions in the sewer. These conversions include main processes in the carbon and the sulfur cycles. Modifications of the original model are made by adding the chemical and biological oxidation reactions of sulfide (Carrera *et al.* 2016). The kinetics of the biological and chemical oxidation of sulfide in the water phase and biofilm are adopted from Nielsen *et al.* (2003) and Nielsen *et al.* (2006). The diagram in Figure 1 describes the key processes in the implemented biochemical model showing the main state variables in the model. The resulting model is in the form of a system of ordinary differential equations (ODEs) that describe the production and consumption rates of the modelled species.

Integration of the biokinetics with the hydrodynamics can either be carried out using a conventional CSTR-in-series approach or by solving the complete advection-diffusion equation for each of the species using the CFD method. Both of these are considered in the present work, where the CSTR model is used for calibration and sensitivity/uncertainty analysis of the model to its input parameters, while the CFD model is used to assess uncertainty with respect to the influence of solids settling and non-homogeneous reactions in the biofilm on the walls.

### CSTR-in-series process model

The process model is developed in MATLAB/Simulink, where the sewer pipe is represented as a series of CSTRs. Sensitivity analysis for the number of the CSTRs was conducted as a preliminary step. It was found that including more than 80 CSTRs resulted in no noticeable change in the concentration profiles, which indicated that this is an adequate representation of the plug flow behaviour. Therefore, all subsequent results implement 80 CSTRs. As the sewage flow rate and characteristics largely vary along day-time and night-time, dynamic inlet conditions are necessary to be considered for better prediction of the dynamic and peak concentration of the species. Species concentrations at the outlet of the series of the CSTR are monitored for the calibration and sensitivity analysis of the model.

#### Calibration and uncertainty/sensitivity analysis

Studies on the uncertainty/sensitivity of the sewer biochemical models are scarce in the literature. Therefore, most of the initial values for the biochemical model parameters used here are adopted from Calabro *et al.* (2009) and Hvitved-Jacobsen *et al.* (2013) as an initial parameter state. Using 24-hour dynamic measurement data, the key parameters of the model are calibrated to obtain a reasonable match between the model prediction and the dynamic measurement data for total COD, soluble COD (sCOD), VFA, and total sulfide. For further quantification of the significance of the role of each of the key parameters in calibrating the model and potential source of model output variance, uncertainty/sensitivity analysis is conducted. Following the work of Sin *et al.* (2011), the analysis is conducted into two steps. First, Monte Carlo method is used to explore the propagation of uncertainty from the model input to the output. Then, the Monte Carlo results are used for analysis by graphical representation and by fitting to multivariate linear functions of the model input using the Standardized Regression Coefficient (SRC) method. This study is concerned with identifying the most influential parameters in the model. However, SRCs could be used to determine the non-influential parameters as well, if the coefficient of determination (R^{2})is higher than 0.7, as in the case of our study (to be shown later) (Cosenza *et al.* 2013). Further study could be conducted using one of the variance-based sensitivity analysis methods such as the Extended-FAST method (Cosenza *et al.* 2013) to explore the higher order effect of the model parameters.

Monte Carlo simulations are performed for the uncertainty/sensitivity analysis of the model parameters using 1,000 simulations in a randomly sampled parameter space. A pre-defined variational range for each parameter, except the parameters pertaining to the anoxic processes, as illustrated in Table 1, was determined. A guideline for specifying the variational range and default values are provided by Brun *et al.* (2002). Parameters are classified into three uncertainty classes (Table 1) with variations of 5, 25 and 50%, around the preliminary calibrated parameter values, where these values correspond to classes 1, 2, and 3, respectively. A parameter sampling matrix, *S*, is created using the Latin Hypercube Sampling (LHS) technique to ensure full coverage of the range for each parameter variation (Helton & Davis 2003; Sin *et al.* 2009). The calibrated parameter values are used as the mean values of the variation range and reference values.

. | Parameter . | Symbol . | Unit . | Calibrated value . | Uncertainty class . |
---|---|---|---|---|---|

Biomass growth in bulk water parameters | |||||

1 | Maximum aerobic specific growth rate | μ_{H,O2} | [d^{−1}] | 7 | 3 |

2 | Yield constant for heterotrophic biomass in water phase | Y_{hw} | [gCOD/gCOD] | 0.55 | 1 |

3 | Maintenance energy requirement rate constant | q_{m} | [d^{−1}] | 1 | 3 |

4 | Saturation constant for readily biodegradable substrate | K_{sw} | [gCOD/m^{3}] | 1 | 3 |

5 | Saturation constant for dissolved oxygen | K_{O} | [gO_{2}/m^{3}] | 0.05 | 3 |

Biomass growth in biofilm parameters | |||||

6 | Half-order aerobic reaction rate constant per unit area for biofilm surface | K_{1/2} | [gO_{2}^{0.5} m^{−0.5}d^{−1}] | 4 | 3 |

7 | Yield constant for heterotrophic biomass in biofilm | Y_{hf} | [gCOD/gCOD] | 0.55 | 1 |

8 | Saturation constant for readily biodegradable substrate in biofilm | K_{sf} | [gCOD/m^{3}] | 5 | 3 |

9 | Saturation constant for dissolved oxygen | K_{O} | [gO_{2}/m^{3}] | 0.05 | 3 |

Particulate hydrolysis parameters | |||||

10 | Rapidly hydrolysis rate constant | K_{h1} | [d^{−1}] | 12 | 3 |

11 | Saturation constant for the rapidly hydrolyzable substrate | K_{X1} | [gCOD/gCOD] | 1.5 | 3 |

12 | Slowly hydrolysis rate constant | K_{h2} | [d^{−1}] | 5 | 3 |

13 | Saturation constant for the slowly hydrolyzable substrate | K_{X2} | [gCOD/gCOD] | 0.5 | 3 |

14 | Efficiency constant for anaerobic hydrolysis | η_{h,ana} | — | 0.18 | 2 |

15 | Relative efficiency constant for the biomass in the biofilm | ɛ | — | 0.15 | 2 |

16 | Biomass concentration in the biofilm | X_{bf} | [gCOD/m^{2}] | 10 | 3 |

Fermentation parameters | |||||

17 | Fermentation rate constant | q_{ferm} | [d^{−1}] | 2 | 3 |

18 | Saturation constant for fermentation | k_{ferm} | [gCOD/m^{3}] | 20 | 3 |

Hydrogen sulfide formation parameters | |||||

19 | Rate constant for sulfide formation | a | [gS^{0.5} m^{−0.5} h^{−1}] | 0.0032 | 3 |

20 | Sulfate saturation constant | K_{so4} | [gS/m^{3}] | 2 | 3 |

Hydrogen sulfide oxidation parameters | |||||

21 | Rate constant for sulfide oxidation in biofilm | K_{s(II)ox,f} | [gS^{0.5}gO_{2}^{−0.5}md^{−1}] | 12 | 3 |

22 | Chemical oxidation reaction order | n_{1} | — | 0.9 | 2 |

23 | Chemical oxidation reaction order | n_{2} | — | 0.15 | 2 |

24 | Rate constant for chemical sulfide oxidation of molecular sulfide | K_{H2S,c} | [(gSm^{−3})^{1−n1}(gO2 m^{−3})^{n2}d^{−1}] | 0.96 | 3 |

25 | Rate constant for chemical sulfide oxidation of ionic sulfide | K_{HS,c} | [(gSm^{−3})^{1−n1}(gO2 m^{−3})^{n2}d^{−1}] | 12 | 3 |

26 | Maximum rate constant for biological sulfide oxidation at the pH_{opt} value | K_{S(II),b,pH,opt} | [(gSm^{−3})^{1−n1}(gO2 m^{−3})^{n2}d^{−1}] | 19.92 | 3 |

27 | Constant for sulfide oxidation rate function of pH level | Ω_{S(II)b} | — | 25.00 | 2 |

. | Parameter . | Symbol . | Unit . | Calibrated value . | Uncertainty class . |
---|---|---|---|---|---|

Biomass growth in bulk water parameters | |||||

1 | Maximum aerobic specific growth rate | μ_{H,O2} | [d^{−1}] | 7 | 3 |

2 | Yield constant for heterotrophic biomass in water phase | Y_{hw} | [gCOD/gCOD] | 0.55 | 1 |

3 | Maintenance energy requirement rate constant | q_{m} | [d^{−1}] | 1 | 3 |

4 | Saturation constant for readily biodegradable substrate | K_{sw} | [gCOD/m^{3}] | 1 | 3 |

5 | Saturation constant for dissolved oxygen | K_{O} | [gO_{2}/m^{3}] | 0.05 | 3 |

Biomass growth in biofilm parameters | |||||

6 | Half-order aerobic reaction rate constant per unit area for biofilm surface | K_{1/2} | [gO_{2}^{0.5} m^{−0.5}d^{−1}] | 4 | 3 |

7 | Yield constant for heterotrophic biomass in biofilm | Y_{hf} | [gCOD/gCOD] | 0.55 | 1 |

8 | Saturation constant for readily biodegradable substrate in biofilm | K_{sf} | [gCOD/m^{3}] | 5 | 3 |

9 | Saturation constant for dissolved oxygen | K_{O} | [gO_{2}/m^{3}] | 0.05 | 3 |

Particulate hydrolysis parameters | |||||

10 | Rapidly hydrolysis rate constant | K_{h1} | [d^{−1}] | 12 | 3 |

11 | Saturation constant for the rapidly hydrolyzable substrate | K_{X1} | [gCOD/gCOD] | 1.5 | 3 |

12 | Slowly hydrolysis rate constant | K_{h2} | [d^{−1}] | 5 | 3 |

13 | Saturation constant for the slowly hydrolyzable substrate | K_{X2} | [gCOD/gCOD] | 0.5 | 3 |

14 | Efficiency constant for anaerobic hydrolysis | η_{h,ana} | — | 0.18 | 2 |

15 | Relative efficiency constant for the biomass in the biofilm | ɛ | — | 0.15 | 2 |

16 | Biomass concentration in the biofilm | X_{bf} | [gCOD/m^{2}] | 10 | 3 |

Fermentation parameters | |||||

17 | Fermentation rate constant | q_{ferm} | [d^{−1}] | 2 | 3 |

18 | Saturation constant for fermentation | k_{ferm} | [gCOD/m^{3}] | 20 | 3 |

Hydrogen sulfide formation parameters | |||||

19 | Rate constant for sulfide formation | a | [gS^{0.5} m^{−0.5} h^{−1}] | 0.0032 | 3 |

20 | Sulfate saturation constant | K_{so4} | [gS/m^{3}] | 2 | 3 |

Hydrogen sulfide oxidation parameters | |||||

21 | Rate constant for sulfide oxidation in biofilm | K_{s(II)ox,f} | [gS^{0.5}gO_{2}^{−0.5}md^{−1}] | 12 | 3 |

22 | Chemical oxidation reaction order | n_{1} | — | 0.9 | 2 |

23 | Chemical oxidation reaction order | n_{2} | — | 0.15 | 2 |

24 | Rate constant for chemical sulfide oxidation of molecular sulfide | K_{H2S,c} | [(gSm^{−3})^{1−n1}(gO2 m^{−3})^{n2}d^{−1}] | 0.96 | 3 |

25 | Rate constant for chemical sulfide oxidation of ionic sulfide | K_{HS,c} | [(gSm^{−3})^{1−n1}(gO2 m^{−3})^{n2}d^{−1}] | 12 | 3 |

26 | Maximum rate constant for biological sulfide oxidation at the pH_{opt} value | K_{S(II),b,pH,opt} | [(gSm^{−3})^{1−n1}(gO2 m^{−3})^{n2}d^{−1}] | 19.92 | 3 |

27 | Constant for sulfide oxidation rate function of pH level | Ω_{S(II)b} | — | 25.00 | 2 |

The values of should be in the range of −1 to 1. For the values of to be valid, the coefficient of determination of the linear regression should be high enough, i.e. , to ensure that the model is adequately linear. The absolute values of are used to determine the sensitivity of the model output to the corresponding input and, consequently, the contribution of that model input to the output variance.

### Distributed parameter CFD model

While the process model divides the geometry of pipe into a number of CSTRs, assuming homogeneity in each CSTR, the geometry is discretized into a number of computational cells in the CFD model. In each of these cells, the governing conservation equations are integrated using the finite volume method to arrive at a coupled set of algebraic equations. The geometry of the computational domain is based on that described in measurement data, where the domain is discretized using ANSYS ICEM meshing software to produce a structured-like mesh with an O-grid profile at the cross-section. Due to the length of the domain, the computational cells are created with a high aspect ratio to reduce the number of the cells, hence the computational cost. Therefore, it is necessary to produce cell faces perpendicular to the flow direction to reduce the error coming from the interpolation of the different variables at the face centroids. A grid independence test was done on three different meshes containing 1.82, 3.17, and 5.88 million cells. The velocity profiles at different sections were compared and it was determined that the difference in the velocity values between the intermediate the finest mesh is less than 5%. Therefore, the intermediate mesh was used as a compromise between the accuracy and the computational cost.

ANSYS FLUENT CFD software was used for the simulation of the hydrodynamics of the flow in the pipe in addition to the settling behaviour of the solids in the flow. Due to the dynamic nature of the sewer system, a transient simulation was performed to simulate the dynamic behaviour of the flow and species concentration at the domain inlet. A transient table of the inlet conditions was defined for the CFD case providing the variation of the inlet flow rate and the concentrations of the simulated species. The biological and chemical reactions in the sewage were simulated along with the transport of each species in the sewage. The simulated species in the sewage were defined based on the integrated biochemical model. In the following sections, the governing equations for the hydrodynamics, modelling of the solids settling and surface reactions in the biofilm model development and the integration between the CFD and the biochemical model are presented.

#### Governing equations

*i*direction. Momentum conservation equation for turbulent flow is described as:

Here, is the mean pressure, is the dynamic viscosity, is Reynolds stress tensor, and is the velocity fluctuation. The Reynolds stress tensor, the final term in Equation (4), is closed using the realizable k-*ɛ* turbulence model.

#### Biochemical model integration with CFD

Here is the concentration of the species, is the diffusivity of the scalar in water, and is the source term of the species. The turbulent diffusivity is included in and is defined using the computed turbulent viscosity assuming a turbulent Schmidt number of 0.7.

##### Production/consumption rate of the wastewater constituents

The source terms, , for the different species is defined by the corresponding conversion rates in the biochemical model and is implemented using user-defined functions. The source term describes the production and consumption rates of each species. These rates were computed using the local values of the dependent species concentrations in each computational cell, thereby accounting for the heterogeneity.

##### Surface reaction modelling

A technique is developed for the surface reactions in the biofilm such that these reactions are modelled as a superficial flux from the wall faces. A user-defined function is used to identify wall faces and, for all cells adjacent to a wall, the flux from the wall face is defined with direction defined based on the net generation/consumption of the corresponding species and magnitude defined by the biochemical model. More details and verification of the technique used for the simulation of the surface reactions are provided in the Supplementary information.

*a*. Special attention is drawn to this parameter since it is directly affected by the homogenization of the surface reactions model as described in the sulfide formation rate:

Here, , , , and are the concentrations of fermentable COD, VFA, soluble oxygen, and rapidly hydrolysable substrate, respectively. is the oxygen saturation constant and the factor is the homogenization factor that converts the areal sulfide formation rate to volumetric form in the process model structure. In the CFD model, this homogenization factor is not required, such that the model uncertainty resulting from the introduction of this parameter can be assessed though the CFD results.

##### Solids settling

*et al.*(1999) to model the setting of the solids in a secondary clarifier. The settling velocity of the solids is a function of the solids concentration, which is evaluated as a sum of the particulate species in each computation cell. The settling velocity is then calculated using the double exponential model (Takács

*et al.*1991) which is defined as:where is a reference settling velocity, is the non-settleable concentration, and and are model parameters. The model parameters for this model are highly dependent on the nature of the solids in the wastewater. For approximation, the data from Patziger & Kiss (2015) measured for primary sludge are used to determine the model parameters for this study. Takács's model was originally proposed for activated sludge. However, it was adapted in the work of Patziger & Kiss (2015) to model the settling of solids in the primary clarifier. The maximal settling velocity in the primary clarifier is much higher than that of activated sludge. This is due to the relatively higher density of primary sludge and the larger particle size that varies in the range of 0.01–0.5 mm (Patziger & Kiss 2015). This results in a settling velocity of 10–11 cm/s at concentration of 100–200 mg/L suspended solids.

## RESULTS AND DISCUSSION

### Monte Carlo simulations

The global uncertainty of the model is assessed using Monte Carlo simulations without assuming any correlation between the input parameters. Time series of the concentrations of soluble COD, VFA and sulfide at the pipe outlet were monitored for each simulation. The upper and lower bands, along with the 10th and 90th percentiles, of the Monte Carlo simulation outputs are shown in Figure 2. A higher band represents the largest difference between an instance of the Monte Carlo simulation output and the mean profile of the all simulations in the positive direction. Conversely, the lower band represents the largest difference in the negative direction. These give an overall indication of the uncertainty level of the model, where a larger difference between the Monte Carlo bands and the mean indicates a higher uncertainty.

### Uncertainty/sensitivity analysis

#### Graphical representation

Figure 3 shows the influence of the key parameters of the model (as defined in Table 1) on the model outputs, with respect to the uncertainty envelope of the Monte Carlo simulations output. The parameters related to the biomass growth in the bulk water, which describe the suspended biomass growth kinetics on the sCOD, show no effect on the model outputs. This can be explained on the basis that the present study involves a pressure main where the dissolved oxygen concentration is low enough that anaerobic conditions prevail. The same conclusion can be made for the biomass growth in the biofilm parameters as well. The hydrolysis processes are more impactful on the model output. The anaerobic hydrolysis rate constant,, and efficiency,, show no direct effect on the sulfide and VFA levels. However, they have a significant effect on the sCOD, since these parameters control the rate of conversion of particulate COD to fermentable substrate. Therefore, by increasing and , the concentration of the sCOD increases. However, and are not the predominant parameters for determining the variation in the sCOD. This can be concluded by comparing with the upper and lower bands of Monte Carlo simulations. Sulfide formation rate, *a*, and saturation constant of rapidly hydrolysable substrate,, play a significant role in determining the correct sCOD level as well.

The average sensitivity measure, , of and on sCOD are 53.8% and 25.4%, respectively. The saturation constant for the rapidly hydrolysable substrate shows a high influence on the sCOD with relative sensitivity of 19.8%. Fermentation rate, , is another key parameter in calibrating the model where its effect is clear on the VFA levels only (Figure 3(b)). The sensitivity measure of the fermentation rate constant on the VFA concentrations is 67.6%, whereas the other parameters do not exceed 13.3% with respect to VFA. The sulfide formation constant, *a*, is a key parameter that has a large influence on all of the model outputs. This is clear in Figure 3, where its value has a significant impact on the concentration of sCOD and VFA, and it is the major parameter that controls the sulfide concentration output as shown in Figure 3(c). The sensitivity measure of the sulfide formation rate constant on the sCOD, VFA, and sulfide are 23.3%, 13.3% and 81.6%, respectively.

From the data presented, it is concluded that the sulfide formation rate constant shows a significant influence on all of the monitored species. It is the predominant parameter affecting the predicted sulfide. Conversely, the fermentation rate constant is the main parameter controlling the concentration of VFA, but it has almost no effect on the other species. This can be explained by the fact that the fermentation process is responsible for converting the fermentable substrate to acetate, which is another form of sCOD. As a result, the overall sCOD does not change as a result of fermentation. In addition, the sulfide formation expression is a function of sCOD, thereby the concentration of the sulfide is not determined by the fermentation process rate. The anaerobic hydrolysis efficiency and hydrolysis rate constant of the rapidly hydrolysable substrate have a predominant influence on sCOD since anaerobic conditions prevail in pressure mains. For sulfide formation, the colloidal COD is considered as a rapidly hydrolysable substrate and is added to sCOD in the sulfide formation rate expression. Therefore, the effect of the anaerobic hydrolysis efficiency does not propagate to the sulfide concentration output of the model. In general, the hydrolysis process parameters are the most influential parameters on the sCOD output.

#### Standardized linear regression (SRC)

Three multivariate linear models are fitted using linear regression (implemented in MATLAB) for the averaged concentrations of sCOD, VFA, and total sulfide at the sewer outlet. The fitted models can reasonably reproduce the time-averaged output of the Monte Carlo simulations as depicted in Figure 4. The figure shows a comparison between the probability distribution of the Monte Carlo simulations and the fitted model outputs where clear overlap is obtained except some differences (blue and orange part of the bars illustrate higher probability by the Monte Carlo simulations and by the fitted models, respectively). This is illustrated by the determination coefficient reported in Table 2 for each of the model outputs where for all the outputs. This indicates that the linear effect of parameter variability could account for the majority of average output variation. The regression coefficients, , can be used to evaluate the contribution of each of the model parameters on the overall variance. Furthermore, the condition of the standardized regression coefficients is generally satisfied, as shown in the last row of Table 2. The parameters that show the greatest influence () on the model output hence highest sensitivity and contribution in the overall variance can be deduced from Table 2. The SRC results support the results obtained from the graphical representation, where the parameters that have the most important influence on the results of the model are the same, namely, , , , , and *a* are the most influential parameters in the model. Further to this, the sign of the standardized regression coefficients indicates the direction of the influence of the parameter on the model output while the magnitude reflects the strength of the parameter effect.

. | Soluble COD . | VFA . | Total sulfide . | |||
---|---|---|---|---|---|---|

Parameter . | 0.974 . | 0.981 . | 0.981 . | |||

R^{2}
. | Estimate . | p-Value
. | Estimate . | p-Value
. | Estimate . | p-Value
. |

μ_{H} | −0.050 | 0.000 | −0.023 | 0.000 | 0.035 | 0.000 |

K_{sw} | −0.001 | 0.859 | 0.000 | 0.913 | −0.003 | 0.086 |

K_{O} | 0.002 | 0.678 | −0.006 | 0.148 | −0.015 | 0.000 |

K_{sf} | −0.012 | 0.018 | −0.005 | 0.234 | −0.001 | 0.672 |

K_{1/2} | −0.007 | 0.176 | 0.001 | 0.739 | 0.004 | 0.025 |

q_{m} | −0.004 | 0.470 | −0.005 | 0.266 | 0.005 | 0.003 |

K_{h1} | 0.773 | 0.000 | 0.038 | 0.000 | −0.001 | 0.578 |

Kx_{1} | − 0.299 | 0.000 | −0.013 | 0.003 | 0.000 | 0.968 |

η_{h,ana} | 0.361 | 0.000 | 0.000 | 0.929 | −0.002 | 0.269 |

ɛ | 0.045 | 0.000 | 0.063 | 0.000 | −0.001 | 0.740 |

a | − 0.340 | 0.000 | − 0.161 | 0.000 | 0.953 | 0.000 |

K_{so4} | 0.035 | 0.000 | 0.019 | 0.000 | −0.091 | 0.000 |

q_{ferm} | 0.005 | 0.323 | 0.970 | 0.000 | 0.001 | 0.726 |

K_{ferm} | 0.006 | 0.274 | −0.067 | 0.000 | 0.001 | 0.446 |

n_{1} | 0.003 | 0.565 | 0.002 | 0.729 | −0.035 | 0.000 |

n_{2} | 0.002 | 0.642 | −0.003 | 0.447 | 0.240 | 0.000 |

K_{s(II)ox} | −0.003 | 0.587 | −0.002 | 0.629 | − 0.151 | 0.000 |

K_{H2S,c} | 0.002 | 0.746 | −0.002 | 0.612 | 0.001 | 0.627 |

K_{HS,c} | 0.001 | 0.920 | 0.003 | 0.481 | −0.001 | 0.494 |

K_{S(II),b,pH,opt} | 0.002 | 0.690 | −0.002 | 0.612 | −0.026 | 0.000 |

Ω_{S(II)b} | −0.001 | 0.804 | 0.007 | 0.109 | 0.000 | 0.973 |

0.977 | 0.977 | 1.001 |

. | Soluble COD . | VFA . | Total sulfide . | |||
---|---|---|---|---|---|---|

Parameter . | 0.974 . | 0.981 . | 0.981 . | |||

R^{2}
. | Estimate . | p-Value
. | Estimate . | p-Value
. | Estimate . | p-Value
. |

μ_{H} | −0.050 | 0.000 | −0.023 | 0.000 | 0.035 | 0.000 |

K_{sw} | −0.001 | 0.859 | 0.000 | 0.913 | −0.003 | 0.086 |

K_{O} | 0.002 | 0.678 | −0.006 | 0.148 | −0.015 | 0.000 |

K_{sf} | −0.012 | 0.018 | −0.005 | 0.234 | −0.001 | 0.672 |

K_{1/2} | −0.007 | 0.176 | 0.001 | 0.739 | 0.004 | 0.025 |

q_{m} | −0.004 | 0.470 | −0.005 | 0.266 | 0.005 | 0.003 |

K_{h1} | 0.773 | 0.000 | 0.038 | 0.000 | −0.001 | 0.578 |

Kx_{1} | − 0.299 | 0.000 | −0.013 | 0.003 | 0.000 | 0.968 |

η_{h,ana} | 0.361 | 0.000 | 0.000 | 0.929 | −0.002 | 0.269 |

ɛ | 0.045 | 0.000 | 0.063 | 0.000 | −0.001 | 0.740 |

a | − 0.340 | 0.000 | − 0.161 | 0.000 | 0.953 | 0.000 |

K_{so4} | 0.035 | 0.000 | 0.019 | 0.000 | −0.091 | 0.000 |

q_{ferm} | 0.005 | 0.323 | 0.970 | 0.000 | 0.001 | 0.726 |

K_{ferm} | 0.006 | 0.274 | −0.067 | 0.000 | 0.001 | 0.446 |

n_{1} | 0.003 | 0.565 | 0.002 | 0.729 | −0.035 | 0.000 |

n_{2} | 0.002 | 0.642 | −0.003 | 0.447 | 0.240 | 0.000 |

K_{s(II)ox} | −0.003 | 0.587 | −0.002 | 0.629 | − 0.151 | 0.000 |

K_{H2S,c} | 0.002 | 0.746 | −0.002 | 0.612 | 0.001 | 0.627 |

K_{HS,c} | 0.001 | 0.920 | 0.003 | 0.481 | −0.001 | 0.494 |

K_{S(II),b,pH,opt} | 0.002 | 0.690 | −0.002 | 0.612 | −0.026 | 0.000 |

Ω_{S(II)b} | −0.001 | 0.804 | 0.007 | 0.109 | 0.000 | 0.973 |

0.977 | 0.977 | 1.001 |

The values in bold are those with the most significant effects on the model outputs.

Improvement of the model accuracy can be obtained by reducing the uncertainty of these parameter values. The uncertainty reduction required for the model parameters, in order to achieve a specific reduction in variance in the model outputs, can be quantified (Sin *et al.* 2011). Therefore, the most influential parameters should be studied in more detail, with more laboratory experiments being conducted to reduce the uncertainty of the overall model. For instance, the sulfide production rate constant, *a*, is emphasized to be the most influential parameter on the model, as it significantly affects all the outputs of the model as illustrated from the graphical and analytical results. Thus, special attention should be drawn to get an accurate value of this parameter.

### CFD model results

First, the influence of the homogenization of the surface reactions in the biofilm is studied using the CFD model. Three simulations of the biochemical conversions along with the flow field in the sewer pipe were conducted using different values of *a*. The upper and lower limits of the *a* variation range used in the Monte Carlo simulation are used along with best fit value obtained from the calibration using Weighted Sum of Squared Errors (WSSE) criteria, which is . The mass-weighted average of the concentrations of the simulated species are monitored at the outlet of pipe and plotted in Figure 5. The 1-D model simulations for the same parameter set are plotted on the same figure to compare between the outputs of the two models. The comparison shows that similar variation of the results of the CFD model is obtained as that of the 1-D CSTR-in-series model structure. It is obvious that the CSTR-in-series model tends to smear out the high dynamic variation in the inlet concentration. This could be explained by the fact that the CFD model is, in general, better in modelling the advection and diffusion influence than the 1-D CSTR-in-series model. However, the mean values of the different concentrations predicted by the CFD model and 1-D CSTR-in-series model are similar. Thus, the detailed modelling of the surface reactions in the case of pipe flow shows that the approximation used in the process model of the homogenization of the surface reactions is acceptable in this case. However, this cannot be generalized for all segments of the sewer network and other surface reactions in the biological treatment stages in WWTP. The CFD model demonstrated here could be used in such cases to determine whether or not the heterogeneous nature of the reactions is important to take into account or whether a homogenized CSTR-in-series model is acceptable. Moreover, the dynamic response of the 1-D CSTR-in-model could be examined using the CFD model to ensure that the variation is in the design limit of the chemical treatment dosing of sewer systems.

Figure 6(a) shows the spatial distribution of the total COD at a pipe cross-section at the middle of the simulated pipe. It is clear that because of the settling of the particulate COD, a higher concentration is noticed at the bottom of the pipe. However, the variation of the total COD is not significant across the cross-section of the pipe, such that the completely mixed assumption of the 1-D CSTR-in-series model is deemed to be acceptable. The same can observed from the soluble COD contours as shown in Figure 6(b), due to the utilization of the soluble COD by the SRB in the biofilm, a higher utilization rate is observed near the wall. Therefore, a lower concentration of the soluble COD is near the wall. However, the variation of the soluble COD concentration across the cross-section of the pipe is negligible. This due to the diffusion of the soluble matter and the hydrolysis process which compensates for the soluble COD consumption by the SRB. Since the production of sulfide is mainly at the wall, slightly higher sulfide concentration was noticed near the wall. However, due to diffusion, this difference between the sulfide concentration at wall and in the bulk water is negligible.

To assess the impact of solids settling on the model outputs, two simulations were performed: one with the solids settling term included in the CFD model and the other without considering this term. Settling of the solids does have noticeable influence on the mass-weighted average of the species concentrations if compared with the simulations that disregards the solids settling (results not shown). This indicates that the solids settling for the low solids concentration in the sewer may not have a significant influence on the overall simulation outputs in straight pipes. Again, this conclusion should not be generalized to all sewers, and this methodology can be considered as a useful tool in assessing the influence of settling for other operating conditions as well. For example, the case of intermittent flow is not considered in this work and could result in solids settling being important.

For illustration of the effect of the solids settling on the profiles of the particulate species along the pipe length, the particulate COD profiles at different distances from the inlet are plotted in Figure 7. From the figure, the concentration gradients in the radial direction are evolved from the uniform distribution assumed at the inlet (zero gradient) to more noticeable values based on the solids concentration. However, due to transient concentration of the solids at the inflow of the sewer pipe, the concentration gradient of the particulate COD is not uniform in the longitudinal direction. It is noticeable that the radial concentration gradient increases with the level of the solids concentration as described in Equation (6).

In conclusion, the influence of the detailed modelling of the solids settling and chemical and biological reactions heterogeneity has a marginal in the case of straight pressure sewer pipes. This could be explained by the fact that the gradients of the different species due to solids setting and heterogeneous model of the sulfide formation are low since the convection of the species in the flow direction is predominant on the distribution of the species concentrations.

## CONCLUSIONS

SRC analysis of outcomes from time-averaged outputs from the biochemical WATS model demonstrated that parameters could account for >90% of the output variance observed in soluble COD, VFA, and total sulfides. Parameters pertaining to hydrolysis, efficiency of the anaerobic processes, and sulfide production rate constant per unit area of biofilm area were the most significant predictors of these outputs. Comparison of 1-D biochemical and multidimensional CFD analysis, including settling, turbulent flow, and biochemical reactions indicated that the 1-D model was an adequate representation of the multi-dimensional model, with increased short-term time-dependent behaviour in the multi-dimensional model. Further studies should be conducted in more complex segments in the sewer system that could be hot spots for hydrogen sulfide production. In such locations, it would be beneficial to assess the accuracy of the process models in predicting the hydrogen sulfide levels in light of the highly heterogeneous reactions. While CFD is not justified by improved predictive capability in output values for the case considered, it provides mechanistic analysis into the way that sedimentation, reaction, and hydraulic processes can interact.

## ACKNOWLEDGEMENTS

This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) [grant numbers RGPIN-2017-04078, CRDPJ-488704-15].

## DATA AVAILABILITY STATEMENT

Data cannot be made publicly available; readers should contact the corresponding author for details. Further information and figures are available in Supplementary Information.