Variability of green infrastructure performance due to climatic regimes across Sweden

In the context of increasing urbanization and global warming, there is a growing interest in the implementation of green infrastructure (GI) across different climates and regions. Identifying an appropriate GI design criteria is essential to ensure that the design is tailored to satisfy local environmental requirements. This article aims to compare the hydrological performance of GI facilities in eleven Swedish cities by isolating the effect of climatic conditions using an identical GI design configuration. Long-term simulations based on 23-years of meteorological time-series were used as inputs for the Storm Water Management Model (SWMM) with Low Impact Development (LID) controls representing two types of facilities: a biofilter cell (BC) and a green roof. (GR). Large differences in potential annual and seasonal runoff retention were found between locations, driven mainly by the extent of winter/spring season, and the distribution of precipitation patterns (for BCs) and the sequence of rainy days-dry periods and evapotranspiration rates (for GRs). Winter/spring and summer demonstrated the highest/lowest differences between the seasons, results that suggest that implications for design might be aligned to the spatio-temporal distribution of precipitation patterns, and runoff regimes generated by snowmelt and rain-on-snow events, in locations where snowmelt represent high portion of runoff generation.


Introduction
Full-scale pilot tests and wide-spread implementation have demonstrated the potential of urban green infrastructure (GI) for coping with increased stormwater runoff caused by urbanization, limiting both runoff flows and volumes (Marsalek, 2006;Hatt et al., 2009;Nawaz et al., 2015;Harper et al., 2015;Pennino et al., 2016;Blecken et al., 2011;Klemm et al., 2017).GI aims to manage stormwater flows and volumes using facilities that mimic natural, pre-development hydrological cycle through a combination of detention, infiltration, and evaporation processes (e.g., biofilter cells and green roofs) (Coffman et al., 2004;Matlock and Morgan, 2011;USEPA, 2016;CIRIA, 2017).Despite the growing trend of using GI practices and the extensive research in this field in subtropical, mediterranean and temperate climates (Ahiablame et al., 2012), there are still significant knowledge gaps concerning their long-term performance, particularly in cool-temperate and cold regions (Nordberg and Thorolfsson, 2004;Fach et al., 2011;Muthanna et al., 2013;Kratky et al., 2017).
Across the world, recommendations and guidelines support and facilitate the design of GI; such guidelines exist at a variety of geographic scales, sometimes applying to a particular local context and at others applying to vast geographic areas (Sage et al., 2015).If the spatial resolution of design guidelines is defined at a national level or to an extensive region, it may result in facilities being designed based on what are effectively different hydrological criteria, making it difficult to predict real performance simply by transferring guideline recommendations from one site to another, or by assuming a linear relationship to scale the design.
Indeed, hydrological field studies of GI have reported large variations in retention capacity of GI, explained not only by design configuration, e.g., geometry, substrate and vegetation (Stovin et al., 2015;Johannessen et al., 2018;Berretta et al., 2018) but also by ambient conditions, which vary by location, such as precipitation patterns and climatic conditions (Stovin et al., 2012;Stovin et al., 2015;Sims et al., 2016).In cold and cool-temperate climates, GI performance varies seasonally due to reduced infiltration capacity assoicated with soil frost (Muthanna et al., 2013;Kratky et al., 2017), reduced or prevented by soil freezing, plant dormancy and snow accumulation and melting, snow cover properties and distribution, movements through snowpack, snow-soil interactions, and routing of generated snowmelt runoff (Pratap Singh, 2001;Lönnqvist et al., 2021a,b), making design in these regions particularly challenging (Caraco and Claytor, 1997;Semadeni-Davies, 2003).Thus, cold climate conditions create an additional complexity in the design of GI facilities, influenced by freeze-thaw mechanisms with respect to infiltration capacities and hydraulic retention (Muthanna et al., 2010;Kratky et al., 2017).
In Sweden, which stretches over temperate and continental climate zones (Cfb, Dfc and Dfb) (Peel et al., 2007), national recommendations for the design of stormwater control measures contain criteria for media composition and facility sizing based on both field monitoring and laboratory experiments (Svenskt Vatten, 2016;Larm and Blecken, 2019).Although a single set of guidelines applies to the country as a whole, it has not yet been confirmed that homogenous performance can be expected for facilities located throughout this large geographic area.At country scale, national depth-duration-frequency (DDF) statistics were long based on the approach of Dahlström (1979).Currently, general recommendations provide functional requirements for dimensioning and design of stormwater systems based on statistics of historical rainfall data and future projections (Svenskt Vatten, 2011b).In contrast, a study conducted by Olsson et al., 2017aOlsson et al., , 2017b) demonstrated a clear regionalization of high-intensity rainfall patterns, clustering the country into four regions: North (N), Central (C), Southwest (SW) and Southeast (SE).Given the regionalization of both rainfall and temperature patterns, as well as the clear link between these factors and GI hydrological performance, adapting regional water quantity management criteria or guidelines to a local scale requires variations with of GI functioning with climate patterns to be considered.The adoption of efficient GI solutions in urban areas located in cold climates (which account for 24.6% of global land according to Peel et al., 2007), strongly depend on the guidance and evidence given to practitioners to account for these management criteria.
In recent years, several physically-based numerical models have been created to predict GI hydrological performances (Broekhuizen et al., 2019).The introduction of computer models into urban drainage analysis allows more accurate computations of urban runoff, suited to the small temporal and spatial scales typical of urban and peri-urban applications (Fletcher et al., 2013;Moghadas et al., 2016).Similarly, the effect of frozen soils and rain-on-snow events on the reduction of infiltration capacities and increasing runoff volumes and its link to the GI performance was described in numerous earlier studies using snowmelt modelling tools (Semadeni-Davies, 2003;Valeo and Ho, 2004;Dobre et al., 2017;Moghadas et al., 2018).The widespread adoption of accurate modelling tools, with the advantage of a multi-scale analysis replication, allows scientific evidence of how GI facilities perform over time and under varying seasonal conditions to be provided (Moghadas et al., 2016a,b).Thus, the effect of freeze-thaw cycles and climate patterns on GI facilities will be investigated using a modelling approach.
Previous investigations have identified geographic variations in expected GI hydrological performance due to differences in rainfall and evapotranspiration patterns in humid/semi-arid continental climates (Stovin et al., 2015;Sims et al., 2016), as well as temperate oceanic and subpolar oceanic climates (Johannessen et al., 2017).However, to date, no studies have considered variations in expected long-term (i.e., annual or pluriannual) GI performance in regions with cold climates, where in addition to differences in precipitation and evapotranspiration, differences in temperature patterns influence cycles of frozen ground, snow cover, rain on snow and snowmelt.These conditions add additional complexity to GI performance: snow forms a stock of water on facilities during freezing conditions, which is later mobilized during warmer dry periods or rain on snow events (Muthanna et al., 2010;Kratky et al., 2017;Hamouz and Muthanna, 2019).In addition, the infiltration capacity of GI facilities has been shown to be reduced when soil is frozen (Moghadas et al., 2016;Paus et al., 2016;Balstad et al., 2017), making overflow more likely.While adapting design to local climatic condition in cold climate areas requires these mechanisms to be considered, accounting for them requires the use of additional models which entail substantial uncertainty.
Hence, in the present work, an innovative approach is presented, which combines deterministic and stochastic modelling tools (machine learning and probabilistic models) to represent variations GI performance across multiple Swedish urban areas influenced by cold-climate hydrological processes using historical meteorological data.The present investigation seeks to (1) provide a baseline for understanding the differences in the hydrological performance of a biofilter cell (BC) and green roof (GR) at regional and large-scale level and (2) evaluate whether variations in performance necessitate different design criteria for GR and BC facilities between regions in Sweden.The assessment of GI hydrological performance across different climate gradients in Sweden was conducted using the Low Impact Development (LID) component of the semi-distributed SWMM model, with a modelling chain accounting for the effects of frozen soils and rain-on-snow events.

