the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Multivariate autoregressive modelling and conditional simulation for temporal uncertainty analysis of an urban water system in Luxembourg
Jairo Arturo TorresMatallana
Ulrich Leopold
Gerard B. M. Heuvelink
Uncertainty is often ignored in urban water systems modelling. Commercial software used in engineering practice often ignores the uncertainties of input variables and their propagation because of a lack of userfriendly implementations. This can have serious consequences, such as the wrong dimensioning of urban drainage systems (UDSs) and the inaccurate estimation of pollution released to the environment. This paper introduces an uncertainty analysis in urban drainage modelling, built on existing methods and applied to a case study in the HauteSûre catchment in Luxembourg. The case study makes use of the EmiStatR model which simulates the volume and substance flows in UDS using simplified representations of the drainage system and processes. A Monte Carlo uncertainty propagation analysis showed that uncertainties in chemical oxygen demand (COD) and ammonium (NH_{4}) loads and concentrations can be large and have a high temporal variability. Furthermore, a stochastic sensitivity analysis that assesses the uncertainty contributions of input variables to the model output response showed that precipitation has the largest contribution to output uncertainty related with water quantity variables, such as volume in the chamber, overflow volume, and flow. Regarding the water quality variables, the input variable related to COD in wastewater has an important contribution to the uncertainty for the COD load (66 %) and COD concentration (62 %). Similarly, the input variable related to NH_{4} in wastewater plays an important role in the contribution of total uncertainty for the NH_{4} load (34 %) and NH_{4} concentration (35 %). The Monte Carlo (MC) simulation procedure used to propagate input uncertainty showed that, among the water quantity output variables, the overflow flow is the most uncertain output variable, with a coefficient of variation (cv) of 1.59. Among water quality variables, the annual average spill COD concentration and the average spill NH_{4} concentration were the most uncertain model outputs (coefficients of variation of 0.99 and 0.82, respectively). Also, low standard errors for the coefficient of variation were obtained for all seven outputs. These were never greater than 0.05, which indicates that the selected MC replication size (1500 simulations) was sufficient. We also evaluated how the uncertainty propagation can more comprehensively explain the impact of water quality indicators for the receiving river. While the mean model water quality outputs for COD and NH_{4} concentrations were slightly above the threshold, the 0.95 quantile was 2.7 times above the mean value for COD concentration and 2.4 times above the mean value for NH_{4}. This implies that there is a considerable probability that these concentrations in the spilled combined sewer overflow (CSO) are substantially larger than the threshold. However, COD and NH_{4} concentration levels of the river water will likely stay below the water quality threshold, due to rapid dilution after CSO spill enters the river.
 Article
(5946 KB) 
Supplement
(885 KB)  BibTeX
 EndNote
