Modelling stormwater runoff , quality , and pollutant loads in a large urban catchment

Identification of stormwater runoff, its pollution load, and their implications for land use is essential in implementing stormwater management strategies. Hydrologic modelling provides an opportunity to assess them at limited data resources. In this study, the stormwater management model SWMM5 is applied for model development for a large basin in Tallinn. A geographic information system tool is used for subcatchment delineation, identification of directly connected impervious areas (DCIAs), and preparation of catchment input parameters. The model is calibrated and verified using sampled storm events to estimate event mean concentrations and annual loads. The predictive capability of the model for quantity is good and for quality moderate. The findings from the model show the percentage of the impervious area in the large catchment to be low at 19.7%. Although DCIAs, in particular roads and roofs, have relatively smaller areas they significantly impact runoff production (up to 75%) and loads (up to 66% total phosphorus and 71% total suspended solids). The first flush at the beginning of runoff is less important in case of a low intensity of rainfall, but heavy rain and snowmelt generate substantial runoff and pollution loads. When grab sampling is applied, it should focus on the medium and large events within 6 hours of storm commencement in order to achieve better mass estimations.


INTRODUCTION
Stormwater runoff has induced flooding and degraded the receiving water bodies in many urbanized areas.Runoff reduction and the control of pollution loads at coastal areas has been a goal for the Baltic Sea states for a number of past decades (HELCOM, 2002(HELCOM, , 2007)).This is reflected in national strategies, acts, and regulations such as, in Estonia, the Estonian Water Act (Riigikogu, 2016), Tallinn stormwater strategy until 2030 (Tallinna Linnavolikogu, 2012), and Tallinn City development plan (Tallinna Linnavolikogu, 2013), among others.The estimation of actual stormwater quality and quantity has been a requirement in evaluating the compliance of stormwater management regulations and in implementing effective control measures, and it is an inherently difficult task when there is a lack of sufficient information (Chiew and McMahon, 1999;Park and Stenstrom, 2008).Although effective, extensive monitoring campaigns are not always feasible due to resource availability and associated uncertainties.In such a situation, stormwater modelling is a helpful tool, which uses limited data resources and can simulate time intervals beyond the monitoring period (Vezzaro and Mikkelsen, 2012).
The constituents of stormwater runoff mostly depend on the rainfall and catchment characteristics (Nazahiyah et al., 2007;Liu et al., 2012Liu et al., , 2013) ) as well as on the upstream human activities (Brezonik and Stadelmann, 2002;Park et al., 2009).Models are primarily underpinned by rainfall (rainfall intensity, rainfall duration, antecedent dry condition) and catchment characteristics such as land use type and surface imperviousness (Park et al., 2009;Liu et al., 2013).However, it is also necessary to consider urban forms such as the road layout and spatial distribution of urban areas, which have significant influences on pollutant generation (Liu et al., 2012(Liu et al., , 2013)).Moreover, mixed land use, typical in large catchment basins, produces a higher variability in stormwater quality, as there is a high interspersion of various land use types (Lee et al., 2009).Therefore, a large catchment behaves in a different way than a small urban catchment in runoff and pollutant loading.
One of the important catchment characteristics is surface imperviousness.A large proportion of impervious surface accelerates the runoff rate and produces shorter times of concentration or lag times as well as an increased peak flow and runoff volume; however, the proximity of impervious surfaces to a drainage system plays a significant role in these interactions (Shuster et al., 2005).Estimation based on total impervious area overestimates the quantity and quality of stormwater (Shuster et al., 2005;Ebrahimian et al., 2016) and a larger area of the catchment is more crucial to add up the amount.Effective impervious area (EIA) is used by many researches instead of total impervious area to effect closer results.For EIA directly connected impervious areas (DCIAs) need to be accounted because the part that is not directly connected to storm drainage has less impact on the output (Lee and Heaney, 2003;Ebrahimian et al., 2016).Identification of the potential DCIAs is essential in implementing the possible solution for reducing runoff and pollution load.
The reviews of the stormwater models by Elliott and Trowsdale (2007) and Jayasooriya and Ng (2014) found that the US Environmental Protection Agency's Storm Water Management Model (EPA SWMM) has an effective performance in simulating stormwater quality and quantity.It is a widely used water quality model that has the capacity for both single-event and continuous simulation in the prediction of flows and pollutant concentrations.Many studies have indicated that SWMM has reasonable accuracy when the model outcomes are calibrated and validated (Jayasooriya and Ng, 2014).There is broad applicability of SWMM for simulating runoff and quality dynamics: hydrology assessment for pre-and post-development condition (Jang et al., 2007), pollutant wash-off water quality analyses (Tsihrintzis and Hamid, 1998;Temprano et al., 2006;Lee et al., 2010), combined sewer overflow modelling and assessment (Zhang and Li, 2015), flood forecasting (Han et al., 2014), and stormwater treatment facilities modelling and assessment (Burszta-Adamiak and Mrowiec, 2013;Zhang et al., 2013).
Many researchers applied the model for assessing runoff and loads with calibration and verification (Tsihrintzis and Hamid, 1998;Temprano et al., 2006;Nazahiyah et al., 2007;Tan et al., 2008;Lee et al., 2010;Chow et al., 2012;Mancipe-Munoz et al., 2014;Rosa et al., 2015).These studies are mostly for small catchment basins.Tan et al. (2008) and Mancipe-Munoz et al. (2014) worked on comparatively large urban catchments but mainly focused on runoff calibration, although they worked also for continuous events.There are few studies on quality dynamics for large catchment basins.In Estonia, stormwater modelling can rarely be found.Hood et al. (2007) used SWMM to estimate the flow and pollution load of a moderate size subcatchment in Tallinn, but they did not evaluate its predictive capabilities adequately.Six small subcatchments with distinctive land use (transportation, residential, and commercial) in Tallinn were analysed for runoff and suspended solids by Koppel et al. (2014).The runoff dynamics of mixed land use can be different from single land use and can provide a resultant runoff coefficient that is different from the coefficient of single land use (Lee et al., 2009).
The studies by Hood et al. (2007) and Koppel et al. (2014) used rainfall data from Harku station, which is almost 5−7 km away from the studied sites, but no stations nearer the sites were available at that time.However, an analysis using distant stations can lead to unconvincing results.There is also room to look at the impact of directly connected impervious areas (DCIAs).Lee and Heaney (2003) modelled the hydrologic performance of DCIAs and reported that DCIAs are the main areas contributing to runoff and have a most pronounced effect on urban hydrology.The DCIA or the connectivity to the urban area at the catchment scale influences the hydrologic response (Yang et al., 2011;Burns et al., 2015).
This study is mainly focused on model development for the Mustoja basin in Tallinn.It aims to identify the sensitive input parameters and potential DCIAs that influence the runoff and loads in a large catchment basin.It will estimate event mean concentrations (EMCs) and mass loads and finally use these results to evaluate the practiced sampling campaigns.

