Accounting for the Spatio-Temporal Variability of Pollutant Processes in Stormwater TSS Modeling Based on Stochastic Approaches

: Stormwater quality modeling remains one of the most challenging issues in urban hydrology today. The processes involved in contaminant generation and transport are very complex, with many associated uncertainties, including uncertainty arising from process variability. In this study, the spatio-temporal variability of build-up/wash-off processes in a heterogeneous urban catchment within the Parisian region is assessed based on three stochastic modeling approaches integrated into the physically based distributed hydrological model, the Urban Runoff Branching Structure (URBS) model. Results demonstrate that accounting for process variability at the scale of a hydrological element is important for analyzing the contamination recorded at the catchment outlet. The intra-event dynamics of total suspended solids (TSS) were most accurately selected for the stochastic exponential SWMM model, as this model succeeded not only in simulating the general trend of TSS concentrations ﬂuctuations but also in replicating multiple peaks observed in pollutographs. The advantage of this approach is that it captures the stochastic nature of the processes with minimal prior knowledge and without extensive calibration, though further enhancement is necessary for it to become a useful tool to support decision making.


Introduction
Stormwater runoff is an important source of contamination in urban environments [1][2][3]. The design and implementation of management tools for stormwater quality control require knowledge of stormwater quality loads and dynamics, which are commonly derived using mathematical models [4].
Multiple types of urban stormwater quality models exist in the literature with varying degrees of complexity, ranging from generic empirical formulations to conceptual and physically-based models [5,6]. The performance of these models has constantly been put into question since, to date, no reliable stormwater quality model exists despite extensive research on the subject [7]. Conceptual models remain fairly common for stormwater quality modeling and have been integrated into software such as SWMM and Flupol [8,9], despite their limitations [10,11].
The URBS model, developed by the French Institute of Science and Technology for Transport, Development and Networks (IFSTTAR) in Nantes, France [29] was selected for the hydrologic modeling component of this study. This model provides a spatially distributed representation of the catchment and a physically-based description of the hydrological dynamics in surface and subsurface areas. Choosing a process-based distributed hydrological model was considered essential for investigating contaminant transport because the model provides a detailed description of spatial variations in land use characteristics on the elemental scale of a cadastral parcel.
Catchment morphology is represented as a set of urban hydrological elements (UHE) connected to a hydrological network routing the flow toward the catchment outlet. The UHE includes a cadastral parcel and its neighboring street segments ( Figure 1). The characteristics of the UHE are determined from common urban databanks consisting of surface area, impervious fraction, vegetation fraction, slope, length, location of connection point to the hydrological network, and the depth of the drainage network at this point. The hydrological network is composed of a series of connected street gutters and sewer conduits characterized by their length, slope and diameter [29].
To account for the heterogeneity of land surface characteristics, three vertical profiles were defined based on the distinction of land uses for each UHE, including roofs, roads and natural soils. The hydrological processes at the UHE scale were modeled for each vertical profile in four reservoirs, representing the interception over the surface, the surface area, the vadose and the saturated zone. The processes include: interception by trees over the surface area, water infiltration into the soil, surface runoff, evaporation of water intercepted at the ground surface, plant transpiration and drainage of soil water. A detailed description of the processes and their representative equations can be found in [29].
First, the runoff generated by each UHE is drained on the surface through a flow path according to the connection point of the first manhole in the sewer network, using a travel time function. Then, the flow is transported through the sewer to the catchment outlet using the Muskingum-Cunge scheme. If necessary, the process representation in this model allows for individual simulation of the runoff derived from each land use.  [29]). In the figure, an urban hydrological element (UHE) is delineated by a dashed bold line, which encompasses a cadastral parcel and its adjacent street surface. The connection point to the hydrological network is represented as (Pc).  [29]). In the figure, an urban hydrological element (UHE) is delineated by a dashed bold line, which encompasses a cadastral parcel and its adjacent street surface. The connection point to the hydrological network is represented as (Pc).
To account for the heterogeneity of land surface characteristics, three vertical profiles were defined based on the distinction of land uses for each UHE, including roofs, roads and natural soils. The hydrological processes at the UHE scale were modeled for each vertical profile in four reservoirs, representing the interception over the surface, the surface area, the vadose and the saturated zone. The processes include: interception by trees over the surface area, water infiltration into the soil, surface runoff, evaporation of water intercepted at the ground surface, plant transpiration and drainage of soil water. A detailed description of the processes and their representative equations can be found in [29].
First, the runoff generated by each UHE is drained on the surface through a flow path according to the connection point of the first manhole in the sewer network, using a travel time function. Then, the flow is transported through the sewer to the catchment outlet using the Muskingum-Cunge scheme. If necessary, the process representation in this model allows for individual simulation of the runoff derived from each land use.