Study area and methods
A rainfall-runoff model with long-term simulations was applied to estimate the volume-reduction performance of a BC and a GR for eleven different urban areas in Sweden.The model aims to isolate the effect of climatic conditions using an identical design configuration and operational conditions.To do so, identical models were forced with historic meteorological data from each location to evaluate runoff volume reduction and its seasonal variation for each site.An uncertainty analysis on model outputs was conducted by identifying the model input design parameters subject to uncertainty and quantifying confidence intervals of the hydrological performance estimated.Local climate indices were defined as a basis for discussing limitations and potential recommendations for the design of infiltration facilities under different Swedish climatic regimes.

Site descriptions
Eleven urban areas across Sweden, representing a span in seasonal temperature and precipitation patterns along a north-south gradient (55-68 • N) and a coastal-inland gradient were selected (Fig. 1, and Table 1, in appendix in the supplementary material).The locations include the most densely populated areas of Sweden and represent a variety of climatic conditions, ranging from temperate to cold climate according to the Köppen-Geiger climate classification (Peel et al., 2007).The northernmost locations (Kiruna, Luleå, Umeå, Östersund and Ornskoldsvik) are classified as sub-arctic climates (Dfc), and urban areas in the southern part represent a humid continental climate (Dfb in Gävle, Stockholm, Eskilstuna, and Örebro) and an oceanic climate (Cfb in Gothenburg and Malmö).Based on the clustering results obtained by Olsson et al., 2017), a division of Sweden into four regions was used for a regional analysis of the urban areas.The regional classification of Sweden corresponds to a subarctic climate (Dfc) for the N and C regions, and both humid continental and oceanic (Dfc and Cfb) for some areas in the SE and SW regions.

Input data
Available historical meteorological data between 1997 and 2017 were retrieved from the national network operated by the Swedish Meteorological and Hydrological Institute (SMHI).The time-series are quality-controlled in the database system of SMHI and have an overall coverage of around 90% of the complete recording period.Values were not available for the exact position of Kiruna, Umeå, Östersund and  Orebro, so the nearest available stations (within a radius of less than 30 km from the urban area) were used (Table 1, annexe 1 supplemental material).The network of automatic stations consists of Geonor weighing gauges, which measure precipitation with a time-step of 15 min.Air temperature and relative humidity are measured, respectively, by a resistance thermometer and a hygrometer, with reported hourly values, and wind is measured with an anemometer which measures the highest instantaneous recorded value during an interval period of 1 h.Daily potential evapotranspiration (PET) was estimated using the FAO Penman-Monteith method.PET was estimated for the period 1997-2019 as a function of daily maximum and minimum temperature ( • C), daily wind speed (km /hr), and dew point ( • C).Due to the absence of daily solar radiation data for the period of analysis, an estimation was derived from air temperature differences (Hargreaves radiation formula), and its relation to the degree of cloud cover in each location (Allen et al., 1998).In addition, soil temperature measurements were collected from the integrated monitoring system in natural ecosystems controlled by the Swedish University of Agricultural Sciences (SLU).Data from Aneboda, Gammtratten, Gårdsjön and Kindla stations were used at shallow, intermediate, and deep soil layers (5-10 cm, 25-30 cm, and 58 cm respectively) with a recording period between 1996 and 2018 for most stations.

Model configuration and facilities design
The semi-distributed Storm Water Management Model (SWMM) was used for long-term simulations of runoff quantity using rainfall/runoff and snowmelt processes with a 15-min time-step.GI facilities are conceptually represented in SWMM by separated horizontal layers (from top to bottom), where water fluxes between the layers are described by flow continuity equations.The design variables that affect the hydrologic performance of LID controls include the properties of the media (soil and gravel) contained within the unit, the vertical depth of its media layers, the hydraulic capacity of any underdrain system used, and the surface area of the unit itself (Rossman and Huber, 2016a).
For BCs, the model represents the facility with three horizontal layers: a surface layer, a soil layer, and a storage layer beneath the soil, and assumes a sealed bioretention facility with an impermeable liner around the base and sides, and a perforated underdrainage pipe (Fig. 2).BCs were modelled as part of the sub-catchment area (3.2% of the total area), and the baseline scenario was defined as a catchment with 100% imperviousness.The design configuration of BCs (Table 1) were based on recommendations for design of infiltration facilities in Sweden (Larm and Blecken, 2019).GRs are represented similarly to BCs, except that they use a drainage mat instead of a gravel aggregate in its storage layer, no underdrain, and no inflow from the contributing catchment area (Fig. 2).GRs were modelled as 100% of the sub-catchment area covered by the facility with dimensions (Table 2) selected according to general guidelines for GR design (Capener et al., 2017).The dynamic-wave routing method for overland flow calculations were used, which solves the complete one-dimensional Saint-Venant flow equations (Rossman and Huber, 2016b).The Green-Ampt infiltration method was used to model infiltration of rainfall into the upper soil zone of both BC and GR facilities.

