## Abstract

This work proposes a Bayesian non-informative reconstruction of virtual state variables in the representation of stormwater total suspended solids pollutographs by the traditional wash-off models, based on 255 rainfall events measured in a 185 ha French urban catchment. Results from event-based analyses revealed the missing representation of an essential process in the traditional rating curve (RC) model (simplest wash-off model) for 56% of the rainfall events. The unsatisfactory performances of the RC model are found to be not necessarily linked to antecedent dry weather conditions, as assumed by a great number of accumulation/wash-off models. Statistical tests suggest that non-representable rainfall events by the RC model are randomly distributed in time. The proposed Bayesian reconstructions of a potential process missed by the RC model exhibit a suitable identifiability at an intra-event scale. However, these reconstructions are not interpretable from the traditional accumulation/wash-off notions, i.e. in terms of a unique state of virtual available mass over the catchment that is decreasing over time, due to their high unrepeatability regarding their shape and their low prediction capacity for other rainfall events.

## INTRODUCTION

During the past 40 years, modelling the dynamics of stormwater total suspended solids (TSS) loads at the outlet of urban catchments has been mainly addressed by the idea of accumulation/wash-off (originally by Sartor *et al.* (1974)). Although this conceptual model does not constitute a rigorous physical description of the system, it is usually considered as a physically based model in the sense that it is based on mass conservation principles and its parameters can be interpretable from a physical point of view (Bonhomme & Petrucci 2017). Indeed, this accumulation/wash-off model has been implemented in a massive amount of literature and a wide variety of urban catchments (size, land use, complexity), for purposes such as real time control, climate change assessment, risk analysis, water management, and in multiple commercial software tools. However, the unsatisfactory performance of this approach is frequently reported, as well as the difficulty of generalizing its results to real world applications, especially as urban catchments are large and complex (e.g. Vaze & Chiew 2003; Deletic *et al.* 2012).

The accumulation/wash-off original model relies on the idea of a virtual ‘mean’ state of TSS mass over the catchment that is accumulated during dry periods and is later on washed off by rainfall and runoff, producing the TSS load at the outlet of the catchment. However, this notion from Sartor *et al.* (1974), followed by various analogue formulations (e.g. Egodawatta *et al.* 2009; Freni *et al.* 2009; Crobedu & Bennis 2011) has been developed under the following experimental/methodological settings: (i) controlled experimental conditions in small scale systems, (ii) limited number of TSS data, (iii) limited number of rainfall events, and (v) limited assessment of uncertainty in data and model parameters. This suggests how uncertainty linked to application of these models outside their limits (e.g. in large and complex urban catchments – see Liu *et al.* (2012)) is expected to be amplified.

Furthermore, the existence of a ‘mean’ state of available TSS mass for large catchments (from 100–200 ha) has been questioned from a low identifiability and a high spatial variability of model parameters (Bonhomme & Petrucci 2017). This hypothesis can be nourished from satisfactory results when applying a simplified rating curve (RC) model (Huber *et al.* 1988) without including a ‘mean’ available mass term (e.g. Kanso *et al.* 2005). Under this perspective, one can ask for the existence of a deterministic global process missed by the RC model, essential to represent the pollutant loads, which is oversimplified or misinterpreted by the accumulation/wash-off idea.

For this purpose, the time variable parameters (TVPs) concept has been introduced in the hydrological and environmental context as a powerful statistical model-based approach to describe unobserved processes or state variables (e.g. Young 2013). The TVP method reconstructs the potential temporal dynamics of the optimal parameters for a given model. TVPs can be used with transfer function models estimated using instrumental variables (e.g. Young 2013), or with general conceptual models within Bayesian frameworks (Renard *et al.* 2010). Indeed, Bayesian total error analysis (BATEA) (Kuczera *et al.* 2006) has emerged as a promising model-based TVP estimation technique for reconstructing unmeasured inputs or state variables, adaptable to non-linear and complex model structures, including flexibility in the error model of the reconstructed state variable (e.g. for rainfall time series reconstruction – see Leonhardt *et al.* (2014)).