Model Implementation
Input files were set up through the processing of the topographic and land use data derived from urban databanks, using two free GIS software: QGIS and OrbisGIS. Analysis of the data shows that the catchment encompasses 274 cadastral parcels, which represent as many individual UHEs ( Figure 2). Among the 274 UHEs, 56 are not directly connected to an adjacent street, because they represent isolated parcels inside a block. Model Implementation Input files were set up through the processing of the topographic and land use data derived from urban databanks, using two free GIS software: QGIS and OrbisGIS. Analysis of the data shows that the catchment encompasses 274 cadastral parcels, which represent as many individual UHEs ( Figure 2). Among the 274 UHEs, 56 are not directly connected to an adjacent street, because they represent isolated parcels inside a block. Roads account for 31% of the total area of the catchment, while roofs account for 26%. Half of the impervious surface consists of roads, 40% are roofs and the remaining fraction is parking lots and sidewalks. The hydrological network consists of 102 segments: 55 of which represent street gutters and 47 represent sewer pipes. Gutters are considered to be small-sized pipes with a 0.25 m diameter. Main sewer pipes have a total length equal to 1165 m, as well as an ellipsoidal form with a 2.3 m height and a 1.3 m width. Minor pipes have a total length of 446 m, with a 0.3 m diameter circular section.
The second step for model implementation is the estimation of model parameters, which include 15 parameters for water budget and two parameters for transfer (Table A3). Determining the parameters of physical models is not a straightforward task since the perfect physical model does not exist and certain parameters need to be calibrated to the specific catchment [30,31]. Parameters that can be determined from field measurements and those that should be calibrated were distinguished. The values of the parameters estimated from the literature and used by [29] for a similar urban catchment located in a suburb of Nantes, France, were kept the same. These included the values of tree interception parameters, the maximum reservoir capacities for impervious land uses, the groundwater drainage coefficient, the root depth and the representative position of the vadose zone. The hydraulic conductivity of the street and the natural soil were estimated based on field measurements for the Le Perreux catchment. Other parameters corresponding to the hydrodynamic characteristics of the soil were calibrated.
Since UHE-based models are not widely studied in the literature, they are commonly calibrated using the simple artificial trial and error method. Parameters were continuously adjusted and the effects of parameter adjustments for each simulation were quickly assessed from simple graphical The second step for model implementation is the estimation of model parameters, which include 15 parameters for water budget and two parameters for transfer (Table A3). Determining the parameters of physical models is not a straightforward task since the perfect physical model does not exist and certain parameters need to be calibrated to the specific catchment [30,31]. Parameters that can be determined from field measurements and those that should be calibrated were distinguished. The values of the parameters estimated from the literature and used by [29] for a similar urban catchment located in a suburb of Nantes, France, were kept the same. These included the values of tree interception parameters, the maximum reservoir capacities for impervious land uses, the groundwater drainage coefficient, the root depth and the representative position of the vadose zone. The hydraulic conductivity of the street and the natural soil were estimated based on field measurements for the Le Perreux catchment. Other parameters corresponding to the hydrodynamic characteristics of the soil were calibrated.
Since UHE-based models are not widely studied in the literature, they are commonly calibrated using the simple artificial trial and error method. Parameters were continuously adjusted and the effects of parameter adjustments for each simulation were quickly assessed from simple graphical analyses and the value derived from the Nash Sutcliffe coefficient.