Urban snowmelt modelling
In cold climates, the distinction between rain and snow is important to consider due to the different contribution to the runoff generation.The type of precipitation is directly determined from the model inputs by using a dividing temperature factor, which is 0 • C (Rossman and Huber, 2016b).The simulations all begin during winter (on January 1); therefore, as an initial condition, a homogenous depth of snow, represented as snow water equivalent (SWE), is applied to the entire area of the catchment, based on historical weather observations (SMHI meteorological network) of mean snow depths during winter periods for each urban area.Initial SWE values were considered by assuming that snow covers 100% of the catchment area and represent the average of total snow depth over the entire catchment (Moghadas et al., 2018).SWMM accounts for snow redistribution by assigning fractions of snow moved out of the catchment, moved from impervious to pervious surfaces in the catchment and converted into immediate melt.Due to the absence of snow management practice statistics for each study location, snow removal parameters in SWMM were based on a survey of Swedish cities (Bengtsson and Westerström, 1992) (see Table 3).
Snowmelt is simulated in SWMM using an energy balance method based on Anderson's U.S National Weather Service procedures (1973) during rainfall periods and a degree-day method during dry weather.The latter uses a linear function of the difference between air temperature, T a , and a base temperature, T base (0 • C), using a melt rate factor  Depending on the location Fraction transferred out to the watershed 10% Fraction converted into immediate melt 25% 1,2 Melt coefficient values follow a sinusoidal curve according to the calendar day (Rossman and Huber, 2016a).
1 Values according to site-specific conditions for each urban area (Wern, 2015).2, 3 Values taken from observed distribution of snow in Luleå, Sweden (Bengtsson and Westerström, 1992).
I. Mantilla et al. (Rossman and Huber, 2016b).Melt coefficient parameters used in the energy balance method are described in Table 2. Previously calibrated melt-rate factors from Luleå, Sweden, which vary throughout the year to account for variations in solar radiation, were chosen to account for properties of undistributed and ploughed/stored snow during the melting period (Valeo and Ho, 2004;Moghadas et al., 2018).The capacity of the snowpacks to hold liquid water during snowmelt was assumed to be 10%, based on field measurements of water movement dynamics on snow packs (Gerdel, 1945), and remained constant during the whole period of simulation.

Accounting for reduced infiltration capacity due to soil frost
To account for reduced infiltration capacity during winter due to soil frost, seasonal periods were defined at each location based on soil temperature and different saturated hydraulic conductivity (Ksat) factors were applied for each season.Ksat values were determined for each season based on observed infiltration rates and via infiltration test (Paus et al., 2016).Seasonal periods were defined based on daily average soil temperature, following the methodology defined by SMHI (Vedin, 1995).

Soil temperature prediction
As soil temperature data was not available for each urban area, an artificial neural network (ANN) model was used to predict soil temperatures at a depth of 10 cm below the ground.The ANN model was applied using antecedent soil temperature data (retrieved from Gammtratten, Kindla, Gårdsjön and Aneboda monitoring stations) with time lags of 1, 2 and 3 days and antecedent air temperatures (retrieved from the SMHI stations for each urban area) with time lags of 1 and 2 days as inputs.The selection of stations to apply the ANN model was based on the closest geographical location to each urban area.A spatial analysis method (Thiessen polygon) was applied by assigning an aerial significance to each station with soil temperature data.
A multilayer perceptron (MLP) model was constructed, and a backpropagation algorithm was used for training of the ANN time series, with a model structure that consists of three layers of neurons: (1) an input layer, (2) an output layer, and (3) and an intermediate layer (Bhagat et al., 2021).The time series were pre-processed prior to the development of the ANN models.Model inputs (soil temperature and air temperature) were segregated into a training set (70%) and a test set (30%) to improve the learning performance of the model by a constant update of weights and biases during a repetitive cycle of data training.Training and test sets of data were normalized into the range [0, 1] before the implementation of ANN models due to the use of a sigmoid as an initiation function, as is described in equation A. Sigmoid functions have a domain of all real numbers, with a common range between 0 and 1, and has a characteristic of monotonically increasing (Witten et al., 2017).
The process of model calibration or dataset training was based on a backpropagation algorithm by adjusting the weights and computing the error between the output and the corresponding target value and propagating this backward through the network to adjust weights and produce the desired output (Abdalla et al., 2021).The optimum number of neurons in the hidden layers was based on a trial error process, starting with a small number of neurons, which gradually was increased until obtaining the lowest forecasting error (mean square error function).The validation process was conducted to evaluate the fit between the ANN model outputs and independent dataset (sample of data used for the training).The applied predictive models were compared to the observed data quantitatively using performance metrics, including the goodness of best-fit (e.g., Nash and Sutcliffe), the absolute error measure (RMSE) and the correlation coefficient (r), and quantitatively, by inspection of output plots (Fig. S3a and Table S3a, annexe 3 supplemental material).

