High-resolution hydrometeorological modelling of the June 2013 flood in southern Alberta, Canada

From June 19 to June 22, 2013, intense rainfall and concurrent snowmelt led to devastating floods in the Canadian Rockies, foothills and downstream areas of southern Alberta and southeastern British Columbia. The complexity of the topography in the mountain headwaters, presence of snow at high elevations and other factors challenged hydrological forecasting of this extreme event. In this study, the ability of the Global Environmental Multi-scale hydrological modelling platform (GEM-Hydro), running at a 1-km grid spacing, to simulate hydrometeorological conditions in several Alberta rivers 15 during this event is assessed. Four quantitative precipitation estimation (QPE) products were generated using the Canadian Precipitation Analysis (CaPA) system by varying (i) station density and (ii) horizontal resolutions (10, 2.5 and 1 km) of the GEM precipitation background. CaPA at 2.5 and 1 km including all available stations in the headwaters provided the most accurate estimation of intensity and total amount of precipitation during the flooding event. Using these products to drive GEM-Hydro, it is shown that QPE accuracy dominates the ability to predict flood volumes. Initial snow conditions also 20 represent a large additional source of uncertainty. Default GEM-Hydro simulations starting with almost no snowpack at high-elevations led to a systematic underestimation of flood volume and peak flow. Gridded estimates of snow water equivalent from the Snow Data Assimilation System (SNODAS) were also considered. They led to contrasting abilities to simulate flood discharge volumes and a consistent overestimation in the headwater catchments, illustrating the strong need for a reference snow product in the mountains of Western Canada. Finally, GEM-Hydro did not predict peak flow timing and 25 hydrograph shape well. Model sensitivity tests show that it could be improved by adjusting the Manning coefficients, suggesting the need to revisit the routing parameters. There may be a need to include water management effects on flood hydrographs as well. These results will guide the development of GEM-Hydro as a hydrological forecasting system in Western Canada. Hydrol. Earth Syst. Sci. Discuss., https://doi.org/10.5194/hess-2019-152 Manuscript under review for journal Hydrol. Earth Syst. Sci. Discussion started: 23 April 2019 c © Author(s) 2019. CC BY 4.0 License.