Stochastic Approach
While stormwater quantity runoff is well understood and sometimes accurately replicated by hydrological models, reliable stormwater quality models are almost nonexistent [7]. In commonly-used quality models, the variability of pollutant following a storm event is not adequately represented. In this context, we developed a new methodology for quality modeling that incorporates the temporal and the spatial variability of pollutant generation and transport. Three stochastic approaches were used to simulate catchment stormwater quality by accounting for the inter-event variability of pollutant loads and the variability of pollutant loads with respect to land use. A stochastic approach assumes that the model inputs are purely random. Although the cause and effect of pollutant generation and transport have been widely investigated, the available data are often insufficient to treat these processes in a rigorous, deterministic manner.
The main difference between a deterministic and a stochastic approach lies in the values associated with the model's parameters. In deterministic approaches, a unique value is assigned to each parameter while in a stochastic approach; the parameters are assigned with a probability distribution from which a value is sampled for each iteration. Thus, the values of the parameters are not the same at each sampling, which will subsequently vary the value of the outcome.
The concept of the proposed methods relies on a spatio-temporal conceptualization of both the event mean concentration (EMC) of corresponding pollutants and the parameters of the quality models intended to simulate the instantaneous variation of pollutant concentration. The EMC and the model parameters are described following the generic form ϕ(p;ev;i;Ω). ϕ represents the spatio-temporal variable quantity, whose possible values are included in the sampling range Ω, which depend on the parcel p, rainfall event ev and land use i. ϕ is a random variable chosen for each land use within every parcel at the beginning of a new rainfall event from the sampling range, with a uniform probability distribution. p is the index of the UHE, while ev represents the rainfall events sample. The value of i comprises two components: the roof and the road. The sampling range Ω comprises the realizations of ϕ within its lower and upper boundaries.
The simulation methodology consists of generating the stochastic population of EMC and quality model parameters, by drawing these variables from a uniform probability distribution, for each land use within each parcel at the beginning of a rainfall event. Depending on the used model (cf. 2.2.2), 1644, 548 or 1096 random combinations of EMCs and water quality parameters were generated at the scale of the urban catchment. Then, the TSS concentrations dynamics were derived at the scale of each UHE, following the equation of the water quality model. The instantaneous TSS loads were then calculated and routed to the catchment outlet, where the final TSS concentrations dynamics are obtained.

The Models
Three TSS models were applied in this approach. The first model is the commonly used exponential wash-off model, SWMM [8]. For this model, ϕ represents three variables which are the EMC and the model parameters C 1 and C 2 (Equation (1)). After the selection of the values of the EMC and the parameters, the initial available mass at the beginning of the rainfall event is calculated using the following Equation (2): where: M ero (t) is the eroded mass at time t (g/m 2 ) during the time step dt; Q(t) is the runoff rate at time t expressed as water depth per time step (mm/h); M build−up (t) is the available mass for erosion (g/m 2 ) at time t; C 1 and C 2 are the wash-off parameters; M build−up (initial) is the available mass at the beginning of the rainfall event (mg); EMC is the event mean concentration (mg/L); V total is the total runoff volume for the corresponding rainfall event (liter); Q(t 1 ) is the runoff rate at the first time step of the rainfall event (mm/h); dt is the time step (minutes); and N is the total time step of the event.
Next, the eroded mass is calculated using Equations (1) and (3): The second model uses a homothetic relationship of ratio K between the concentrations of TSS and the runoff. For this model, ϕ represents one variable: EMC. The EMC is first drawn from Ω and the ratio K is then calculated as a function of the EMC and the instantaneous flow using Equation (4): where K is the homothety ratio (mg·s/L 2 ); EMC is the event mean concentration (mg/L); V total is the total runoff volume for the corresponding rainfall event (L); Q(i) is the instantaneous flow at time i (L/s); dt is the time step; and N is the total time step of the event.
The pollutograph is then obtained following Equation (5): The third model is based on the equations of M(V) curves. M(V) curves indicate the distribution of pollutant mass vs runoff volume [32]. For this model, ϕ integrates two variables: EMC and the model parameter b which is the exponent of the M(V) curve. Equation (6) is applied to simulate the pollutograph for this model: where Cumulative eroded mass (t) is the cumulative eroded mass at time t (mg); EMC is the event mean concentration (mg/L); V total is the total runoff volume for the corresponding rainfall event (L); Cumulative runo f f volume (t) is the cumulative runoff volume at time t (liter); and b is the M(V) curve coefficient. The instantaneous eroded mass M eroded (t) is then calculated by subtracting the cumulative eroded mass at time step t and the cumulative eroded mass at time t − 1.
It should be noted that for all modeling approaches, the accumulated mass is considered equal to the eroded mass and it is calculated by the using of the EMC. This assumption is very important as it partially accounts for the variability in processes related to particle size distribution, where the entire fraction of build-up susceptible to displacement is being transported by the runoff.