Definition of seasons -Ksat factors
The seasonal periods were defined based on the predicted soil temperature time-series for each location in the study.The boundary between the seasons was drawn when the parameter maintained a trend of continued falling or rising during a consecutive number of days (five to seven days depending on the season).To define the transition between the seasons, a constant trend (increase or decrease) in the temperature and number of consecutive days was used according to the SMHI seasonal definition.In the meteorological division of seasons, winter is defined as the period when the daily average temperature permanently is 0 • C or lower and summer as if it is permanently 10.0 • C or higher.Consequently, it is spring and autumn, respectively, when the average temperature of the day is greater than 0.0 • C but less than 10.0 • C. Since the arrival of spring is calculated over a period of seven days, while the arrival of summer is only estimated over five days in a row, the summer could theoretically arrive earlier than spring.To avoid this, a requirement that summer cannot arrive earlier than spring was added (Johannessen, 1970;Vedin, 1995).
To account for seasonal variations in the infiltration capacity of the facilities as a function of changes in soil temperature (i.e.soil freezing), hydraulic conductivity (Ksat) was varied based on previously reported experimental data (Paus et al., 2016;Balstad et al., 2017).Thus, summer months represent a scenario with natural conditions (100% of full soil infiltration capacity), and the remaining seasons, correspond to intermediate (spring, 80%, and autumn, 60%) and low infiltration capacities (winter period, 40%).Fig. 3 shows the monthly mean Ksat (mm/hr) values used to account for seasonal variability by applying factors to the LID facilities in SWMM.There is a reduction of Ksat during winter (i.e., December to February) as low soil temperatures and/or frost reduce the infiltration capacity in the soil media during the colder months).Depending on the soil moisture content during freezing, a formation of porous/concrete frost may reduce the soil permeability.Ksat values for the summer months represent the infiltration capacity of the soil media in natural state.

Data analysis and interpretation
To calculate hydrological performance the results of a scenario of a 100% impervious catchment were compared with/without the implementation of the facilities (a BC and GR).The period of 1997-2019 was used to evaluate the performance in terms of percent runoff volume reduction (ROR) at urban and regional scales.The hydrologic performance is estimated according to Eq.A at different time scales, including annual (ROR AN ) and long-term volumetric retention across the full multi-year simulation period (ROR LT ).In addition, seasonal volume reductions were calculated to highlight the dependency of performance on air temperature and freeze-thaw cycles (respectively ROR sp , ROR su , ROR au , and ROR wi for spring, summer, autumn, and winter).To evaluate the hypothesis of significantly different performance between the urban areas, a non-parametric statistical hypothesis test (paired Wilcoxon test with a threshold of P = 0.05) was conducted based on the results for ROR AN and ROR LT .

Uncertainty analysis
Uncertainties in volume reduction were quantified using a Monte Carlo (MCM) technique for both BCs and GRs to quantify and analyse model output variability caused by the assumed input parameters (criteria for the design configuration of the facilities).The analysis was conducted for each urban area including 1000 iterations for each year of the period covered , ending up in 23,000 simulations in total; afterward the results were aggregated at a regional level.
The most influential LID parameters for the water balance calculations identified from literature were selected for uncertainty analysis and combinations built by sampling random values for each parameter from uniform distributions were conducted (Table 4.).The maximum melt coefficient, fraction free water capacity (FWF) and snow water equivalent (SWE) were selected to account the uncertainty in snowmelt runoff processes (Semadeni-Davies, 2003;Moghadas et al., 2018).In addition, snowmelt is simulated in SWMM with a variant of degree-day melt-rate coefficients depending on land-use characteristics, timing in the melt season, and location.The relevant values were chosen to account for different properties of undisturbed and ploughed/stored snow during the melting period (Bengtsson and Westerström, 1992).
The potential to reduce runoff volumes by GRs are highly influenced by the inter-event conditions (initial soil saturation) and their capacity of storage restoration (Johannessen et al., 2017).Therefore, soil thickness (t) and porosity (p) were selected as the most influential parameters in the runoff retention capacity.(Leimgruber et al., 2018).These findings have an agreement with results obtained by Krebs et al. (2016), who identified the porosity as the most sensitive parameter in terms of total water holding capacity.Furthermore, GR facilities were found sensitive to drainage mat (vf), hydraulic conductivity (k), and conductivity slope, which were thus also considered for the uncertainty analysis (Leimgruber et al., 2018).

Comparison of climates across locations
Monthly precipitation values for cities located in the N, C and SW regions range from 20 mm in the driest months (February to April) to > 65 mm in the wettest months (August-October) (Fig. 4a).The SW region has the highest precipitation volumes for each month, with volumes that range from 43 mm (in April) to 87 mm (in October).In general, the largest amount of precipitation in all four regions typically falls during summer, with a dominant effect of convective and orographic rain processes occurring due to the action of high summer temperatures (Svenskt Vatten, 2011a); lower precipitation volumes fall during the spring season.Frequent precipitation events are more present in the SE, with a higher number of rainy days during the year.
Mean temperatures in the four regions vary along a north-south gradient, with the maximum temperatures of over 18.5 • C in the SW region (for July), and the minimum in February (− 9 • C) in the N region.In general, for all regions, the month with highest mean temperatures corresponds to July (>12 • C), and the lowest in January (<− 1 • C); during winter, temperatures are on average 10-15 • lower in the N compared to the SW region.The evapotranspiration rates have a normal distribution across the year, with higher daily PET values during the summer months caused by the higher solar radiation (Fig. 4b).Daily PET rates are highest in the SE region for the winter and spring seasons (>1 mm/day from March) and lowest for the N region in almost all months (except for July).During summer months (June and July), PET values exceeded 3.5 mm/day in all regions.
Finally, the greatest snow depths occurred during January to April for all the regions.On average, the greatest snow depth occurs in the N region (0.38-0.5 m during January to April), followed by the C region (0.18-0.22 m in January to March).The lowest snow depths were found in the SW region and the average snow depth was zero between June and September for all the regions.The N region has the longest snow cover of the regions (non-zero average snow depths from October to May), followed by the Central region.Table 5 provides climatic indices at annual and seasonal basis for all the urban areas.
On the other hand, the transitions between seasons are significantly faster in northern Sweden where the difference in solar radiation changes rapidly.Spring arrives later in the N and C regions, but no significant differences were found for the arrival of summer, except for Malmo where the soil temperature exceeded 10 • C for more than three months.The definition of seasonal periods was conducted based on the soil temperatures estimated from the ANN models.The periods for the winter, spring, summer, and autumn are described in Table 4.The shortest winter period was found in the SE and SW regions, with Goteborg and Orebro as the urban areas with fewest days with negative soil temperatures.