Combined sewer systems are important components of urban water infrastructure. These systems are typically found in old and large cities (Baker, 2009; Litrico and Fromion, 2009) and are designed to transport the water generated and accumulated in an urban catchment to the receiving water body. During normal conditions, all water is transported to the treatment facility before it is released to the environment. This is the socalled throttled outflow or passforward flow (Hager, 2010). However, during extreme conditions with heavy precipitation, the combined sewer overflow (CSO) discharges excess water directly to nearby streams, rivers, lakes, or other water bodies (Baker, 2009). The CSO contains polluted water and solid matter (Hager, 2010), which, when released to the environment, can have a damaging impact on the water quality status of the receiving waters (BachmannMachnik et al., 2018; Gasperi et al., 2012). CSO pollutant load emissions are of similar or greater magnitude than the emissions from wastewater treatment plants (Gasperi et al., 2012; BachmannMachnik et al., 2018). CSO discharge impacts are mainly high peak flows, high organic loads from single events, which can lead to oxygen depletion, and ecotoxic concentrations of ammonia (NH_{3}) (Miskewitz and Uchrin, 2013; BachmannMachnik et al., 2018). To reduce pollution in receiving waters it is important to minimise CSO load and concentration.
One of the main variables is the chemical oxygen demand (COD), which is an indicator of organic compounds in water. It is used to measure the effluent quality (Viana da Silva et al., 2011). High levels of COD are correlated with a decrease in the amount of dissolved oxygen (DO) available for aquatic organisms. A depletion in the DO concentration in the water column from near 9 mg L^{−1} (the maximum solubility of oxygen in estuarine water on an average summer day) to below 2 mg L^{−1} is referred to as hypoxia. If hypoxic conditions are reached, the health of the ecosystem is affected, causing physiological stress, and even death, to aquatic organisms (on Environmental and atural Resources – CENR, 2003). Ammonium (NH_{4}) is another important variable and is an indicator of nitrogen compounds in water. Concentrations of NH_{4} in water and wastewater are relevant because high levels of nitrogen in receiving waters can cause eutrophication and, therefore, excessive growth of algae and other microorganisms, resulting in oxygendissolved depletion and fish toxicity (Huang et al., 2010).
To better assess environmental impacts, numerical models are applied in urban hydrology to simulate CSO emissions into the environment. It is recommended, however, that such modelling approaches consider the inherent uncertainty associated with the system representation and the approximation of the model to the reality (Hutton et al., 2011). Moreover, the model inputs are also not free of errors, and associated uncertainties will also propagate to the model output (Heuvelink, 1998).
A total of five approaches for representing the presence or absence of uncertainty and how it is represented in the context of urban water systems are often distinguished (Walker et al., 2003; Refsgaard et al., 2007; van der Keur et al., 2008; Bach et al., 2014) as follows: (1) determinism, (2) statistical uncertainty, (3) scenario uncertainty, (4) recognised ignorance, and (5) total (unrecognised) ignorance. Following van der Keur et al. (2008), determinism applies when we have knowledge with absolute certainty about the system under analysis. This is the ideal world case, which is not realistic for urban hydrology systems. The statistical approach is useful when it is possible to describe uncertainty in statistical terms, i.e. when uncertainty can be characterised by probability distribution functions (pdf's). The scenario approach, in contrast, applies when quantitative probabilities cannot be determined and, instead, qualitative measures of uncertainty are used. It is used when possible outcomes of uncertain inputs are known but the probabilities of these outcomes are not (Brown, 2004). There is also no claim that the list of possible outcomes (scenarios) is exhaustive. Recognised ignorance occurs when there is an awareness of lack of knowledge but without any further possibility to process and address the recognised uncertainty. This is the case for very complex functional or inherently unidentifiable relationships, when, for example, predictions are infeasible due to the chaotic behaviour of the system or when our understanding of the system behaviour is too limited (van der Keur et al., 2008). This is common in social systems where the behaviour of humans and groups of humans may often unpredictable. Finally, total ignorance is the state of complete lack of awareness about imperfect knowledge (van der Keur et al., 2008). It is the opposite of determinism and reflects a state in which we do not know that we do not know (Walker et al., 2003). Among the approaches described above, in this paper we will use the statistical approach to characterise and propagate uncertainties.
A total of three main sources of uncertainty in the context of performance evaluation analysis and design of urban water infrastructure and urban drainage modelling are identified (Walker et al., 2003; Neumann, 2007; Deletic et al., 2012). First, model input uncertainty is related to errors in input data, i.e. in driving forces such as precipitation. Second, parameter uncertainty is related to the uncertainty regarding the (calibrated) parameters of the model. Third, model structural uncertainty relates to uncertainty due to model conceptualisation and simplification. For instance, an urban drainage model might ignore certain subprocesses, such as evaporation or chemical transformation, or might simplify a nonlinear relation between model variables to a linear relation. These types of uncertainties are not captured in model input and model parameter uncertainty and are represented by model structural uncertainty. The focus of this work is on the propagation of model input uncertainty.
Regarding methods for uncertainty propagation analysis, a distinction can be made between analytical methods, such as the Taylor series method (Heuvelink, 1998), and numerical techniques, such as Monte Carlo (MC) simulation. Numerical techniques are more flexible and hence more convenient for analysing uncertainty propagation with complex models (Zoppou, 2001). MC simulations are computationally demanding, especially in the case of complex models, but they can still be used if there are sufficient computational resources (Bastin et al., 2013), among others, because it can greatly benefit from parallel computing.
Although uncertainty propagation analysis has been applied extensively in hydrologic modelling (e.g. Beven and Binley, 1992; Kuczera and Parent, 1998; Hutton et al., 2011; Vrugt et al., 2003a, b; Vrugt and Robinson, 2007; Renard et al., 2010; Datta, 2011), the number of applications of longterm simulations in urban drainage modelling is limited and typically does not consider the influence of the temporal and spatial correlation in the analysis of the propagation of input uncertainty. Temporal correlation occurs in uncertain dynamic variables, such as precipitation and COD of household wastewater, because values of these variables over short time lags will be more similar than over large time lags. The same concept applies to variables that are spatially distributed (Webster and Oliver, 2007). It is important to take the temporal (and spatial) correlation of uncertain inputs into account because this may have a major influence on the outcomes of an uncertainty analysis (Heuvelink, 1998). In this paper, we perform a temporal uncertainty propagation analysis in urban water modelling using a MC simulation. As a case study, we use the simplified model EmiStatR (TorresMatallana et al., 2018b) to predict wastewater volume, COD, and NH_{4} concentrations in CSOs for three urban–rural subcatchments of the HauteSûre catchment in the northwest of Luxembourg.
The objectives of this study are to (1) select and characterise the main sources of input uncertainty accounting for the temporal auto and crosscorrelation within EmiStatR; (2) propagate input uncertainty through EmiStatR, taking into account the temporal auto and crosscorrelation of uncertain dynamic inputs; and (3) quantify and assess the contributions of each uncertainty source to model output uncertainty dynamically (over time) for the Luxembourg case study.
2.1 The EmiStatR model
EmiStatR is used to simulate CSO flows and water quality concentrations. Details regarding the conceptual and mathematical model are provided in TorresMatallana et al. (2018b). A list of the EmiStatR inputs and outputs is provided in Appendix A. The main components of the EmiStatR model are (1) dry weather flow (DWF), including infiltration flow (IF), (2) pollution of DWF, (3) rain weather flow (RWF), (4) pollution of RWF, (5) combined sewer flow (CSF) and pollution, and (6) combined sewer overflow (CSO) and pollution. Figure 1 illustrates the scheme of the sewer system analysed.
Basically, the total dry weather flow, Q_{DWF} (L s^{−1}), is calculated as follows:
where ${Q}_{{\mathrm{DWF}}_{t}}$ (L s^{−1}) is the dry weather flow at time t, and ${Q}_{{\mathrm{s}}_{t}}$ (L s^{−1}) is the dry weather flow of the residential sewage in the catchment at time t, calculated as 86 400^{−1} pe_{t} qs_{t} (where $\mathrm{86}\phantom{\rule{0.125em}{0ex}}\mathrm{400}=\mathrm{24}\times \mathrm{60}\times \mathrm{60}$ is a measurement unit conversion factor), with pe_{t} (PE) the population equivalents of the connected CSO structure at time t, and qs_{t} $\left(\mathrm{L}\phantom{\rule{0.125em}{0ex}}{\text{PE}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\text{d}}^{\mathrm{1}}\right)$ is the individual water consumption of households at time t. ${Q}_{{\mathrm{f}}_{t}}$ (L s^{−1}) is the infiltration flow at time t that enters the pipes from groundwater flow through cracks and joints, calculated as ${A}_{\mathrm{imp}}\cdot {q}_{{\mathrm{f}}_{t}}$, where A_{imp} (ha) is the impervious area of the catchment, and ${q}_{{\mathrm{f}}_{t}}$ $\left(\mathrm{L}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\text{ha}}^{\mathrm{1}}\right)$ is the infiltration water inflow flux (specific infiltration discharge from groundwater flow) at time t. Variables qs_{t} and pe_{t} are dynamic and can be defined as time series with daily, weekly, and seasonal patterns.
The contribution of rainwater to the combined sewage flow, Q_{r} (m^{3} s^{−1}), is derived from precipitation as follows:
where $\frac{\mathrm{1}}{\mathrm{6}}$ is a factor for units conversion, ${P}_{t{t}_{\mathrm{fS}}}$ is precipitation at time t−t_{fS} (mm min^{−1}), t_{fS} is a delay in time response related to flow time in the sewer system, A_{imp} is the impervious area of the catchment (ha), A_{total} is the total area of the catchment (ha), C_{imp} is the runoff coefficient for impervious areas (–), and C_{per} is the runoff coefficient for pervious areas (–). From ${Q}_{{\mathrm{r}}_{t}}$, the CSO volume calculation is based on the exceeding volume stored in the combined sewer overflow chamber (CSOC). The CSO volume depends on four CSOC stages, namely (1) filling up, (2) CSO spill volume, (3) stagnation, and (4) emptying. The sum of the total dry weather flow, ${Q}_{{\mathrm{DWF}}_{t}}$, and the rain water flow, ${Q}_{{\mathrm{r}}_{t}}$, is called combined sewer flow at time t, ${Q}_{{\mathrm{CSF}}_{t}}$.
The COD load, B_{COD, Sv} (g), in the spill overflow volume is calculated as a function of the spill overflow volume at time t, ${V}_{{\mathrm{Sv}}_{t}}$ (m^{3}), a combined sewer mixing ratio at time t, ${\mathrm{cs}}_{{\mathrm{mr}}_{t}}$ (–), the mean dry weather pollutant concentration at time t, ${C}_{{\mathrm{COD}}_{t}}$ (mg L^{−1}), and the concentration due to rainwater pollution at time t, ${\mathrm{COD}}_{{\mathrm{r}}_{t}}$ (mg L^{−1}) as follows:
The variable ${V}_{{\mathrm{Sv}}_{t}}$ depends directly on the water volume in the CSO chamber at time t, ${V}_{{\mathrm{Chamber}}_{t}}$ (m^{3}). It is computed as follows:
where ${V}_{{\mathrm{r}}_{t}}$ is the rain weather volume at time t accumulated during a time interval Δt (min), ${V}_{{\mathrm{dw}}_{t}}$ (m^{3}) is the total dry weather volume (amount of dry weather water in combined sewage flow) at time t, ${V}_{{\mathrm{d}}_{t}}$ is the volume of throttled outflow to the wastewater treatment plant (WWTP) at time t (m^{3}), V (m^{3}) is the CSOC volume, and ϵ is a numerical precision term set equal to 10^{−5} (m^{3}). While V_{Sv}, cs_{mr}, and C_{COD} are dynamic, ${\mathrm{COD}}_{{\mathrm{r}}_{t}}$ can either be dynamic or assumed constant if the pollution concentration is assumed constant in time. ${C}_{{\mathrm{COD}}_{t}}$ (mg L^{−1}) is calculated as follows (TorresMatallana et al., 2018b):
where C_{COD, S} is the COD sewage pollution per capita (PE) load per day $\left(\mathrm{g}\phantom{\rule{0.125em}{0ex}}{\text{PE}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\text{d}}^{\mathrm{1}}\right)$. Similar equations, as above, apply to the second water pollution indicator NH_{4}.
2.2 Sewer system in the HauteSûre catchment
The study area is composed of three subcatchments of the HauteSûre catchment in the northwest of Luxembourg. The combined sewer system drains three villages, namely Goesdorf (GOE), Kaundorf (KAU), and NocherRoute (NOR). The local sewer system downstream of each village has a CSO structure to store pollutant peaks in the first flush of combined sewage flows. Figure 2 depicts the location of the CSO structures and the delineation of the subcatchments. The main land use types in the villages are residential, smaller industries, and farms. Outside the villages, forest and agricultural arable and grassland are the dominating land uses. The receiving water bodies of the CSO structures are tributaries of the river Sûre (or Sauer in German).
2.3 Input data
The input variables of the EmiStatR model are shown in Table S1 in the Supplement. Following TorresMatallana et al. (2018b), seven input variables were calibrated, namely water consumption (qs_{t}), infiltration flow (${q}_{{\mathrm{f}}_{t}}$), flow time structure equivalent to the time of concentration to the combined sewer overflow structure (t_{fS}), runoff coefficient for impervious area (C_{imp}), runoff coefficient for pervious area (C_{per}), orifice coefficient of discharge (C_{d}), and the initial water level (Lev_{ini}). The main objective of the calibration process is to appropriately represent the water volume in the CSOC.
The observed precipitation (P_{t}) is a 1 year time series for 2010 at a 10 min time interval, measured at stations EschsurSûre and Dahl (Fig. 2). The variable water consumption (qs_{t}) is also dynamic and represented as a time series with a daily pattern according to factors proposed in the German ATVDVWKA 134 guideline (Evers et al., 2000).
The hydraulic variable measured is the water level in the CSOC, namely Lev (m). The temporal resolution of the measurements of Lev is 30 s. Regarding the wastewater quality (WWQ) characterisation, values of C_{COD, S} and ${C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{S}}$ in the wastewater were derived from DWF measurements at Goesdorf, Kaundorf, and NocherRoute. A total of 91 2 h composite samples were taken and measured in the laboratory to determine the concentrations of COD (mg L^{−1}) and NH_{4} (mg L^{−1}). This led to the identification of seven at Goesdorf on 4 May 2011, 48 between 19 June and 21 July 2010 at Kaundorf, and 36 between 9 March and 2 August 2011 at NocherRoute. The variables COD_{f} and ${{\mathrm{NH}}_{\mathrm{4}}}_{\mathrm{f}}$ were set to zero because the pollution contribution of the infiltration water is negligible in the study area. The contribution of ammonium from rainwater ${{\mathrm{NH}}_{\mathrm{4}}}_{\mathrm{r}}$ was assumed constant and set to 2.00 mg L^{−1}, while COD_{r} was equal to zero. Table S1 summarises the base values of the general input variables, and Table S2 presents the general characteristics of each CSO structure for each village and the base values of input variables. These base values were used when running EmiStatR in the deterministic mode (see Sect. 3.1). Some of the variables were calibrated based on observations in the CSOC to simulate water level and concentrations and loads of pollutants spilled in the CSO to the stream, river, or lake. These variables are water consumption (qs), infiltration flow (q_{f}), time flow (t_{fS}), runoff coefficient for the impervious area (C_{imp}), runoff coefficient for the pervious area (C_{per}), orifice coefficient of discharge (C_{d}), and initial water level (Lev_{ini}).
2.4 Selection of model input for uncertainty quantification
Following recommendations from Nol et al. (2010), not all model inputs were taken into account in the uncertainty propagation analysis. Only inputs that are very uncertain and to which the model output is very sensitive were included because these are the ones that have the largest contribution to output uncertainty (Heuvelink, 1998; Sect. 4.4). The level of uncertainty of the inputs was defined by expert judgement and similar case studies in the literature. A quick scan was used to determine the model sensitivity to each of the model inputs by running EmiStatR in deterministic mode with input base values given in Table S1. The level of model sensitivity was defined by analysing the mathematical model structure and components of the model, expert judgement, and simulations with EmiStatR. Inputs that rank high on both the level of uncertainty and on model sensitivity were selected and included in the uncertainty propagation analysis.
2.5 Uncertainty quantification of selected model input
Because we used a statistical approach, probability distribution functions (pdf's) are the basis for representing the uncertainties of the selected model inputs. This constitutes the most difficult step of an uncertainty propagation analysis and is done in different ways for constants and dynamic variables, as explained in the following subsections.
2.5.1 Uncertain constants
Following Heuvelink et al. (2007), an uncertain continuous numerical constant C can be characterised by its marginal (cumulative) pdf (mpdf) as follows:
Usually, a parametric approach can be taken, meaning that a common shape for F_{C} is chosen (e.g., normal, lognormal, exponential, or uniform) so that the mpdf is reduced to a number of parameters. In this study, the input variables that are in this category are water consumption (qs), infiltration inflow (q_{f}), total area (A_{total}), impervious area (A_{imp}), the runoff coefficients for impervious area (C_{imp}) and pervious area (C_{per}), population equivalents (pe), flow time structure (t_{fS}), and initial water level (Lev_{ini}).
2.5.2 Univariate autoregressive modelling
Dynamic uncertain inputs may be temporally autocorrelated. This may dramatically influence the outcome of an uncertainty propagation analysis and must therefore be accounted for. One way of doing this is by assuming an autoregressive order one (AR(1)) model as follows:
where y_{t} is the uncertain input at time t, μ is its mean, ϕ is the autoregressive parameter ($\mathrm{0}\le \mathit{\varphi}<\mathrm{1}$), and w_{t} is a Gaussian white noise time series with mean zero and variance ${\mathit{\sigma}}_{w}^{\mathrm{2}}$. The initial value y_{0} is taken from a normal distribution with mean μ and variance σ^{2}. The parameters of the model can be estimated based on observations, or in the absence of observations, suitable values are taken based on expert judgement or literature reference values. Note that the effect of the initial condition usually fades out quickly and hence is not of great concern.
The implementation of the AR(1) model in R was done via the R function arima.sim of the R base package stats (RCoreTeam and contributors worldwide, 2017), both for model calibration and simulation.
2.5.3 Multivariate autoregressive modelling
In the case of multiple uncertain dynamic inputs, crosscorrelation between these inputs may also need to be included. For example, C_{COD, S} and ${C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{S}}$ and their uncertainties are likely correlated. This can be done using a multivariate AR(1) model (Luetkepohl, 2005), which is a natural extension of the univariate AR(1) model, as follows:
where Y(t) is a vector of inputs at time t, A is a square matrix with parameters that define how the variables at time t+1 depend on those at time t, μ is now a vector of means, and ε(t) is a vector of zeromean normally distributed white noise processes. We further assume that the variance–covariance matrix C of ε(t) is time invariant. The initial value Y_{0} is assumed normally distributed and uncorrelated (Λ is a diagonal matrix). In order to estimate the vector μ and matrices A and C, a sample of the variables of interest is needed. Parameter estimation is done by means of the R package mAr (Barbosa, 2015).
2.5.4 Input precipitation model
In case precipitation is selected as an uncertain input to be included in the uncertainty analysis, it too must be characterised by a pdf. Since precipitation, however, is not normally distributed and has many zeros, we cannot make the Gaussian assumption, and hence we cannot use the approach described in Sect. 2.5.2 to model its dynamic behaviour and uncertainty. In addition, we usually have precipitation measurements nearby so we need to condition the simulations to these measurements. Recall from Sect. 2.3 that, in the case study, precipitation data are recorded at stations EschsurSûre and Dahl.
TorresMatallana et al. (2017) present a model to simulate precipitation inside a target catchment given a known precipitation time series in a nearby location outside the catchment, while accounting for the uncertainty that is introduced due to spatial variation in precipitation. The method used for input precipitation uncertainty characterisation is essentially the same as the application of a Kalman filter/smoother (Kalman, 1960; Webster and Heuvelink, 2006). Calibration of the model requires precipitation time series at two locations near the catchment of interest. Once the model is calibrated, it is used to simulate precipitation inside the target catchment from a single precipitation time series nearby the catchment. Details of the calibration and conditional simulation are presented in Sect. S4.
2.6 Uncertainty analysis
We used a MC simulation (Hammersley and Handscomb, 1964; Kalos and Whitlock, 2008) to analyse how input uncertainty propagates through the EmiStatR model because it is flexible and straightforward to implement. It is also feasible in our case study because EmiStatR is a relatively simple model that does not involve a long computation time.
2.6.1 Monte Carlo simulation
The MC method runs the EmiStatR model repeatedly, each time using different model input values sampled from their pdf. The method thus consists of the following steps:

Repeat n times:

Generate a set of realisations of the uncertain model inputs at 10 min resolution.

Run the model at 10 min resolution and store output for this set of realisations. Later, in order to compute the summary statistics, a temporal aggregation of the model output to 1 h intervals is done.


Compute and store sample statistics from the n model outputs.
Here, n is the number of MC runs, i.e. the MC sample size. Common sample statistics that measure the uncertainty are the standard deviation and quantiles of the distribution of MC outputs, such as the difference between the 0.95 and 0.05 quantile, which can be easily calculated from the n MC outputs.
Sampling from the pdf of uncertain inputs was done using simple random sampling.
2.6.2 Monte Carlo output summary
Proper presentation of MC outputs is important for achieving the most from the experiment. Therefore, summary statistics are one important way to summarise the MC outputs. Commonly, a MC study yields n model outputs, which are stored in the MC result matrix X in Boos and Osborne (2015). From this matrix, various statistics can be computed. Basic summary statistics include the mean μ_{MC}, the standard deviation (σ_{MC}), and the variance ${\mathit{\sigma}}_{\mathrm{MC}}^{\mathrm{2}}$. From these, we can compute the coefficient of variation CV_{MC} (σ_{MC}∕μ_{MC}), which is a dimensionless expression of relative uncertainty. The coefficient of variation is a standardised measure of the spread of a sampling distribution, which is useful because it allows a direct comparison of the variation in samples with different units or with very different means (Marwick and Krishnamoorthy, 2019). We computed estimates and standard errors for these statistics and also for the interquartile range (IQR_{MC}), 0.005 (ζ_{0.005}) and 0.995 (ζ_{0.995}) quantiles, and the 99 % width of the prediction band (ζ_{w, 0.99}).
2.6.3 Bootstrap computation for Monte Carlo summary
Following Boos and Osborne (2015), who argue that “Good statistical practice dictates that summaries in MC studies should always be accompanied by standard errors,” we used the bootstrap method to compute the standard errors of all MC statistics. These tools are particularly relevant in a case without analytic solutions (Boos, 2003). According to Boos and Osborne (2015, p. 228) standard errors for MC output statistics are often not computed, which is an additional computational step on top of the overall analysis. Standard errors are straightforward to compute for simple statistics, such as the sample mean over the replications of the MC output, but are more difficult to compute for more complex statistics, such as medians, sample variances, and the classical Pearson measures of skewness and kurtosis. Therefore, to avoid burdensome computations, we opted to compute the standard errors by the bootstrap method. We briefly explain the bootstrap method below. For a detailed explanation, we refer to Efron (1979).
To compute the bootstrap variance of estimators, we follow the logic given by Boos and Osborne (2015). From a MC sample ${Y}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots}\phantom{\rule{0.125em}{0ex}},{Y}_{n}$, we draw a random sample of size n, with replacement ${Y}_{\mathrm{1}}^{*},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots}\phantom{\rule{0.125em}{0ex}},\phantom{\rule{0.125em}{0ex}}{Y}_{n}^{*}$, and compute an estimator $\widehat{\mathit{\theta}}$ of the MC statistic θ from this resample. We independently repeat this process B times, resulting in a sample of estimators ${\widehat{\mathit{\theta}}}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots}\phantom{\rule{0.125em}{0ex}},{\widehat{\mathit{\theta}}}_{B}$. Then the bootstrap variance estimate, ${\widehat{V}}_{B}$, is the sample variance of this sample of estimators as follows:
where $\overline{\widehat{\mathit{\theta}}}$ is the mean of the sample of estimators. The MC standard error, se, is simply the square root of the bootstrap variance.
We implemented, in stUPscales (TorresMatallana et al., 2019), specific routines for computing, by means of the bootstrap method, the MC estimators and their standard error for all MC statistics, where the variance of the model output is the most important. We compared our results with the results obtained using the Monte.Carlo.se R package (Boos et al., 2019).
2.6.4 Contributions of input variables to total uncertainty
A number of m+1 MC analyses are needed to compute the contributions of input variables to total uncertainty, where m is the number of model input variables selected for uncertainty quantification. The first MC analysis, MC_{tot}, is done to compute the total output uncertainty by stochastically varying all input variables. The uncertainty associated with the first variable x_{1} is quantified by a second MC analysis MC_{1} in which only x_{1} is equal to its deterministic value, while the other input variables vary stochastically. Similarly, the other MC simulations MC_{2}, MC_{3}, … , and MC_{m} are used to quantify the uncertainty for the variables x_{2}, x_{3}, … , and x_{m}.
To quantify the contributions of individual input variables to the total uncertainty of the model inputs, the stochastic sensitivity S_{i} for each uncertain input x_{i} is computed. It is defined as follows:
The index represents the main effect contribution of each input factor to the total variance of the output. The larger the index, the more important the input uncertainty. We computed stochastic sensitivity indices per time step and aggregated contributions for the whole year. For plotting purpose, we aggregated the outputs from a 10 min time step to hourly time steps. The aggregation was done for each individual MC run before the contributions were computed.
2.7 Water quality impact
The results of the MC uncertainty propagation were also compared with the water quality standards. Standards are introduced to evaluate the impact of emissions of COD and NH_{4} in CSOs in the receiving water. However, as Toffol (2006) recognises, although there are European emission standards for wastewater treatment plant effluent, standards for combined sewer overflow are not so clear. According to Steinel and Margane (2011), the European Water Framework Directive (WFD) is mainly concerned with the natural state of waters. Therefore, emission standards for effluent discharge are not set. The European Union Council Directive 91/271/EEC (1991) sets standards for COD and total nitrogen; hence, similar values have been adopted in many European member states. For more details about guidelines and design procedures in Europe, see Blumensaat et al. (2012). We assessed the emissions according to the German guideline ATVDVWKA 128 (1992), which is the standard for the dimensioning and design of storm water structures in combined sewers and commonly used in Luxembourg. The Austrian ÖWAVRB 19 (2007) is also taken into account because it provides key reference guidelines for the design of urban water infrastructure in central Europe. Three main indicators are taken into account, namely hydraulic impact, COD concentration, and acute ammonium toxicity.
2.7.1 Hydraulic impact
According to the Austrian guidelines, and as summarised by Kleidorfer and Rauch (2011), the evaluation of the hydraulic impact is given by the following:
where $\mathrm{0.1}\le {f}_{\mathrm{h}}\le \mathrm{0.5}$, Q_{1} (L s^{−1}) is the maximum sewer overflow discharge with return period of 1 year, and Q_{r1} (L s^{−1}) is the maximum water discharge in the river with return period of once per year. The factor f_{h} is taken as 0.1 in more sensitive streams, whereas it is 0.5 for streams with more stable bed and higher recolonisation potential Toffol (2006). Time series of daily values recorded in 2006 to 2013 of the river Sûre at Heiderscheidergrund were used to compute the daily flow expected, with a return period of once per year (1.01 years), Q_{r1}.
We used the German guideline ATVDVWKA 128 (1992), which computes the throttle discharge at CSOs, Q_{t, CSO} (L s^{−1}) as follows:
where $\mathrm{7.5}\le {f}_{t}\le \mathrm{15}$ and A_{imp} (ha) are the impervious areas connected to the combined sewer system. For the overflow flow MC output mean, the 0.95 and 0.995 quantile, we computed the exceedance percentage over the thresholds, which is calculated as the proportion of time steps exceeding the number of total time steps in the year (8759 time steps at 1 h).
2.7.2 COD concentration
Steinel and Margane (2011, Table 14) present the effluent standards for discharging into freshwater adopted in selected European countries. A COD concentration of 125 mg L^{−1} is reported for European Union countries. Austria has stricter rules, with a standard of 90 mg L^{−1} for populations between 50 and 500 inhabitants and 75 mg L^{−1} for populations greater than 500 inhabitants. The Goesdorf population was 1025 inhabitants by 2001 and was 1297 inhabitants by 2011 (Statec, 2020). For NH_{4} a similar approach was used.
2.7.3 Acute ammonium toxicity
Following Kleidorfer and Rauch (2011), who argue that “the ammonia (NH_{3}) concentration depends on the ammonium (NH_{4}) concentration and on the dissociation equilibrium between NH_{3} and NH_{4} (which is influenced by temperature and pHvalue)”. According to Kleidorfer and Rauch (2011), the Austrian guideline ÖWAVRB 19 (2007) establishes a maximum value of 2.5 mg L^{−1} for the ammonium (NH_{4}) concentration calculated for a 1 h duration for salmonid streams. For cyprinid streams, a maximum value of 5.0 mg L^{−1} is recommended.
3.1 Selection of model inputs for uncertainty quantification
In this section, we assess the degree of uncertainty and sensitivity for all input variables, following the procedure described in Sect. 2.4. We summarise the results in Tables 1 and 2. To better support our decisions, we also include a graphical assessment of the degree of uncertainty and sensitivity of each input, as in TscheiknerGratl et al. (2017, see Fig. S1).
3.1.1 Wastewater
Water consumption, qs, is a fairly uncertain input variable, and the model output is sensitive to this variable. Volume and flow of CSO are sensitive to changes in qs. Regarding water quality output, the total load of NH_{4} is very sensitive to changes in qs. Pollution of sewage as COD load per capita per day, C_{COD, S}, is the first selected input variable for propagation of uncertainty, due to the fact that it is both a very uncertain input variable and the model output (average and 99.9 percentile overflow COD concentration) is very sensitive to it. Pollution of sewage as NH_{4} load per capita per day, ${C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{S}}$, is also included in the uncertainty propagation analysis. It is a very uncertain input variable, and the model output (overflow load and concentrations of NH_{4}) is very sensitive to it. The variables C_{COD, S} and ${C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{S}}$ are very uncertain because these are correlated to the temporal and spatial pattern of water consumption, which has a daily, weekly, and seasonal temporal variability.
3.1.2 Infiltration water
Inflow of infiltration water, q_{f}, is a very uncertain input variable because this inflow depends on the number of anomalies in the pipes (cracks or wrong connections) that allows infiltrations flowing into and out of the system. The distribution of these anomalies has a strong random component, and hence, q_{f} is very uncertain, and the model output is sensitive to it. Although this is a very uncertain input, the quickscan analysis showed that model output sensitivity is not very high, as indicated in Table 1. For this reason, we did not include this variable in the uncertainty propagation analysis.
Pollution of infiltration water as COD load per capita per day, COD_{f}, and pollution of infiltration water as NH_{4} load per capita per day, ${{\mathrm{NH}}_{\mathrm{4}}}_{\mathrm{f}}$, are not uncertain because in the HauteSûre study area the values of these variables are negligible.
3.1.3 Rainwater
Precipitation, P, is the main driving force of the model, and given the spatial variability of the rain fields, this input is considered very uncertain. The model output, additionally, is very sensitive to it. As a consequence, this input variable is treated as the third input variable in the uncertainty propagation analysis. Pollution of runoff as COD concentration, COD_{r}, is the fourth input variable considered in the uncertainty propagation, given that it is a very uncertain and very sensitive input variable, particularly to the load and concentration of COD in the overflow.
Pollution in rainwater as NH_{4} concentration, ${{\mathrm{NH}}_{\mathrm{4}}}_{\mathrm{r}}$, is considered fairly uncertain. The model output (overflow load and concentration of NH_{4}) is very sensitive to it. Although model output is very sensitive to this model input variable, model input uncertainty is not very high, as indicated in Table 1. For this reason, it was not included in the uncertainty propagation analysis.
3.1.4 Subcatchment
The model is very sensitive to the total area A_{total} and to the runoff coefficient for the pervious area (C_{per}) and sensitive to the impervious area A_{imp} and to the runoff coefficient for impervious area (C_{imp}). However, we did not include A_{total} in the uncertainty analysis because it can be fairly accurately derived from spatial databases, and hence, their uncertainty is not large. Although model output is very sensitive to the input variable C_{per}, the uncertainty about this variable is not very high, as indicated in Table 1. The reason behind this is that C_{per} can be derived fairly accurately from geographic information system (GIS) products, such as land use and soil type maps. Therefore, we did not include this variable in the uncertainty propagation analysis. The population equivalents, pe, is a sensitive variable but not very uncertain. Hence, this variable was not included in the uncertainty analysis. The theoretical largest flow time in the catchment t_{fS} is not uncertain and not sensitive.
3.1.5 CSO structure
Although model output is very sensitive to maximum throttled outflow Q_{d, max} and volume V, these are not included in the uncertainty analysis because their values are accurately known. The same is true for the variables curve level–volume lev2vol, orifice diameter D_{d}, and discharge coefficient C_{d}. These variables are accurately known and, therefore, not considered as being uncertain variables. The initial water level in the chamber Lev_{ini} is very uncertain, but the model output is not sensitive to this variable. Therefore, Lev_{ini} was not included in the uncertainty analysis.
3.2 Uncertainty quantification of selected model input
After ranking all inputs on levels of uncertainty and model sensitivity, we selected four input variables to be included in the uncertainty analysis. These are C_{COD, S}, ${C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{S}}$, COD_{r}, and P (Table 2).
3.2.1 Sewage per capita COD and Ammonium
The fit of pdf's for the two uncertain inputs, C_{COD, S} and ${C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{S}}$, was based on measurements under dry weather flow conditions. Measurement campaigns were done in Goesdorf from 28 April to 24 June 2011, in Kaundorf from 22 June to 18 August in 2010 and from 20 July to 5 August in 2011, and in NocherRoute from 18 November 2010 to 27 April 2011. Samples of COD and NH_{4} in mg L^{−1} (91 in total for each variable) were analysed. An average wastewater amount was calculated for Goesdorf (153 L PE^{−1} d^{−1}), Kaundorf (112 L PE^{−1} d^{−1}), and NocherRoute (94.3 L PE^{−1} d^{−1}). Table 3 presents summary statistics of the dry weather flow measurements of COD and NH_{4} and the corresponding value of C_{COD, S} and ${C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{S}}$. COD is converted to C_{COD, S} by means of a simple conversion from mg L^{−1} to g PE^{−1} d^{−1}, by multiplying COD by the measured per capita flow (112 L PE^{−1} d^{−1}) and dividing by 1000. NH_{4} was converted to ${C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{S}}$ in a similar way.
Closer inspection showed that C_{COD, S} and ${C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{S}}$ observations are best characterised by a lognormal distribution (Sect. S5). Since C_{COD, S} and ${C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{S}}$ are dynamic and crosscorrelated, we calibrated a bivariate AR(1) model with state vector $\mathit{Y}=\left[\text{log}\right({C}_{\mathrm{COD},\phantom{\rule{0.125em}{0ex}}\mathrm{S}})$ $\text{log}\left({C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{S}}\right){]}^{T}$. The estimated parameters of the model, using the methodology described in Sect. 2.5.3, are as follows:
The defined multivariate autoregressive model also captures the dynamic behaviour, temporal correlation, and crosscorrelation of the input variables, deriving the probability distributions of C_{COD, S} and ${C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{S}}$ from measurements in the HauteSûre catchment, which agreed well with values reported in the literature (Katukiza et al., 2014; Heip et al., 1997).
3.2.2 Runoff COD concentration
Regarding COD_{r}, due to the fact that no field measurements were available, expert judgement and reference values from the literature were the basis for characterising the pdf of this input variable. The variable was assumed to be lognormally distributed with a mean value of 71 mg L^{−1}. Although, House et al. (1993) and Welker (2008) reported a higher value, namely 107 mg L^{−1} for COD_{r}, we selected a lower value due to the specific characteristics of the CSO system in the HauteSûre catchment. The value of 150 mg L^{−1} as standard deviation of COD_{r} leads to a coefficient of variation (SD ⋅ mean^{−1}) equal to 2.11, which is greater than the coefficient of variation for C_{COD, S} (0.84). We allow the standard deviation of COD_{r} to be greater than the standard deviation of C_{COD, S} because COD measurements in rain water are very uncertain.
3.2.3 Input precipitation model
Precipitation and its associated uncertainty was modelled as an autoregressive model conditioned to the observed precipitation at a nearby measurement station. We assumed a multivariate lognormal distribution and included a temporal correlation of the simulated time series. Calibration of the precipitation model is done with the mAr package, as explained in Sect 2.5.4, and using a 10 min precipitation time series of stations EschsurSûre and Dahl for 2010. Upon calibration of the multivariate autoregressive model, we proceeded with the conditional simulation of Y_{c} (Sect. S4.2, Eq. 5). For this, we computed the parameters of the model, as shown in Eq. (14). The model parameters are given by (TorresMatallana et al., 2017) the following:
Next, we generated conditional simulations of the 10 min precipitation for 2010 for each subcatchment using the approach described in Sect. 2.5.4. Note that this involves simulating logtransformed precipitation, which can easily be transformed to precipitation data using the antilog. The simulation procedure was repeated as many times as simulated precipitation time series were required for the MC uncertainty propagation analysis.
The simulated precipitation time series captured the main statistics of the observed time series well. The reader can find evidence for this in the Supplement (Table S4 and Fig. S2). Despite the satisfactory performance of the proposed method, some cases showed an overestimation of the simulated precipitation, mainly due to high values of the ratio of the multiplicative factor δ(t). This behaviour was also recognised by McMillan et al. (2011), who stated that the multiplicative factor used in their study “does not capture the distribution tails, especially during heavy precipitation where input errors would have important consequences for runoff prediction”.
3.3 Uncertainty analysis
Model output sensitivity and the degree of uncertainties evaluation of each model input helped to define the four input variables included in the uncertainty analysis, namely C_{COD, S}, ${C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{S}}$, COD_{r}, and P. In this section, we present the results of the uncertainty propagation for these four selected input variables to the model output, both for water quantity (volume in the combined sewer overflow chamber, CSOC, and overflow volume and flow) and for water quality (loads and concentrations of chemical oxygen demand, COD, and ammonium, NH_{4}).
3.3.1 Monte Carlo output and uncertainty quantification
The seven output variables from EmiStatR were analysed by MC input uncertainty propagation. A detailed description of the Monte Carlo simulation size and timing is presented in the Supplement (Sect. S6).
Figure 3 illustrates the uncertainty propagation outcomes for the first Monte Carlo simulation, in which all input variables vary stochastically. The MC simulations were performed for the entire year of 2010 at 10 min time steps, which were aggregated to hourly time steps in the figure. The aggregation function was used for precipitation, CSO chamber volume, CSO spill volume, and loads was the sum, whereas for CSO spill flow and concentrations the aggregation function was the mean. The top of Fig. 3 shows input precipitation as main driving input. For illustration purposes, two events with a 2 d duration each are shown. The first event occurred in spring (May 2010) and the second in autumn (September 2010), and this shows that the uncertainty is high when there is a high precipitation event. The more intense the precipitation input, as seen in Fig. 3 in the top left (May event), the greater the uncertainty bandwidth for overflow flow and for COD and NH_{4} loads (Fig. 3e–h) and concentrations (Fig. 3i–l). The MCestimated statistics and the standard errors (se) are presented in Table 4. The table shows the uncertainty quantification of outputs obtained from the MC uncertainty propagation for the first MC simulation (all selected input variables are uncertain).
Table 4 shows the standard deviation (SD) and the coefficient of variation (cv) for the seven output variables considered in the uncertainty propagation. For the volume in the CSO chamber, V_{Chamber}, the annual mean standard deviation, σ_{MC}, (8.27 m^{3}) is lower than the mean, μ_{MC}, (92.51 m^{3}). This goes along with an annual mean coefficient of variation (CV_{MC}) of 0.100. A (CV_{MC}) greater than one means large uncertainty. The overflow spill volume, V_{Sv}, had a coefficient of variation of 0.070, while it was 1.585 for the overflow flow, Q_{Sv}. This shows that the relative uncertainty of the overflow flow is very large. Regarding the overflow COD load, the annual mean (1.18 kg) is similar as the annual mean standard deviation (1.27 kg). Similar behaviour was observed for the overflow COD concentration, which had an annual mean value of 170 mg L^{−1} and a standard deviation of 162 mg L^{−1}. For overflow NH_{4} load and overflow NH_{4} concentration, the annual mean also had the same order of magnitude as the annual mean standard deviation. Overflow COD and NH_{4} loads had a coefficient of variation of 0.087 and 0.075, respectively, whereas the coefficient of variation for concentrations were 0.988 and 0.815, respectively. This suggests that overflow concentrations are more uncertain.
Low standard errors (se) for the coefficient of variation were obtained for all seven outputs. These were never greater than 0.05, which indicates that the selected MC replication size (1500 for mc1) is a suitable value. This holds for all output statistics because, in all cases, the standard error is small in comparison to the estimated value.
3.3.2 Contributions of input variables to total uncertainty
The contributions of input variables to the total uncertainty of the model inputs were also computed using the procedure described in Sect. 2.6.4. A total of four MC simulations with a total of 6000 runs were performed to estimate S_{i} (Eq. 10). Afterwards, four contributions were evaluated per time step and aggregated for the whole year. Following Eq. (10), the per time step contribution of input variables to output variables in terms of percentage of variance, stochastic sensitivity S_{i} of the input variables C_{COD, S}, ${C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{S}}$, COD_{r}, and P were calculated. An example of the contributions analysis per time step is presented in Fig. 4. Here we remark that a high uncertainty over time is shown mainly for the spring event.
The aggregatedovertime contributions of input variables to output variables in terms of percentage of variance and stochastic sensitivity S_{i} of the input variables were also calculated (Table 5). Note that P is the only source of uncertainty for V_{Chamber} and V_{Sv}, while uncertainty in NH_{4} inputs only propagates to NH_{4} outputs, which is similar for COD (Fig. 4).
We found, as expected, that precipitation, P, is the only source of uncertainty from all uncertain input considered for water quantity output variables V_{Chamber} and V_{Sv}. Regarding average values for the whole year, for the water quality output variables B_{COD, Sv} and ${C}_{\mathrm{COD},\phantom{\rule{0.125em}{0ex}}\mathrm{Sv},\phantom{\rule{0.125em}{0ex}}\mathrm{av}}$, C_{COD, S} has the largest contribution to the output variance, which is about 66 % for B_{COD, Sv} and about 62 % for ${C}_{\mathrm{COD},\phantom{\rule{0.125em}{0ex}}\mathrm{Sv},\phantom{\rule{0.125em}{0ex}}\mathrm{av}}$. The second variable that contributes to the uncertainty of these COD output variables is P, with about 3 % for B_{COD, Sv} and 9 % for ${C}_{\mathrm{COD},\phantom{\rule{0.125em}{0ex}}\mathrm{Sv},\phantom{\rule{0.125em}{0ex}}\mathrm{av}}$. Similarly, the input variable ${C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{S}}$ plays an important role in the contribution of total uncertainty for ${B}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{Sv}}$ (on average about 34 % of the variance for the whole year) and ${C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{Sv},\phantom{\rule{0.125em}{0ex}}\mathrm{av}}$ (about 35 %). Equally contributing to the uncertainty of these NH_{4} output variables is P with about 66 % for ${B}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{Sv}}$ and 65 % for ${C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{Sv},\phantom{\rule{0.125em}{0ex}}\mathrm{av}}$. From these results, we can infer that precipitation is a main source of uncertainty for all six outputs considered.
3.4 Uncertainty and water quality impact
Quantification and assessment of the water quality impact is an important step after the uncertainty propagation. As described in Sect. 2.7, the assessment of water quality standards was done by taking into account the reference thresholds recommended in the European Union guidelines for COD and the German and Austrian guidelines for hydraulic impact and acute ammonium toxicity.
3.4.1 Hydraulic impact
From the time series of daily values for 2006 to 2013 of the river Sûre, a daily flow expected with a return period of once per year (1.01 years), Q_{r1} of 16 m^{3} s^{−1} was computed at Heiderscheidergrund, which corresponds with the entire catchment area of the HauteSûre storm water system (182.1 ha). Therefore, we estimated the river daily flow in the Goesdorf CSO structure as a proportion to 30 ha, which is equal to 2.6 m^{3} s^{−1}. Following Eq. (11), the maximum sewer overflow discharge with return period of 1 year, Q_{1} can have a value between 0.26 and 1.32 m^{3} s^{−1}. Accordingly, with the German guideline ATVDVWKA 128 (1992; Eq. 12), two additional thresholds are defined for the maximum sewer overflow discharge with return period of 1 year for the Goesdorf catchment (A_{imp}=5.0 ha). Q_{1} is expected to vary between 37.5 and 75.0 L s^{−1}. We contrasted these values with those obtained from the uncertainty analysis. From Table 4, we obtained a 1 h mean value for the overflow spill flow, Q_{Sv}, of 55.5 L s^{−1}, 90 % prediction band width of 164.1 L s^{−1}, and standard deviation of 64.5 L s^{−1}. Figure 3a and b present the overflow spill flow for the two periods chosen for illustration. The upper dashed red line indicates the 75 L s^{−1} threshold, and the lower dotted red line indicates the 37.5 L s^{−1} threshold. Table 6 (top) shows the exceedance percentage of overflow spill flow over the 37.5 and 75.0 L s^{−1} thresholds for the mean, 0.95 quantile, and 0.995 quantile. We found a 0.49 % exceedance of the mean value over the 37.5 L s^{−1} threshold and about 1.7 % for the quantiles. As expected, slightly lower percentages were found for the 75.0 L s^{−1} threshold.
3.4.2 COD concentration
A reference COD concentration emission in CSOs was presented in Sect. 2.7.2. For the European Union, a value of 125 mg L^{−1} is used. We obtained a 1 h average spill COD concentration with a mean of 170 mg L^{−1}, standard deviation of 162 mg L^{−1}, and a 90 % prediction band width of 450 mg L^{−1}. Figure 3i and j present the average spill COD concentration. The upper dashed red line indicates the 125 mg L^{−1} threshold and the lower dotted red line the 90 mg L^{−1} threshold. The mean COD concentration in the overflow volume was higher than the thresholds. However, note that when entering the river system it will quickly be diluted, suggesting that the negative impact on the environment will be dampened by the receiving water body.
Table 6 (centre row) shows the exceedance percentage of overflow COD concentration over the 90 and 125 mg L^{−1} thresholds for the mean, 0.95 quantile, and 0.995 quantile. We found a 1.62 % exceedance of the mean value over the 90 mg L^{−1} threshold and about 1.8 % for the quantiles. Slightly lower percentages were found for the 125 mg L^{−1} threshold for the mean value (1.03 %). For the quantiles equal values were found as for the 90 mg L^{−1} threshold.
3.4.3 Acute ammonium toxicity
We compared the acute ammonium toxicity reference values presented in Sect. 2.7.3 (2.5 mg L^{−1} for the ammonium concentration calculated for a 1 h duration for salmonid streams, and for cyprinid streams, a maximum value of 5.0 mg L^{−1}) was calculated with the values we found for ammonium. An average spill NH_{4} concentration with a mean of 7.19 mg L^{−1}, standard deviation of 5.66 mg L^{−1}, and 90 % prediction band width of 16.64 mg L^{−1} was obtained. Figure 3k and l show the average spill NH_{4} concentration for the two periods chosen for illustration. The ammonium (NH_{4}) concentrations in the overflow flow are higher than the reference values, which are given for concentrations in the river.
Table 6 (bottom row) shows the exceedance percentage of overflow NH_{4} concentration over the 2.5 and 5.0 mg L^{−1} thresholds for the mean, 0.95 quantile, and 0.995 quantile. We found a 1.8 % exceedance of the mean and quantile values over the 2.5 and 5.0 mg L^{−1} thresholds. A slightly lower percentage (1.1 %) was found for the 5.0 mg L^{−1} threshold with regards to mean value.
This study aimed to select and characterise the main sources of input uncertainty in urban water systems, while accounting for temporal auto and crosscorrelation of uncertain model inputs, by propagating input uncertainty through the EmiStatR model and quantifying and assessing the contributions of each uncertainty source to model output uncertainty dynamically (over time). In the following discussion, we start with the uncertainty and water quality impact of the model outputs to the environment, in relation to the uncertainty analysis. Next, we discuss the accuracy of Monte Carlo analysis, followed by a discussion of other sources of uncertainty. Finally, we highlight some limitations and possible solutions to the approach used in this work.
4.1 Uncertainty and water quality impact
Next, we discuss how the uncertainty propagation analysis done gives additional insight regarding hydraulics, COD concentration, and the acute ammonium toxicity impact on water quality over the river Sûre due to the CSO discharges under study. After doing the uncertainty propagation analysis, we not only have the predictions of model outputs, but we also know how uncertain these are. An added value arises when we take into account the uncertainty information. For the case of the overflow spill flow, the expected model output (mean of 55.5 L s^{−1}) is below the environmental threshold of 75 L s^{−1}, but the 0.95 quantile (164.1 L s^{−1}) is very much above the threshold. This indicates that there is a considerable chance of it being above the threshold.
Regarding water quality outputs, although the mean model output for COD and NH_{4} concentrations is fairly above that of the thresholds, the 0.95 quantile is 2.7 times above the mean value for COD concentration and 2.4 times above the mean value for NH_{4}. Also, here we can conclude that we are not certain that we are below the threshold because there is a considerable probability that the true values are above, even though the expected value is below the threshold.
We were able to compute the water quantity and quality at the CSO outlet to the river. We found that water quality (COD and NH_{4}) were sometimes above the environmental threshold. Even if the expected value was below the threshold, there could still be a considerable probability that the quality was above the threshold because of the large uncertainty. Therefore, policy and decisionmakers and water managers need to be aware of this because, whenever concentrations are above the threshold, this may harm the environment. Nevertheless it is worth noting that we computed the concentration in the outlet of the CSO. When this spilled water enters the river, it will quickly mix with the much cleaner river water and concentrations will drop quickly, so it is only a local problem. How local it is and how the river water quality is distributed in space and time is not an easy problem to solve and requires the use of hydrological and hydraulic river models, e.g. SIMBA (IFAK, 2007) or MIKE 11 (DHI, 2017). Those models have been well developed and, for some of them, uncertainty analyses have also been done (Beven and Binley, 1992; Refsgaard, 1997; Beven and Freer, 2001; Vrugt et al., 2003a, 2008; Beven et al., 2010; AndrésDoménech et al., 2010; Beven, 2012; JervesCobo et al., 2020; Yu et al., 2020), but obviously such uncertainty analyses can only be done if the inputs and the uncertainty associated with these inputs to these models are known. One of these inputs is the inlet from CSO. That is where our paper makes a very valuable contribution because our work has quantified water quantity and quality of CSO structures, including uncertainty, and that is exactly what these river models need to be able to do an uncertainty propagation analysis. We also recognise other attempts at quantity (e.g. Sriwastava et al., 2018) and quality, especially with measurements taken at CSOs, which demonstrate that the measured water quality at the WWTP influent is expected to render a low representativity of the conditions at the CSOs (e.g. Brombach et al., 2005; DiazFierros T et al., 2002). We present some comparisons with these studies in the following lines.
Sriwastava et al. (2018) apply an uncertainty propagation to a complex hydrodynamic model to quantify uncertainty in sewer overflow volume. They used MC for the uncertainty propagation and Latin hypercube sampling (LHS) as an efficient sampling scheme. Although LHS ensures a full coverage of the sample space and provides a faster convergence than simple random sampling, the LHS application in the case of dynamic model inputs (e.g. precipitation, COD, and NH_{4} inputs) is not trivial, and its implementation is more complex than in the case of sampling from static variables (i.e. uncertain constants). In our study, we sampled the time series of dynamic inputs using an implementation in stUPscales (TorresMatallana et al., 2019, 2018a).
DiazFierros T et al. (2002), in a study in the city of Santiago de Compostela (northwestern Spain; population about 100 000 inhabitants), where a combined sewer system feeds to a grossly undersized wastewater treatment plant, reported an event mean concentration (DiazFierros T et al., 2002; Table 4) for the output variables ${C}_{\mathrm{COD},\phantom{\rule{0.125em}{0ex}}\mathrm{Sv},\phantom{\rule{0.125em}{0ex}}\mathrm{av}}$ and ${C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{Sv},\phantom{\rule{0.125em}{0ex}}\mathrm{av}}$ of 329.1 and 8.7 mg L^{−1}, respectively. These values are larger than those found by Brombach et al. (2005) and are more in agreement with our findings, especially for the case of ${C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{Sv},\phantom{\rule{0.125em}{0ex}}\mathrm{av}}$. DiazFierros T et al. (2002) reported values of ${C}_{\mathrm{COD},\phantom{\rule{0.125em}{0ex}}\mathrm{Sv},\phantom{\rule{0.125em}{0ex}}\mathrm{av}}$ as high as 1073 mg L^{−1}, which agrees with the righthand tail of the distribution obtained in our study (i.e. a 0.995 quantile of 909.7 mg L^{−1}). Similarly, for the case of ${C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{Sv},\phantom{\rule{0.125em}{0ex}}\mathrm{av}}$, DiazFierros T et al. (2002) reported values as high as 32.5 mg L^{−1}, which is comparable with the 0.995 quantile (29.20 mg L^{−1}) found in our study.
It is worth noting that, regarding measurements taken at CSOs, the measured water quality at the WWTP influent is expected to render a low representativity of the conditions at the CSOs, as reported by DiazFierros T et al. (2002) and Brombach et al. (2005). Thus, when comparing model outputs with independent measurements, one should bear in mind that discrepancies between measured and predicted values are not only caused by errors in model inputs, model parameters, and model structure but are also the result of errors in the water quality measurements.
4.2 Accuracy of Monte Carlo analysis
Regarding the Monte Carlo replication size for uncertainty propagation, we presented the results in Fig. S4 for three output variables and three replications size 250, 1000, and the selected 1500 (Nash–Sutcliffe efficiency, NSE, closer to 1.0 for most of the output variables). We compute replications for 50, 100, and from 250 to 2000 at steps of 250 replications and the comparison of two equal MC runs (MC1 and MC2) with different seed for the pseudorandom number generator. The results suggest that the output variables related to COD (load and concentration) have a larger dispersion when we compare MC1 and MC2 for the same replications size. This is also reflected in the larger standard errors reported in Table 4 for, for example, the overflow COD load. Nevertheless, 1500 runs are a feasible MC replication size for running a relatively simple and fast model as EmiStatR (7.29 min in average execution time, using parallel computing and 50 cores for a time series with 4464 time steps). For a more complex full hydrodynamic model with a high computational burden, 1500 replications repeated four times to compute contributions may be not possible. Therefore, we suggest checking the intermediate results of the MC convergence test. We found that for quantity variables as the spill overflow volume and quality variables as the overflow NH_{4} load in 250 replications (7.10 min in average execution time using parallel computing and three cores for a time series with 4464 time steps) per individual MC execution was enough, which makes an execution of this kind of uncertainty propagation more feasible.
Figures 3 and 4 show that there is a large uncertainty for the early May event and smaller uncertainty for the September event. This is due mainly to the presence of a large dry period before the spill event in May, i.e. a shorter dry period preceding the spill flow leads to a lower uncertainty. This finding suggest also an interaction between the antecedent dry period and the concentration of pollutants.
4.3 Other sources of uncertainty
In this work, we only looked at input uncertainty and not at parameter and model structural uncertainty. Further research can be done on those topics. Neumann (2007) address how uncertainty ranges for parameters of fullscale systems are obtained and how model structure uncertainty manifests itself and can be quantified for a performance evaluation and the design of urban water infrastructure. MorenoRodenas et al. (2019) also studied and depicted how a model parameter is an important source of uncertainty. They emphasised that “still, uncertainty analysis is seldom applied in practice, and the relative contribution of the individual model elements is poorly understood.”. Also, they highlighted that, after inferring the river process parameters with system measurements of flow and dissolved oxygen, combined sewer overflow pollution loads became the dominant uncertainty source along with rainfall variability. These findings agreed with our results.
BachmannMachnik et al. (2018) recognised that the most important parameters causing uncertainties in the sewer system model are the connected area and the storm water runoff quality. Our analysis confirms these findings, specifically regarding the storm water runoff quality. In our study, the input variable runoff COD was an important source of uncertainty with relation to the annual mean overflow COD released from the CSO.
4.4 Limitations and possible improvements
Despite the extensive temporal uncertainty propagation analysis, the approach also has some limitations which we present hereafter, addressing possible solutions in future work.