Boundaries of the Sampling Ranges
The boundaries of the sampling ranges for the EMC and model parameters are based on a firm hypothesis. The values of the EMC are closely related to the type of land use, while wash-off parameter ranges are very wide and highly dependent on the calibration period. SWMM wash-off parameters only have physical meaning in relation to the characteristics of the studied catchment [33]. Even when parameters are optimized, a high level of uncertainty is associated with their calibrated values [13], revealing important challenges in defining sampling ranges for water quality parameters. The sampling ranges of the EMC for roads were defined by measurements collected at the inlet of the main boulevard within the catchment [11], while intervals for roof surfaces were determined by available data provided in the literature [3]. The boundaries of the sampling ranges for SWMM wash-off parameters C 1 and C 2 correspond to the lowest and highest values obtained from the best cases investigated in an earlier study conducted within the same catchment [11]. As for the M(V) curve coefficient b, the boundaries are selected to replicate first flush and uniformly distributed events, since the occurrence of last flush events within the region was rare during the monitored period. The sampling ranges for the wash-off parameters are identical for all land uses, while the EMC is distinguished for each land use. With the exception of the EMC for roofs, the sampling ranges for the wash-off parameters and the EMC for roads are determined by prior knowledge acquired on the site. This method faces difficulties in setting appropriate values for sampling ranges.
A summary of the modeling approaches and the corresponding inputs and parameters is shown in Table 1.

Catchment and Data Description
The studied catchment is located in the residential district of Le Perreux sur Marne in Val de Marne, a suburb east of Paris ( Figure 3). The catchment area is~12 ha and has an average slope of 2.6%. The site consists mainly of residential houses and commercial shops spread along both sides of the highly trafficked (>30,000 vehicles/day) Boulevard Alsace Lorraine, a main road which runs across the catchment. Impervious surfaces account for nearly 70% of the total catchment area. The catchment is drained via a separate sewer network that routes the flow to the outlet, located at the north-eastern edge of the catchment. The flow is collected from the surface by way of manholes before it is transported into the main pipes.
Water quantity and quality data were collected at the catchment outlet, which is equipped with a Nivus flowmeter recording flow measurements at a 2-min time step. Water quality data, including measurements of turbidity, was obtained using a multi-parameter probe DS5 OTT at a 1-min time step. The relationship between the turbidity and the concentration of TSS was established based on laboratory analysis conducted on several rainfall samples. The turbidity measurements were transformed into TSS concentrations using the following linear equation with a correlation factor R 2 > 0.98: [TSS] = 0.9 × Turbidity (7) where [TSS] is the total suspended solids concentration (mg/L); and Turbidity is the turbidity measurement in NTU.

Catchment and Data Description
The studied catchment is located in the residential district of Le Perreux sur Marne in Val de Marne, a suburb east of Paris ( Figure 3). The catchment area is ~12 ha and has an average slope of 2.6%. The site consists mainly of residential houses and commercial shops spread along both sides of the highly trafficked (>30,000 vehicles/day) Boulevard Alsace Lorraine, a main road which runs across the catchment. Impervious surfaces account for nearly 70% of the total catchment area. The catchment is drained via a separate sewer network that routes the flow to the outlet, located at the north-eastern edge of the catchment. The flow is collected from the surface by way of manholes before it is transported into the main pipes. Water quantity and quality data were collected at the catchment outlet, which is equipped with a Nivus flowmeter recording flow measurements at a 2-min time step. Water quality data, including measurements of turbidity, was obtained using a multi-parameter probe DS5 OTT at a 1-min time step. The relationship between the turbidity and the concentration of TSS was established based on The precipitation data was collected at a meteorological station installed on the rooftop of a building near the studied catchment. The station consisted of a tipping bucket rain gauge, which has a 0.1 mm resolution.
GIS data required for this study were obtained in collaboration with the National Institute of Geographic and Forest Information (IGN), France. An important dataset providing a detailed description of the parcels, buildings, vegetation and roads is available for all the French territories. The sewer network database was obtained from the local authority represented by the General Council of the Department of Val de Marne. GIS data were used to estimate the different morphological characteristics needed to describe the urban hydrological elements and to build the runoff branching network.