This paper aims to revisit the accumulation/wash-off idea by using Bayesian TVP reconstructions, scrutinizing for evidence of a deterministic global process and its interpretability as a ‘mean’ state of available pollutant mass being washed by rainfall. Furthermore, this accumulation/wash-off model structure is sought to be reformulated by analysing the inter-event repeatability of the reconstructed state variables. The experimental/methodological settings aimed at overcoming the limitations of previous studies by considering: (i) an urban catchment of 185 ha, (ii) TSS load measurements with a temporal resolution of 2 minutes (with data uncertainty), (iii) 255 rainfall events (from 2004 to 2011), (iv) the state-of-art calibration and parameter uncertainties assessment (Bayesian methods, e.g. Kuczera *et al.* (2006)).

## MATERIALS AND METHODS

### Accumulation/wash-off model

*et al.*(1974) describes the stormwater TSS load (kg) at the outlet of an urban catchment as the product of two factors: a sediment mass

*M*(

*t*) (kg) and a removal factor

*W*(

*t*) (-). In this representation,

*M*(

*t*) is a decaying state variable of a virtual available stock mass which in the end limits the load production given by a wash-off process

*W*(

*t*). Sartor

*et al.*(1974) describe the two terms as follows, including a set of calibration parameters: where

*M*

_{0}(kg) is the initial condition from which

*M*starts to decay (

*M*(0) =

*M*

_{0}, with

*t*

*=*0 the beginning of a rainfall event), and

*a*and

*r*express the influence of the stormwater flow

*Q*(

*t*) (m

^{3}/s) on the removed TSS mass. Several different formulations for

*M*and

*W*during rainfall events can be found in the literature and some examples are given in Table 1.

Formulations | Reference | ||
---|---|---|---|

M (kg) | W (-) | Parameters | |

M_{0}, a, r | Sartor et al. (1974) | ||

M_{0}, C1, C2 | Freni et al. (2009) | ||

M_{0}, C1, C2 | Egodawatta et al. (2009) | ||

1 | M_{0}, r | Huber et al. (1988) (RC model) |

Formulations | Reference | ||
---|---|---|---|

M (kg) | W (-) | Parameters | |

M_{0}, a, r | Sartor et al. (1974) | ||

M_{0}, C1, C2 | Freni et al. (2009) | ||

M_{0}, C1, C2 | Egodawatta et al. (2009) | ||

1 | M_{0}, r | Huber et al. (1988) (RC model) |

*et al.*1988), where

*M*is considered constant (

*M*

*=**M*

_{0}) and the wash-off term shows a non-linear relation to the flow:

In this case, the calibration parameters are reduced to two (*M*_{0} and *r*).

As a generality, it can be seen that for most model structures 0 < *M*(*t*) ≤ *M*_{0} and d*M*/d*t* ≤ 0 for all *t* during rainfall events, *M* being described as either a time variable or a constant process (e.g. Table 1). Indeed, the RC model in Equation (2) can be directly linked with any time variable formulation of *M* (as in Equation (1)), by making *M*_{0} a TVP. Therefore, this work proposes to undertake Bayesian reconstructions of the ‘virtual’ state variable *M* by adopting the RC traditional model in Equation (2) and replacing *M*_{0} by a TVP, (*t*), keeping *r* as a calibration parameter (formulation F1). The (*t*) estimation in F1 can be directly compared to a time constant or variable traditional *M* formulation (e.g. Table 1). As a complementary analysis, a second formulation (F2) is explored, in which the RC model is modified with (*t*) (-) as the TVP to be reconstructed and *M*_{0} as a calibration parameter.

### Experimental setup