Incorporation of the spatial distribution of model inputs. Specifically for precipitation,Breinholt et al. (2012) stated that, due to a poor representation of the spatial precipitation that is measured by point gauges and the complexity of the sewer systems, large output uncertainty can be expected. We also infer that we obtained a large output uncertainty due to neglecting the inherent spatial variability in precipitation. Therefore, we suggest that further research is needed to account for spatial variability in precipitation that can bring light an understanding of how this variability impacts the output uncertainty and thus quantify it properly. This issue should also be related to the problem of change of support. When modelling precipitation, we also ignored the support effect, i.e. we ignored that the subcatchment area is much greater than a point. Future research may address this issue of changeofsupport. Studies that tackled this issue are available, e.g. Leopold et al. (2006); Wadoux et al. (2017); Cecinati et al. (2018).

Linkage of submodels and uncertainty compensation effect. TscheiknerGratl et al. (2019) addressed the question as to whether there is an increase in uncertainty by linking integrated models or whether a compensation effect could take place by which overall uncertainty in key water quality parameters decreases. Some further insight into this topic could be obtained by quantifying uncertainties at submodel level and analysing whether the uncertainty at submodel level is greater or smaller than at the overall level. With our implementation, this is not a difficult task because EmiStatR has a stringent modular design in which it is easy to analyse outputs and their uncertainties at submodel level.