Introduction
From June 19 to 22, 2013, heavy rainfall and snowfall (at high elevations) with local amounts exceeding 300 mm fell over a broad region of the Canadian Rocky Mountains of southeastern British Columbia and southwestern Alberta, and the foothills and adjacent plains of southern Alberta, Canada (Milrad et al., 2015;Liu et al., 2016;Kochtubajda et al., 2016;. This heavy precipitation resulted from moisture convergence from the Pacific and from the Canadian Prairies and 5 US Great Plains through evapotranspiration (Milrad et al., 2015;Li et al., 2017). At high elevations, rain and then snow fell on a deeper than normal, late-lying snowpack leading to rain-on-snowmelt which enhanced runoff generation (Fang and Pomeroy 2016;Pomeroy et al., 2016a). This heavy rainfall combined with a rapidly melting alpine snowpack triggered severe flooding in the Oldman, Bow and Red Deer river basins of southern Alberta and the Elk River of southern British Columbia and impacted many communities, including Calgary, the largest city in Alberta . The floods 10 led to the evacuation of 100,000 people and caused five fatalities and over $6 billion of damage to infrastructure such as roads, railways, bridges, parks and homes .
This extreme weather event has been investigated in details in the last few years from the perspectives of atmospheric and climate modelling. Li et al. (2017) showed that the Weather Research and Forecast (WRF) atmospheric model at convectionpermitting resolution (3-km) demonstrated reasonable skill in simulating the spatial and temporal evolution of precipitation 15 patterns before and during the flood. They investigated with the model the dynamical features that led to the generation of heavy rainfall. Milrad et al. (2017) carried out sensitivity experiments with the WRF model at the same resolution as Li et al. (2017) to determine the impacts of local topography on the extreme rainfall event. Teufel et al. (2017) presented complementary results using the Canadian Regional Climate Model at lower resolution (0.11 ∘ ) than WRF. In particular, they showed that wet antecedent soil moisture conditions in the Canadian Prairies and US Great Plains impacted precipitation 20 patterns, intensities and amounts through enhanced evapotranspiration in these upwind regions.
Despite its severe hydrological consequences, this extreme weather event has received little attention from a hydrological modelling point of view. Teufel et al. (2017) discussed the impact of initial snow and soil conditions on the amplitude and timing of streamflows. However, their study only used daily discharges from one hydrometric station located in Calgary along the Bow River for model evaluation, and this station was heavily impacted by reservoir regulation making it unsuited 25 for this purpose. Li et al. (2017) and Milrad et al. (2017) did not use their high-resolution WRF simulations to generate hydrological simulations of the flood. At the local scale, Fang and Pomeroy (2016) studied in detail the impact of antecedent soil moisture and snowpack conditions on runoff generation during the flood for a small (9.4 km ! ) and well-instrumented mountain research basin. Pomeroy et al. (2016a) used the data from the same basin to better understand and simulate the complex rain-on-snowmelt runoff generation during the event. They showed that successful modelling of the event at small 30 scales required detailed consideration of precipitation phase, rain-on-snow processes, slope/aspect impact on snowpack energetics, forest canopy-hydrological process interactions, antecedent soil moisture, sub-surface storage and routing of runoff water through both sub-surface and overland flow pathways. However, the modelling studies mentioned here did not Hydrol. Earth Syst. Sci. Discuss., https://doi.org /10.5194/hess-2019-152 Manuscript under review for journal Hydrol. Earth Syst. Sci. Discussion started: 23 April 2019 c Author(s) 2019. CC BY 4.0 License. evaluate the potential and limitations of large-scale hydrological models to simulate the dynamics of the flood at hourly time scales over large river basins. Furthermore, within the context of the development of a national hydrological forecasting system for Canada (GWF project, 2018), there is a need to better assess the ability of the hydrological systems currently used in Canada to simulate such extreme events.
In the recent years, Environment and Climate Change Canada (ECCC) has developed the Global Environmental Multi-scale 5 hydrological modelling platform (GEM-Hydro, Gaborit et al. 2017). GEM-Hydro combines the GEM-Surf surface prediction system (Bernier et al., 2011) including the SVS (Soil, Vegetation and Snow) land surface scheme Husain et al., 2016) and the WATROUTE routing scheme (Kouwen, 2010). GEM-Hydro is designed "to perform reliable streamflow simulations with a distributed model over large and partly ungauged basins" (Gaborit et al., 2017). The model follows the same approach as the Weather Research and Forecasting Model hydrological modelling system (WRF-10 Hydro, Gochis et al., 2015) and applies a land surface scheme to calculate hydrological processes, with runoff converted to streamflow using hydrological routing. Such an approach was initially developed for Numerical Weather Prediction (NWP) despite limitations in the representation of spatial variability of hydrological processes in such schemes (Clark et al., 2015).
WRF-Hydro serves as the core model for the United States National Water Model (Maidment, 2017). GEM-Hydro has been evaluated in details in the Great Lakes region (Gaborit et al., 2017) and is a component of the Water Cycle Prediction 15 System which provides daily operational forecasts of streamflow and lakes level in this region (Durnford et al., 2018).
However, the skill of GEM-Hydro during extreme rainfall and snow melting events in mountainous terrain is still largely unknown.
The primary objective of this study is to evaluate the ability of GEM-Hydro to reproduce the observed hydrology in hindcast mode for the June 2013 flood in southern Alberta. Another main objective is to assess the sensitivity of hydrological 20 simulations to the main factors controlling flood dynamics, such as precipitation forcing accuracy, initial snowpack conditions and river routing parameters. For this purpose, a specific configuration of the Canadian NWP model GEM (Côté et al., 1998;Girard et al., 2014) was deployed over the region to produce atmospheric forcing at resolutions ranging from 1 to 10 km to drive GEM-Hydro. In particular, the precipitation simulated by GEM was combined with precipitation gauge measurements from different observation networks using the Canadian Precipitation Analysis (CaPA) system (Mahfouf et 25 al., 2007;Lespinas et al., 2015;Fortin et al., 2018) to generate several quantitative precipitation estimation (QPE) products for the flood. These products were used to assess how QPE accuracy influences the uncertainty in simulating streamflow from the mountain headwaters downriver to cities that were impacted by the flood such as Calgary and High River. The role of initial snowpack conditions at high elevations prior to the flood and uncertainties in river routing parameters were also considered. 30 The paper is organized as follows. Section 2.1 presents the study area and the different datasets used in our study. It also describes the GEM-Hydro modelling platform and the configurations of the different experiments carried out with this model in our study. Section 3 evaluates the different hydrological inputs used in this study (initial snowpack conditions, QPE Hydrol. Earth Syst. Sci. Discuss., https://doi.org/10.5194/hess-2019-152 Manuscript under review for journal Hydrol. Earth Syst. Sci. Discussion started: 23 April 2019 c Author(s) 2019. CC BY 4.0 License. products) and the resulting hydrological simulations. Section 4 contains a discussion of the main results of this study.
Finally, concluding remarks are presented in Section 5.

Study area and data
The study area covers the three main river basins in southern Alberta that were strongly impacted by the June 2013 flood: the 5 Red Deer, Bow and Oldman river basins (Fig. 1). These rivers drain the Rocky Mountains and their foothills, flowing eastward towards the Canadian Prairies and eventually joining to form the South Saskatchewan River. The black box in Fig.   1 highlights the location of the Elbow and Highwood rivers -both spilled over their banks during the flood, affecting the cities of Calgary and High River.
Precipitation data were obtained from four different networks (Fig. 1). Automatic synoptic stations (SYNOP network) 10 maintained by Environment and Climate Change Canada and other organizations such as Alberta Agriculture and Forestry provide a good coverage over the prairies and the forested foothills but only a sparse coverage of the Canadian Rockies, particularly in Banff National Park where the stations are mainly located in the valleys (see for example the upper Bow river basin, Fig. 1). Therefore, to fill in data gaps, precipitation data were taken from Alberta Environment and Parks because of their good coverage of the higher elevations of the Canadian Rockies headwaters. Data from high elevation mountain 15 stations of the Canadian Rockies Hydrological Observatory (CRHO, University of Saskatchewan, https://www.usask.ca/hydrology/CRHO.php) were also used in the upper Bow River Basin. Finally, stations from the American Cooperative Observer Network (COOP) reporting 6-h precipitation amounts were included. Theses stations are mainly located in the US and on the western side of the study area and were referred as stations from the SHEF network (Standard Hydrometeorological Exchange Format) in Lespinas et al. (2015). They were not included in the 6-h operational 20 precipitation analysis at the time of the flood.
Information on snow conditions before the flooding event and during the winter of 2012/2013 were obtained from 11 automatic snow pillows from Alberta Environment and Parks located in southern Alberta. These stations measure hourly snow water equivalent (SWE). Additionally, outputs from the SNOw Data Assimilation System (SNODAS) from the US National Operational Hydrologic Remote Sensing Center (NOHRSC) were included in this study. SNODAS estimates 25 various snow properties (including SWE and snow depth) by merging satellite, airborne, and ground-based snow data with a numerical simulation of snowcover (Barrett, 2003). SNODAS data are available at 1-km spatial resolution and 24-hour temporal resolution and cover continental US as well as part of Canada (up to 54 ∘ N). In southern Alberta, the snow pillows from Alberta Environment and Parks are included in SNODAS since they provide relevant information for mountain snow conditions, including those affecting the Columbia River flowing from Canada to the US. The ECCC snow analysis 30 (Brasnett, 1999) was not considered in our study since its spatial resolution (10 km) is too coarse to accurately represent snow conditions in the complex topography of the Canadian Rockies. Moreover, the ECCC snow analysis only contains information on snow depth and therefore lacks the crucial information that SWE is for hydrology.
Hourly streamflow time series were obtained from the Water Survey of Canada for stations located in the Jumpingpound Creek, Elbow River and Highwood River basins (Fig. 2). These rivers were selected for the evaluation of GEM-Hydro simulations since they were strongly impacted by the flood and are only slightly affected by regulation. Two stations were 5 destroyed during the flooding event and removed from the study. Overall, 10 stations were kept for the evaluation of hydrological simulations (Fig. 2).

The GEM-Hydro modelling platform
GEM-Hydro is a distributed runoff-modelling platform developed at ECCC for hydrological modelling across Canada (Gaborit et al., 2017). It is composed of two components: (i) the GEM-Surf surface prediction system (Bernier et al. 10 2011) and (ii) the WATROUTE routing scheme (Kouwen, 2010). GEM-Hydro has been extensively evaluated over the Great Lakes (Gaborit et al., 2017) and is a component of the Water Cycle Prediction System running operationally at ECCC over the Great Lakes and St. Lawrence River watershed (Durnford et al., 2018).
GEM-Surf simulates the evolution of six types of surfaces: land, glacier, inland water, sea, ice (over sea or lake) and urban areas. Over land, it uses the SVS (Soil, Vegetation and Snow) land surface scheme Husain et al., 2016). 15 SVS relies on a multiple energy budget approach to solve the energy and mass balance at the earth surface. The land tile is divided into three main covers: bare ground, low vegetation and high vegetation. Then, SVS solves independent energy budgets using a force restore approach for bare ground, vegetation, snow covering bare ground and low vegetation, and snow below high vegetation . SVS uses a single-layer energy balance scheme to simulate the snowpack evolution. Hydrological processes are represented in SVS using a multi-layer scheme and assuming Darcian flow between 20 the soil layers . SVS simulates surface runoff using lateral flow for sloping surfaces (Soulis et al., 2011) and base flow at the bottom layer. WATROUTE routes surface runoff, lateral flow and base flow produced by SVS to the basin outlet. In WATROUTE, runoff and lateral flow directly feed the streams while base flow is provided to a lower zone storage (LZS) compartment, representing surficial aquifers, which releases water to the streams.

Model configuration and experiments 25
In this study, GEM-Hydro was used in a stand-alone (offline) mode to produce a hindcast of the hydrology of the June 2013 flood. The model was driven by a new set of atmospheric data combining short-term numerical weather forecast and precipitation analysis. These atmospheric forcings were produced at different resolutions to assess the sensitivity of GEM-Hydro to the resolution of forcing data. The impact of initial snow conditions on streamflow simulation was also assessed. In the following, the main features of the new atmospheric forcing data are detailed, as well as the different experiments carried 30 out with GEM-Hydro.

Atmospheric simulations
The latest operational version of the GEM NWP model used at ECCC (Côté et al., 1998;Girard et al., 2014;Caron et al., 2015;Milbrandt et al., 2016) was used to recreate the atmospheric conditions during the flooding event. GEM was configured with three one-way nested domains with 10, 2.5 and 1-km grid spacing centred over the study area in southern Alberta. This nested approach has been applied in previous studies at ECCC (e.g. Vionnet et al., 2015;Leroyer et al., 2018). 5 Two domains at the kilometre scale (2.5 and 1 km) were used since the previous studies of Li et al. (2017) and Milrad et al. (2017) showed that models at convection-permitting resolution performed best for this event. Details of the simulation domains are listed in Table 1. Figure 3 shows the extents of the 2.5 and 1-km grids over western Canada. Table 1 Morrison and Milbrandt, 2015).
Land surface characteristics were specified at each grid cell of each domain using different geophysical datasets described in Table 2.