Overall and annual performance of BCs
Overall, the minimum and maximum ROR LT ranged from 67% in Kiruna to 90% in Östersund.Kiruna is the northernmost city considered in this study (located in the N region); while it has one of the lowest mean annual precipitations, which could be expected to lead to better performance, it also has one of the lowest potential evapotranspiration and the greatest number of winter days, during which infiltration capacity is assumed to be reduced.Gothenburg showed the second-lowest ROR LT (77%), being the urban area with most days with precipitation and the highest mean annual precipitation; among the key factors affecting the BCs hydrological performance, initial conditions (i.e., initial soil moisture) have been demonstrated to have a strong effect on BC retention capacities (Stovin et al., 2012;Jiang et al., 2019), suggesting that locations with a high frequency of rainfall events (more evenly distributed rainfall patterns) are likely to show lower ROR.LT.
Similarly, low temperatures during cold seasons have previously been associated with lower rain garden hydrological performance (Muthanna, 2010;Suihko, 2016), a result attributed mainly to lower infiltration capacities caused by soil frost.Malmo, which presented a ROR LT of 83%, is an urban area characterized by high-intensity rainfall events, but with fewer annual days with precipitation.It is somewhat surprising that Ostersund (located in the C region) had the greatest ROR LT, as it is an urban area with one of the coldest climates, with mild precipitation and evapotranspiration patterns.The remaining cities all fit into a relatively small range (80-85%) of ROR LT , with cities towards the southern part of the country as the urban areas with the highest performances (Stockholm, Eskilstuna, and Malmo).The analysis on an annual basis showed high inter-annual variability for the 23-year time series in the different urban areas, mainly due to differences in yearly climatic conditions.The lowest ROR A in any city was found in Kiruna (35%), results which agreed with patterns of longterm performance (ROR LT).Similarly, minimum ROR AN were found in years with considerable high mean annual precipitation (>670 mm) and mild potential evapotranspiration values (<430 mm); for example, the minimum ROR AN in Orebro and Gothenburg corresponds to a year with an annual precipitation of 1115 and 1386 mm (the highest values during 23-year period of analysis) and an annual evapotranspiration of 423 and 292 mm, respectively.On the other hand, all urban areas had a maximum ROR A exceeding 90%.Overall, locations with monthly average evapotranspiration >90 mm and few days with precipitation (annual days with precipitation <100) showed higher ROR LT.
In addition, Kiruna and Ostersund were distinguished by having the most and least dispersion in annual performance (σ = 12%, and σ = 4.9%).Four cities located in the same region (SW and C region) showed similar variability of ROR A : Gothenburg and Malmo (σ = 6.9% and 7.3%) and Ornskoldsvik and Gavle (σ = 9.1% and 9.3%).Besides, Kiruna, Gothenburg and Ostersund showed significant differences in terms of ROR AN with all the urban areas (Wilcoxon paired test).At a regional level, no significant differences were found between par of cities in the SE region (Eskilstuna -Stockholm, Orebro -Stockholm) or between a pair of cities located in the C region (Ornskoldsvik -Gavle).

Seasonal BC hydrological performance
Overall, the lowest seasonal ROR s were found in Kiruna during the winter season (35%), primarily driven by low soil temperatures, the extent of winter days (period with minimum Ksat values, and therefore a reduced BC infiltration capacity), and almost negligible PET rates (Fig. 5).The second-lowest ROR s were observed in Gavle during the spring season (46%), which could be expected due to the high number of freeze-thaw cycles, and therefore, an additional contribution of runoff due to snowmelt.Indeed, the joint occurrence of rain and snowmelt have been found to generate on average 5 times greater runoff volumes compared to the summer season in the city of Kiruna (Moghadas et al., 2018), conditions that create an extra runoff volume to be retained by the BC.Urban areas in the southern part of Sweden (located in the SE and SW region) had the lowest ROR s during the spring season (53% for Gothenburg and 45% for Stockholm); However, Orebro and Eskilstuna showed lower performance during the summer season (69% and 67%).
A higher number of significant differences (see Table 6) were observed between cities during the spring season (34), a period where the effect of rain-on-snow events and snowmelt regimes have been found previously to have a high impact on the BC hydrological performance (Géhéniau et al., 2015;Moghadas et al., 2016), as compared with other seasons (27,21,27 for Winter, Summer, Autumn, respectively), The fewest significant differences were observed during the summer season, a period when infiltration rates correspond to the natural conditions in the BC soil media, and the climatic conditions in terms of MMP and MPET are more similar between the urban areas.Moreover, the ROR s in Kiruna was found to have the most significant differences with respect to remaining urban areas (p-values < 0.05), followed by Ostersund by having significant differences during all seasons with Gavle, Stockholm, and Orebro.At a regional level, no significant differences were found between the four seasons in the SW and SE regions.Overall, performance did not appear to be grouped by the regionalization identified by Olsson et al. (2019) for high-intensity precipitation patterns.

