Development of the Surface Urban Energy and Water balance Scheme (SUEWS) for cold climate cities

. The Surface Urban Energy and Water Balance Scheme (SUEWS) is developed to include snow. The processes addressed include accumulation of snow on the different urban surface types: snow albedo and density aging, snow melting and re-freezing of meltwater. Individual model parameters are assessed and independently evaluated using long-term observations in the two cold climate cities of Helsinki and Montreal. Eddy covariance sensible and latent heat ﬂuxes and snow depth observations are available for two sites in Montreal and one in Helsinki. Surface runoff from two catchments (24 and 45 ha) in Helsinki and snow properties (albedo and density) from two sites in Montreal are also analysed. As multiple observation sites with different land-cover characteristics are available in both cities, model development is conducted independent of evaluation. The developed model simulates snowmelt related runoff well (within 19 % and 3 % for the two catchments in Helsinki when there is snow on the ground), with the springtime peak estimated correctly. However, the observed runoff peaks tend to be smoother than the simulated ones, likely due to the water holding capacity of the catchments and the missing time lag between the catchment and the observation point in the model. For all three sites the model simulates the timing of the snow accumulation and melt events well, but underestimates the total snow depth by 18–20 % in Helsinki and 29– 33 % in Montreal. The model is able to reproduce the diurnal pattern of net radiation and turbulent ﬂuxes of sensible and latent heat during cold snow, melting snow and snow-free periods. The largest model uncertainties are related to the timing of the melting period and the parameterization of the snowmelt. The results show that the enhanced model can simulate correctly the exchange of energy and water in cold climate cities at sites with varying surface cover.