Description of the studied catchment area
The study was conducted in the Mustoja basin, a large catchment basin in Tallinn City (Fig. 1).Tallinn is the capital of Estonia where approximately 32% of the population is centred.It is in the north-eastern part of Europe on the shores of the Gulf of Finland.It has a humid continental climate with moderately warm summers and cold snowy winters.Normally, average monthly temperatures range from -5.7 °C in February to 16.3 °C in July (EWS, 2015).Total annual rainfall in Tallinn is 704 mm and average monthly rainfall varies from 32 mm in April to 86 mm in August.Snow cover usually lasts from mid-December to late March.
The study catchment comprises approximately 10.24 km 2 (30% of the Tallinn area), which for the most part spreads over the Kristiine and Mustamäe districts in the upstream side.Runoff mainly flows through underground pipe networks and ditches to the downstream natural channel where three pipes under Marja, Haabersti, and Mustjõe streets intersect.Finally this water is discharged to Kopli Bay and the Baltic Sea.The land use within the catchment is mainly residential with detached houses and apartment buildings in the upstream side.Industrial and commercial areas are dominant features in the downstream side within the catchment.
In this study, EPA SWMM version 5.1 (SWMM5) was used to simulate stormwater quantity and quality.It is a comprehensive hydrological and water quality model used for single event or long-term events in urban areas (Rossman, 2010).The model is common throughout the world for planning, analysis, and design related to stormwater drainage systems in urban areas.
The US EPA provides this model and related graphical user interface free of charge.The SWMM conceptual model comprises four major environmental compartments: (i) the atmosphere compartment, e.g.rain gauge objects, accounts for precipitation and pollutants from the air; (ii) the land surface compartment, which is represented through subcatchment objects, models areas that receive precipitation and generate runoff; (iii) the groundwater compartment, e.g.aquifer objects, receives infiltration from the land surface and provides input to the transport compartment; and (iv) the transport compartment, whose routes flow from runoff source areas through a network of pipes, channels, etc. (Rossman, 2010).SWMM offers a selection of three different built-in infiltration models and flow routing methods.The infiltration loss on the pervious area was estimated using Horton's equation (Horton, 1933) and the runoff transport was computed using the dynamic wave routing method under the complete Saint-Venant equations that consider the backwater effects with 30 s time steps.