Accounting for crosscorrelation between the inputs precipitation and runoff COD concentration. It is worth noting that we did not include correlation between COD_{r} and P. Including such a correlation would yield a more realistic model of the uncertainty because these variables are known to have a strong correlation. It is highly recommended that correlations between COD_{r} and precipitation be included because loads in chemical oxygen demand are correlated with the overland flow due to precipitation, which may transport distributed pollutants to the sewer system. Also, the inputs C_{COD, S} and ${C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{S}}$ can be related with a daily curve that reflects the pattern of consumption in the household, like the German ATVDVWKA 134 curve. We used the latest version of EmiStatR (version 1.2.2.0), which considers these kinds of patterns.

Absence of high frequency water quality observations to compare with model outputs and uncertainty prediction bands. In order to gain an understanding of the temporal dynamics of nutrients (nitrogen, N, and phosphorus, P), Yu et al. (2020) applied highfrequency monitoring in a groundwaterfed, lowlying urban polder in Amsterdam (the Netherlands). They argued that although spatial and temporal concentration patterns from discrete sampling campaigns of water quality parameters, such as E. coli, showed a clear dilution pattern, the temporal patterns of N and P were still poorly understood, given their reactive nature and more complex biogeochemistry. Therefore, highfrequency measurement is a key factor in understanding these temporal dynamics and patterns.