The proposed methods are tested with information about 255 rainfall events from the Chassieu urban catchment (Lyon, France), measured between 2004 and 2011. This is one of the experimental sites of the OTHU project (Field Observatory for Urban Hydrology – www.othu.org). The basin is a 185 ha industrial area drained by a separate storm sewer system, with imperviousness and runoff coefficients of about 0.72 and 0.43 respectively. Flow rate *Q* (L/s) and *TSS* concentrations (mg/L) at the outlet of the catchment are measured in a 1.6 m circular concrete pipe with a 2 minute time-step resolution. *Q* is calculated from water depth with a relative standard uncertainty from 15% to 25%. Regarding *TSS* concentrations, the drainage water is pumped into an off-line monitoring flume in a shelter, where turbidity measurements (NTU) are converted into *TSS* concentration values by local calibration functions (see Sun *et al.* 2015). The standard uncertainties in the estimated *TSS* concentrations ranged from 11% to 30%. The uncertainties of the *TSS* load (kg) are propagated by the law of propagation of uncertainties (LPU) (ISO 2009), obtaining relative standard uncertainties from 15% to 50%. Standard uncertainties become higher for higher values of the measurands (*Q*, *TSS* concentration, *load*).

### Bayesian reconstruction of virtual state variables by time variable parameters

The reconstruction of a virtual state variable by TVPs consists in solving a calibration problem where a TVP is represented as a time series of parameters rather than a singular time-constant parameter. In principle, every time-step of a TVP can be considered as an independent parameter in the inference scheme, with a length equal to the length of the output series load. However, the dimensionality of this problem will be massive, risking over-parametrization and delivering incorrect estimations (Vrugt *et al.* 2009). To avoid this over-parametrization in the inference, the TVP time series has usually a coarser temporal resolution than the output data. For this purpose, a time window strategy is adopted, reducing the dimensionality of the problem by dividing the TVP reconstruction into equally spaced time windows (see e.g. Sun & Bertrand-Krajewski 2013). For this case, a further test with non-equally spaced time windows based in the cumulative volume reported results analogous to the equally spaced time windows approach.

*t*) for F1 or (

*t*) for F2) is proposed to be estimated jointly with the other set of regular calibration parameters θ by means of a Bayesian inference scheme. Therefore, the joint posterior probability density function (PDF) of θ and

*TVP*are calculated by Equation (3) as

*p*(θ,

*TVP*/

*Q*(

_{obs}*t*),

*load*(

_{obs}*t*)), given the input

*Q*(

_{obs}*t*) and output

*load*(

_{obs}*t*) data with some prior knowledge about θ and

*TVP*(from BATEA in Kuczera

*et al.*(2006)). The values of θ and

*TVP*values with the maximum likelihood in the estimation of the joint posterior PDF

*p*(θ,

*TVP*/

*Q*(

_{obs}*t*),

*load*(

_{obs}*t*)) are assumed be the ‘optimal’ parameters θ

*and*

_{opt}*TVP*((

_{opt}*t*)

*for F1 or (*

_{opt}*t*)

*for F2). where*

_{opt}*n*is the total observed load values in

*load*(

_{obs}*t*). The

*load*(

_{sim}*Q*(

_{obs}*t*), θ,

*TVP*) is the simulated load by the input flow rate series

*Q*(

_{obs}*t*) and a set of θ and

*TVP*.

*p*(θ) and

*p*(

*TVP*) represent a prior belief about the probability that a candidate set of θ and

*TVP*values are ‘true’ (assumed as a non-informative uniform distribution for all cases). Equation (3) allows the explicit separation of the model error of

*TVP*(second Pi product) from the error model of the output

*load*(first Pi product). With the purpose of finding a TVP estimation ‘as constant as possible in time’ (and therefore less informative), the error model of

*TVP*(second Pi product) is assumed to be proportional to the TVP's own variance

*Var*(

*TVP*). This constraint is aimed at avoiding ‘curve fitting’ the TVP with