Model development
Model development can be divided into two steps.In the first step, runoff modelling is constructed, which requires information on catchment characteristics, conveyance system, rainfall, and infiltration.In the second step, the quality model is developed using the runoff model, which requires build-up and wash-off components.
Catchment characteristics are catchment area (A), catchment width (W), average slope (S o ), percentage of the impervious area (% imp), surface depression storage, and surface roughness.Surface depression storage includes impervious (D imp ) and pervious depression storage (D per ) while surface roughness includes impervious (N imp ) and pervious surface roughness (N per ).
Subcatchments were delineated using the drainage networks and surface slope.The drainage network was provided by AS Tallinna Vesi in digital format, which includes the location/elevation of stormwater pipes and manholes, pipe diameter, pipe material, and year of installation (Table 1).
The Mustoja catchment basin is mainly served with a separate storm drainage system of approximately 51 km of pipes with the diameter varying from 0.15 to 2 m and 4 km of drain ditches constructed within 79 years.The surface slope was determined using a 1-m resolution digital elevation map (DEM) provided by the Estonian Land Board.A geographic information system (GIS) toolbox was used to prepare the catchment properties.The Mustoja basin was divided into 378 subcatchments (see in Fig. 1) ranging in size from 0.06 to 23.8 ha with the slope varying in the range 0-11.1%.The width of each catchment was computed applying the longest flow path method within the catchment.The width ranged from 25.8 to 3002.8 m.
For estimating the percentage of the impervious area, the land use map for the year 2014 was obtained from the Estonian Land Board.According to the Estonian National Topographic Database, the Mustoja basin had impervious surfaces of almost 50%, most of which are roads, fields of production sites, and buildings (Table 2).Pervious surfaces were mainly formed by the green areas, which made up a quarter of the catchment area, and by forest and private yards with 1/5 of the catchment area.The total impervious area (TIA) does not contribute to the actual urban runoff.Instead, EIA  or the portion of TIA that is hydraulically connected to the storm sewer system is an important parameter in determining this runoff (Lee and Heaney, 2003;Ebrahimian et al., 2016).The DCIAs were identified and separated from the unconnected parts following the procedure explained by Lee and Heaney (2003).Field investigation, areal maps, online maps, and storm pipe details were used, and imperviousness was spatially analysed using GIS.Numerous land uses were grouped together to form a simple classification including residential, commercial, industrial, forest, water bodies, roads, and roofs (Table 3).Roads and roofs are potential urban land uses for runoff and pollutants (Ballo et al., 2009).Approximately 55.7% of the areas were determined as residential (R), 7% as industrial (I), 10.2% as roads (Rd), 23% as residential roofs (Rr), 3.9% as commercial roofs (Cr), and 0.2% as forest and water (Table 3).Of the total catchment area DCIAs were found to form nearly 27% (278 ha) of which TIA made up 79% (217 ha) and EIA 56% (156 ha).We estimated EIA because all the runoff does not enter the inlets despite the fact that DCIA only represents the land use connected to the drainage system.Nevertheless, when the directly not connected impervious area (DNCIA) that is not connected to the drainage system was considered, TIA was 44% (448.9 ha) and EIA was 19% (198.4 ha).The major part of EIA was DCIA roads, followed by DCIA commercial roof, industrial areas, and residential areas.

Rainfall details
Stormwater runoff estimation in the Tallinn region has usually been based on the Tallinn-Harku station, about 6 km from the pilot basin, where the meteorological observations started in 1805.The rainfall varies in time, space, and altitude among stations even within a catchment.Rainfall measuring equipment must be close to or within the catchment in order to take account of the climatic factors, characterization and modelling of the stormwater drainage system in a statistically reliable way (Löwe et al., 2014).Since 2013, AS Tallinna Vesi has installed nine other rainfall measuring stations.Within the Mustoja basin, in the downstream area near Saarma Street, the Department of Environmental Engineering of Tallinn University of Technology installed a tipping bucket rain gauge 'Saarma' in May 2014.The influential stations were identified as Tondi 90 and Saarma (see Fig. 1) applying the Thiessen polygon method (Fetter, 2001) using a GIS toolbox.The Tondi 90 station was within the basin in the upstream area near Tondi Street.
One minute rainfall data from these two stations were used for the simulation so that the stations close to the subcatchments could feed them the corresponding data.4).Over 2015, 144 to 156 events occurred with average rainfall intensity 0.1 to 4.8 mm/h, producing 523 to 550 mm of rainfall.When the intensities were averaged over the years, the intensity amounted to 0.5-0.7 mm/h.Most of the storms (90%) during these two years were below 1.4 mm/h.An extreme event during the record period reached the peak of 19.4 mm/h on 22 Aug. 2014.The average total event rains were small in range, 1.7−2.2mm for a duration of 6.1−8.1 h and inter-event time of 52.8−59.8h.