Implications for BC design
The introduction of seasonal Ksat values (constant infiltration capacities for each season) showed that when initial infiltration rates are sufficiently high, it is possible to obtain satisfactory hydrological performances even during cold months.A potential recommendation is to design BCs with a minimum Ksat value of 10 mm/h during the cold months to reduce the effect of frozen soils as was pointed out previously by Paus et al. (2016); Sondre et al. (2018).This value is considerably lower than infiltration rates found during the summer period, as findings obtained from a survey of BCs in Sweden with Ksat values ranging between 300 and 3620 mm/h (Beryani et al., 2021).Similarly, as there are no Swedish design standards for Ksat, often initial infiltration rates of 12-300 mm/h is recommended according to Fassman et al. (2013).Moreover, the design of BCs for cold climates requires a trade-off between utilizing soil material with high infiltration rates to minimize the effect of the frozen soils, and still maintaining sufficient water quality improvement.Besides, an investigation of the metal removal effectiveness in replicate laboratory BC columns (at 2 • C, 8 • C, and 20 • C) showed that materials with high infiltration rates (up to 1000 mm/h) can also provide significant treatment of stormwater (Blecken et al., 2011).
Because coarse soil with low clay content drains more efficiently than finer soils, it is less likely to be saturated at the onset of freezing weather, making the formation of concrete frost in the BC, which is associated with the highest loss of infiltration capacity, less probable (Caraco and Claytor, 1997;LeFevre et al., 2009;Moghadas et al., 2016).Along these lines, the Minnesota Stormwater Manual (Gulliver and Anderson, 2008) recommends BC soil media containing less than 5% clay.In the present study, soil properties (e.g., porosity and field capacity) representing sand-based filter materials, were found to provide considerable runoff volume reductions even in the winter season across the urban areas.On the other hand, 3% of the catchment area relative to the facility area was used according to dimensioning values recommended by Larm and Blecken (2019); larger catchment areas would result in higher hydraulic loads and likely lower performance.Therefore, dimensioning criteria in cold regions should account for seasonal variations in infiltration rates, resulting in higher required surface areas compared to facilities designed for warm climates (Kratky et al., 2017).
Although retention capacities were found with a high inter-annual variability due to annual variations in climatic conditions, further key recommendations might be focused on four aspects: 1) determination of the required area of the facility to catchment area ratio (BCA/CA), by including lower infiltration rates during the cold season, and its implication on the facility area to achieve a certain runoff retention target.; 2) adequate selection of soil composition and determination whether local soils have sufficient infiltration capacity during cold conditions.If this is not the case, engineered bioretention media should replace local soils.In addition, the use of soils with high permeability (sand composition) and minor additions of topsoil to support vegetation growth.Since the critical factor influencing infiltration into frozen soil is the soil water content, soils with high hydraulic conductivity and low field (both of which will increase the probability of having relatively low soil moisture when frost occurs) are recommended (Zhao and Gray, 1999;Moghadas et al., 2016a,b).

Overall and annual performance of GRs
Overall, the highest and lowest ROR LT were found in the wettest and driest locations, respectively, with a performance that ranged between 24% in Gothenburg to 57% in Eskilstuna.Results showed that ROR LT was lower for locations with higher precipitation volumes and shorter drying periods, outcomes that agreed with previous studies where drier locations exhibited better retention capacities (Sims et al., 2016).Accordingly, Eskilstuna is characterized by a moderate temperate climate, with mild monthly precipitation and evapotranspiration volumes, and longer drying periods (low number of days with precipitation), conditions that favoured restauration of the GR retention capacity.Similarly, Johannessen et al. (2018) highlighted the sensitivity of GRs to antecedent dry weather periods (ADWP), and consequently, locations with higher potential evapotranspiration values provide better retention performances.On the other hand, Gothenburg is the urban area considered in this study with the highest precipitation volumes and number of days with rainfall, conditions that lead to higher initial soil moistures in the GR substrate.Likewise, previous studies have found a positive correlation between heavy rainfall events, and lower retention capacities (Nawaz et al., 2015), results that agreed with lower ROR LT in Gothenburg and Malmo.Similarly, Ornskoldsvik, Ostersund, and Gavle (urban areas located in the C region) were found with very similar ROR LT (ranging from 48 to 49%).The remaining cities all fit into a relatively small range (40-45%) of ROR LT .
In general, urban areas located in the C region (Ornskoldsvik, Ostersund and Gavle) exhibited higher ROR AN , being locations characterized by mild precipitation volumes and a number of days with precipitation (<570 mm and 150 days respectively), and mild MPET rates (87 mm).On the other hand, the minimum ROR AN were found in Kiruna and Gothenburg (11% and 12%), with a clear trend of lowest performances in urban areas located in the N and SW regions (Umeå, Luleå, and Malmo), with ROR AN that ranged between 14% and 18%.Similarly, locations towards the north of Sweden showed higher dispersion of ROR AN (13% for Kiruna, Luleå, and Umeå; 12% for Ostersund).Eskilstuna showed the highest number of significant differences in ROR AN with respect to the other urban areas (9 paired significant differences, Table 7).In addition, significant differences were observed between cities within regions of similar high-intensity rainfall (i.e., Kiruna -Luleå, Luleå-Umeå, in the N region; Orebro -Eskilstuna, Stockholm -Orebro and Eskilstuna -Stockholm in the SW region; Gothenburg -Malmo, in the SE region).However, no significant differences were found in the urban areas located in the C region.

Seasonal performance
Overall, most cities showed lower seasonal performance in autumn and winter as compared with spring and summer (Fig. 6), with the exception of Gothenburg and Malmo, which showed lower performance in spring.Previous studies of seasonal variations in runoff retention have shown that GR hydrological performance is reduced in winter and increased in summer, due to seasonal variations in potential evapotranspiration values and the sequence of rainy and dry periods (Mentens et al., 2006;Uhl and Schiedt, 2008;Nawaz et al., 2015).However, Stovin et al. (2012) found almost equal retention volumes during the

Table 6
Significant differences between the urban areas for RORs for BCs and GRs.W, SP, S AND AU represent annual, winter, spring, summer, and autumn seasons.spring and summer seasons, explained mainly by larger rainfall depths that occurred during the warm months (June to August).The combined effect of low evapotranspiration rates during the cold months (negligible values for urban areas in the N region), and mild precipitation volumes (mean monthly precipitation >54 mm for urban areas in the SW region), derived in minimum ROR s during the winter season in Gothenburg, Kiruna,Lulea,and Umeå (19%,27%,34%,and 27%).On the contrary, locations with high precipitation coupled with high potential evapotranspiration rates are likely to have the higher ROR s ; for example, Eskilstuna showed the maximum performance in summer (63%), a season where the potential evapotranspiration rates are at their highest.Overall, ROR s for the remaining cities vary in a relatively large range (23%-62% for maximum) of ROR s .
On the other hand, the greatest number of significant differences in ROR s were observed during the winter season in terms of ROR s (40, 32, 30, 36 for winter, spring, summer, autumn, respectively).This is in accordance with previous results, which showed that while the amount of rainfall greatly impacts GR retention, seasonal variations of climate factors (e.g., air temperature, and wet hours) that affect evapotranspiration rates also have considerable impact on GR retention (Sims et al., 2016).Consequently, Gothenburg was found as the urban area with most significant differences in ROR s with respect to other locations (37 paired significant differences, see Table 7), with clear seasonal differences (for the four seasons) with all the cities except Kiruna and Umeå (cities located in the N region).