*load*(

_{obs}*t*).

Both error models, for the *load* and TVP estimations (first and second Pi product respectively), are assumed to be independent and normally distributed, with the error variances and , respectively. is considered heteroscedastic, being equal to the square of the standard uncertainty of each observed value *load _{obs}*(

*t*) (e.g. Sun & Bertrand-Krajewski 2013). On the other hand, is considered to be homoscedastic and is estimated as another parameter in the set θ (e.g. Sage

*et al.*2015), expressed as (kg) for F1 and (-) for F2. The same Equation (3) is used for calibrating the RC model, by omitting the second Pi product term, delivering

*p*(θ/

*Q*,

*load*) as a posterior probabilistic characterization of θ (with θ

*also for the optimal parameters).*

_{opt}The DREAM algorithm (Vrugt 2016) is applied to solve Equation (3) and to evaluate the three formulations (the traditional RC, the F1 and F2 formulations) for the different θ and *TVP* variants (Table 2). The logarithmic form of Equation (3) is implemented to ensure numerical stability and a maximum number of simulations of 6 × 10^{5} is used with a maximum of 30 parallel Markov Chains to reach convergence (Gelman and Rubin (GR) convergence criteria >1.2, see Gelman & Rubin (1992)). The TVP reconstructions are undertaken with a resolution of 12 time windows (i.e. equivalent to 12 parameters, see Table 2), balancing between the convergence of the algorithm (GR > 1.2) and capturing the global dynamics of the TVPs (see application in Sun & Bertrand-Krajewski (2013)). Table 2 also summarizes the parameters' minimum and maximum values search ranges in the prior uniform distributions *p*(θ) and *p*(*TVP*) in Equation (3). These ranges are defined for *r* and *M*_{0} based on the literature (Kanso *et al.* 2005) and are equally adopted for each window of (*t*) and (*t*). The maximum of the error variances and are defined, respectively, as four times the standard deviation of *M*_{0} and *r* in the RC calibration.

Formulation | Parameters θ | TVP | Number of parameters θ + TVP | Minimum/maximum search range values |
---|---|---|---|---|

RC model | r (-), M_{0} (kg) | – | 2 | r (0, 5); M_{0} (0, 8 × 10^{5}) |

F1 | r (-), (kg) | (t) (kg) | 2 + 12 = 14 | r (0, 5); (t) (0, 8 × 10^{5}) for each t; (0, 2.4 × 10^{5}) |

F2 | M_{0} (kg), (-) | (t) (-) | 2 + 12 = 14 | M_{0} (0, 8 × 10^{5}); (t) (0, 5) for each t; (0, 1.5) |

Formulation | Parameters θ | TVP | Number of parameters θ + TVP | Minimum/maximum search range values |
---|---|---|---|---|

RC model | r (-), M_{0} (kg) | – | 2 | r (0, 5); M_{0} (0, 8 × 10^{5}) |

F1 | r (-), (kg) | (t) (kg) | 2 + 12 = 14 | r (0, 5); (t) (0, 8 × 10^{5}) for each t; (0, 2.4 × 10^{5}) |

F2 | M_{0} (kg), (-) | (t) (-) | 2 + 12 = 14 | M_{0} (0, 8 × 10^{5}); (t) (0, 5) for each t; (0, 1.5) |

The parameters θ and the reconstruction of virtual state variables *TVP* are estimated for each individual rainfall event (event-based calibration). This eliminated the need for a ‘dry build up’ model (e.g. Freni *et al.* 2009) as the initial sediment mass (*M*_{0}) is directly estimated for each rainfall event.

For each *i*th calibration event, the results for the different formulations proposed in Table 2 (posterior probabilities and optimal parameter sets) provide the basis for revisiting some of the traditional accumulation wash-off ideas. This is performed by looking at:

- (i)
a relationship between the Nash–Sutcliffe efficiency (NS

) (comparing_{i}*load*(_{sim}*t*)and_{i}*load*(_{obs}*t*)) obtained from the RC model event-based calibration, as an evidence of a process potentially missed by this model structure, and the antecedent dry weather period of the rainfall event (ADWP_{i}). This linkage is evaluated by considering further characteristics of the rainfall events (mean and maximum flow rate and the date of the event) in addition to the ADWP, into a principal component analysis (PCA) (_{i}*relation of RC model performance to inter-event characteristics*). - (ii)
the NS

(comparing_{i}*load*(_{sim}*t*)and_{i}*load*(_{obs}*t*)) for the event-based reconstructions_{i}*TVP*in the F1 and F2 formulations, besides the well-posedness and identifiability of the inferred joint posterior PDFs,_{i}*p*(θ,*TVP*/*Q*,*load*)(_{i}*intra-event identifiability of a potential process missed by RC model*). - (iii)
similarities in the shape or dynamics of the

*TVP*_{opt}_{i}estimation with other reconstructions, obtained from different rainfall events (*inter-event identifiability from repeatability*). - (iv)
the capacity of a given set of θ

_{opti}and*TVP*_{opt}_{i}to represent the load in another rainfall event, evaluated by the NS (*inter-event transferability*, see Bardossy & Singh (2008)). - (v)
further hypotheses about a potential missing process based on physical knowledge about the system and the obtained results (

*interpretability*).

## RESULTS AND DISCUSSION

The available data and simulation results from a local event-based calibration with the RC model and the F1 and F2 formulations applied to event *i* = 16 are shown as an example in Figure 1 (event from 23/09/2004 22:00 to 24/09/2004 07:00). The input hydrograph *Q _{obs}*(

*t*)

_{16}is shown in Figure 1(a) and the observed TSS

*load*(

_{obs}*t*)

_{16}in Figure 1(b) (measurements in solid black lines with corresponding 95% coverage intervals), jointly with the optimal TSS load simulation

*load*(

_{sim}*t*)

_{16}for RC (bottom solid line, in green), F1 (next to top solid line, in blue) and F2 (next to bottom solid line, in red)

*.*Figure 2 shows a summary of the

*p*(θ,

*TVP*/

*Q*,

*load*)

*inference for F1 (left) and F2 (right), including the reconstruction of the optimal TVP (*

_{16}*t*)

*and (*

_{16}*t*)

*, with the optimal values (*

_{16}*t*)

_{opt}_{16}and (

*t*)

_{opt}_{16}and the 95% coverage intervals (Figure 2(a) and 2(b), respectively). The correlation matrix of the parameters θ for the formulations F1 and F2 is also estimated for the example event 16 (Figure 2(c) and 2(d), respectively).

The three event-based Bayesian formulations proposed in Table 2 are undertaken, as for the example event 16 shown above, with each of the 255 rainfall events of the dataset. This, given that the aim of the present study is to analyse and identify potential key processes missed by a model structure (RC model) and not to carry out a traditional calibration/verification process. The results are discussed in the following sections.

### Relation of RC model performance to inter-event characteristics

The RC model local calibrations report adjustments between the simulated and measured loads of NS* _{i}* < 0.8 for 142 of the 255 events (56%). Therefore, the RC model can be considered by itself as unsatisfactory by local calibrations, suggesting the lack of a key process in this model structure when applied to a large (and complex) urban catchment. Furthermore, a PCA shows no relation between the representability of a given event by the RC model (black NS

*<0.8 and grey NS*

_{i}*> 0.8) and different physical characteristics of the event (maximum flow rate (m*

_{i}^{3}/s), mean flow rate (m

^{3}/s), ADWP (days), date of the event (hours)) (data scaled by a

*z*normalization with zero mean and unitary standard deviation, see Kreyszig (1979)) (Figure 3). Thus, if the unsatisfactory performances of the RC model (events with NS < 0.8 in black) are due to the lack of a potential deterministic process, this process seems not to be dependent on the ADWP, as assumed by a great number of accumulation/wash-off models (no separation between black and grey events in Figure 3). In addition, a Wald Wolfowitz test (Wald & Wolfowitz 1943) suggested that the distribution of these ‘black’ or ‘grey’ events in time, as a time series vector of binary states, follows a random sequence (

*p*-value < 0.05).