Storm events and sampling details
Three storm events were used for model calibration and four events for validation.The calibrated events were during 9 September 2014 (event 1), 22−23 September 2014 (event 2), and 4−5 December 2015 (event 3).The rainfall characteristics (intensity, duration, antecedent dry days (ADD), and total rain) during these events are presented in Table 5. Event 1 and event 2 were captured after a long dry period of 12.1 and 12.6 days; the former had intense rain with 2.2 mm/h of rain of short duration (2.3 h) and the latter had medium rain with 0.8 mm/h of a long duration (26 h).Event 3 was recorded during wet weather conditions with 1.5 ADD and was below 1.4 mm/h; this falls between the medium and 90 percentile across all recorded years.Event 1 did not have quality data; therefore, it was used only for runoff calibration.
The automatic flow measurement unit installed by AS Tallinna Vesi provided flow rate for each 10 min for event 1 and event 2. At the same time, Tallinn University of Technology installed a water level measurement gauge to form a discharge-rating curve.The time interval sampling approach called sampling approach 1 (SA1) was used to collect samples during the events.Grab samples were also taken two times a week and a total of 104 samples were taken during the period from 6/11/2014 to 16/12/2015, which is called random sampling approach or sampling approach 2 (SA2).On each occasion, the flow was measured along with instantaneous water level using a flow tracker.Samples were analysed by competent water quality laboratories certified by the Estonian Accreditation Centre (EAK), which consistently followed the standard EVS-EN ISO/IEC 17025:2006(EVS, 2015).Many parameters, e.g.temperature, dissolved oxygen, conductivity, pH, biological oxygen demand (BOD), total suspended solids (TSS), ammonia, total nitrogen (TN), total phosphorus (TP), chloride, oils, microbiological parameters, and heavy metals, were measured.Heavy metal samples were taken every week, and a total of 30 samples were determined for analysing the content of Cd, Cr, Cu, Pb, Ni, Fe, and Zn.In this study, TSS and TP parameters with a strong linear relationship (regression coefficient 0.72−0.86)were simulated and the Water Quality Laboratory at Tallinn University of Technology was involved in analysing these parameters.
Other input parameters for catchment properties were adopted from the range provided in SWMM User's manual (Huber et al., 1988), books (Bedient and Huber, 1988;Wanielista, 1990), and published papers (Temprano et al., 2006;Nazahiyah et al., 2007), e.g.0.3 mm for D imp , 2.5 mm for D per , 25% for zero depression storage, 0.01 for surface roughness coefficient for overland flow on impervious portion, 0.3 for surface roughness coefficient for overland flow on pervious portion, etc. (Table 6).The infiltration parameters (maximum infiltration rate, minimum infiltration rate, and decay constant) were often difficult to identify and were adopted from the range provided by Bedient and Huber (1988) and Nazahiyah et al. (2007) based on the soil type used.Pollutant buildup depends on the land use, dirt and dust accumulation, and ADD, whereas pollutant wash-off depends on land use, runoff rate, and removal efficiency.The build-up and wash-off parameters provided by Chow et al. (2012) for different land uses were similar to the values used by Koppel et al. (2014) in the local context, and those values were adopted to the associated land use types.

Sensitivity and model performance test
The robustness of the results from the model was analysed for sensitivity by varying the values of the input parameters.For example, % imp and width were changed within ±10% and other parameters within the range as presented in Table 6.Sensitivity coefficient (Sc), which measures the effect of change in one factor on another (James and Burges, 1982;Chow et al., 2012), was the indicator for this analysis.The model goodness of fit was tested with four types of indicators as used by Chow et al. (2012) and Koppel et al. (2014): correlation coefficient, relative error (RE), normalized objective function (NOF), and Nash-Sutcliffe coefficient (NSC).Correlation coefficient (CC) is a general statistical measure of the linear relationship between observed and predicted values.Generally, the best calibration requires CC to be as close to 1 as possible.The RE for flow and peak flow measures the average tendency of the simulated data to be larger or smaller than observed (Chow et al., 2012).The optimal value of RE is 0.0, but within ±10% can be acceptable with low magnitude values indicating accurate model simulation.Let O denote the observed value, S the simulated value, O i the ith observed and S i the ith simulated value, and n the total number of observations.Then RE is calculated as follows: NOF, which is the root mean square deviation (RMSD) normalized with the mean of the observed values, is expressed as follows: The optimal value of NOF is 0 but 0 to 1 can be acceptable for calibration (Kornecki et al., 1999).NSC is calculated by Eq. ( 3) (Nash and Sutcliffe, 1970). ) NSC ranges between 0 and 1.0, with NSC = 1 being the optimal value.Values between 0.0 and 1.0 are generally viewed as acceptable levels of performance.

Results of sensitivity analysis
In the studied basin, the sensitive input parameters were found to be % imp, W, D imp , and N imp (Table 7).Among them, % imp has a significant influence on runoff flow rate at Sc = 0.82 and on peak flow at Sc = 0.68.Catchment width has a moderate influence with Sc of 0.44 for flow rate and 0.36 for peak flow.Both D imp and N imp have a negative coefficient, which indicates that the output values will increase with a decrease in these input parameters.In this large basin, they are not highly significant being in the range from -0.019 to -0.029 for D imp and from -0.008 to -0.018 for N imp .However, D imp is comparatively more influential for initial peak flow, having Sc of -0.029, which indicates that the depression storages regulate the first peak in the hydrograph.Tan et al. (2008) and Chow et al. (2012) found a similar influence for % imp.However, the width of their subcatchments was more sensitive to peak flow than runoff depth, which was not the case in this study.Instead, the influence there was nearly equal on peak flow and flow rate.Unlike in their study, in our case the basin was large and the sensitivity to peak flow was not strong due to width and N imp .