Precipitation analysis
To test their influence on the hydrology of the June 2013 flooding event, as simulated by GEM-Hydro, new QPE products were generated. They relied upon the Canadian Precipitation Analysis (CaPA) system (Mahfouf et al., 2007;Lespinas et al., 2015;Fortin et al., 2018). CaPA combines precipitation observations with a background field obtained from a short-term meteorological forecast using optimal interpolation to produce a QPE on a regular grid. Precipitation observations consisted 25 of rain gauges and radar QPEs . The current operational version of CaPA at ECCC produces 24-h and 6-h QPE at 10-km resolution over a domain covering most of North America and at 2.5-km resolution over a domain covering most of Canada (Fortin et al., 2018). km. These additional stations substantially increase the density of stations in the Rockies and their foothills (Fig. 1). Radar data were not included in the different QPE products since they tend to underestimate precipitation amounts in this mountainous region for this event (Kochtubajda et al., 2016). Table 3

Surface and hydrological modelling
GEM-Hydro was used to simulate the evolution of the surface and hydrological conditions during and after the flooding event (from 18 June 2013 12 UTC to 26 June 2013 12 UTC). GEM-Hydro simulations combine successive integration of GEM-Surf (including SVS) and WATROUTE. In our study, GEM-Surf ran on the same grid and used the same geophysical 15 fields as GEM 1 km (Tables 1 and 2). Seven layers down to 1.4 m depth were used to represent the vertical layering of soil, following Gaborit et al. (2017). The routing with WATROUTE was implemented over a 1-km grid covering the Red Deer, Bow and Oldman river basins (Fig. 1). Flow directions were derived from the HydroSHEDs database (Lehner et al., 2008).
The default (uncalibrated) version of GEM-Hydro was used in the study. Parameters in SVS were the same as in the version of SVS used for NWP Husain et al., 2016). Parameters in WATROUTE were derived using the standard 20 procedures implemented in the operational hydrological system of ECCC (Durnford et al., 2018). In particular, two Manning's roughness parameters are used for flood plain flow, n f , and channel flow, n c , based on values derived from the literature (Chow, 1959). n f varies spatially with vegetation type (0.035 for low vegetation and 0.075 for high vegetation) and seasonally to represent temporal changes in vegetation. For channel flow, the value of n c is affected by the roughness of the channel surface on the bottom and sides as well as the roughness of the underside of the ice cover overlying the river when 25 ice is present. It varies with drainage area between 0.03 and 0.05 to represent the overall decrease of Manning coefficient from the headwaters to the river mouth (Decharme et al. 2010). These variations are effective at large scale and over our region of interest, n c is almost constant, equals to 0.044. Calibrating GEM-Hydro is a challenging task due to model computation time (Gaborit et al., 2017). The main purpose of GEM-Hydro as used by ECCC is to produce hydrological simulations over large river basins including ungauged basins. Therefore, the model has to be able to produce reasonable 30 simulation accuracy with default parameters, because it is difficult to obtain maximized local performances when calibrating a hydrological model over large area. Calibrating it locally on each gauge from upstream to downstream is not feasible in practice due to model computation time, and would violate parameter consistency in space as well as lack robustness in time (Gaborit et al., 2017). For these reasons, GEM-Hydro was not calibrated in this study.
Eight simulations of the flood were carried out with GEM-Hydro. They differ in terms of resolution of the atmospheric driving data, precipitation forcing and initial conditions. Table 4 describes their main characteristics. The atmospheric forcings required for GEM-Surf include shortwave and longwave irradiation at the surface, surface pressure, air temperature, 5 specific humidity, precipitation, and wind. Continuous atmospheric forcings from 18 to 22 June were obtained from successive GEM forecasts at different resolutions (Table 1). For the period of 22 nd to 26 th June, only operational forecasts of the RDPS at 10-km resolution were used as forcings. In all instances, forecast hours 6 to 11 were used to avoid spin-up errors. For the forecasts at 10 and 2.5 km, surface pressure, air temperature, humidity, and precipitation phase were downscaled to correct for elevation differences between the grid of the forecasts and the 1-km grid used by GEM-Surf as in 10 Bernier et al. (2011). The different precipitation datasets produced with CaPA were used as precipitation forcing data (Table 3). Six-hourly precipitation outputs from CaPA were then disaggregated into hourly accumulations in accordance with the temporal precipitation structure of the GEM forecast.
Initial conditions for the surface variables on 18 June 2013 12 UTC were taken from two different 1-km GEM-Surf experiments running. Both experiments were initialized on June 1 st , 2012 using RDPS 10-km surface fields interpolated to 15 the 1-km grid and then integrated for more than one year until 18 June 2013 12 UTC to allow a sufficient spin-up of surface variables, especially soil moisture. The hourly meteorological forcings for this one-year spin-up period were taken from successive RDPS forecasts at 10-km resolution issued 4 times per day (00, 06, 12 and 18 UTC). Similarly to all other experiments, 6-11 forecast hours were used. The meteorological forcing was then downscaled to the 1-km GEM-Surf grid using the method described in Bernier et al. (2011). Following Carrera et al. (2010), precipitation was taken from the 20 RDPS from 1 November 2012 to 30 April 2013 and from the operational version of CaPA at 10 km for the rest of the period.
CaPA underestimates solid precipitation amounts (Schirmer and Jamieson, 2015), since the precipitation gauges used in CaPA are prone to wind-induced undercatch of winter precipitation (Fortin et al., 2018).
In the first experiment, OPL, simulated snow conditions evolved in response to the atmospheric forcing without constraint from observations. Using a similar setup, Carrera et al. (2010) and Separovic et al. (2014) have shown that GEM-Surf tends 25 to underestimate snow water equivalent on the ground (SWE) since the 10-km forcing from RDPS does not resolve local topographic snowfall enhancement. For this reason, a second experiment, SND, used estimated SNODAS SWE to correct GEM-Surf modelled SWE over the region (Sect. 2.1). The correction was carried out on May 1 st 2013, close to peak SWE to adjust for over-winter bias in accumulated precipitation. On that date, simulated snow depth in GEM-Surf was used to estimate SWE using simulated snow density, and this SWE was adjusted to match SNODAS SWE. After insertion on May 30 1 st , the SND simulation kept running until 18 June 2013 12 UTC using the same atmospheric forcing as simulation OPL. The  Very little snow remained in the OPL simulation, which is inconsistent with the presence of high mountain snow reported prior to the flood (Liu et al., 2016;Pomeroy et al., 2016a,b). On the other hand, the SND simulation predicted the persistence of snowpacks at high elevations prior to the event, mainly in the Upper Bow River Basin (detailed topographical heights are 15 given in Fig. 1) which is consistent with observations in the region.  (Milrad et al., 2015) and to evaluate atmospheric models (Liu et al., 2016;Milrad et al., 2017;Teufel et al., 2017). Stations from the CRHO, COOP and ABE networks were not included in this reference precipitation analysis and were used on Fig. 5b to provide an independent evaluation of CaPA 10 km Def over the Rockies and their foothills. 25