Evaluation Criterion of the Model Performance
The performance of both hydrological and water quality models was evaluated based on observational and statistical comparisons of simulations and measurements. A general overview of model behavior, in terms of dynamic replication of hydrographs and pollutographs, was first assessed by assessing the graphical plots of flows and TSS concentrations. The goodness of fit was then evaluated using multiple statistical criteria that estimate a normalized error between simulated and observed variables ( Table 2). The chosen mathematical efficiency criteria for model evaluations are the Nash Sutcliffe coefficient (C Nash ), the determination coefficient (C R 2 ), and the root mean square error (C RMSE ).
Choosing appropriate criteria for evaluating goodness of fit is very critical and highly dependent on the objectives of the modeler [34]. As the purpose of the hydrological model is to give accurate estimates of the instantaneous flow rate and the total runoff volume, the C Nash and C R 2 were chosen for its evaluation. The Nash coefficient is widely used as an indicator of the overall fit of hydrographs [35]. However, since the water quality modeling approach applied in this study aims to evaluate the variability of the processes on the pollutographs rather than fitting the measurements to the simulations, the evaluation of modeling performance in terms of Nash is not appropriate. Instead, the determination coefficient was used to assess whether the variability follows the same trend as the measurements. Therefore, the C R 2 was chosen to evaluate the water quality model, and was coupled with the C RMSE to have a general overview of the residuals between the measurements and the simulations. Table 2. Evaluation criteria of the performance of hydrological and water quality models. Sim t is the simulated value at time t; Obs t is the observed value at time t; Sim is the average of the simulated values; Obs is the average of the observed values; t 1 is the beginning of the simulated event; t n is the end of the simulated event; n is the total duration of the simulated event.

Hydrological Model Water Quality Model
Nash Sutcliffe

Model Application
Calibration of the hydrological model was carried out on 30 rainfall events that occurred between 8 October 2014 and 1 January 2015 (Table A1). These events represented typical rainfalls occurring within this area, characterized by a short duration and low intensity. Maximum rainfall intensities varied between 1.1 and 42.3 mm/h, while precipitation depth ranged between 2 and 22.1 mm. The duration of rainfall events varied from 52 min to 12 h, while the average antecedent dry weather period was equal to 2.5 days. Validation of the hydrological model was carried out on four rainfall events throughout the period of 31 March 2015 to 27 April 2015. The rainfall events had precipitation depths ranging from 2 to 6.7 mm. Maximum rainfall intensities varied between 2.3 and 43.3 mm/h. Hydrological simulations were compared to the measurements in terms of runoff volume and instantaneous flows.
As for water quality modeling, nine rainfall events occurring between 8 October 2014 and 1 January 2015 were simulated for each approach and their characteristics are summarized in Table A2. For each rainfall event, ten stochastic simulations were performed, since the variability of the results following each simulation was low and given the limitation of the computing time. The total washed-off load and TSS dynamics were determined for each simulation, and the evaluation criteria were calculated. This was followed by a statistical analysis to determine the minimum, maximum and mean values of the evaluation criteria.

Water Flow Simulations
The evaluation criteria for water flow simulations are summarized in (Table 3). The calibration results demonstrate a strong level of agreement between simulated runoff volumes and the observations with C Nash and C R 2 ≥ 0.9.
The model gave a satisfactory calibration performance, succeeding in replicating the instantaneous flow wherein the C Nash and C R 2 were higher than 0.7. Visual inspection of the simulated and observed hydrographs showed that the model provided a good replication of the overall trend in runoff generation (Figure 4). A minor time offset is observed for the second event on 8 October and the event on 13 December. The measured peak flows are slightly ahead of simulated peaks, possibly due to the simple routing scheme used to describe the transfer function. Using more robust numerical schemes based on Barre de St Venant equations might give a better performance. However, this effect is minimal in all cases and does not occur in all events.  The model was also validated on four events throughout the period from 31 March 2015 to 27 April 2015. Results indicate that the dispersion of the observed runoff volumes and flow are accurately explained by the simulations. However, the CNash for both runoff volumes and instantaneous flow are slightly lower than those calculated for the calibration period with values equal to 0.57 and 0.58, respectively. This is attributed to the parameters, which may not be suitable for representing the hydrological processes for this period.