Calibration and verification results
The model is evaluated for the predictive capabilities through four indicators using Eqs (1)−(3).The results in Table 8 and Table 9 show that the CC, NOF, and NSC are within the acceptable range for runoff quantity for all events, indicating the modelled and measured runoff to be in a good relationship.The model for runoff simulation is acceptable at three important indicators: CC is 0.87−1 i.e. close to 1; NOF is 0.1−0.3,i.e. between 0 to 1; and NSC is 0.3-1, i.e. close to 1.However, RE is beyond ±10% at flow error less than 20% and peak error is between −27.4% and +21.2%.Similarly, the indicators such as CC and NOF for quality simulation (Table 9) suggest that the predictive capability of the model is in an acceptable range being 0.4−1 for CC and 0.1−0.5 for NOF.Here NSC and RE are not in a good range; NSC is poor for TSS but it is good for TP at 0.4−0.7,while RE can go beyond the acceptable range for both quality parameters.Overall, the model is acceptable for runoff quantity, but it provides less accurate estimation for quality performance.Calibration through increasing the number of events can reduce the error to some extent.The calibrated results for runoff quantity are presented in Table 6 and for quality in Table 10.These are the best-fit values for the input parameters after testing the goodness of fit.The output fits at 0.9 factor of % imp (Table 6); as a result the ultimate value is 19.7%,where the runoff coefficients for different land use vary as follows: for R from 0.05 to 0.2, Rr from 0.15 to 0.20, RC from 0.2 to 0.26, I from 0.09 to 0.45, Cr from 0 to 0.56, and Rd from 0.03 to 0.57, depending on the proximity of the connection to the drainage system.Temprano et al. (2006) in residential areas of Spain and Hood et al. (2007) in the highly residential Tallinn area found lower % imp values of 15.9% and 7.2%, respectively.A relatively large pervious area within the subcatchments and the effect of mixed land use have resulted in lower imperviousness.The depression storage used in this model is the average of dry and wet weather, which is determined as 0.7 mm for the basin.It is higher than the values obtained by Hood et al. (2007) and Koppel et al. (2014), where they proposed 0.15 to 0.29 mm for the Tallinn area.Depression storage for impervious area is sensitive to the initial peak flow, and it varied from 0.3 to 1 in this study depending on the weather conditions.Chow et al. (2012) provided this storage for residential area 0.2 in wet weather and 0.8 in dry weather, while it can increase to 0.75 in wet and 1.05 in dry conditions for commercial area.Temprano et al. (2006) suggested for residential area a higher value of 2.5 mm.Therefore, the value obtained for this basin can be justified.Manning's roughness for impervious area is 0.0135, since the impervious surface characteristics within this area are concrete/asphalt paving and/or a gravelled surface (Huber et al., 1988).A similar value is also used by Koppel et al. (2014) and Chow et al. (2012) in their studies.
Build-up and wash-off parameters determine the amount of pollutants discharged to drains.The more impervious the area, the more wash-off coefficients were found.The build-up is slightly higher in the commercial area than in the residential area and the industrial area has a slightly higher build-up than the commercial area (Table 10).Our analysis is different from the previous studies listed in Table 10 in the aspect of land use separation.In our study, roads and roofs were separated to investigate their effect on build-up and wash-off components.Wash-off coefficients for DCIA roads and roofs are higher than for DNCIA, which indicates that they are crucial for pollutant load.Maximum build-up coefficients are relatively small compared to findings by Temprano et al. (2006) and Hood et al. (2007) and the model fit values are similar to the findings of Chow et al. (2012), suggesting that there are cleaning activities on the basin upstream.Street cleaning is active in the basin as it is one of the action plans of the stormwater strategy in Tallinn (Tallinna Linnavolikogu, 2012).In the model, the cleaning efficiency used is 30% to 50% for TSS and up to 90% for TP in DCIA roads at an interval of 7 to 14 days.
Four storm events − 08/06/2014, 06/11/2014, 21/05/2015, and 06/08/2015 − were used for model verification.Evaluation was performed by comparing the observed and simulated runoff volume and peak flow for quantity and observed/simulated event load and peak load for TSS and TP using their RE (Table 11).The storm events for this verification have few observations and are not suitable for use with other indicators.A single event of 06/08/2015 was recorded to verify the quality performance.For volume RE% is ±10% ranging from −9.6% to 6% and for peak flow RE% is nearly ±10% in a range from −11.0% to 7.8%.Such RE% ranges indicate that the model can sufficiently predict stormwater runoff.However, the quality prediction is moderate for TSS and weak for TP having RE% beyond ±10%.Therefore, the verification results also suggest that more events of quality observations are needed to calibrate and verify quality performance.Leecaster et al. (2002) proposed seven storms or ~50% of the storms to get annual concentrations within 10% uncertainty whereas Bertrand-Krajewski (2007) recommended at least ten storms to calibrate regression models with a smaller confidence interval.The higher the number of storms captured and calibrated, the lesser model error can be expected.