Quantitative precipitation estimation
Cumulative precipitation was underestimated at most of the stations, in particular in the upper Elbow and Highwood river basins.
The introduction of CRHO, COOP and ABE stations in the precipitation analysis (CaPA 10 km New) modified the pattern of cumulative precipitation during the flood ( Fig. 6a and d), even as the precipitation background at 10-km grid spacing remained the same between the two simulations (Table 3) Fig. 6. Additional details associated with the influence of the topography on the precipitation background are seen at 2.5 km (Fig. 6c) and at 1 km (Fig. 6e) in areas of complex terrain, especially in the Bow River Basin. 5 Figure 7 shows an evaluation of the different QPE products against CRHO, COOP and ABE station observations. For CaPA 10 km Def , which does not assimilate these 3 networks, this is an independent evaluation. For the other datasets, a "leaveone-out" cross validation method was used, where stations are removed one by one and an analysis value is estimated at their location using stations in the vicinity, following Lespinas et al. (2015). This analysis value was then compared to that observed. Figure 7a shows that CaPA 10 km Def captured the general tendency of observed cumulative precipitation but 10 missed the highest precipitation amounts, in particular in the Bow River Basin (Fig. 5b). This led to a large negative bias in precipitation. Performances were strongly improved with CaPA 10 km New at the same horizontal resolution (Fig. 7b).
Cumulative precipitation greater than 200 mm was better captured but CaPA 10 km New still had a negative bias, which was removed with CaPA 2.5 km and CaPA 1.0 km. Overall, CaPA 2.5 km showed the best performance with a lower bias and RMSE, and higher correlation coefficient (Fig. 7c). Compared to CaPA 2.5 km, CaPA 1.0 km tended to overestimate 15 cumulative precipitation greater than 200 mm for several stations of the Bow River Basin, which led to a slight positive bias.
Improved performance in estimating cumulative precipitation during the flood is not sufficient alone for accurate hydrological simulations, and the precipitation forcing should also reproduce the 6-h precipitation amounts well. Figure 8 uses quantile-quantile (QQ) plots to compare the distribution of observed and analysed precipitation for three datasets (CaPA 10 km New, CaPA 2.5 km and CaPA 1.0 km) and for the three main river basins in southern Alberta. As in Fig. 7, analysed 20 precipitation is derived from the leave-one out method. The concordance correlation coefficient (Lawrence and Lin 1989), , is used to estimate the agreement between analyzed and observed 6-h precipitation amounts. ranges from -1 to 1. 1 (-1) indicates perfect concordance (discordance) whereas 0 indicates no correlation. Over the Bow River Basin, all datasets reproduced the distribution of 6-h precipitation amounts well, except for amounts above 60 mm. For this basin, the best performance was obtained with CaPA 1.0 km (especially over 40 mm). The same conclusion can be made for the Oldman 25 River Basin where CaPA 10 km New and CaPA 2.5 km underestimated the 6-h precipitation depths greater than 35 mm.
Finally, the three version of CaPA has similar performance over the Red Deer River Basin, where no values greater than 40 mm were observed. Overall, Fig. 8 illustrates the added value of precipitation analysis at the kilometre scale for better capturing the intensity of extreme mountain precipitation that may be important for hydrological simulations.