Introduction
Today more than half of world's population resides in urban areas, and this fraction is expected to increase in the next decades (Martine and Marshall, 2007). Thus, the ability to understand and forecast the urban climate is crucial for sustainable urban planning and our quality of life. The exchanges of heat and water between the surface and the atmosphere are of great importance to urban climate studies. These exchanges describe the surface forcing in numerical weather prediction, air quality and climate models.
In urban areas several land surface models, with different complexity, simulate these energy exchanges, but none of the models consistently outperforms the others . The latent heat flux is commonly underestimated and sometimes even ignored, which further increases the direct heat emissions to the atmosphere. Furthermore, most of these models only concentrate on the surface-atmosphere interactions without any connection to the water cycles in urban areas. Similarly, several hydrological models for simulating urban drainage and the surface runoff in urban areas have been developed (Mitchell et al., , 2008Easton et Published by Copernicus Publications on behalf of the European Geosciences Union. Jacobson, 2011), but these do not typically consider the full energy balance.
Both in land surface and hydrological model studies, urban areas located in cold climates have been little studied despite their particular sensitivity to regional and global climate change. Thus appropriate, robust, well-tested modelling tools are needed. Modelling studies of cold cities are focused on a few sites mainly in North America (e.g. Valeo and Ho, 2004;Lemonsu et al., 2010;Leroyer et al., 2010) and Scandinavia (e.g. Semádeni-Davies et al., 1998). These emphasize the need for the correct description of snow cover in hydrological models. Snow affects surface energy partitioning via albedo and snowmelt, re-freezing and the phase-changerelated energy fluxes. The energy required for snowmelt can be of the same magnitude as the sensible and latent heat fluxes . Snow impacts water availability and its melt may cause springtime floods in urban areas (Semádeni-Davies and Bengtsson, 1998). To keep cities operational, snow is often redistributed within neighbourhoods and/or is transported away (Semádeni-Davies and Bengtsson, 1998Bengtsson, , 1999, which impacts both the energy and water cycles. The lack of observational data in urban areas with continuous winter snow cover makes the determination of model parameters and flux evaluation challenging. Surfaceatmosphere exchange of sensible and latent heat can be measured directly using the eddy covariance technique, but these observations are relatively rare, especially in cold climate cities. Notable exceptions include the work of Lemonsu et al. (2008), Vesala et al. (2008), Bergeron and Strachan (2012) and Nordbo et al. (2012a, b). These studies have found a strong seasonality in the energy exchanges and a need for the correct estimation of anthropogenic heat emissions from building sources, notably heating in winter. Similarly, the few hydrological studies have shown strong seasonality in stormwater runoff and differences in the amount of the snowmelt when compared to natural environments (Bengtsson and Westerström, 1992;Semádeni-Davies and Bengtsson, 1998;Valtanen et al., 2013).
The purpose of this study is to develop a model that can correctly simulate both the energy and water balances in cold climate cities. The model developed is included in the Surface Urban Energy and Water Balance Scheme (SUEWS, Järvi et al., 2011) with particular attention to the accumulation and melting of snow. The development and independent evaluation of the model uses several years of data collected in Helsinki (60 • N, 24 • E) and Montreal (45 • N, 73 • W). These include turbulent fluxes of heat and water measured with the eddy covariance technique, stormwater runoff and snow properties. In addition to snow related processes, the parameterization of the leaf area index has been improved to be more applicable for cold climate cities (Appendix A).

The Surface Urban Energy and Water balance Scheme (SUEWS)
The Surface Urban Energy and Water Balance Scheme SUEWS  simulates the urban energy and water balance components on a local or neighbourhood scale using hourly meteorological forcing data. These data inputs are kept to a minimum to enhance the flexibility of the model and commonly include: measured solar radiation (probably the least frequently measured), air temperature, relative humidity, surface air pressure, wind speed and precipitation. In addition, SUEWS requires information about the characteristics of the area to be simulated, such as surface cover fractions of paved surfaces, buildings, evergreen trees/shrubs, deciduous trees/shrubs, irrigated and non-irrigated grass, water, population density and building and tree heights. Rates of evaporation/interception for a single layer for each of the surface types are calculated and below each surface type, except water, there is a single soil layer. At each time step (5 min to 1 h), the moisture state of each surface and soil type is calculated. Horizontal water movements at the surface and in the soil are incorporated. Latent heat flux is calculated with a modified Penman-Monteith equation and sensible heat flux as a residual from the available energy minus the latent heat. The model contains several sub-models, for example, for net all-wave radiation (NARP, Offerle et al., 2003;Loridan et al., 2011), storage heat fluxes (Grimmond et al., 1991), anthropogenic heat fluxes and external irrigation.

New developments
The new version of SUEWS presented here incorporates a parameterization for snow cover. Previously, snow cover was a required input that was assumed to cover the whole grid area and only directly impacted the radiation. Now, accumulation and melting of snow are estimated, with impact to net all-wave radiation, evaporation and other water balance components included. For each surface type, the energy and water balances are calculated separately for snow-free and snow-covered areas and the model outputs are weighted according to their respective fractions. The energy and water flow calculations in the snow-free surface types follow those in the original version of the model . Here we present the equations related to the snow covered surface which is treated as a single snow layer.
The energy balance of the snow covered surface modified for urban areas can be written as (e.g. Oke, 1987;Cline, 1997) where Q M is the latent heat storage change caused by melting or freezing, Q s,I is the change in the storage heat of Geosci. Model Dev., 7, 1691-1711, 2014 www.geosci-model-dev.net/7/1691/2014/ the snow, Q * is the net all-wave radiation, Q F is the anthropogenic heat flux, Q H and Q E are the turbulent sensible and latent heat fluxes, Q P is the heat released by liquid precipitation on the snow, Q g is the heat exchange between the snow and the soil below and Q A is the net advective heat flux. Snowmelt occurs if the net energy input to the snow is positive (i.e. right-hand side of the Eq. (1) > 0). Q F is calculated based on cooling and heating degree days . Advection occurs at a number of scales. The micro-scale (or sub-grid-scale) advection is not resolved in the model, but rather embedded within the coefficients obtained using model optimization. The inter-grid advection is assumed to be negligible. This is consistent with the eddy covariance fluxes used to assess the model. To resolve advection at this scale would require the model to be embedded in a mesoscale model. The ground heat flux Q g is not separately resolved and is assumed to be included within the parameterization coefficients. The link to the snow mass balance is through Q E or evaporation (E): where P is precipitation (snowfall, rain), F is water that freezes on a snow-free surface, R is the runoff from the snowpack, T R is the transport of snow from the study area (e.g. via snow clearing) and S WE is the change in (liquid and solid phase) snow water equivalent (S WE ).

Surface albedo
Snow affects Q * by modifying the albedo of the surface and thus the reflected short-wave radiation, and the upwelling long-wave radiation as the surface temperature of snow and snow-free surface are different. The snow albedo (α s ) varies with snow age for each time step ( t), based on whether it is the "cold snow period" when melting does not occur (Baker et al., 1990): or the "warm snow period" when snowmelt occurs (Verseghy, 1991): For simplicity, the warm snow period is defined as the time when air temperature (T a ) is above 0 • C. α min s is the minimum snow albedo, τ d is a period of 1 day (86 400 s), and τ a and τ f are time constants related to the snow aging. After new snowfall, when S WE exceeds 2 mm (Koivusalo and Kokkonen, 2002), the snow albedo is reset to its maximum (fresh snow) value (α max s ). The upward long-wave radiation uses a constant snow emissivity.

Snow heat storage
The net heat storage in the snow can be considered for describing the convergence or divergence of sensible heat fluxes within the snowpack volume. This is calculated using the objective hysteresis model (OHM; Grimmond et al., 1991): where a 1 , a 2 and a 3 parameters are set by the model user. The first term describes the direct heating by radiation, the second term the hysteresis of the warming and cooling phases and the third the time lag. Q s,I is negative when the snowpack loses energy and the snowpack cools increasing the "cold content" of the snow (energy needed to heat the snow to 0 • C), and positive when the snow is heated towards 0 • C and the cold content is filled. Cold content is the total energy needed before the melting of snow can start (Bengtsson, 1982).

Energy for melting and freezing
There are two main approaches to estimate the snowmelt and refreezing of the meltwater (M) and the related energy (e.g. Martinec, 1989;Tobin et al., 2013): (1) the energy balance method, where M is calculated as a residual from the other energy balance components and (2) the degree-day method where M is calculated using daily or hourly air temperatures and possibly solar radiation. Although the first is more physically based it requires more input variables, whereas the latter uses more readily available variables. Comparisons of the two methods have found insignificant differences in the melted water calculated (Kustas et al., 1994;Debele et al., 2010). However, the site-specific degree-day parameters need to be assessed (Bentsson, 1984). In SUEWS, the second approach is used via a radiationtemperature index for each surface type i (Kustas et al., 1994;Semadeni-Davies et al., 2001;Tobin et al., 2013). Snowmelt induced runoff is delayed by the re-freezing of melted water (Bengtsson, 1982), particularly in spring, when the diurnal variations in air and snow surface temperatures are large. Daytime melt-water refreezes after sunset, releasing energy. Traditionally, the degree-day methods have utilized a daily time step, but in urban areas this shows poor performance (Bengtsson, 1984). Therefore, an hourly time step is utilized here. Melting and freezing occur as a function of air temperature (T a ) and Q * under three conditions: with factors for radiation melt a r (mm W −1 h −1 ), temperature melt a t and freezing a f (mm • C −1 h −1 ) which are typically linearly related with a f ≤ a t (Tobin et al., 2013). M i cannot be larger than the amount of solid snow in the pack and the amount of freezing water cannot exceed the amount of water in the snow. The energy consumed in melting and re-freezing is where ρ w is the water density at 0 • C (kg m −3 ) and L f is the latent heat of fusion at 0 • C (J kg −1 ).
Besides re-freezing of melted water, the snowmelt runoff from the snowpack is delayed by the amount of water the snow can hold (Bengtsson, 1982;Semádeni-Davies and Bengtsson, 1998). In SUEWS, this liquid water retention capacity (C R ) is calculated as a function of snow density (ρ s , kg m −3 ) (Anderson, 1976;Jin et al., 1999): where C R min and C R max are the minimum and maximum capacities and ρ e is a threshold density set to 200 kg m −3 . With time, the snow density changes (Verseghy, 1991): to a maximum snow density ρ max s with a time constant τ r . τ h is the seconds in an hour (3600 s h −1 ). After snowfall, ρ s is calculated as the weighted average of the fresh (ρ min s ) and previous snow densities.

Heat release by rain on snow
A rain-on-snow event provides heat, when the precipitation temperature is above the liquid/solid threshold (T lim ) : where c w is the specific heat capacity of water (J kg −1 K −1 ) and P i is the precipitation on ith surface (in m s −1 ). Here, it is assumed that the temperature of the precipitation is at the air temperature . Rain stays as a liquid and is routed to meltwater store.

Latent heat flux and evaporation
To calculate the latent heat flux (Q E ), a modified Penman-Monteith equation is used with a negligible surface resistance for the snow covered surfaces and an available energy that is constrained by snowmelt and re-freezing of the meltwater: where s is the slope of the saturation vapour pressure curve over ice (Pa • C −1 ) calculated according to Lowe (1977), γ is the psychrometric constant (Pa • C −1 ), c p is the heat capacity of air (J kg −1 K −1 ), ρ is the density of air (kg m −3 ), V is the vapour pressure deficit (Pa) and r a is the aerodynamic resistance (s m −1 ). To calculate the r a for the snow surface, the roughness length for heat and water vapour (z 0v , m) is calculated using (Voogt and Grimmond, 2000) z ov = z 0m exp(−20), where z 0m is the roughness length for momentum (m).

Change in snow water equivalent
For the water mass balance calculations, the model adopts a 5 min time step in order to respond to precipitation and snowmelt events. When the surface is completely covered by snow, the snow water equivalent of the ith surface type (S WE,i ) is calculated: If melt occurs (M i > 0) the water is held in the snowpack until the liquid water holding capacity C R i is exceeded. The excess water goes directly to runoff (R i ). If the surface is partially covered with snow, the excess water is added to the snow-free surface storage (S i ) and the snow-free surface equations are used . If a negative S WE,i occurs, the calculated evaporation is assumed to be too large and is reduced by an equivalent amount (constrained by E i ).
Snow from paved and built surfaces (T R,i ) can be transported out from the study area. The amount removed is calculated as amount of excess snow above a defined threshold (S WE,Lim ). This behaviour is neighbourhood specific (based on, for example, city or neighbourhood ordinances, snow clearance priorities). The S WE is assumed to be reduced to the S WE,Lim at the next site-specific snow clearing time period. People are also assumed to redistribute snow (e.g. paths are cleared and the snow is piled elsewhere) within the study area, and this is considered via depletion curves (Eq. 15a-c).
The snowpack starts to form when the surface temperature T s < 0 • C and under two conditions: when solid precipitation occurs and/or when water on a snow-free surface freezes (F i ). The snow depth s d (mm) is:

Surface fraction of snow
One of the most important factors controlling the energy balance and snowmelt is the patchiness of snow (Swenson and Lawrence, 2012). This is particularly important in urban areas, where snow clearing from streets and roofs takes place regularly . During the melt period, surface-type-specific depletion curves are used to approximate the fraction of snow cover (f s,i ) as a function of Geosci. Model Dev., 7, 1691-1711, 2014 www.geosci-model-dev.net/7/1691/2014/ S WE (e.g. Ek et al., 2003;Valeo and Ho, 2004). These are a function of surface-specific maximum snow water equivalent S max WE that control the initiation of snow patchiness (Swenson and Lawrence, 2012). For vegetated surfaces, the Swenson and Lawrence (2012) form of the function is used with coefficients estimated using the data for vegetated surfaces from Ek et al. (2003): As this function was developed for climate models, its application to smaller scales does require caution. For paved and built surfaces, the equations were derived from Valeo and Ho's (2004) data: The forms of the depletion curves are shown in Appendix B. The different curves between vegetation and impervious surfaces are used as human activities redistribute snow. For example, large roadside snow piles are created that melt slowly through the spring. In contrast, during the accumulation period snow is assumed to fall evenly on all surfaces.

Measurement sites and measurements
The model is applied in two cities that typically have extended periods of snow cover: Helsinki and Montreal. As multiple observation sites with different land-cover characteristics are available in both cities, model development is conducted independent of evaluation.

Helsinki, Finland
Meteorological and hydrological observations from three areas of Helsinki are used ( Fig. 1). At the Kumpula (Ku, SMEAR III) site, both meteorological forcing and evaluation data are measured (Järvi et al., 2009a). In addition, the observed runoff from Pasila (Pa) and Pihlajamäki (Pi) catchments are used for model development and evaluation. Kumpula is located 4 km north-east of the Helsinki city centre in a suburban area and 3.8 km from Pa and Pi (Fig. 1). Both Pa and the built sector of Ku (Ku1) have large areas of impervious surfaces (62 %). At both sites, the buildings are mostly offices with mean heights of 15 m (Pa) and 11 m (Ku1). Pasila has pedestrian areas at two heights with extensive concrete surfaces creating a complex morphology. The other two sectors around the SMEAR III flux tower (Ku2, Ku3) and the Pi catchment are more vegetated (Table 1). Pihlajamäki, with 34 % impervious surfaces, is a typical suburban area in Helsinki with multi-family block houses.
Tower based eddy covariance (EC) sensible and latent heat fluxes measured at 31 m, with an ultrasonic anemometer (Metek, USA-1) and a closed-path infrared gas analyser (LI-7000, Li-COR Biosciences, Lincoln, NE, USA) at Ku are used. Post-processing of the 10 Hz data use commonly accepted procedures described in detail in Järvi et al. (2009b) and Nordbo et al. (2012a).
Tower-top air temperature (platinum resistant thermometer, Pt-100, "in-house"), wind speed (Thies Clima 2.1x, Goettingen, Germany) and incoming and outgoing short-and long-wave radiation (CNR1, Kipp&Zonen, Delft, Netherlands) are used to force and test the model. Other forcing data measured on a nearby roof (24 m a.g.l.) include air pressure (Vaisala DPA500, Vaisala Oyj, Vantaa, Finland), relative humidity (Vaisala HMP243), and precipitation (rain gauge, Plu-vio2, Ott Messtechnik GmbH, Germany). Snow depth measured next to the tower by the Finnish Meteorological Institute is used in the model evaluation.
Runoff was monitored at 1 min intervals using an OCM Pro CF flow meter (Nivus GmbH, Eppingen, Germany) mounted in the two catchment storm flow discharge pipes from September 2010 to 30 April 2011 (See Table 2 for data availability). In Pi, excess pipe flow was observed causing  (2012) runoff at unexpected times. Because of the water quality observed, this is thought to be associated with pipe leakage in household water systems. From the beginning of September, the excess pipe flow observed was 0.0038 m 3 s −1 which had increased to 0.0125 m 3 s −1 by the end of the measurement campaign. This pipe flow was removed from the analysis when assessing the runoff as pipe leakage is not modelled currently.

Montreal, Canada
Two residential areas with impervious cover of 71 % (Rl, Rosemont-La-Petite-Patrie borough) and 49 % (Pr, Pierrefonds-Roxboro borough) 18 km apart were modelled (see Bergeron and Strachan (2012) for map). The more densely populated Rl has two to three storey buildings, whereas the suburban Pr is a single family house residential area (Table 1). At both sites, a tower mounted (26 m a.g.l.) sonic anemometer (CSAT3, Campbell Scientific Canada Corp., Edmonton, AB, Canada) and an open-path infrared gas analyser (LI-7500, Li-COR Biosciences, Lincoln, NE, USA) provided the 20 Hz data that are post-processed to EC fluxes of sensible and latent heat (Bergeron and Strachan, 2012). Forcing data of air temperature and relative humidity (HMP45C-212 at Rl, HMP45C at Pr, Campbell Scientific Canada Corp.), pressure (barometric pressure sensor, RM Young Model 61205V, RM Young Company, Michigan, USA) and radiation (CNR1) are from the EC tower at 25 m a.g.l. Snow depths were monitored in a back lawn of Pr and on a roof at Rl with snow ranging sensors (SR5, Campbell Scientific Canada Corp.). Snow properties, including snow density and albedo, were regularly (weekly: 2007/2008 winter or twice every month: 2008/2009 winter) observed for undisturbed snow, sidewalks, lawns and rooftops. Observations from Coteau-du-Lac (35 km south-west from Pr) and Pierre Elliot Trudeau Airport data (7 km from Pr and 16 km from Rl) (National Climate Data and Information Archive of Canada, 2013) are used to create a precipitation data set with separation of snow/rain.

Model runs
In Helsinki, SUEWS was run for Ku for 3 years (January 2010 till December 2012) and for the two catchments, for 16 months (January 2010-April 2011). At all three sites, the first 7 months are a spin-up period for the model that is neither used in model development nor testing. The spin-up time allows the model to become independent of the initial conditions set by the user. Even in urban areas, soil moisture initial state has a large impact on urban land surface model performance (Best and Grimmond, 2013). The remainder of the periods are used for model development and evaluation. In Ku, data prior to 2012 are used to develop and adjust model parameters: Q * , upward shortwave radiation and upward and downward long-wave radiation are used to adjust Geosci. Model Dev., 7, 1691-1711, 2014 www.geosci-model-dev.net/7/1691/2014/  (Pa) is used to constrain the temperature and radiation melt rates (Eq. 6), retention capacity of the snow (Eq. 8) and the limit for the liquid precipitation. Q * and its components, Q H and Q E , the snow depth from Ku in 2012 and the runoff from the medium-intensity catchment are used to independently evaluate the model. The meteorological data measured at the Ku site are used to force the model for all three Helsinki sites. The data are gap filled using the procedures described in . Due to the very different characteristics surrounding the Kumpula tower, the model is run for the three surface cover areas within a 1 km radius circle. The flux time series evaluated against observations are combined from the surface cover areas (Ku1-3) based on the prevailing wind direction.
In Montreal, only the first of the 22 months (December 2007 till September 2009) is used as a spin-up. The short spin-up time is chosen to allow two snowmelt periods in model development and testing. The remainder of the suburban data set (Q * , Q H , Q E , snow depth, snow density and albedo) is used for the model development: snow density and albedo are used to determine the shape of the snow aging curves (Eqs. 3, 4 and 9); Q * , the surface and snow albedo; and Q H and Q E , the other snow related parameterizations. The urban site observations are used for independent evaluation of the model. The model domain is a 1 km radius circle around the flux tower.
The results are analysed by considering snow-free, cold snow and melting snow periods. For snow-free periods, the simulated snow depth is zero, whereas the cold snow and melting snow periods are separated by the air temperature 0 • C.

Evaluation statistics
Several statistics are used to evaluate the model performance (e.g. Daley, 1991). Linear regression is used to describe the linear dependence between the independent variable, in this case the observations (X Obs ), and the dependent variable, the model output (X Mod ). The slope (S) relative to 1, and intercept (I ) relative to zero, provide information on the model performance. Further, goodness of fit is evaluated using the root mean square error (RMSE): where N is the number of data points. Like the intercept in the linear regression, the RMSE has the units of the variables being evaluated and it depends on the magnitude of the mean variables. Therefore, it is useful to normalize the RMSE (nRMSE) relative to the range of values observed : When comparing the performance of the model to simulate different variables, the RMSE can also be normalized with the standard deviation of the observations σ Obs (Taylor, 2001): In addition, the mean bias error (MBE) between the modelled and observed time series is considered: where the overbar indicates an average. Ideally, the RMSE, nRMSE, sRMSE and MBE would all be zero.

General weather conditions
The weather conditions during the modelled period for Helsinki and Montreal are shown in Figs. 2 and 3, respectively. Daytime solar radiation exhibits a strong seasonal pattern, with the 15 • latitudinal difference causing more rapid changes and stronger amplitude in Helsinki than in Montreal. In summer, K ↓ reaches 970 W m −2 in Montreal, whereas in Helsinki, the maxima remain below 830 W m −2 . In winter, the solar radiation in Helsinki is very low (< 120 W m −2 ), whereas Montreal peaks are below 400 W m −2 . Despite the difference in K ↓, air temperatures are fairly similar in both cities. Daily maxima mean temperatures are around 26 • C in summer, while the minimum daily mean temperature in winter in Helsinki is −20 • C and in Montreal −23 • C. In both cities precipitation is quite evenly distributed throughout the year.
During the 3 years of measurements, the daily snow depths in Helsinki are all below 0.8 m, with a longer snow period in winter 2010/2011 than 2011/2012. The timing of snowpack formation depends strongly on the year. In 2010, it was initiated in November, whereas in the following winter this was delayed until January 2012. This will have a large influence not only on both natural energy and water exchanges, but also urban activities. In Montreal, snowpack depth and timing has large variability between years; for example, a 1 m snow pack is observed in March 2008 with melting in late April, compared to only 0.6 m the next year, with snow melting by the end of March.

Snow properties
The time constants for describing the aging of snow, the minimum and maximum snow albedo, and density were determined by optimization using observations undertaken at the suburban site in Montreal (Pr). The observed snow properties are treated as averages from the measured surface types in order for them to be compatible with the scale of the simulations. To evaluate the snow albedo, the observed reflected shortwave radiation (K ↑) in Helsinki in 2011 is also used. To assess the radiative exchanges, SUEWS is run using the radiation measurements source area (99 % field of view) estimated as a 31 m radius circle around the 31 m tall measurement tower (Nordbo et al., 2012a). The surface cover characteristics are different for this area than those within the turbulent flux source area; with 49 % paved surface, 4 % buildings, 3 % deciduous trees/shrubs and 44 % grass.
Geosci. Model Dev., 7, 1691-1711, 2014 www.geosci-model-dev.net/7/1691/2014/ Comparison of these observations with the Lemonsu et al. (2010) (hereafter Le10) aging functions used for the suburban site in Montreal shows that modifications to the coefficients are needed for both snow albedo and density (Fig. 4). The Le10 maximum density of 350 kg m −3 is too small for the current observations. Thus, the maximum snow density is set to 400 kg m −3 ; the minimum value ρ min s is kept at 100 kg m −3 . In addition, the time constant τ r is decreased to 0.043. After these modifications, the simulated snow density follows the behaviour of the median observations well (Fig. 4a). Similarly from the observations, the minimum (α min s ) and maximum (α max s ) snow albedo are set to 0.18 and 0.85, respectively, which differ from Le10 (α min n = 0.15-0.30 across the different surface cover types). In Lemonsu et al. (2010) snow albedo aging time constants (τ f = 0.174, τ a = 0.008) could not be fully evaluated due to a lack of data. However, τ a compared to our observations is too small. When increased to 0.018, τ f decreases to 0.11. Again, good correspondence between the observed snow albedo and model output (Fig. 4b), and between the observed and modelled K ↑ in Helsinki in 2012 (not shown) are seen. During the cold snow period, the linear fit statistics are S = 0.68 ± 0.02 W −1 m 2 , I = 0.27 ± 0.47 W m −2 (RMSE = 11.3 W m −2 , N = 2232), and during the warm snow period S = 0.50 ± 0.01 W −1 m 2 and I = 0.85 ± 0.47 W m −2 (RMSE = 11.4 W m −2 , N = 604). One likely reason for the poorer model performance during the warm snow period is the sensitivity of the albedo to the fraction of snow covered surface. In the model, the fraction of snow is parameterized based on the maximum S WE , but it is likely that this is site dependent at a neighbourhood scale due to redistribution and transport of snow. However, as the other net all-wave radiation components are larger in magnitude than K ↑, the negative bias during the melting period is likely to have small impact on the available energy.

Melt and freezing factors
The freeze and melt factors (a r and a t ), representative for the neighbourhood scale, are optimized using runoff from Pa and snow depth from Ku (Helsinki). SUEWS was run using a r values between 0.0008 and 0.002 mm W −1 h −1 using 0.0001 mm W −1 h −1 resolution, and a t between 0.05 and 0.15 mm • C −1 h −1 with 0.01 mm • C −1 h −1 resolution. The 146 modelled combinations were analysed with respect to the amount of meltwater accumulated during the snow covered period and the timing for complete snowmelt ( Table 3). The smallest differences compared to the observations are obtained with a r = 0.0016 mm W −1 h −1 and a t = 0.12 mm • C −1 h −1 . Thus, these are used in the model runs. Values are slightly larger, but of the same order of magnitude, than those obtained for hourly factors at an Arctic watershed in Alaska (a r = 0.001 mm W −1 h −1 and a t = 0.095 mm • C −1 h −1 ; Kane et al., 1997). Unfortunately, no hourly values for urban areas were found in the literature. However, using these factors the daily melt rates are the same order of magnitude as those that have been typically reported for urban areas (Bengtsson and Semádeni-Davies, 2011).

Snow storage heat
To determine the storage heat flux coefficients a 1 , a 2 and a 3 for snow (Eq. 5), shallow water values were used as an initial basis with a 1 = 0.50, a 2 = 0.21 and a 3 = −39.1 (Souch et al., 1998). Given the assumption that the snow heat capacity is around half that of water (Rogers and Yau, 1996), a 1 is set to 0.25 for snow. The other two coefficients (a 2 and a 3 ) are assessed relative to their effect on the sensible heat flux by running SUEWS for Pr over a range of values during the snow covered period. Pr was chosen to optimize Q s,I due to its more homogeneous surface characteristics when compared to other sites. The RMSE between the observed and modelled Q H varies between 44.4 and 49.0 W m −2 and MBE between −23.1 and 28.0 W m −2 , when a 2 varies between 0 and 1.2 and a 3 between −60 and 0 ( Table 3). The optimal result with minimum RMSE = 48.2 W m −2 and MBE = 0.19 W m −2 is obtained with a 2 = 0.60 and a 3 = −30. Thus, these coefficients are used in the model to calculate the snow storage heat or the internal energy of the snow. The values imply a smaller slope or fraction of radiative energy entering/leaving (a 1 ), a greater hysteresis (a 2 ) and a similar phase or time lag (a 3 ) for snow relative to water. Heuristically this appears appropriate.   Table 4 lists the parameters, both for the snow covered and snow-free surface, used in the model runs. The snow parameters are optimized (Sect. 3.1), whereas the limit values for the snow transport (S WE,Lim ) were estimated based on maximum mass allowances of water at the Pa catchment. The same values were used at all sites as no additional information was available. Sensitivity analyses (Sect. 3.4) suggest the model is fairly insensitive to S WE,Lim despite the expectation for values to be neighbourhood specific. Figure 5 shows the daily observed and modelled runoff from the two catchments in Helsinki. The grey line separates the non-snow and snow related runoff as the continuous winter snow cover formed on 18 November 2010. At both sites, the model simulates the snowmelt induced runoff well, reproducing both the spring melt peak and the recession in April. When the model is run treating the catchments as a whole, it tends to overestimate the runoff peaks and to be more peaked than observations (Fig. 5), which have smaller but longer runoff peaks. Partly this can be explained by the absence of time lags for the water to move from the most distant points (hydrologically and hydraulically) because the catchment is modelled as one unit (in the current setup). However, in terms of hourly performance, the correlation between the observed (R obs ) and modelled (R mod ) runoff is good with S = 1.20 (mm h −1 ) −1 and I = −0.01 mm h −1 (RMSE = 0.14 mm h −1 , r = 0.75) in Pa, and S = 1.24 (mm h −1 ) −1 and I = 0.02 mm h −1 (RMSE = 0.16 mm h −1 , r = 0.60) in Pi. The coefficients are calculated for periods when both R mod and R obs are nonzero (675 and 760 h in Pa and Pi, respectively). In Pa, the model underestimates the cumulative runoff over the snow covered periods by 3 % as R mod = 82 mm and R obs = 85 mm (Fig. 6). In Pi, the cumulative runoff is overestimated by 19 % as R mod = 97 mm and R obs = 81 mm. Before the continuous snow cover, the model performs similarly at both catchments. Notably, the model overestimates runoff at Pi with high intensity precipitation. The overestimation is seen in the linear correlation between R obs and R mod as S = 1.29 (mm h −1 ) −1 and I = 0.04 mm h −1 (RMSE = 0.20 mm h −1 , r = 0.68, N = 668) as well as in the modelled cumulative runoff, which is 47 % higher than the observed in Pi (27 and 42 mm, respectively) (not shown). In Pa, S = 1.26 (mm h −1 ) −1 (RMSE = 0.20 mm h −1 , r = 0.90, N = 743), and the cumulative runoff is underestimated by 4 % by the model (R obs = 84 mm and R mod = 88 mm). Some of these differences are caused by the forcing precipitation and other meteorological variables being from the flux site Ku. This would particularly affect the model performance during convective precipitation, which accounts for 88 % of the precipitation events between April and September (Punkka and Bister, 2005).

Snow depth
The model calculates snow water equivalent (S WE ) and snow depth (s d ) separately for each surface type. Due to snow removal and the different surface characteristics, s d behaves differently for the vegetated and built surfaces. This can be seen when the modelled s d for each surface (paved, building, grass and tree) is plotted with the observations for Helsinki (Fig. 7) and Montreal (Fig. 8). Unfortunately, the observations are each representative of individual point and surface types, whereas the model values are for the different surface Geosci. Model Dev., 7, 1691-1711, 2014 www.geosci-model-dev.net/7/1691/2014/  types at the neighbourhood scale. Thus, some differences between the modelled and observed s d are expected. In Helsinki, the point observations are made in an open space that corresponds most appropriately to the modelled grass surface. Data for 2011 and 2012 are plotted separately as the first year is used in the model parameter determination, whereas the latter is an independent data set (Fig. 7). In both years, the model reproduces the accumulation of snow and melt events well, but underestimates the snow depth by 84 mm on average compared to the observations. The measured maximum snow depth in 2011 is 720 mm, whereas the modelled snow depth above grass is 587 mm. Similarly, for 2012 the observed depth is 630 mm and the modelled value is 504 mm. This underestimation is caused by either the underestimation of modelled S WE or by overestimation of snow density as the snow depth is a function of these two variables (Eq. 14). The model starts to accumulate snow 4 days later than the observations in January 2012, but later in the year the observed and modelled snow cover appear on the same date (29 November). In 2011, the snowmelt is observed to have ended on 15 April, 5 days earlier than simulated, whereas  in 2012 the snowmelt finished on 12 April, 1 day later than modelled.
For Montreal, s d is calculated separately for the suburban (Pr) and urban (Rl) sites for January 2008-April 2009. In Pr, the observations are made on a lawn corresponding to the modelled grass surface (Fig. 8a). The model follows the accumulation and melt events well, but like Helsinki, snow depths are underestimated, especially in the 2007-2008 winter. The maximum observed s d is 1020 mm, while 670 mm was modelled for grass. In winter 2008-2009, the maximum observed and simulated snow depths are 665 and 460 mm. In 2008, the modelled snow starts to accumulate on the same day (8 December) as observed. Snowmelt is completed on 20 April 2008, which is 3 days after the modelled date. In 2009, the modelled snowmelt finishes 1 day before observed (30 March).
In Rl, s d is observed on a building roof. This has both lower snow amounts and earlier melt compared to the lawn observations in Montreal (Fig. 8b). The model simulates this behaviour well, but again underpredicts the depths. The observed s d maxima are 390 and 415 mm for the two winters, while 301 and 285 mm are modelled, respectively. Accumulation of snow takes place on the correct day in Rl, and the snowmelts on the same day as observed (26 March) in 2008, and 9 days later (7 March) in 2009 than observed.
The underestimation of s d is also impacted by the precipitation measurements, as the difference between observed and modelled values begins during the snow accumulation period. Precipitation measurements are known to underestimate snowfall, especially due to wind effects (Goodison et al., 1998;Savina et al., 2012).

Turbulent and radiative energy fluxes
The simulated Q * , Q H and Q E are assessed for snow-free, cold snow and warm snow periods (Table 5, Fig. 9), with the diurnal behaviour of both the observed and modelled fluxes for the independent data sets in Helsinki and Montreal considered (Figs. 10 and 11).
Generally, the best simulated flux of the three is Q * independent of whether there is snow on the ground or not. For the cold and warm snow periods, the RMSE varies between 27-41 W m −2 (nRMSE = 0.037-0.061) and 31-41 W m −2 (nRMSE = 0.041-0.057) across the sites. At all sites, Q * is underestimated in cold snow conditions with MBE between −34 and −13.5 W m −2 . Mostly this underestimation is related to the downward long-wave radiation that is calculated Geosci. Model Dev., 7, 1691-1711, 2014 www.geosci-model-dev.net/7/1691/2014/  from air temperature and relative humidity (Loridan et al., 2011 -who suggest that use of cloud cover data as input with this technique is better). This parameterization works less well in cold conditions than above 0 • C. In Montreal, the warm snow underestimation is even larger (MBE = −37 to −25 W m −2 ), compared to Helsinki where the underestimation decreases to −5 to −3 W m −2 . Especially during the warm snow periods, the fraction of snow cover plays an important part in the model performance. It affects both the snow albedo and outgoing long-wave radiation via surface temperature. The best performance for Q * is under snow-free conditions, when the MBE is between −10 and 8 W m −2 and the nRMSE is clearly lower than for the periods with snow cover (Fig. 9a).
The scatter in the model performance is larger for Q E than for the other energy balance components, with cold snow periods having the best, and warm snow periods the poorest model performance (Fig. 9c). The RMSE during cold snow varies between 9-12 W m −2 (nRMSE = 0.038-0.059), and for warm snow between 22-34 W m −2 (nRMSE = 0.071-0.170). The increase in RMSE during warm snow periods is understandable as the energy consumed in melting snow and freezing meltwater is higher and thus errors in the degree-day-method propagate more easily to Q E (as well as to Q H ). During melting periods there can be advection from snow-free surfaces to the snowpack altering the energy balance as specified in Eq. (1) (Bengtsson and Semádeni-Davies, 2011). MBE varies between −11 and 4 W m −2 when there is snow on ground. During snow-free periods, the model underestimates Q E at all sites with MBEs between −27 and −11 W m −2 , RMSEs between 26 and 35 W m −2 and nRMSE between 0.058 and 0.066.
In SUEWS, Q H is calculated as a residual from other energy fluxes; therefore, the net error accumulates in Q H . Despite this, the model is able to simulate its behaviour well. When there is snow on ground, the RMSE varies between 23 and 50 W m −2 and nRMSE between 0.067 and 0.118. During the cold snow periods, the simulated heat fluxes are slightly better than during warm snow periods, similar to Q E . The model overestimates Q H during snow cover, except in Rl during cold snow periods, and MBEs vary between −25 and 9 W m −2 . In summer, the performance of the model in simulating Q H improves following the performance of Q * and Q E . The model performance for the energy fluxes is more dependent on the period of analysis than the site where it is run. An exception to this is Q H at R1, where the model overestimates and shifts the diurnal peak flux earlier compared to the observations (Fig. 11). This appears whether there is snow on ground or not, suggesting that this is caused by the snowfree storage heat flux which is underestimated by the model or anthropogenic heat flux, which is overestimated. RMSEs obtained for the warm snow periods in Pr are higher than Le10 obtained for the same suburban area in Montreal using the Town Energy Balance model in spring 2005. However, direct comparisons are difficult as in their 1 month of observations snow cover is present only on some days compared to the longer time period evaluated here.

Energy balance of urban snow covered surface
Snow cover and the related energy storage and the energy related to phase change alter the surface energy balance. The components at the most built-up site Rl are evaluated here (Fig. 12). During cold snow periods, the daytime energy balance is dominated by the net all-wave radiation (Q * ) and the sensible heat flux (Q H ), reaching 119 W m −2 and 113 W m −2 . Q H is fuelled by both Q * and Q F (reaching 46 W m −2 ), and it accounts on average for 68 % of the daytime (10:00-15:00 LT) available energy. The dominance of Q * and Q H are typical also for natural cold snow packs (e.g. Oke, 1987). Only 12 % of Q * + Q F is dissipated by evaporation, whereas the storage fractions are 10 and 8 % at the snow and snow-free surfaces, respectively. At night, on the other hand, the urban surface loses long-wave radiation causing the internal energy of the snow to decrease, that is, the cold content of the snow increases. At the same time, the snow-free surface loses some energy (around 10 W m −2 ) and both Q F and Q E remain positive (by more than 10 W m −2 ), with Q H less than 5 W m −2 .
During warm snow periods, Q * is clearly the most important component of the surface energy balance reaching 200 W m −2 in daytime. Now the daytime Q F is slightly smaller than during the cold periods reaching 35 W m −2 . Most of the energy, but clearly less than during the cold snow period, is partitioned to Q H (46 %), with the second largest contribution going to snow-free surface heat storage (29 %). Evaporation consumes 17 % of the energy, and only 4 % and 3 % is stored in the snow and consumed by snowmelt. The largest Q H and Q S are consistent with the observations obtained by Le10 at the suburban site, although they documented a larger contribution going to snow related processes than to evaporation. Moreover, the modelled fractions during the snow covered periods are of the same order of magnitude as obtained for observations at the same site (Bergeron and Strachan, 2012).
When the ground is free of snow, most energy (Q * + Q F ) again goes to Q H (188 W m −2 , 45 %) followed by the storage heat flux (175 W m −2 , 42 %). Due to the high impervious nature of the surface, daytime Q E reaches 50 W m −2 , which is only 12 % of the available energy. The resulting daytime Bowen ratio (Q H /Q E ) is 3.7, which corresponds well with the expected relationship of the Bowen ratio with the site's vegetation fraction (Loridan and Grimmond, 2012).

Model sensitivity
To better understand the impact of both the optimized values and those estimated (Table 4) without detailed observation on the model performance, sensitivity analyses were undertaken. The analysis included the power of the vegetation depletion curve (Eq. 15a), limit for the transport of snow from paved and building surfaces (S WE,Lim ), snow heat storage (a 2 , a 3 ) and the meltwater coefficients (a r , a t ). SUEWS was run using the three independent datasets (Ku in 2012, Rl and Pi) changing the parameters by ± 30 % using a 10 % step. The results were compared to hourly measured Q H , Q E and runoff and the RMSE determined for each site and variable (Fig. 13). The other parameters were held constant during each set of analyses.
The impact of the coefficients on heat storage in snow pack is shown only for Q H as their effect on Q E and runoff is small or non-existent. S WE,Lim and the meltwater coefficients have the largest impact on the heat fluxes at both sites, whereas the smallest effect is for the power used in the depletion curve.    indicating that the model is fairly insensitive to changes in the studied parameters.

Conclusions
The Surface Urban Energy and Water Balance Scheme (SUEWS) is developed to simulate the energy and water The sensitivity are to changes in: the power in the depletion curve (Eq. 15a) (C1), S WE,Lim (C2), meltwater coefficients (a r , a t ) (C3) and storage heat flux coefficients (a 2, a 3 ) (C4). The final model values are indicated (*). For other details see text.
balances in cold climate cities with special attention on the simulation of snow cover. The new model considers the accumulation of snow, snow properties including snow water equivalent, snow depth, snow density and albedo and snowmelt and refreezing of meltwater based on an hourly degree-day method. The development and independent evaluation is undertaken using observations from three sites in Helsinki and two sites in Montreal. Each of these sites varies in terms of surface cover characteristics. In Helsinki, the observations include stormwater runoff from two catchments and turbulent fluxes of sensible and latent heat from one site. In Montreal, the observations include snow properties as well as the turbulent fluxes at both sites.
The model developments include an improved description of vegetation phenology (LAI) in cold climate cities. The leaf-off period based on daily air temperature was accelerated using a combination of daily air temperatures and day length. Updated aging functions for snow density and albedo in urban areas were developed based on snow observations in Montreal; an improved equation for the degree-day method was used to calculate snowmelt and freezing of the melted water; and new parameter values were developed to calculate the snow storage heat flux using the objective hysteresis model (OHM).
The enhanced model can correctly simulate the winter and springtime melt-related runoff, but the runoff peaks tend to be sharper than the observations partly due to the absence of time lag to let the water flow to the observation point at the catchment discharge point. Despite this, the modelled cumulative runoff during the snow covered periods corresponds well with the observations. The formation and melting of the snowpack is simulated well both in Helsinki and Montreal, but the snow depth is underestimated either due to overestimation of the snow density or underestimation of snow water equivalent. Following the hydrological variables, the net radiation and turbulent sensible and latent heat fluxes also are modelled well. The model simulates their diurnal behaviour throughout the year, but the largest uncertainties occur during the snowmelt period at all sites. This is related to the uncertainties in determining the snow covered surface fractions, as well as the propagating uncertainties from the calculation of melt and freezing related energies based on the degree-day method.
The model can correctly simulate the energy and water cycles in cold climate cities and can potentially be used independently for urban planning purposes or nested to a mesoscale or global scale atmospheric model. However, some of the parameterizations are still city and site dependent; more observations from cold climate cities are needed to create more generalized formulations.