Stormwater pollution load: TSS and TP
The simulated event-based concentrations and loads are presented in Table 12 and in Fig. 3.The EMCs of TSS in event 2 and event 3 exceed the national stormwater limit value of 40 mg/L (Vabariigi valitsus, 2013).The literature review by Göbel et al. (2007) showed that EMCs variation is large and depends on urban forms or land use.In their study the TSS loadings ranged from 0.2 to 937 mg/L with roofs (13−120 mg/L) and highdensity traffic areas (99−937 mg/L) having high TSS discharges.Loadings of TP were in the range from 0.01 to 0.5 mg/L; 0.06−0.5 mg/L came from roofs and 0.23−0.34mg/L from traffic areas.In five small catchments with a total area of 116 hectares in a Polish city the concentrations of TSS varied in the range 1.8-736 mg/L (mean ~31 mg/L) (Barałkiewicz et al., 2014).Compared to Estonian regulations, the concentrations of TP in our study were below the national stormwater limit value of 1 mg/L (Vabariigi valitsus, 2013) but higher than the limit of poor status of river quality levels of class IV (i.e.0.1 mg/L in annex 4, Regulation No. 59) (Keskkonnaminister, 2010).A study of three stormwater catchments in Paris resulted in the TP content range of 0.3−3.52 mg/L (Zgheib et al., 2012) and in a Polish city the TP content varied from 0.02 to 0.57 mg/L (mean ~0.17 mg/L) (Barałkiewicz et al., 2014).The EMC results of the Mustoja basin fell in the range provided by different studies, although simulated TP loads are higher than the usual range.Illegal discharges into the sewerage system can be suspected in the basin, which probably attributed to the higher TP loads.The peak concentrations were higher than the national limits: ~3 to 7 times as high as the limit for TSS and 1.4 to 2.7 times the national TP limits.The effect on runoff can be observed for up to 19.7 h for event 1 where a single rainfall occurred.The baseflow in this event has a great influence on increasing the volume of runoff (Fig. 3), accounting for approximately 67% of the total runoff.Therefore, the duration of runoff is crucial when estimating event total volume.The stormflow volume showed an increasing tendency proportional to total rainfall whereas stormflow mass load is likely to increase in proportion to rainfall intensity.In all three events, stormflow was more polluted than baseflow, making up nearly 90% of the total load for TSS and TP, which do not respond to the volume of runoff whether it is higher or lower than the baseflow.When simulating continuous rainfall for the years 2014 and 2015, the evaporation loss was also considered according to the daily average temperature obtained from the Harku station.Water from the surface is continuously vaporized during the dry period.Annual evaporation was estimated as 0.6 ML (0.014%) in both years.The annual outfall loadings were found to be 97.8 tonne TSS and 1.5 tonne TP from 4400 ML of runoff in 2014 and 110.7 tonne TSS and 1.7 tonne TP from 4500 ML of runoff in 2015 (Fig. 4).The simulated SWMM results were compared with the results from three monitoring programmes.The first one is the monitoring programme during 2012−2014 conducted by the Department of Environmental Engineering of Tallinn University of Technology (DEE) and AS Tallinna Vesi at the outlet of the Mustoja basin approximately 500 m downstream from the studied outlet.This programme was commissioned by Tallinn Environmental Board and the samples were taken six times a year.The same methodology was used in the second monitoring programme but it was conducted by the Estonian, Latvian, and Lithuanian Environment Group (ELLE) in 2015.The third monitoring programme was SA2 conducted by DEE, which was different in methodology in terms of sampling interval as the samples were taken twice a week.The total annual runoff was calculated based on daily average flow, and the total annual loads were calculated based on mean flow and mean concentrations of pollutants (Fig. 4).
In the simulation, the system rainfall was a combination of rainfalls from the Tondi 90 and Saarma stations while DEE, ELLE, and SA2 used rainfall from the Harku station (Fig. 4).Compared to rainfall from the Harku station, system rainfall is by 25% lower for 2014 and 11% lower for 2015.However, the modelled runoffs differ from the DEE and SA2 runoffs to some extent.The modelled runoff is nearly 10% higher than DEE in 2014 and 8% lower than SA2 in 2015.On the contrary, the ELLE runoff is significantly higher at 54% than the simulated runoff, although the annual rainfall is just 11% higher.One reason could be the method of calculation and another reason could be errors in measured values, which can make a difference in the annual runoff and loads.In the ELLE measurement, higher measured flow rates and concentrations have resulted in higher runoff and loads.ELLE TSS is enormously high being about 400% the simulated runoff.However, TSS loads from DEE and SA2 are 35% and 60% lower than the SWMM results for 2014 and 2015.There are quite significant differences in TP loads: 60% between DEE and SWMM, 59% between SA2 and SWMM, and 43% between ELLE and SWMM.These deviations can be attributed to the load calculation method and further to the mean concentration and mean flow because the error in the means could magnify the estimates.Overall, the modelled runoff and loads are higher than estimations from DEE and SA2, but they are significantly lower than ELLE measurements.