Hydrological simulations 30
Initial snow and soil conditions from the OPL and SND simulations were used in combination with the atmospheric forcing and the QPE products at different resolutions to provide initial and boundary conditions for eight hydrological simulations of  (Table 4). These simulations were evaluated using discharge measurements from stations on Jumpingpound Creek, Elbow River and Highwood River (Fig. 2). Figure 9 shows the percent bias (PBIAS) of the simulated streamflow volumes for all eight hydrological simulations of the flood, against observed discharge volumes from 20 June 0 UTC to 25 June 0 UTC for different hydrometric stations in our catchments of interest. PBIAS is used to assess the simulation's overall water budget fit during the flooding event: a negative (positive) value denotes a general 5 tendency to under-(over-)estimate streamflow (discharge) volume during the flooding event. Figure 9 shows a strong underestimation of flood discharge volume for all three gauging locations (station 05BH015 for Jumpingpound Creek, station 05BJ010 for the Elbow River and station 05BL024 for the Highwood River) in simulation 10km_Def_OPL using initial conditions from the OPL run. This underestimation was also observed for all other hydrometric stations shown in Fig. 9 and is consistent with the general underestimation of precipitation in CaPA 10 km Def over the 10 foothills and the mountainous part of these basins ( Fig. 5b and 7a). Including additional stations and keeping the same resolution in the precipitation analysis (experiment 10km_New_OPL) leads to a reduction of the negative bias in flood discharge volume. For example, at the outlet of the Elbow River (station 05BJ010), PBIAS decreased from -56% with 10km_Def_OPL to -24% with 10km_New_OPL. Increasing the horizontal resolution of the precipitation analysis and of the atmospheric forcing from 10 to 2.5 km (experiments 2.5km_OPL) brought additional improvements and contributed to 15 reducing the negative PBIAS for discharge at these three river gauging stations. Finally, the 1.0km_OPL simulation gave similar results with to that of 2.5km_OPL for these hydrometric stations The four simulations using initial conditions from the OPL run all show an underestimation of the flood discharge volume for gauged sites in three river basins, especially for the Elbow River (-20% in experiment 1.0km_OPL) and the Highwood River (-28% in experiment 1.0km_OPL). This may be explained by the underestimation of initial SWE in the basin 20 headwaters in the OPL simulation as shown in Fig. 4. For this reason, the SND simulation, which incorporates SNODAS SWE to estimate peak snow accumulation, was considered as an alternative to obtain initial snow and soil conditions. Figure 9 shows that this substitution has no impact on hydrological simulations for Jumpingpound Creek as no snow was present in this drainage basin on 18 June 2013 in either simulation (Fig. 4). In contrast, the SND simulation with SNODASinformed initial conditions improved the estimation of flood volume for the Highwood River (station 05BL024) in all 25 simulations with a decrease in the negative bias and nearly unbiased estimates from the 2.5km_SND and 1.0km_SND simulations. For the Elbow River, results are different and using SND simulations led to an overestimation of the flood volume in all the experiments using the CaPA analysis with additional stations. It should be noted that the SNODAS analysis overestimated SWE close to peak accumulation at the snow pillow located in the headwaters of the Highwood River Basin (Fig. 4). This overestimation of flood discharge volume when using run SND was also evident when simulations were 30 compared to observations from the hydrometric stations located in the headwaters of the river basins, station 05BJ004 for the Elbow River and station 05BL019 for the Highwood River. The streamflow discharge hydrographs during the flood are shown in Fig. 10 for two stations on the Elbow River and in Fig. 11 for two stations on the Highwood River. For each river basin, one station located in the headwater of the catchment and one near the outlet were chosen for analysis. These figures compare the evolution of basin-averaged SWE, basinaveraged cumulated precipitation and streamflow for four GEM-Hydro simulations: 10km_Def_OPL, 10km_Def_SND, 1.0km_OPL, 1.0km_SND. The filled area on the hydrographs (Fig. 10d, f and Fig. 11d, f) illustrates the impact of using 5 initial conditions from run SND instead of run OPL for a given meteorological forcing. For all stations, the 10km_Def_OPL and 10km_Def_SND simulations failed to capture streamflow dynamics and strongly underestimated the peak flow.
Improved results were obtained using the precipitation analysis and the meteorological forcing at 1 km. For the Elbow River ( Fig. 10), simulated peak flow was still underestimated with the 1.0km_SND simulation providing the best results in terms of magnitude. It also gives the best results in terms of timing but still remains late compared to the hydrometric observations 10 (5h and 6h at station 05BJ004 and 05BJ010, respectively). While very little snow is present in simulations using the OPL simulation for initial conditions, the decrease in basin-averaged snow mass between 18 and 25 June 2013 reached -62% in the 1.0km_SND simulation and -56% in the 10km_Def_SND simulation at station 05BJ004. For the Highwood River ( Fig. 11), the 1.0km_SND simulation overestimated the peak flow at both stations. This overestimation reached 59 % at station 05BL019 (Fig. 11f). The basin draining to this station experienced both intense precipitation and rapid snowmelt on 15 20 June 2013. The timing of the peak flow was delayed by 11 and 14 hours in 1.0km_SND and 1.0km_OPL simulations respectively.
The delay in simulated peak flow is found at the outlets of the Elbow River and the Highwood River for all combinations of initial snow conditions and QPE. This suggests that uncertainties in modelling of river routing with WATROUTE potentially affects the simulations. In particular, erosion and changes in river channels were reported for many rivers impacted by the 20 June 2013 flood (Pomeroy et al., 2016b) and a shift from sub-surface to overland flow was noted (Pomeroy et al., 2016a;Fang and Pomeroy, 2016), suggesting that the default routing parameters used in WATROUTE may not be suitable for this extreme event. To test this assumption, additional routing experiments were made using the 1.0km_SND and 1.0km_OPL simulations that benefited from the most accurate precipitation forcing ( Fig. 7 and 8). The Manning's coefficients over the GEM-Hydro simulation domain were successively decreased by 10, 20, 30 and 40% and the timing of the simulated peak 25 flow compared to discharge observations. Figure 12 compared the hydrographs obtained in the default configuration with those where Manning coefficients were reduced by 30 % because this gave the best peak flow simulations. For the Elbow River ( Fig. 12a and b), adjusting the Manning coefficients allowed GEM-Hydro to better capture the rapid increase in streamflow observed at both stations. It also improved the timing of the peak flow, for example the delays were reduced from + 10 (+ 5) hours to +1 (-2) hours in the 1.0km_OPL (1.0km_SND) simulations at station 05BJ010. The magnitude of 30 the peak flow was only slightly impacted by decreasing Manning's coefficient. A similar improvement in the simulated timing of the peak flow was found for the Highwood River (Fig. 12d)  by 33% compared to 14% in the default configuration. An early and overestimated peak flow was also found in the headwaters of the Highwood River Basin (Fig. 12c) with the 1.0km_SND simulation when using updated Manning's coefficients. This shows that simple, spatially uniform adjustment of Manning's coefficients is not sufficient to capture the peak flow timing and magnitudes observed at various hydrometric stations. This is likely due to the substantial spatial and temporal variability in runoff, storage and channel processes (Fang and Pomeroy, 2016) and that such a hydraulic calibration 5 is necessarily combined with other sources of uncertainties in the timing and magnitude of intense precipitation and snowmelt in the land surface scheme, and in the other hydraulic parameters of the routine scheme such as channel width and geometry.