### Intra-event identifiability of a potential process missed by the RC model

Although there is an intrinsic correlation between RC model parameters (Kanso *et al.* 2005), the intra-event identifiability of the *TVP* jointly inferenced with θ gives promising results. The F1 or F2 formulations achieve greater NS values in all the events than the RC model, increasing the values of NS* _{i}* > 0.8 for 60% of cases in which RC reported NS

*< 0.8. In the example shown in Figure 1(b), the NS for the RC model (bottom simulation line, in green) is 0.65, while the F1 results (next to top simulation line, in blue) and F2 (next to bottom simulation line, in red) show NS above 0.8. However, this improvement in the NS values for F1 or F2 formulations can be simply a numerical effect resulting from increasing the number of parameters. Therefore, further well-posedness and identifiability analyses were undertaken regarding the Bayesian inferences for F1 and F2. For example, in the case of θ in F1 and F2, the couples*

_{i}*r*;

_{i}*(for F1) and*

_{i}*M*;

_{0 i}*(for F2) exhibited an appropriate parametric identifiability with unimodal posterior PDFs, with spurious parametric correlations of*

_{i}*Rho*(

*r*;

_{i}*) < 0.3 (F1) and*

_{i}*Rho*(

*M*;

_{0 i}*) < 0.11 (F2) for more than 90% of the events. For the example event 16, the parametric correlations are of*

_{i}*Rho*(

*r*

_{16};

_{16}) = 0.05 and

*Rho*(

*M*

_{0 16};

_{16}) = 0.12 (Figure 2(c) and 2(d), respectively). In addition, undesired correlations between the measured

*load*(

_{obs}*t*)

*and the estimated TVP (*

_{i}*t*)

*or (*

_{opt i}*t*)

*are above 0.6 and 0.5, respectively, for only 25% of the events. These results bring evidence that the*

_{opt i}*TVP*reconstructions might be contributing with additional information (also supported by NS values higher than for RC), without mimicking the measured

_{opt}_{i}*load*dynamic.

### Inter-event identifiability from repeatability

A functional clustering by *k*-centre method is applied to identify groups of the *TVP _{opt}* time-varying curves with similar shapes among the different rainfall events (from Chiou & Li 2007). A number of

*k*groups is defined in order to visualize

*k*different potential ‘repeated’ behaviours in the set of optimal

*TVP*curves (F1 or F2). For interpretability, each of the

_{opt}*TVP*

_{opt}_{i}curves is standardized by the normalization

*z*(

*TVP*

_{opt}_{i}) in the clustering curves analysis (Kreyszig 1979). For the F1 formulation, applying the Chiou & Li (2007) method with

*k*= 2 groups allows separating the (

*t*)

*into two ‘similar’ clusters, shown in Figure 4(a) (light and dark coloured groups of curves, corresponding to 57% and 43% of the events). The mean curve for each of these groups of (*

_{opt}*t*)

*curves is shown in Figure 4(b), along with their 95% coverage interval.*