Absence of a joint spatiotemporal uncertainty analysis. According to Zhou et al. (2020), the limitations in algorithms for classic uncertainty estimates is the reason that only the uncertainty in one dimension (either temporal variability or spatial heterogeneity) is considered, whereas the variation in the other dimension is dismissed, resulting in an incomplete assessment of the uncertainties. Zhou et al. (2020) also showed that classic metrics underestimate the uncertainty through averaging, which means a loss of information in the variation across spatiotemporal scales. To handle this limitation, suitable methods are the 3D variance partitioning for a new uncertainty estimation in both spatiotemporal scales (Zhou et al., 2020) or spatiotemporal geostatistics (Gräler et al., 2016).

Uncertainty analysis with complex models. In this research, we were able to conduct a comprehensive Monte Carlo uncertainty propagation analysis, which required a large number of Monte Carlo runs. This was possible because we used a strongly simplified urban water system model, EmiStatR. For more complex models that take much more computing time, the application of a Monte Carlo uncertainty propagation analysis is more challenging. However, given sufficient, resources it is possible because each model run can be run independently, and hence, the analysis is extremely suitable for parallelisation and cloud computing. In particular, the use of graphics processing units (GPUs) for heavy computation is promising. Some recent examples that demonstrate the potential of GPUs for this purpose are Eränen et al. (2014), Sten et al. (2016), and Sandric et al. (2019). Sriwastava et al. (2018) applied uncertainty propagation to a complex hydrodynamic model by selecting a small subset of dominant input/model parameters that explain most of the model output variance.