Discussion
This paper has examined the ability of the GEM-Hydro hydrological modelling platform to simulate streamflows during the 10 June 2013 Flood in southern Alberta. This is the first time GEM-Hydro has been evaluated in this region of Canada, which is characterized by a transition from headwaters in the complex topography of the Canadian Rockies towards downstream plains. The simulation of this extreme event offers insights into the key factors governing the prediction skill of flooding events in this region and the associated modelling challenges.

Quantitative Precipitation Estimation 15
This study used four QPE products to assess how QPE accuracy influences GEM-Hydro simulations. The QPE products were generated using the Canadian Precipitation Analysis (CaPA) system with varying station density and different horizontal resolutions of the GEM precipitation background. The results of this study show that the QPE product that most resembles the CaPA product that was generated in operations at ECCC at the time of the flood (10-km grid spacing, SYNOP stations only) led to a strong underestimation of precipitation in the mountainous parts of the region and a systematic 20 underestimation of the flood volume in the Elbow and Highwood River catchments. This version of CaPA was used to evaluate precipitation simulated by atmospheric models during the June 2013 Flood in several previous studies (Liu et al., 2016;Teufel et al., 2017). Including additional stations from university and provincial networks led to significant improvements in the 10-km precipitation analysis for the headwaters that received the largest amount of precipitation during the flood with a positive impact on the hydrological response simulated by GEM-Hydro. This illustrates the continuous need 25 to include stations from additional networks in operational precipitation analysis systems used for hydrological modelling as already discussed in several studies (Lespinas et al., 2015;Soci et al., 2016;Fortin et al., 2018). This need is critical in the mountainous regions of Western Canada where the network of stations used in the operational version of CaPA (Fortin et al., 2018) is scarce with stations mainly located at low elevation valley bottoms. Radar data were not considered in our study since Kochtubajda et al. (2016) showed that these data were affected by attenuation during the heavy precipitation periods of 30 the flooding event and produced rate and accumulation lower than gauges estimates in mountainous areas. Increasing the resolution of the GEM precipitation background from 10 km to kilometre scales (2.5 and 1 km) improved the accuracy of the final QPE products in terms of cumulative precipitation as well as the intensity of extreme precipitation during the flooding event. Streamflow simulated by GEM-Hydro benefited from these improvements in forcing data. The largest difference was found in moving from 10-km QPE to the 2.5 km QPE which invoked precipitation forcing data from a convection-permitting NWP system. Li et al. (2017) and Milrad et al. (2017) have shown that atmospheric models at such 5 resolution were required to capture the spatial and temporal variability of extreme precipitation during the June 2013 flooding event. This has also been shown in other mountainous regions (e.g. Richard et al., 2007;Weusthoff et al., 2010;Rasmussen et al., 2011). The results shown here confirm that improving the realism of the precipitation model in complex terrain also improves the precipitation analysis. This shows the potential to improve hydrological forecasting by using the latest operational version of CaPA running at 2.5 km over Canada (Fortin et al., 2018). On the other hand, only marginal 10 improvements in the intensity of extreme precipitation and slight degradation of cumulative precipitation predictions were found when increasing the resolution of the precipitation model from 2.5 km to 1 km. Theses findings are in agreement with Pontoppidan et al. (2017) who showed only minor improvements from 3 km to 1 km in atmospheric simulations of an intense precipitation event in Western Norway. Further evaluations are required in the mountains of Western Canada to clearly assess the benefits and limitations of a 1-km precipitation analysis. In this context, further developments are also 15 required in CaPA to better account for the effect of complex topography in the optimal interpolation scheme (Daly, et al., 1994;Gottardi et al., 2012).