Main contributing impervious surfaces
In the Mustoja basin, the directly connected impervious area or DCIA was determined to be 26.8% of the total land use area (see Table 3).Such percentage of DCIAs was found to be quite effective for the runoff and quality output at 80.1% for peak flow, 75.1% for total runoff volume, 70.5% for TSS, and 66.1% for TP (Table 13).The road area and commercial roof occupy nearly 77% of the total EIA.Thus the DCIA roads and roofs have a higher contribution to the runoff and pollution load, even though the overall runoff coefficient was found to be less at 0.18.Most of the areas in the basin are either pervious or not connected directly to storm drainage.These areas constitute about 73.2% of the total area.

First flush effect
It would be interesting to examine the first flush phenomenon in the basin if it exists, in order to control the initial contaminated portion through isolation or diverting the stormwater from the road or roof surfaces to treatment facilities.The presence of the first flush phenomenon was assessed by plotting the cumulative fraction of the pollutant load against the cumulative fraction of the runoff volume (Bertrand-Krajewski et al., 1998) as in Fig. 5.The higher loading during storm runoff suggests the presence of first flush (Lee et al., 2004).To account for the intensity of the first flush, the pollutant load swept along by the first 30% of the volume was measured (Temprano et al., 2006;Nazahiyah et al., 2007).In Fig. 5, four events (event 2 to event 5) were used for TSS and TP loadings.One event was during the snowmelt period (event 4: 26/01/2016 to 28/01/2016) and another immediately after the snowmelt (event 5: 29/01/2016 to 01/02/2016).Event 2 and event 3 had no influence from snowmelt, as the former had antecedent dry condition of 12.6 days and the latter had antecedent wet conditions of 1.5 days.The pollutant loadings were 38%, 28%, 50%, and 39% of TSS and 45%, 36%, 45%, and 33% of TP swept by 30% of the runoff volume in events 2, 3, 4, and 5, respectively.The degree of the first flush for event 2 is not high, being 38% for TSS and 45% for TP, but the deviation from the diagonal line in Fig. 5 is evident after 40% of the runoff volume.It suggests that the flushing of the pollutant load is higher at a later stage within the last 60% of the runoff.The number of ADD was almost 13 for this event but the intensity of the rainfall played a greater role than ADD, as a similar result was obtained for TP in event 2 but the deviation is smaller.The first flushes at the initial portion of runoff for events 2, 3, and 5 are not significant.
The maximum intensities of rainfall within the 30% of runoff for these events are identified as small at 0.7 mm/h, 1.7 mm/h, and 0.9 mm/h.Instead, the latter part within the last 60% of the volume is more important for the pollutant load.Similar findings were obtained by McCarthy (2009), who found the first flush at the end of the event.Moreover, in a relatively pervious area, Maestre and Pitt (2005) did not observe any first flush.However, Lee et al. (2004) suggested that a seasonal first flush can occur in most cases.Our study also showed that the snowmelt event had a greater influence on the first flush at 50% for TSS and 45% for TP, indicating the effect of a seasonal first flush.Before event 4, there was an extensive period of temperatures below zero, and pollutants accumulated with the snow packing, which was washed off after the temperature became positive during this event.

Comparison of loadings from four sampling approaches
Finally, the sampling approaches are evaluated for the mass estimations, which provide information for choosing the appropriate sampling option.Four sampling approaches: time weighted sampling (SA1), random grab sampling (SA2), grab sampling within 6 h irrespective of storm size (SA3), and grab sampling within 6 h of medium and large storms (SA4) were analysed for annual average flow, annual average pollutants concentration, and annual pollutants load.Samples from SA3 and SA4 are formed from the data set of SA2 using the recorded time of sampling and rainfall time.Volume weighted mass load estimation from time weighted sampling is considered as the base load because the estimation has a smaller error (Leecaster et al., 2002) compared to other sampling approaches.SA2 is a random sampling that does not respond to the time of corresponding storms.SA3 does not take account of the influencing storm size but is taken within 6 h of the commencement of rain.SA4 is the sampling approach recommended for grab sampling by Lee et al. (2007) and Khan et al. (2006), in which grab samples were taken within 6 h of medium and large storm events.The volume of total runoff is calculated based on the runoff coefficient, catchment area, and runoff depth.Comparison of these four approaches with simulated output is presented in Fig. 6.The sampling method SA1 with volume weighted estimation is close to the simulated flow rate and annual loads, which indicates that time weighted sampling with volume weighted calculation as specified by Leecaster et al. (2002) provides a smaller error.It is considered as the actual load in the analysis of sampling approaches.The uncertainty of the model in quality estimation has probably built such an error and difference in TSS and TP from SA1.In the case of flow rates and TSS loads obtained from the sampling approach SA4, they are close to the estimations from SA1 and SWMM simulation at around 350 L/s and 60 t/yr, although TP loads have double the deviation.Random sampling or SA2 has often estimated low values of less than half the flow and 2/3 of loads.This is due to the fact that small rainfall or baseflows were mainly sampled in this approach.SA3, which is taken irrespective of the storm size as usual practice, has all of its mean outputs in the middle range of SA2 and SA4.SA4 with mean estimation, on the other hand, is nearer to the actual flow, TSS, and TP.This approach has limitations in identifying medium and large storms.It is difficult to predict the storm size before the rainfall ends.Alternatively, SA3 samples can be used to overcome the limitation after the storm details are retrieved.Therefore, SA3 is a practical and suitable sampling method, which assists in determining the samples of medium and large storms within 6 h after the start of rainfall, although it requires rainfall information and a greater number of samples than SA4.