Convolution of precipitation to runoff. Equation (2) indicates that the runoff that enters the combined sewage in EmiStatR has the same shape as the precipitation, with an offset. Thus, for precipitation, no convolution is applied. In a new release of EmiStatR, we will incorporate convolution of precipitation when transforming precipitation to runoff.
In this final section, we conclude by highlighting the importance of temporal uncertainty propagation analysis and the selection and characterisation of uncertain model inputs impacting model sensitivity. We also point out that uncertainty propagation analysis helps to identify the sources that contribute the most and can provide better evidence for the impact assessment of pollutant release from sewer systems to the environment, in particular to the receiving waters.

Uncertainty analysis is important because it quantifies the accuracy of model outputs and quantifies the uncertainty source contributions. The latter provides essential information to help take informed decisions about how to improve the accuracy of the model output. But MC uncertainty analysis is only possible if it is computationally feasible. We used a simplified urban water system model with the capability to apply it for the minimisation of transient pollution from urban wastewater systems in parallel mode, which minimises model running time, allowing for uncertainty propagation, longterm simulations, and the evaluation of complex scenarios. These capabilities are crucial also for, for example, realtime control applications, where simplified models of fast running times are desirable.

Input variables that were very uncertain, and for which model output was very sensitive, were selected to be included in the uncertainty propagation analysis. We found four main input variables to be analysed, namely (1) precipitation, P, (2) chemical oxygen demand sewage pollution per capita load per day, C_{COD, S}, (3) ammonium pollution per capita load per day, ${C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{S}}$, and (4) chemical oxygen demand COD_{r} concentration.

Selected input variables for uncertainty propagation can be characterised in terms of the input uncertainty in four specific cases, depending on the type of input variable. (i) Uncertain constant inputs, characterised by their marginal (cumulative) pdf, e.g. water consumption, infiltration flow, impervious area and runoff coefficients; (ii) temporally autocorrelated dynamic uncertain inputs, characterised by univariate time series autoregressive modelling, e.g. COD_{r}; (iii) temporally crosscorrelated multiple dynamic uncertain inputs, characterised by multivariate time series modelling, considering cross and nocorrelations among variables, e.g. C_{COD, S} and ${C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{S}}$; and (iv) rain gauge input precipitation, characterised by autoregressive model conditioned to the observed precipitation (P).

Model input uncertainty propagation through the simplified, combined sewer overflow model (EmiStatR) helped us to understand how uncertainty propagates and how large the uncertainty of EmiStatR outputs is in a case study. Three output variables were considered for water quantity and four variables for water quality. The Monte Carlo uncertainty propagation analysis showed that, among the water quantity output variables, the overflow flow, Q_{Sv}, is the more uncertain output variable and has a large coefficient of variation (cv of 1.585). Among water quality variables, the annual average spill COD concentration, ${C}_{\mathrm{COD},\phantom{\rule{0.125em}{0ex}}\mathrm{Sv},\phantom{\rule{0.125em}{0ex}}\mathrm{av}}$, and the average spill NH_{4} concentration, ${C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{Sv},\phantom{\rule{0.125em}{0ex}}\mathrm{av}}$, were found to have large uncertainty (coefficients of variation of 0.988 and 0.815, respectively). Also, low standard errors (se) for the coefficient of variation were obtained for all seven outputs. They were never greater than 0.05, which indicated that the selected MC replication size (1500 simulations) was a suitable value.

Main sources of uncertainty model outputs. For water quantity outputs it was precipitation, while for COD, water quality outputs were P, C_{COD, S}, and COD_{r}, and for NH_{4}, outputs were P and ${C}_{{\mathrm{NH}}_{\mathrm{4}},\phantom{\rule{0.125em}{0ex}}\mathrm{S}}$.