Initial snow conditions
The June 2013 Flood was characterized by rain falling on a late-lying snowpack at high elevations inducing rain-onsnowmelt and enhancing runoff generation such that there was more runoff than storm precipitation (Fang and Pomeroy, 20 2016). The results shown here confirm the impact of initial snow conditions on simulated flood volume and dynamics, in agreement with previous studies (Pomeroy et al., 2016a;Teufel et al. 2017) and illustrate the challenge of obtaining realistic information on snow conditions in the Canadian Rockies. Indeed, the first set of hydrological simulations used initial snow conditions taken from an open-loop simulation of GEM-Surf and presented a systematic underestimation of flood volume in GEM-Hydro simulations with all QPE products. This is mainly due to the absence of snow in the high elevation headwaters 25 at the time of the flood, which resulted from a strong and systematic underestimation of winter snow accumulation. These results are consistent with the earlier findings of Carrera et al. (2010) and Separovic et al. (2014) and are mainly explained by an underestimation of mountain winter precipitation simulated by GEM at 10-km grid spacing (Schirmer and Jamieson, 2015). This illustrates the challenge of obtaining accurate snow simulations in mountainous terrain when relying on precipitation forcing made of successive forecasts (e.g. Bellaire et al., 2011;Vionnet et al. 2016). In addition, GEM-Surf at 30 1-km grid spacing does not simulate key processes such as wind-induced and gravitational snow transport (e.g. Bernhardt and Schulz, 2010;Vionnet et al., 2014;Musselman et al., 2015) or the effect of slope on incoming radiation (e.g. Oliphant et Hydrol. Earth Syst. Sci. Discuss., https://doi.org/10.5194/hess-2019-152 Manuscript under review for journal Hydrol. Earth Syst. Sci. Discussion started: 23 April 2019 c Author(s) 2019. CC BY 4.0 License. al., 2003) that influence winter snow accumulation and runoff generation in mountainous catchment (DeBeer and Pomeroy, 2017;Brauchli et al., 2017).
The insertion of SWE data from the SNODAS operational snow analysis system (Barrett, 2003) near the time of peak snow accumulation was tested to correct the negative bias in winter snow accumulation and to obtain more realistic snow and soil initial condition at high elevation. GEM-Hydro simulations including SNODAS SWE showed contrasting results. In 5 particular, the overestimation of flood volume found at the stations located in the headwaters of the Elbow and Highwood River basins suggests an overestimation of basin-average SWE prior to the flood. This is in agreement with the overestimation of SWE close to peak accumulation reported using snow pillows observations at high elevations in these basins. Teufel et al. (2017)

Model uncertainties and limitations
The analysis of the hourly flood hydrograph simulations highlights the crucial roles of precipitation forcing, initial snow conditions and hydraulic parameter selection in accurate flood simulations. Using the 1-km QPE as driving data, a reduction of Manning's equation coefficients by 30% compared to the default values in GEM-Hydro was required to capture the timing of peak flow at the outlets of the Elbow River and Highwood River. For these two rivers, the timing of the peak flow 20 was more sensitive than the magnitude to adjustment of Manning's equation coefficients. In a recent study, Lin et al. (2018) demonstrated that the flood prediction skill of the US National Water Model (Maidment, 2017) is affected by flow velocity in the river network. The results shown here suggest that the set of Manning's coefficients used in the default version of GEM-Hydro would benefit from specific selection to individual rivers in southern Alberta. The channel area in WATROUTE, that controls whether flow is contained in the main channel or the overbank part, may need to be selected by 25 river as well. Calibration these parameters based on historical data may not be suitable for an extreme event such as the June 2013 flood when major changes occurred in the geometry of the riverbeds (Pomeroy et al., 2016b), but parameter selection from hydraulic characteristics measured by the Water Survey of Canada and from remote sensing and field surveys is an alternative parameterisation methodology used at multiple scales in the region .
In addition to uncertainties in the routing scheme parameters, GEM-Hydro has several limitations in its SVS land surface 30 scheme for application to mountain and cold-region hydrology. SVS uses a 1-layer snowpack scheme to simulate the snow Hydrol. Earth Syst. Sci. Discuss., https://doi.org/10.5194/hess-2019-152 Manuscript under review for journal Hydrol. Earth Syst. Sci. Discussion started: 23 April 2019 c Author(s) 2019. CC BY 4.0 License. cover evolution over bare ground and low vegetation and below high-vegetation. Due to simplification in the representation of water flow through flat and sloping snowpacks, internal snowpack coupled phase change and energy exchange and the crude turbulent transfer scheme employed, such a scheme will be challenged by rain-on-snow events (Würzer et al. 2016) and by the occurrence of preferential flow pathways, layering and refreezing in snowpacks (Leroux and Pomeroy, 2018;. SVS snow neglects surface layer snow energetics, canopy snow interception, melt, unloading and sublimation, and 5 blowing snow transport and sublimation processes which are necessary for successful simulations of snowpacks in this region (Pomeroy et al., 2008;Fang et al., 2013;Pomeroy et al., 2016a;2016c). SVS considers all topography to be flat meaning that its snowpack energetics in mountains will vary from those observed not only in magnitude but in sign (Pomeroy et al., 2003;Ellis et al., 2013). In addition, SVS does not simulate soil freezing and its strong impact on infiltration to frozen soils (e.g. Zhao and Gray, 1999;Gray et al., 2001) and on hillslope hydrological processes (Carey et al., 10 2007). Pomeroy et al. (2016b) and Teufel et al. (2017) suggested that the presence of frozen soil at high elevations and in valley bottoms modified snowmelt infiltration and contributed to the runoff generation during the June 2013 flooding event.
To overcome some of these limitations, ECCC is currently working on the development of a new version of SVS that will solve the heat diffusion equation in the soil including freeze/thaw processes and include a more advanced multi-layer snowpack scheme (Decharme et al., 2016). 15

Conclusions
This study aims to evaluate the ability of the GEM-Hydro modelling platform to simulate (in hindcast mode) hydrometeorological conditions during the June 2013 flood in southern Alberta and to identify the main factors important for accurately reproducing the flood dynamics. Toward this goal, a 1-km resolution configuration of the model was used in a way similar to the standard coarser resolution configurations used at ECCC. The GEM atmospheric model was used in 20 combination with the CaPA precipitation analysis system to generate several QPE products for the flood with varying station densities (especially in the mountainous part of the catchments) and different horizontal resolutions (10, 2.5 and 1 km). In addition to the QPE products, two sets of initial snow conditions were also considered to assess their impact on hydrological simulations in the Elbow and Highwood river basins, which were strongly impacted during the flooding event. snow accumulation provided initial conditions with substantial snowpacks and coverage at high-elevation. The use of SNOWDAS SWE led to contrasting abilities to simulate flood discharge volumes and a consistent overestimation in the mountainous part of the catchments when combined with the most accurate QPE. These results confirm the importance of accurate initial snow conditions for accurate spring and early summer flood simulation in this region and shows the need for more accurate snow products in the Canadian Rockies. 5 • The model suffers from poor prediction skill for the timing of the peak flow with a systematic delay. By conducting a series of model sensitivity tests, it was shown that improved flood peak timing can be obtained by tuning Mannings's coefficient in the routing scheme. However, uniform tuning of this parameter was not successful in estimating both peak flow timing and magnitude. This clearly highlights large uncertainties in routing parameters during extreme events and the need for revised and specific routing parameters for the rivers of this region. 10 This study is the first evaluation of the flood prediction skill of GEM-Hydro in mountainous terrain and downstream. The hydrological simulations analysed here were driven by QPE products obtained with a precipitation analysis system. The results are promising but further evaluation is required for such extreme events with GEM Hydro driven by deterministic or ensemble precipitation forecast from NWP systems. GEM-Hydro as a hydrological forecasting system is currently being implemented by ECCC in Western Canada. Improvements of the forecasting system are can be expected with future 15 development of a new reference snow product for the Canadian Rockies and the development of a new version of the land surface scheme SVS to better simulate the snow and frozen ground processes impacting cold regions hydrology.

Data availability
The CaPA and GEM data at different resolutions are publicly available on the Federated Research Data Repository (Vionnet et al., 2019). 20

Authors contribution
VV, VF and JP designed the study. VV run the numerical experiments with GEM, CaPA, and GEM-Hydro and was responsible for the preparation of the manuscript. GR designed the CaPA experiments. VF, NG and EG contributed to the preparation of the experiments with GEM-Hydro and the analysis of the results. MA provided the latest version of the landsurface scheme SVS and contributed to the preparation of configuration of SVS used in this paper. All authors contributed to 25 the preparation of the manuscript.