Implications for design
On the other hand, the results highlighted the sensitivity of GR retention capacity and inter-event conditions (e.g., restoring storage capacity through evaporative fluxes, and the sequence between rainy days and dry periods).In general, Eskilstuna showed the highest hydrological performance on an annual and seasonal basis, mainly explained by mild precipitation and PET rates in all months, and a shorter a winter period where PET values are negligible.Therefore, findings suggest the importance of the adequate selection of plant species with higher water use that enhances the restoration of the storage capacity between storm events, as was pointed out previously by Johannessen et al. (2017).Besides, cold climates represent a challenge in the selection of plants species due to a difficulty in ensuring the establishment of dense, and persistent vegetation cover and the adaptability of species to local ambient conditions (Lönnqvist et al., 2021).Similarly, the length of the growing season in each urban area and the impact of cold months on plant dormancy are determining factors in the associated functionality of vegetation in GR (Caraco and Claytor, 1997).

Table 7
Significant differences between the urban areas for RORs for GRs.W, SP, S AND AU represent annual, winter, spring, summer, and autumn seasons.In addition, the study emphasized the important link between precipitation characteristics and required storage capacities.Findings showed minimum annual performances in Gothenburg and Malmo, locations characterized by higher mean annual precipitation and precipitation intensity.The extent to which large events can be retained is limited mainly by the declining storage capacity as the storm depth increases (Stovin et al., 2012).Therefore, locations with climates where the antecedent moisture conditions are favourable will have higher retention for all storm sizes (Sims et al., 2016).In addition to increased GR retention (Sims et al., 2016;Johannessen et al., 2018), the combined effect is positively correlated with vascular plant cover, and plant species diversity, as factors that can enhance the plant development and the GRs general multifunctionality (Lönnqvist et al., 2021).
In addition, key factors that enhance the sustained development of vegetation cover and promote greater survival include the adequate selection of the substrate depth (layer thickness) and its water-holding capacity (Dunnett and Nolan, 2004;Durhman et al., 2007;Dunnett et al., 2008).Similarly, the use of native species with the capacity to establish themselves quickly and tolerate extreme environment conditions (Sedum plants) should be considered in areas with short vegetation seasons and cool temperatures such as Scandinavian climates (Lönnqvist et al., 2021).Finally, GR maintenance may have considerable impacts on the vegetation development and dynamics associated to seasonal changes in cold locations.

Uncertainty analysis 3.4.1. Uncertainty analysis for BCs and GRs: long-term performance
Considering the uncertainties associated with parameter estimation, a 95% confidence interval (CI) of BC volume reduction of each urban area (see Table 7 for all CI) confirmed that Kiruna had the lowest longterm performance of the cities considered, as well as that winter/spring and summer, showed the lowest and highest seasonal performances, respectively (see Fig. 7).Considering all urban areas, the average relative uncertainties were ± 14%, with the highest relative uncertainties in Örnskoldsvik and Kiruna ( ± 20 and 24% respectively), and the lowest in Malmo and Gothenburg ( ± 5 and 7%).On the contrary, conclusions as to the urban area with the highest performance from the deterministic approach were not robust when considering parameter uncertainty, as Malmo rather than Ostersund had the best performance ROR LT when accounting for uncertainties.
Given the level of uncertainty, Kiruna and Gavle stood out as the two locations with clearly lower ROR LT than the other locations.In general, the highest relative uncertainties were found in the N and C region (i.e., urban areas with extended winter periods, and larger snow volumes), which is an indication that the parameters used in the snow melt model (e.g., melt coefficients, snow depth equivalent) are a major source of uncertainty in the modelling chain.Previous studies have found melt coefficients and snow water equivalent (SWE) to be the most influential factors for snowmelt and runoff generation (Semadeni-Davies, 2003).In general, all locations in the SE region showed similar ROR LT CI, with similar relative uncertainties that lie within the range of ± 12%.Overall, the results obtained showed a lower performance for the C and N region (mean ROR = 74% and 76% respectively), followed by the SE and SW regions (mean ROR = 83% and 84%).
Relative uncertainties in the long-term performance of GR were much higher for all the urban areas than those obtained for BCs, with the highest relative uncertainty found in Kiruna ( ± 45%), and those for the remaining cities fitting within the range of ± 27 to ± 36%.This finding indicates that GR performance is more sensitive to both design criteria (i.e., soil thickness, porosity) and snowmelt parameters (SWE and melt coefficients).In particular, Gothenburg and Kiruna have the lowest ROR LT , (51% for both locations, see Fig. 8), confirming that performances decreased significantly for locations with higher precipitation volumes, shorter drying periods, and minimum PET rates.On the contrary, Eskilstuna exhibited the highest performance across the locations (78%), but a high relative uncertainty of ± 38%.Considering that northern locations have both a higher fraction of precipitation in the form of snow, these results again indicate that the estimation of snowmelt model uncertainties leads to higher uncertainties in locations where snowmelt represents a higher fraction of the runoff generation (Kiruna, Luleå, and Umeå), with relative uncertainties between ± 36 to 48%.