CONCLUSIONS AND RECOMMENDATIONS
The model was developed for a large basin in Tallinn and was calibrated and verified for the observed storm events.Additionally, sensitivity analyses for finding sensitive parameters, imperviousness identification, first flush study to search for a possibility of controlling the initial portion of runoff, and evaluation of sampling approaches were carried out.The following conclusions were drawn. Sensitivity analysis showed the model to be sensitive to the percentage of the impervious area for predicting both flow rate and peak flow.Impervious depression storage regulates the initial peak flow.Impervious surface roughness and width of catchment have weak connections to the model predictions. For the studied large basin, the percentage of the impervious area was found to be 19.7%where the runoff coefficients for different land use vary from 0.05 to 0.57 depending on the proximity of the connection to the drainage system.The overall imperviousness is relatively low, depicting a large basin with mixed land use, which has a high impact on reducing the runoff coefficient.The average depression storage for dry and wet weather was used and the effective value was found to be 0.7 mm. The duration of runoff is crucial when estimating event total volume.Stormflow volume has an increasing tendency with total rainfall, whereas stormflow mass load is related to intensity.In the events, pollution is attributed more to stormflow than baseflow. The DCIA is an important factor for impervious estimation.Roads and roofs, which are directly connected to storm drainage, are crucial elements of DCIA and these impervious areas can contribute up to 75% of runoff and 66−71% of load. During low-intensity rainfalls, the first flush was found to be less effective.The first flush can occur later during 60% of the runoff volume.Therefore, the implementation of treatment facilities to control initial runoff and pollution load may not be effective and more research on the first flush of large basins is required.It was found that the first flush effect during the snowmelt period existed, and its impact cannot be ignored.It will be interesting to simulate snowmelt once there are sufficient input parameters and observations.
 A reliable sampling approach at limited resources is grab sampling within 6 h of storm events, as it provides results close to the simulated results and those obtained by the time weighted sampling method.This approach should focus on capturing medium and large storm events.Overall, the model development provided information about catchment's stormwater dynamics even at limited resources.The model developed in SWMM had a good performance for quantity estimation, but the quality results were of moderate accuracy.Calibration and verification using more storm events will increase the accuracy.Its applicability will increase even during wintertime once snow period simulation is included.Nevertheless, this study has found crucial impervious areas, percentage of the impervious area, runoff, EMCs, and pollution loads, and evaluated the sampling approaches used.The information will be helpful in planning pollution reduction strategies, implementing pollution control facilities, and designing monitoring programmes in the field of stormwater management.
The average monthly temperatures (in Harku station) and rainfalls recorded in three stations (Tondi 90, Saarma, and Harku) during Jan. 2014 to Feb. 2016 are presented in Fig. 2. August 2014 and July 2015 were the wettest months while February to April were dry months in 2014 and 2015, and October 2015 was the driest month.Compared with the normal rainfall pattern in the Harku station over 1981−2010 (EWS, 2015), these two years were dry years with approximately 34% and 26% less rain in 2014 and 2015, respectively.The storm events recorded from May 2014 to Dec. 2015 in the Tondi 90 and Saarma stations were similar in rainfall characteristics (Table

Fig. 2 .
Fig. 2. Monthly temperatures and monthly rainfalls in three rain gauge stations and long-term average rain.

Fig. 5 .
Fig. 5. First flush phenomenon for TSS and TP.The data above the diagonal line indicate higher loading during storm runoff.

Fig. 6 .
Fig. 6.Comparison of simulated annual average flow and pollutants loads with four different sampling campaigns.

Table 1 .
Pipe network details in the Mustoja catchment basin

Table 2 .
Land use details within the Mustoja catchment

Table 3 .
Simplified land use classified into DCIA and DNCIA of TIA and EIA

Table 4 .
Rainfall characteristics in 2014 and 2015

Table 5 .
Characteristics of sampled storm events

Table 6 .
Calibrated values of input parameters for simulating runoff quantity

Table 7 .
Sensitivity coefficient for input parameters

Table 8 .
Goodness of fit for runoff quantity simulation

Table 9 .
Goodness of fit for quality simulation

Table 10 .
Calibrated build-up and wash-off input parameters

Table 11 .
Verification of model performance

Table 12 .
Event-based runoff and concentrations of TSS and TP Fig. 3. Event-based baseflow and stormflow loadings.

Table 13 .
Contribution of DCIA to impervious land use