Uncertainty propagation analysis can explain the impact of water quality indicators to the receiving river for the Luxembourg case study more comprehensively. Although the mean model water quality outputs for COD and NH_{4} concentrations are fairly above the thresholds, the 0,95 quantile is 2.7 times above the mean value for COD concentration and 2.4 times above the mean value for NH_{4}. We conclude that we are not certain that environmental thresholds are not exceeded because there is a considerable probability that values are above them, even though the expected value is below the thresholds. This is valid for concentrations in the spilled CSO; therefore, it is important to highlight that the results confirmed our hypothesis that annual mean COD and NH_{4} river concentrations are lower than the released CSO concentrations due to dilution and are thus compliant with the water quality thresholds given by the guidelines consulted.
The code scripts and data sets related to Figs. 3, 4, S3, and S4 are available from Zenodo (https://doi.org/10.5281/zenodo.3928079; TorresMatallana, 2020) and GitHub (https://github.com/ArturoTorres/temporal_uncertainty_paper_reproducible.git, last access: 29 December 2020).
The supplement related to this article is available online at: https://doi.org/10.5194/hess251932021supplement.
JATM is the main author of the text in this article. UL and GBMH contributed to this article with their statistical, geostatistical, and programming knowledge and reviewed and edited the text. JATM developed all R code scripts for the computations and Monte Carlo simulations and performed the simulations and analysis in collaboration with UL and GBMH.
The authors declare that they have no conflict of interest.
We thank the Luxembourgish Administration des Services Techniques de l'Agriculture (ASTA) for the rainfall time series and the Observatory for Climate and Environment (OCE) of LIST for the technical support. We thank Kai Klepiszewski for his advice and involvement during early stages of this research and the two anonymous reviewers for the constructive and valuable comments that helped us to improve this paper.
This research received funding from the European Union's Seventh Framework Programme for research, technological development, and demonstration (Quantifying Uncertainty in Integrated Catchment Studies, QUICS, project; grant no. 607000), and the Luxembourg Institute of Science and Technology (LIST).
This paper was edited by Nadav Peleg and reviewed by two anonymous referees.
AndrésDoménech, I., Múnera, J. C., Francés, F., and Marco, J. B.: Coupling urban eventbased and catchment continuous modelling for combined sewer overflow river impact assessment, Hydrol. Earth Syst. Sci., 14, 2057–2072, https://doi.org/10.5194/hess1420572010, 2010. a, b
Bach, P. M., Rauch, W., Mikkelsen, P. S., McCarthy, D. T., and Deletic, A.: A critical review of integrated urban water modelling – Urban drainage and beyond, Environ. Model. Softw., 54, 88–107, 2014. a
BachmannMachnik, A., Meyer, D., Waldhoff, A., Fuchs, S., and Dittmer, U.: Integrating retention soil filters into urban hydrologic models – Relevant processes and important parameters, J. Hydrol., 559, 442–453, https://doi.org/10.1016/j.jhydrol.2018.02.046, 2018. a, b, c, d
Baker, L. A. (Ed.): The Water Environment of Cities, Springer Science+Business Media, LLC, 233 Spring Street, New York, NY 10013, USA, 2009. a, b
Barbosa, S. M.: Package ”mAr”: Multivariate AutoRegressive analysis, The Comprehensive R Archive Network, CRAN, 1.12 edn., 2015. a
Bastin, L., Cornford, D., Jones, R., Heuvelink, G. B. M., Pebesma, E., Stasch, C., Nativi, S., Mazzetti, P., and Williams, M.: Managing uncertainty in integrated environmental modelling: The UncertWeb framework, Environ. Model. Softw., 39, 116–134, 2013. a
Beven, K. and Binley, A.: The future of distributed models: model calibration and uncertainty prediction, Hydrol. Process., 6, 279–98, 1992. a, b
Beven, K. and Freer, J.: Equifinality, data assimilation, and uncertainty estimation in mechanistic modelling of complex environmental systems using the GLUE methodology, J. Hydrol., 249, 11–29, https://doi.org/10.1016/S00221694(01)004218, 2001. a
Beven, K., Leedal, D., and Alcock, R.: Uncertainty and Good Practice in Hydrological Prediction, VATTEN, 66, 159–163, 2010. a
Beven, K. J.: RainfallRunoff Modelling: The Primer, WileyBlackwell, second edn., Lancaster University, UK, 2012. a
Blumensaat, F., Staufer, P., Heusch, S., Reussner, F., Schütze, M., Seiffert, S., Gruber, G., Zawilski, M., and Rieckermann, J.: Water qualitybased assessment of urban drainage impacts in Europe where do we stand today?, Water Sci. Technol., 66, 304–313, https://doi.org/10.2166/wst.2012.178, 2012. a
Boos, D., Matthew, K., and Osborne, J.: Monte.Carlo.se: Monte Carlo Standard Errors, available at: https://CRAN.Rproject.org/package=Monte.Carlo.se (last access: 29 December 2020), R package version 0.1.0, 2019. a
Boos, D. D.: Introduction to the Bootstrap and World, Stat. Sci., 18, 168–174, 2003. a
Boos, D. D. and Osborne, J. A.: Assessing Variability of Complex Descriptive Statistics in Monte Carlo Studies Using Resampling Methods, Int. Stat. Rev., 83, 228–238, https://doi.org/10.1111/insr.12087, 2015. a, b, c, d
Breinholt, A., Moller, J. K., Madsen, H., and Mikkelsen, P. S.: A formal statistical approach to representing uncertainty in rainfall–runoff modelling with focus on residual analysis and probabilistic output evaluation – Distinguishing simulation and prediction, J. Hydrol., 472–473, 36–52, 2012. a
Brombach, H., Weiss, G., and Fuchs, S.: A new database on urban runoff pollution: comparison of separate and combined sewer systems, Water Sci. Technol., 51, 119–128, https://doi.org/10.2166/wst.2005.0039, 2005. a, b, c
Brown, J. D.: Knowledge, uncertainty and physical geography: Towards the development of methodologies for questioning belief, T. I. Brit. Geogr., 29, 367–381, https://doi.org/10.1111/j.00202754.2004.00342.x, 2004. a
Cecinati, F., MorenoRódenas, A. M., RicoRamirez, M. A., ten Veldhuis, M. C., and Langeveld, J. G.: Considering rain gauge uncertainty using kriging for uncertain data, Atmosphere, 9, 1–17, https://doi.org/10.3390/atmos9110446, 2018. a
Committee on Environmental and Natural Resources – CENR: An Assessment of Coastal Hypoxia and Eutrophication in U.S. Waters, Tech. rep., NSTC, National Science and Technology Council, 725 17th Street, Washington, D.C., USA, 2003. a
Datta, A. R.: Evaluation of Implicit and Explicit Methods of Uncertainty Analysis on a Hydrological Modeling, Ph.D. thesis, University of Windsor, Canada, 2011. a
Deletic, A., Dotto, C., McCarthy, D., Kleidorfer, M., Freni, G., Mannina, G., Uhl, M., Henrichs, M., Fletcher, T., Rauch, W., BertrandKrajewski, J., and Tait, S.: Assessing uncertainties in urban drainage models, Phys. Chem. Earth., 4244, 3–10, https://doi.org/10.1016/j.pce.2011.04.007, 2012. a
DHI: MIKE11, A modeling system for rivers and channels, Reference Manual, DHI Water and Environment, Danish Hydraulic Institute, DHI, Hórsholm, Denmark, 2017. a
DiazFierros, T., F., Puerta, J., Suarez, J., and DiazFierros V., F.: Contaminant loads of CSOs at the wastewater treatment plant of a city in NW Spain, Urban Water, 4, 291–299, https://doi.org/10.1016/S14620758(02)000201, 2002. a, b, c, d, e, f
Efron, B.: Bootstrap methods: Another look at the Jackknife, Ann. Stat., 7, 1–26, 1979. a
Eränen, D., Oksanen, J., Westerholm, J., and Sarjakoski, T.: A full graphics processing unit implementation of uncertaintyaware drainage basin delineation, Comput. Geosci., 73, 48–60, https://doi.org/10.1016/j.cageo.2014.08.012, 2014. a
Evers, P., Heinz, H., Hanitsch, P. H., Koch, G., Naupold, L., Tochtermann, W., Tornow, M., Zander, B., Mahret, H., and Warnow, D.: ATVDVWKA 134E: Planning and Construction of Wastewater Pumping Stations, Tech. rep., DWA, Germany, 2000. a
Gasperi, J., Zgheib, S., Cladière, M., Rocher, V., Moilleron, R., and Chebbo, G.: Priority pollutants in urban stormwater: Part 2 – Case of combined sewers, Water Res., 46, 6693–6703, https://doi.org/10.1016/j.watres.2011.09.041, 2012. a, b
Gräler, B., Pebesma, E., and Heuvelink, G.: SpatioTemporal Interpolation using {gstat}, R J., 8, 204–218, 2016. a
Hager, W. H.: Wastewater hydraulics, second edn., https://doi.org/10.1080/00221686.2011.614723, SpringerVerlag, Berlin, Heidelberg, Germany, 2010. a, b
Hammersley, J. and Handscomb, D.: Monte Carlo Methods, Methuen & Co Ltd, London, 1964. a
Heip, L., Assel, J. V., and Swartentenbroekx, P.: Sewer flow quality modelling, Water Sci. Technol., 36, 177–184, 1997. a
Heuvelink, G. B. M.: Error Propagation in Environmental Modelling with GIS, Reserach Monographs in GIS, CRC Press Taylor & Francis Group, Taylor & Francis Ltd. 11 New Fetter Lane, London EC4P 4EE, UK, 1998. a, b, c, d
Heuvelink, G. B. M., Brown, J. D., and van Loon, E. E.: A probabilistic framework for representing and simulating uncertain environmental variables, Int. J. Geogr. Inf. Sci., 21, 497–513, https://doi.org/10.1080/13658810601063951, 2007. a
House, M. A., Ellis, J. B., Herricks, E. E., HvitvedJacobsen, T., Seager, J., Lijklema, L., Aalderink, H., and Clifforde, I. T.: Urban Drainage – Impacts on Receiving Water Quality, Water Sci. Technol., 27, 117–158, https://doi.org/10.2166/wst.1993.0293, 1993. a
Huang, H., Xiao, X., Yan, B., and Yang, L.: Ammonium removal from aqueous solutions by using natural Chinese (Chende) zeolite as adsorbent, J. Hazard. Mater., 175, 247–252, https://doi.org/10.1016/j.jhazmat.2009.09.156, 2010. a
Hutton, C., VamvakeridouLyroudia, L., Kapelan, Z., and Savic, D.: Uncertainty Quantification and Reduction in Urban Water Systems ( UWS ) Modelling: Evaluation Report, Tech. rep., The European Commission, 2011. a, b
IFAK: SIMBA (Simulation of Biological Wastewater Systems): Manual and Reference, Tech. rep., Institut für Automation und Kommunikation e. V., Magdeburg, Germany, 2007. a
JervesCobo, R., Benedetti, L., Amerlinck, Y., Lock, K., De Mulder, C., Van Butsel, J., Cisneros, F., Goethals, P., and Nopens, I.: Integrated ecological modelling for evidencebased determination of water management interventions in urbanized river basins: Case study in the Cuenca River basin (Ecuador), Sci. Total Environ., 709, 1–18, https://doi.org/10.1016/j.scitotenv.2019.136067, 2020. a
Kalman, R. E.: A new approach to linear filtering and prediction problems, Transactions of the American Society of Mechanical Engineers: Journal of Basic Engineering, 82D, 35–45, 1960. a
Kalos, M. H. and Whitlock, P. A.: Monte Carlo Methods, 2 edn., WILEYVCH Verlag GmbH & Co. KGaA, Weinheim, Germany, 2008. a
Katukiza, A. Y., Ronteltap, M., Niwagaba, C. B., Kansiime, F., and Lens, P. N. L.: Grey water characterisation and pollutant loads in an urban slum, Int. J. Environ. Sci. Technol., 12, 423436, https://doi.org/10.1007/s1376201304515, 2014. a
Kleidorfer, M. and Rauch, W.: An application of Austrian legal requirements for CSO emissions, Water Sci. Technol., 64, 1081, https://doi.org/10.2166/wst.2011.560, 2011. a, b, c
Kuczera, G. and Parent, E.: Monte Carlo assessment of parameter uncertainty in conceptual catchment models: The Metropolis algorithm, J. Hydrol., 211, 69–85, https://doi.org/10.1016/S00221694(98)00198X, 1998. a
Leopold, U., Heuvelink, G. B. M., Tiktak, A., Finke, P. A., and Schoumans, O.: Accounting for change of support in spatial accuracy assessment of modelled soil mineral phosphorous concentration, Geoderma, 130, 368–386, 2006. a
Litrico, X. and Fromion, V.: Modeling and Control of Hydroystems, SpringerVerlag, London, United Kingdom, https://doi.org/10.1017/CBO9781107415324.004, 2009. a
Luetkepohl, H.: New Introduction to Multiple Time Series Analysis, SpringerVerlag, Berlin, Heidelberg, Germany, 2005. a
Marwick, B. and Krishnamoorthy, K.: cvequality: Tests for the Equality of Coefficients of Variation from Multiple Groups, available at: https://github.com/benmarwick/cvequality (last access: 29 December 2020), r package version 0.2.0, 2019. a
McMillan, H., Jackson, B., Clark, M., Kavetski, D., and Woods, R.: Rainfall uncertainty in hydrological modelling: An evaluation of multiplicative error models, J. Hydrol., 400, 8394, https://doi.org/10.1016/j.jhydrol.2011.01.026, 2011. a
Miskewitz, R. and Uchrin, C.: InStream Dissolved Oxygen Impacts and Sediment Oxygen Demand Resulting from Combined Sewer Overflow Discharges, J. Environ. Eng., 139, 1307–1313, https://doi.org/10.1061/(ASCE)EE.19437870.0000739, 2013. a
MorenoRodenas, A. M., TscheiknerGratl, F., Langeveld, J. G., and Clemens, F. H.: Uncertainty analysis in a largescale water quality integrated catchment modelling study, Water Res., 158, 46–60, https://doi.org/10.1016/j.watres.2019.04.016, 2019. a
Neumann, M. B.: Uncertainty Analysis for Performance Evaluation and Design of Urban Water Infrastructure, Ph.D. thesis, Swiss Federal Institute of Technology, ETH Zurich, Switzerland, 2007. a, b
Nol, L., Heuvelink, G. B. M., A. Veldkamp, de Vries, W., and Kros, J.: Uncertainty propagation analysis of an N2O emission model at the plot and landscape scale, Geoderma, 159, 9–23, 2010. a
RCoreTeam and contributors worldwide: The R Stats Package, The R Project for Statistical Computing, 3.5.0 edn., available at: https://stat.ethz.ch/Rmanual/Rdevel/library/stats/html/00Index.html (last access: 29 December 2020), 2017. a
Refsgaard, J. C.: Parameterisation, calibration and validation of distributed hydrological models, J. Hydrol., 198, 69–97, https://doi.org/10.1016/S00221694(96)03329X, 1997. a
Refsgaard, J. C., van der Sluijs, J. P., Højberg, A. L., and Vanrolleghem, P. A.: Uncertainty in the environmental modelling process  A framework and guidance, Environ. Model. Softw., 22, 1543–1556, https://doi.org/10.1016/j.envsoft.2007.02.004, 2007. a
Renard, B., Kavetski, D., Kuczera, G., Thyer, M., and Franks, S. W.: Understanding predictive uncertainty in hydrologic modeling: The challenge of identifying input and structural errors, Water Resour. Res., 46, 1–22, https://doi.org/10.1029/2009wr008328, 2010. a
Sandric, I., Ionita, C., Chitu, Z., Dardala, M., Irimia, R., and Furtuna, F. T.: Using CUDA to accelerate uncertainty propagation modelling for landslide susceptibility assessment, Environ. Model. Softw., 115, 176–186, https://doi.org/10.1016/j.envsoft.2019.02.016, 2019. a
Sriwastava, A. K., Tait, S., Schellart, A., Kroll, S., Dorpe, M. V., Assel, J. V., and Shucksmith, J.: Quantifying Uncertainty in Simulation of Sewer Overflow Volume, J. Environ. Eng., 144, 04018050, https://doi.org/10.1061/(ASCE)EE.19437870.0001392, 2018. a, b, c
Statec: Statistics Portal of the GransDuchy of Luxembourg, available at: https://statistiques.public.lu, last access: 29 December 2020. a
Steinel, A. and Margane, A.: Best management practice guideline for wastewater facilities in karstic areas of Lebanon with special respect to the protection of ground and surface waters, Tech. Rep. 2, Federal Ministry for Economic Cooperation and Development, Bundesanstalt für Geowissenschaften und Rohstoffe, BGR, Hannover, Germany, 2011. a, b
Sten, J., Lilja, H., Hyväluoma, J., Westerholm, J., and Aspnäs, M.: Parallel flow accumulation algorithms for graphical processing units with application to RUSLE model, Comput. Geosci., 89, 88–95, https://doi.org/10.1016/j.cageo.2016.01.006, 2016. a
Toffol, S. D.: Sewer system performance assessment – an indicators based methodology, Ph.D. thesis, Universität Innsbruck, Innsbruck, Austria, 2006. a, b
TorresMatallana, J., Leopold, U., and Heuvelink, G.: stUPscales: an Rpackage for spatiotemporal Uncertainty Propagation across multiple scales with examples in urban water modelling, Water, 10, 1–30, https://doi.org/10.3390/w10070837, 2018a. a
TorresMatallana, J., Leopold, U., and Heuvelink, G.: stUPscales: SpatioTemporal Uncertainty Propagation Across Multiple Scales, available at: https://CRAN.Rproject.org/package=stUPscales (last access: 29 December 2020), R package version 1.0.5.0, 2019. a, b
TorresMatallana, J. A., Leopold, U., and Heuvelink, G. B. M.: Multivariate autoregressive modelling and conditional simulation of precipitation time series for urban water models, European Water, 57, 299–306, 2017. a, b
TorresMatallana, J. A., Klepiszewski, K., Leopold, U., and Heuvelink, G.: EmiStatR: a simplified and scalable urban water quality model for simulation of combined sewer overflows, Water, 10, 1–24, https://doi.org/10.3390/w10060782, 2018b. a, b, c, d
TorresMatallana, J. A.: R code and data to reproduce figures from the “Multivariate autoregressive modelling and conditional simulation for temporal uncertainty analysis of an urban water system in Luxembourg” paper, (Version v1.2) [Data set], Zenodo, https://doi.org/10.5281/zenodo.3928079, 2020. a
TscheiknerGratl, F., Lepot, M., MorenoRodenas, A., and Schellart, A.: QUICS Deliverable 6.7: A Framework for the application of uncertainty analysis, Tech. rep., Delft University of Technology, University of Sheffield and CH2M, Zenodo, https://doi.org/10.5281/zenodo.1240926, 2017. a
TscheiknerGratl, F., Bellos, V., Schellart, A., MorenoRodenas, A., Muthusamy, M., Langeveld, J., Clemens, F., Benedetti, L., RicoRamirez, M. A., de Carvalho, R. F., Breuer, L., Shucksmith, J., Heuvelink, G. B., and Tait, S.: Recent insights on uncertainties present in integrated catchment water quality modelling, Water Res., 150, 368–379, https://doi.org/10.1016/j.watres.2018.11.079, 2019. a
van der Keur, P., Henriksen, H. J., Refsgaard, J. C., Brugnach, M., PahlWostl, C., Dewulf, A., and Buiteveld, H.: Identification of major sources of uncertainty in current IWRM practice. Illustrated for the Rhine Basin, Water Resour. Manage., 22, 1677–1708, https://doi.org/10.1007/s1126900892486, 2008. a, b, c, d
Viana da Silva, A. M. E., Bettencourt da Silva, R. J. N., and Camões, M. F. G. F. C.: Optimization of the determination of chemical oxygen demand in wastewaters, Anal. Chim. Acta, 699, 161–169, https://doi.org/10.1016/j.aca.2011.05.026, 2011. a
Vrugt, J. A. and Robinson, B. A.: Improved evolutionary optimization from genetically adaptive multimethod search, P. Natl. Acad. Sci. USA, 104, 708–711, https://doi.org/10.1073/pnas.0610471104, 2007. a
Vrugt, J. A., Gupta, H. V., Bastidas, L. A., Bouten, W., and Sorooshian, S.: Effective and efficient algorithm for multiobjective optimization of hydrologic models, Water Resour. Res., 39, SWC 51–SWC 519, https://doi.org/10.1029/2002wr001746, 2003a. a, b
Vrugt, J. A., Gupta, H. V., Bouten, W., and Sorooshian, S.: A Shuffled Complex Evolution Metropolis algorithm for optimization and uncertainty assessment of hydrologic model parameters, Water Resour. Res., 39, SWC 11–SWC 114, 2003b. a
Vrugt, J. A., ter Braak, C. J. F., Clark, M. P., Hyman, J. M., and Robinson, B. A.: Treatment of input uncertainty in hydrologic modeling: Doing hydrology backward with Markov chain Monte Carlo simulation, Water Resour. Res., 44, 1–15, https://doi.org/10.1029/2007WR006720, 2008. a
Wadoux, A.C., Brus, D., RicoRamirez, M., and Heuvelink, G.: Sampling design optimisation for rainfall prediction using a nonstationary geostatistical model, Adv. Water Res., 107, 126–138, https://doi.org/10.1016/j.advwatres.2017.06.005, 2017. a
Walker, W., Harremoës, P., Rotmans, J., van der Sluijs, J., van Asselt, M., Janssen, P., and Krayer von Krauss, M.: Defining Uncertainty: A Conceptual Basis for Uncertainty Management in ModelBased Decision Support, Integrated Assessment, 4, 5–17, https://doi.org/10.1076/iaij.4.1.5.16466, 2003. a, b, c
Webster, R. and Heuvelink, G. B. M.: The Kalman filter for the pedologist's tool kit, Eur. J. Soil Sci., 57, 758773, https://doi.org/10.1111/j.13652389.2006.00879.x, 2006. a
Webster, R. and Oliver, M.: Geostatistics for Environmental Scientists, 2nd edn., Wiley, Hoboken, New Jersey, USA, 2007. a
Welker, A.: Emissions of pollutant loads from combined sewer systems and separate sewer systems – Which sewer system is better, in: 11th International Conference on Urban Drainage, edited by: ICUD, Edinburgh, Scotland, UK, IAHR/IWA Joint Committee on Urban Drainage, 2008. a
Yu, L., Rozemeijer, J. C., Broers, H. P., van Breukelen, B. M., Middelburg, J. J., Ouboter, M., and van der Velde, Y.: Drivers of nitrogen and phosphorus dynamics in a groundwaterfed urban catchment revealed by high frequency monitoring, Hydrol. Earth Syst. Sci. Discuss., https://doi.org/10.5194/hess202034, in review, 2020. a, b
Zhou, X., Polcher, J., Yang, T., and Huang, C.S.: A new uncertainty estimation approach with multiple datasets and implementation for various precipitation products, Hydrol. Earth Syst. Sci., 24, 2061–2081, https://doi.org/10.5194/hess2420612020, 2020. a, b, c
Zoppou, C.: Review of urban storm water models, Environ. Model. Softw., 16, 195–231, https://doi.org/10.1016/S13648152(00)000840, 2001. a