Limitation of the approach
As with all modelling studies, the results of the present study consider several assumptions that affect the expected behavior of a theoretical model.These assumptions include the introduction of constant Ksat values for each seasonal period, in order to account for changes in the runoff infiltrated into the facilities (caused by frozen soils).Moreover, the analysis of changes in infiltration capacities involves a complex process driven by thermal and hydrophysical properties of soils, and moisture regimes (Maidment, 1996).Therefore, the inclusion of soil water content could have strengthened the relationship between lower infiltration capacities and soil freezing, as it is a critical factor affecting the thawing process, as was demonstrated by Moghadas et al. (2016).In addition, the present study does not account for the influence of bulk densities, organic matter content, and vegetation cover, as factors that influence the infiltration capacity in the soil media (Paus et al., 2016).Therefore, Ksat values used in the present study represent conservative infiltration rates, and therefore, a simplification of all the processes involved in its seasonal variation.
Similarly, the definition of seasonal periods was conducted independently for the eleven urban areas based on soil temperature thresholds (maintained trend of continued falling or rising soil temperatures during consecutive days).However, an annual variability of the boundary between the seasons was observed, which reduces the robustness of the definition of only one seasonal period for the 23-years of analysis.Results obtained could be affected due to a low resolution of point measurements (in this case, urban areas within the region) for each climatological region.More robust results could be strengthened by the use of additional data from more stations in each region.
On the other hand, previous studies pointed out a lower accuracy for SWMM (rainfall/runoff processes) for short, intense rainfall events following a long ADWP (Hamouz and Muthanna, 2019).The explanation assumes that the storage capacity is regenerated during dry periods and does not consider the influence of additional climatic effects (condensation of atmospheric vapor).In addition, in the evaluation of ROR performance of GRs during the winter periods, some portion of snowfall did not directly contribute to runoff volumes but remained on the roofs in the form of snow, a condition that might reduce the runoff generation during the winter period.The results highlight the need for further investigation on the adequate design configuration of BCs and GRs in a cold climate by providing scientific evidence for winter hydrology and a better understanding of snowmelt processes, as was pointed out previously by Moghadas et al. (2016a,b).Therefore, finding a modelling procedure introducing an optimal level of complexity remains a challenge, which could be achieved by understanding the model's assumptions and limitations (Bergstrom, 1991).

Conclusions
Significant variations in the GI retention performance on an annual and seasonal basis implies that it may be to ensure that the GI design configuration is tailored to satisfy local environmental requirements.In other words, assuming a linear relationship to scale the design by applying findings from one site to another has to be considered carefully.However, high inter-annual variability was also found in all the urban areas, showing that retention capacities vary greatly not only between locations but also in a given location according to annual changes in climatic conditions.Therefore, despite the existence of statistically significant difference, the operational significance of the observed differences may be limited.
Outcomes obtained highlight the relationship between the extent of the winter period (lower infiltration capacities) and negligible evaporation rates during the cold months, as factors that lead to lower annual performances for urban areas located towards the north of the country.Similarly, the high frequency of rainy days, and shorter periods for restoring storage capacity through evaporative fluxes, lead to lower performances in the southern locations.For instance, the analysis for GR found minimum ROR AN in Kiruna and Gothenburg (11% and 12%), being two urban areas with considerable differences in their geographic latitude, implying that similar performances are caused by different factors affecting the model water balance.
In general, winter/spring and summer demonstrated the highest/ lowest differences between the seasons.A statistical test used to determine if there is a significant difference between the urban areas under study showed two clear outliers.In the analysis corresponding to BCs, Kiruna was found to have significant differences with all the urban areas on an annual and seasonal basis, and Gothenburg, showed significant differences with all the urban areas in the four seasons, with exception of Kiruna and Umea.For the analysis of GRs, Eskilstuna and Gothenburg were found to be the urban areas with the most significant differences with respect to the other locations.Special considerations must be focused on the effect of rain-on-snow events, in locations where snowmelt represent a high portion of runoff generation (urban areas located in the N and C region).
While no robust results were obtained to suggest the adoption of regional design criteria for BCs, the performance of GR showed a clear difference in its seasonal distribution, as these facilities are more sensitive to changes in the spatial distribution of precipitation patterns (frontal for coastal zones and convective and orographic for inland zones).In general, results did not show a robust fit with the regionalization of high-intensity rainfall patterns for Sweden (Olsson et al., 2017b).
The results from this study highlight the importance of experimental data collection for green infrastructures over long periods in a variety of climatic regimes, and geographic locations.It also demonstrates that modelling can be a useful tool for apprehending the variability in green infrastructure performance that can be expected due to differences in climatic conditions between locations or in a given location from year to year, which may complement the results from field studies by identifying whether they were collected in favourable or in unfavourable conditions compared with other locations and/or time periods.This indicates the importance of the continuous monitoring of variables (i.e., maximum, and minimum melt coefficients, snow distribution within the catchment area, soil temperatures) in cold climate regions that require supplemented design considerations that address urban winter hydrology, where repeating cycles of frozen ground, snow cover, rain on snow and snowmelt are experienced.

Fig. 1 .
Fig.1.Overview of the observational data sets used in this study, and Swedish regional areas.

Fig. 5 .
Fig. 5. Comparison between urban areas in terms of ROR for BCs.(Red dots represent mean values for each urban area).

I
.Mantilla et al.

Fig. 6 .
Fig. 6.Comparison between regions in terms of runoff volume reduction for GRs [%].Red dots represent mean values ROR for each region.

Fig. 7 .
Fig. 7. Density plots showing distributions of long-term performance for GRs.Red line represents mean ROR LT .

Fig. 8 .
Fig. 8. Density plots showing distributions of long-term performance for GRs.Red line represents mean ROR LT .

Table 1
Hydrological and LID parameters used in SWMM.

Table 2
Snowpack parameters used in the degree-day method in SWMM.

Table 3
Snowpack parameters used in the degree-day method in SWMM.

Table 4
Definition of seasons based on soil temperature.

Table 5
Summary of climatic indices for the period 1997-2019.
Definition of annual, seasonal, and monthly climate indices (P = precipitation, T = temperature).aAnnual number of days with precipitation Average number of days where P > 0.2 mm.b Annual number of winter days Annual number of days given T < 0 • C (days).c Number of freeze and thaw cycles A cycle is defined as one or more days of T < 0 • C followed by a period with T ≥ 0 • C. d Annual number of summer days Annual number of days given T > 10 • C (days).e Total runoff volume reduction (ROR LT-BC ) Long-term runoff volume reduction (%) (25% and 75% quantiles) for BCs.f Total runoff volume reduction (ROR LT-GR ) Long-term runoff volume reduction (%) (25% and 75% quantiles) for GRs.I. Mantilla et al.