TSS Loads
The simulated eroded masses for each event varied within the same ranges for all modeling approaches ( Figure 5), and their performance in terms of predicting washed-off loads were equal. This was expected since the EMC for each simulated event for all models were generated from the same ranges.
The result shows important deviations in the estimates of TSS loads from the measurements for the majority of rainfall events. The inter-event variability fluctuates in narrow ranges for measured loads, whereas it varies in larger ranges for the simulations. The measurements demonstrate that the yield of TSS loads at the catchment outlet is limited and tends towards a stable value. However, the simulations show higher variability from one event to another where no systematic trend in the model estimates is observed.
The models underestimated the measurements for three events, while overestimating them for the rest. This result effectively highlights that the ranges of the EMC might have been overestimated. The ranges of the EMC for the roads were defined based on the minimal and maximal EMC values The model was also validated on four events throughout the period from 31 March 2015 to 27 April 2015. Results indicate that the dispersion of the observed runoff volumes and flow are accurately explained by the simulations. However, the C Nash for both runoff volumes and instantaneous flow are slightly lower than those calculated for the calibration period with values equal to 0.57 and 0.58, respectively. This is attributed to the parameters, which may not be suitable for representing the hydrological processes for this period.

TSS Loads
The simulated eroded masses for each event varied within the same ranges for all modeling approaches ( Figure 5), and their performance in terms of predicting washed-off loads were equal. This was expected since the EMC for each simulated event for all models were generated from the same ranges.
The result shows important deviations in the estimates of TSS loads from the measurements for the majority of rainfall events. The inter-event variability fluctuates in narrow ranges for measured loads, whereas it varies in larger ranges for the simulations. The measurements demonstrate that the yield of TSS loads at the catchment outlet is limited and tends towards a stable value. However, the simulations show higher variability from one event to another where no systematic trend in the model estimates is observed.
The models underestimated the measurements for three events, while overestimating them for the rest. This result effectively highlights that the ranges of the EMC might have been overestimated. The ranges of the EMC for the roads were defined based on the minimal and maximal EMC values calculated over twelve samples collected at the boulevard Alsace Lorraine's inlet during the experimental campaign. Previous analysis of these values showed that they are higher than EMCs calculated for other road catchments due to the site's high traffic density [11]. However, some of the road segments in the catchment are not necessarily frequented at the same density as the boulevard; therefore, their loadings of TSS might be lower. Another reason for the high EMC, might be related to the trapping efficiency in gully pots. The trapping efficiency depends greatly on the size of sediment and the water flow rate. While coarse materials can be retained efficiently, the removal efficiency of fine particles < 100 µm (which represent the highest fraction of particles transported in stormwater runoff [36]) by road side gullies in France does not exceed 15% [37].
The highest deviation is observed for the event of 15 November 2014, where modeling results suggest the probability of mobilizing higher loads of TSS. This is likely due to the high rainfall depth of this event, which suggests the shift of sediment loading towards high values. The simulated washed-off loads are driven by the runoff volume, however in real cases, the mobilized fraction of sediments during a rainfall event is mainly controlled by its particle size and the kinetic energy of the rainfall drop, which is represented by the rainfall intensity. Fine particles are the most susceptible to transport, and intense rainfalls generate higher contaminant yields while also enhancing the mobilization of coarser particles [23,36]. Introducing particle size more explicitly into the model and accounting for the impact of rainfall intensity in the sampling procedure can improve the estimated TSS loads. Rainfall profiles can be defined to classify the low and high intensity rainfalls and then higher EMC will be affected by the more intense events. Different classes of particle size can then be introduced into the model, with new modeling assumptions, such as: (i) total mobilization of fine particles could be assumed, and (ii) a percentage representing the fraction susceptible to mobilization for coarser particles could be defined based on rainfall intensity. efficiency of fine particles < 100 µm (which represent the highest fraction of particles transported in stormwater runoff [36]) by road side gullies in France does not exceed 15% [37]. The highest deviation is observed for the event of 15 November 2014, where modeling results suggest the probability of mobilizing higher loads of TSS. This is likely due to the high rainfall depth of this event, which suggests the shift of sediment loading towards high values. The simulated washed-off loads are driven by the runoff volume, however in real cases, the mobilized fraction of sediments during a rainfall event is mainly controlled by its particle size and the kinetic energy of the rainfall drop, which is represented by the rainfall intensity. Fine particles are the most susceptible to transport, and intense rainfalls generate higher contaminant yields while also enhancing the mobilization of coarser particles [23,36]. Introducing particle size more explicitly into the model and accounting for the impact of rainfall intensity in the sampling procedure can improve the estimated TSS loads. Rainfall profiles can be defined to classify the low and high intensity rainfalls and then higher EMC will be affected by the more intense events. Different classes of particle size can then be introduced into the model, with new modeling assumptions, such as: (i) total mobilization of fine particles could be assumed, and (ii) a percentage representing the fraction susceptible to mobilization for coarser particles could be defined based on rainfall intensity. The evaluation criteria obtained by comparing simulated and observed TSS loads are similar for the three models. The determination coefficients varied around 0.5, demonstrating a direct, yet fairly weak proportionality between the observed and simulated TSS loads. The values of the RMSE coefficients showed that the discrepancy between the observations and simulations for all three Figure 5. Boxplots of simulated wash-off load obtained with the three modeling approaches. The central red mark is the median, the edges are the 25th and the 75th percentiles and the whiskers extend to the extreme values that are not considered as outliers. The measured wash-off load is presented as blue circles. The evaluation criteria obtained by comparing simulated and observed TSS loads are similar for the three models. The determination coefficients varied around 0.5, demonstrating a direct, yet fairly weak proportionality between the observed and simulated TSS loads. The values of the RMSE coefficients showed that the discrepancy between the observations and simulations for all three models were low, as they fluctuated around 0.3. The performance of the third approach is slightly better than the other two.