_{opt}The traditional interpretation of *M* can be preliminarily associated with the dark coloured group in Figure 4(b), as the mean curve shows an apparent decaying trend. However, this hypothesis would lack of an immediate physical interpretation for the remaining 57% of events (light coloured group). Furthermore, the trends of the mean curves for the dark coloured group (negative trend) and light coloured group (positive trend) show to be statistically insignificant (*p*-value < 0.05 for a trend test), due to the strong variability of the curves inside each light or dark coloured group. This effect of randomness is stronger when the variability of (t)_{i} (estimated by calculating the 95% coverage intervals from the posterior distributions) is considered. Similar conclusions are drawn for clustering the (*t*)_{opti} curves into more *k* > 2 groups. Furthermore, the same analysis applied for clustering the (*t*)_{opti} curves (F2 formulation) with *k* ≥ 2 clusters does not show any significant average trend or *‘*repeated’ behaviours in any of the clustered groups of curves (Figure 4(c) and 4(d) for *k* = 2 groups in light and dark coloured). These results reveal the difficulty in identifying or characterizing a unique virtual process potentially missed by the RC regarding an inter-event scale, given the lack of repeatability of the shapes of the *TVP _{opt i}* curves.

### Inter-event transferability

The transferability of TVPs is analysed, investigating how a given *TVP _{opt}*

_{i}time series is able to reproduce another event from the dataset. Results for (

*t*) or (

*t*) are analogous; therefore discussion focuses on (

*t*) estimations. Each of the 29 most transferable estimations of (

*t*)

*is able to explain at least 30 rainfall events (NS > 0.8). On the other hand, for the optimal local estimations θ*

_{opt i}*of RC, 60 estimations are able to explain at least 30 events (NS > 0.8). The flatter curves from (*

_{opt i}*t*)

*tend to be more transferable, as they resemble the constant values in the RC formulation. These results stress the low transferability of the potential missing processes (F1 or F2) to further rainfall events.*

_{opt i}### Interpretability

These results bring evidence of a potential missing process in the RC model when applied to a large urban catchment, which is not necessarily linked to the ADWP, as assumed by a great number of accumulation/wash-off model structures. Although both processes ((*t*) or (*t*)) are suitable candidates to explain RC obstacles from an intra-event analysis, there is no evidence that F1 or F2 is a more valid approach than the other. Indeed, the F1 or F2 formulations show the same explanatory capacity, in the sense that none of them performs better. This highlights the challenge in terms of identifiability and unicity of a potential process missed by the RC model structure, which is hardly identifiable from an inter-event analysis.

## CONCLUSIONS

This work suggests the missing representation of an essential process in the simplified RC model, when reproducing the TSS load pollutograph from flow rate observations in 255 rain events measured at the output of a 185 ha urban catchment. The results indicate that the high unrepeatability of this missing process makes it hardly interpretable in terms of a virtual unique state of available mass in the catchment that is decreasing over time and whose initial value is dependent on dry weather conditions, as assumed by a great number of traditional models. Indeed, this study shows how high-time resolution quality measurements can provide a support to revisit and question existing models, and to potentially allow the development of new stormwater quality model formulations. For example, future research can be oriented to understanding the nature of this potentially missed process by considering the obtained reconstructions as realizations of a stochastic governing process (e.g. Del Giudice *et al.* (2016) with the rainfall potential approach). This stochastic process (if one can be defined) might report an average dynamic and variability (or occurrence) consistent with the obtained reconstructions, within an intra- and inter-event scale, respectively. Furthermore, a deeper understanding of the physical processes inside the sewer system (e.g. TSS resuspension, transfer) could be achieved, for example by including online TSS measurements at different points inside the catchment, nourishing new formulations for this potentially missed process. On the other hand, one should bear in mind that the obtained results are dependent on the studied hydrosystem. Therefore, the proposed methodologies are recommended to be applied in different urban catchments with similar data, in order to provide more generality to the conclusions as presented here.

## ACKNOWLEDGEMENTS

The authors are grateful to the OTHU project in Lyon, France (www.othu.org) for providing the dataset, to COLCIENCIAS and the French–Danish Research Collaboration Program (financed by French Institute in Denmark) for the financial support.

## REFERENCES

*ISO/IEC Guide 98-1:2009(E) Uncertainty of Measurement – Part 1: Introduction to the Expression of the Uncertainty in Measurement*. ISO, Geneva, Switzerland