TSS Dynamics
The dynamics of TSS replicated by the models showed that the performance of the models was not consistent for all events or for all models (Table 4). For some events, high values of determination coefficient are obtained, while for others, the variability in simulations and observations may be very different. The variability in the simulated pollutographs is generally low, since the range of variations of the RMSE coefficients is narrow. The pollutographs of the events for which the models had their best performance with respect to the C R 2 are shown in Figure 6. Visual assessment of model outputs reveals an important advantage of the applied approaches: the ability to replicate multiple peaks of TSS concentrations. This was a forgone result for the homothetic model, but not for the others. Observation of the pollutographs shows that, depending on the model and the rainfall, the simulated beam does not frame the entire observations with a varying level of coverage. However, the width of the beam is narrow, revealing small uncertainties associated with simulated TSS concentrations. Accounting for the variability of pollutant generation and transport mechanisms at the elemental scale of the UHE leads to small uncertainties regarding the estimate of the TSS concentrations dynamics at the outlet. The stochastic representation of the processes thus succeeded in assessing the spatial propagation of uncertainty throughout the catchment area. The stochastic SWMM approach gave the best performance among the developed water quality models, while the homothety and the M(V) curve approaches failed to perform accurately. All three approaches gave the same results with respect to the RMSE evaluation criteria, where the calculated CRMSE all fell within the same range. However, the determination coefficients show the limitations of the model in coping with the intra-event variability of TSS concentrations, particularly for App-2 and -3. Visual observation of the pollutographs shows that the simulated dynamics with the stochastic The stochastic SWMM approach gave the best performance among the developed water quality models, while the homothety and the M(V) curve approaches failed to perform accurately. All three approaches gave the same results with respect to the RMSE evaluation criteria, where the calculated C RMSE all fell within the same range. However, the determination coefficients show the limitations of the model in coping with the intra-event variability of TSS concentrations, particularly for App-2 and -3.
Visual observation of the pollutographs shows that the simulated dynamics with the stochastic SWMM model were in accordance with the general decreasing tendency of TSS concentrations. However, coverage of the pollutographs regarding simulations was variable. For the first event, the model failed to replicate the initial rise in the concentrations, but succeeded in replicating the second peak, as well as the small rise at the end of the event. For the other two events, the simulations replicated the highest peak occurring at the beginning of the rainfall event, but overestimated the concentrations, since the beam was higher than the measurement. Therefore, the decrease in concentrations was accurately replicated for the second event, while this was underestimated for the steady phase recorded for the third event. The pollutographs simulated with the homothetic model did not succeed in replicating the initial peaks, but were able to better cope with steady phase fluctuations, with an underestimation for the first event and an overestimation recorded for the second and third. The M(V) curve model completely failed in simulating the dynamics of TSS. The obtained pollutographs were constant during the majority of the event, with the appearance of a small peak at the end of two of the events. The most acceptable simulation is for the event of 15 November, as a small peak of TSS appears at the beginning of the pollutographs, which rapidly decreases into the constant phase.
The first peak fluctuations in the observed TSS concentrations dynamics was not replicated in either approach. While the peaks were mostly overestimated by the SWMM model, they were underestimated by the other two models. This results in misleading estimates of the total TSS load discharged, which are critical for the integration of further downstream models, thus affecting the choice of stormwater pollution mitigation strategy.
The simulated TSS concentrations are dependent on the EMC values and the wash-off parameters. The values of these variables are controlled by the probability distribution and the boundaries of the sampling space. The first implementation of this approach assumed that the variables controlling water quality are generated following a uniform probability distribution. Even though there is no straightforward justification to choose the probability distribution to fit the EMC, the lognormal distribution has always been preferred by researchers to describe the randomness of the EMC [38,39]. Assuming the log-normality of the EMC sample in the proposed stochastic approach incorporates more knowledge into the wash-off process, this enhances the selection of the EMC. As for the boundaries of the sampling ranges, they were based on literature and prior knowledge of the basin. However, the latter information was not available for ungauged catchments, which gave the modeler more freedom in setting range limits. In this case, the sampling ranges must be enlarged to encompass the most possibilities. This requires the generation of an important set of parameters, which is not possible without the implementation of an automatic sampling technique. For this, using a Monte Carlo for EMC and parameters iterations could be explored.

Spatial Variability of TSS Loads
The spatial variability of the simulated eroded TSS loads (mean value) using the first approach (App-1) for the UHE, roads, alleys and roofs is presented in Figure 7 for the rainfall event on 15 November, which had the best performance.
The TSS washed-off loads were first calculated for each land use at the scale of the UHEs, using the total runoff volume and the total EMC drawn from the uniform probability distribution within the predefined ranges. Then, the total TSS loads for each UHE were calculated as the sum of TSS loads generated by each land use. The values of TSS loads were then averaged for the ten simulations that were performed. This calculation was made possible given the configuration of the URBS model, which allows inclusion of the runoff volume not only at the catchment's outlet but also at the outlet of each UHE, while distinguishing the contribution of each land-use. The TSS washed-off loads were first calculated for each land use at the scale of the UHEs, using the total runoff volume and the total EMC drawn from the uniform probability distribution within The spatial patterns show distinct yielding of TSS loads for each UHE and land use within the UHE, confirming the capacity of the stochastic approach to account for the variability within every element. The standard deviation of the loads calculated for the UHE is 35 g/m 2 , while it is 20 g/m 2 and 6 g/m 2 , respectively for the roads and the roofs. On the UHE scale, higher loads were derived from parcels connected to the boulevard. The loads generated by the roads mostly ranged between 1.5 and 3 g/m 2 , regardless of the road fraction. However, assuming that all roads are subject to the same traffic density, larger areas must have a larger contaminant yield, just as areas in closer proximity to traffic must also have a higher contamination than isolated areas [40]. As such, the distance of the cadastral parcel from the adjacent road and the fraction of the adjacent road connected to the parcel should be incorporated into the sampling of the EMC, in order to give a more explicit account of the effect of traffic on pollutant generation. The relationship between these parameters and the EMC could be represented by accounting for an exponentially decreasing function with the distance to the street.

Conclusions
In this study, a stochastic approach for stormwater quality modeling was developed. The approach was integrated into the physically-based distributed hydrological model (URBS), and accounted for the spatial and temporal variability of the processes governing pollutant production and transport in an urban catchment. The methodology adopted consisted of generating the EMC and corresponding wash-off parameters at the beginning of each rainfall event for each land use type, using a uniform probability distribution. Three wash-off models were proposed, based on SWMM exponential equations, a homothetic relation between the concentration of TSS and the flow and finally, the M(V) relationship.
The interest of the developed approach compared with other modeling approaches is that it incorporates the inherent temporal and spatial variability of the build-up/wash-off processes, which is rarely considered in typical uncertainty assessment of stormwater quality modeling. Accounting for temporal variability makes it possible to distinguish the TSS yielding from storm events with different intensities and durations while accounting for the spatial variability allows the explicit incorporation of the effect of the catchment morphology on TSS yielding, which differs greatly between roofs and roads.
The results show that the proposed approach would benefit from further enhancements before being considered as a reliable simulation tool for TSS dynamics and loads. The systematic overestimated TSS loads by all three models would have a critical impact on stormwater pollution modeling, and on the development of mitigation tools. Only the exponential wash-off model led to intra-event fluctuations of TSS dynamics resembling the observed trend.
Author Contributions: This study was designed by S.A.A., F.R., C.B. and G.C. The manuscript was prepared by S.A.A. and revised by F.R., C.B. and G.C.
Funding: This study is supported by OPUR observatory. The authors gratefully acknowledge OPUR research program for financial support and the ANR project TrafiPollu for providing the data.

Conflicts of Interest:
The authors declare no conflict of interest.