Start Submission

Reading: Hydrological modeling of freshwater discharge into Hudson Bay using HYPE


A- A+
Alt. Display

Research Article

Hydrological modeling of freshwater discharge into Hudson Bay using HYPE


Tricia A. Stadnyk ,

Department of Geography, University of Calgary, Calgary, AB; Environmental Science and Engineering Program, University of Northern British Columbia, Prince George, BC, CA
X close

Matthew K. MacDonald,

Ontario Power Generation, Toronto, ON, CA
X close

Andrew Tefs,

Department of Geography, University of Calgary, Calgary, AB, CA
X close

Stephen J. Déry,

Environmental Science and Engineering Program, University of Northern British Columbia, Prince George, BC, CA
X close

Kristina Koenig,

Manitoba Hydro, Winnipeg, Manitoba, CA
X close

David Gustafsson,

Swedish Meteorological and Hydrological Institute, Norrköping, SE
X close

Kristina Isberg,

Swedish Meteorological and Hydrological Institute, Norrköping, KE
X close

Berit Arheimer

Swedish Meteorological and Hydrological Institute, Norrköping, SE
X close


This study details the enhancement and calibration of the Arctic implementation of the HYdrological Predictions for the Environment (HYPE) hydrological model established for the BaySys group of projects to produce freshwater discharge scenarios for the Hudson Bay Drainage Basin (HBDB). The challenge in producing estimates of freshwater discharge for the HBDB is that it spans over a third of Canada’s continental landmass and is 40% ungauged. Scenarios for BaySys require the separation between human and climate interactions, specifically the separation of regulated river discharge from a natural, climate-driven response. We present three key improvements to the modelling system required to support the identification of natural from anthropogenic impacts: representation of prairie disconnected landscapes (i.e., non-contributing areas), a method to generalize lake storage-discharge parameters across large regions, and frozen soil modifications. Additionally, a unique approach to account for irregular hydrometric gauge density across the basins during model calibration is presented that avoids overfitting parameters to the densely gauged southern regions. We summarize our methodologies used to facilitate improved separation of human and climate driven impacts to streamflow within the basin and outline the baseline discharge simulations used for the BaySys group of projects. Challenges remain for modeling the most northern reaches of the basin, and in the lake-dominated watersheds. The techniques presented in this work, particularly the lake and flow signature clusters, may be applied to other high latitude, ungauged Arctic basins. Discharge simulations are subsequently used as input data for oceanographic, biogeochemical, and ecosystem studies across the HBDB.

Knowledge Domain: Ecology
How to Cite: Stadnyk, T.A., MacDonald, M.K., Tefs, A., Déry, S.J., Koenig, K., Gustafsson, D., Isberg, K. and Arheimer, B., 2020. Hydrological modeling of freshwater discharge into Hudson Bay using HYPE. Elem Sci Anth, 8(1), p.43. DOI:
 Published on 19 Aug 2020
 Accepted on 29 Jun 2020            Submitted on 09 Jan 2020
Domain Editor-in-Chief: Steven Allison; University of California Irvine, US
Associate Editor: Julian D. Olden; School of Aquatic & Fishery Sciences, University of Washington, US

1. Introduction

The Hudson Bay Drainage Basin (HBDB) comprises over a third of the Canadian landmass, contains important hydroelectric infrastructure and agricultural lands, and accounts for over a fifth of freshwater exports into the pan-arctic ocean system via the Hudson Bay Complex (HBC comprised of Hudson Bay, James Bay, Ungava Bay, Foxe Basin, and Hudson Strait; McClelland et al., 2006). The HBC itself is important to Arctic marine wildlife, and as a shipping route for Canada. It is a large region of primary production, influenced mainly by the timing of freshwater inputs (largely meltwater). Changing freshwater inputs are essential to the annual formation, decay and break up of sea ice, with seasonal ice cover being important for the HBC pelagic life (Manak and Mysak, 1989; Saucier et al. 1998, 2004). HBDB terrestrial discharge influences Arctic Ocean circulation, sea ice dynamics, biological and biogeochemical processes within the HBC. Due to the importance of terrestrial freshwater fluxes into the HBC for bay-wide circulation and ice-formation, continental-scale hydrologic modeling is key to the combined modeling/observation BaySys group of projects (Barber, 2014).

The importance of the HBDB necessitates a comprehensive hydrological monitoring and prediction framework. Although hydrometric gauges remain relatively dense in southern, lower-latitude regions of the basin, higher-latitude counterparts remain poorly gauged. Additionally, the number of hydrometric gauges has declined across Canada in recent decades (Coulibaly et al., 2013; Mlynowski et al., 2011). As of 2013, 40% of the HBDB is ungauged and 27% has never been gauged. Hydrological modeling is required to fill spatial and temporal gaps in the observational record, and for making long-term streamflow projections.

Over the same period this basin has undergone significant climate change impacting the distribution and variability in hydrology of the basin (MacDonald et al. 2018), there has been increasing regulation imposing a similar trend towards flattening of the average annual hydrograph (Déry et al., 2018). This creates a need to better understand the individual roles of climate change and anthropogenic alteration of the water budget of the basin to properly identify and mitigate contemporary change, but more importantly, for accurate projection of future changes in freshwater discharge. This is a key objective of the BaySys project, and a requirement of any modeling system developed to examine freshwater discharge to Hudson Bay (Barber, 2014). Added to this need is the diverse and complex hydrologic interactions within the basin, which need must be accurately defined.

Across the vast and diverse landscape of the HBDB, water is stored variably in lakes and wetlands, and as snow and ice. These storage dynamics provide several hydrological modeling challenges, such as accurately representing both natural and anthropogenically controled lake storage-discharge dynamics and prairie pothole non-contributing areas (NCAs), as well as the changing role and distribution of frozen soils across the basin.

The Global Lakes and Wetlands Database (Lehner et al., 2011) shows that the HBDB has amongst the greatest concentrations of lakes in the world, with approximately 10% of the HBDB surface area covered by lakes. Lakes modify streamflow by storing and releasing large volumes of water, buffering extreme floods, and increasing the potential for evaporative loss. The multitude of lakes across the HBDB conflicts with the parameterizations found in most hydrological models, which account only for the natural storage-discharge dynamics of the largest lakes (e.g., HYPE (HYdrological Predictions for the Environment; Lindström, 2016), Hydrotel (Fortin et al., 2001) and WATFLOOD (Kouwen, 1988, 2016)). Even more uncommon are parameterizations of storage-discharge relationships in small, internal and regulated lakes (i.e., model subgrid). Several models permit a direct representation of subgrid, internal lakes, yet parameterizing these lakes in numerous ungauged natural basins, and in managed landscapes or subbasins, remains a challenge across the HBDB.

Much of the Nelson River Basin (western HBDB) consists of semi-arid and low relief prairie terrain that generates relatively low runoff volumes. Ubiquitous pothole depressions (e.g., wetlands, sloughs and small natural lakes) are often the terminus of this limited runoff, resulting in numerous internally closed basins (Spence and Woo, 2003). Since prairie streamflow generation depends largely on the interconnectivity of potholes and their connections to river networks, which is also often regulated or controlled through diversions and culverts, streamflow volumes can vary widely from year-to-year (Stichling and Blackwell, 1957). Approximately 8% of the HBDB contains prairie terrain that does not contribute directly to streamflow for a 2-year return period (Prairie Farm Rehabilitation Administration [PFRA], 1983), which poses a significant modeling challenge, particularly because of the variability in connectivity at differing return periods. Early approaches ignored the dynamic nature of the Prairie Pothole Region (PPR) and simply masked out (i.e., removed) published NCAs from the stream network (PFRA, 1983; Wood et al., 1992). More recently, the storage and dynamic interactions between adjacent pothole wetlands have been successfully simulated at small scales using high resolution topographic data (Shook and Pomeroy, 2011; Shook et al., 2013; Huang et al., 2013). Success at modeling the dynamics of the NCA region at relatively coarse resolutions has been achieved using calibrated probability distribution functions of dynamic contributing area models (Muhammad et al., 2019; Mekonnen et al., 2014; Mengistu and Spence, 2016), or with the aid of high-resolution elevation data and a US national wetland inventory (Evenson et al., 2016), whilst others ignore the NCAs altogether at the expense of model accuracy and representativeness.

Prediction of spring melt volumes is additionally complicated by the extensive permafrost and frozen soils underlying the HBDB. During spring melt, frozen soils reduce infiltration into the subsurface, causing relatively more overland flow (Hayashi, 2013). Thawing upper soil layers contribute to runoff as spring warming progresses. Continuous and discontinuous permafrost, common at latitudes above 51°N, reduce groundwater-surface water interactions (Woo and Winter, 1993). Several hydrological models have been improved by including soil freezing processes (Cherkauer and Lettenmaier, 2003; Pomeroy et al., 2007; Wang et al., 2010). Suppressing infiltration under frozen conditions influences runoff timing but can have mixed results on model skill, suggesting that, although frozen soils effectively restrict infiltration at small scales, their parameterization is not always appropriate at larger scales due to the heterogeneous nature of soil properties and flow through macropores (Pitman et al., 1999; Slater et al., 1998).

Given the importance of describing changes in freshwater discharge to Hudson Bay for the BaySys project, and the necessity to separate the relative contributions of anthropogenic and climate driven influences, we develop a hydrologic model of the HBDB. Here we present the development of that hydrologic model, and unique modifications made to account for numerous regulated and natural lakes, prairie potholes or ineffectively drained regions, and underlying frozen soils across the domain. Careful consideration of hydrometric gauge density and proximity is also a consideration within the expansive HBDB domain. Our objective is to summarize enhancements made to the Arctic implementation of the Hydrological Predictions for the Environment (A-HYPE) hydrological model to handle the spatially heterogeneous (natural and anthropogenically controlled) lake-storage discharge relationships, Prairie NCAs and frozen soils. We present an innovative calibration strategy to account for uneven gauge density, and apply it to simulate discharge throughout the HBDB, both in gauged and ungauged watersheds. This work provides the foundation for the BaySys project by producing a reliable hydrological model that is used to provide historic and future projections of runoff and streamflow from the HBDB terrestrial landmass, into the HBC.

2. Hudson Bay Drainage Basin

The HBDB (Figure 1) makes up the landmass between three North American continental divides: the Laurentian to the south, Arctic to the north and Pacific to the west. Runoff from the Canadian provinces of Alberta, Saskatchewan, Manitoba, Ontario and Québec, the Northwest Territories and Nunavut, and four American States (Montana, North Dakota, South Dakota, and Minnesota), all eventually discharge into the HBC. The basin spans 26° latitude, 54° longitude, and rises from sea level to over 3200 m on the eastern slopes of the western Rocky Mountain Ranges.

Figure 1 

Major rivers, lakes, reservoirs and non-contributing area (PFRA, 1983) of the Hudson Bay Drainage Basin. DOI:

The landscape of the HBDB spans seven ecozones: from the southwestern Montane Cordillera, across the semiarid Prairies, to the cooler and wetter Boreal and Hudson Plains, through the Taiga Shield to the Southern and Northern Arctic (Canadian Council of Ecological Areas, 2004). Permafrost is continuous north of the Cape Henrietta Maria area (boundary between Hudson and James Bays), but transitions to sporadic and isolated moving further south toward the Hudson Bay lowlands. A small portion of the basin, in the Canadian Rockies, contains glaciers. A large portion of the HBDB, particularly along the Precambrian Shield and Boreal Forest region, is covered by open surface water in the form of lakes and wetlands.

The HBDB is composed of several regional river systems, which are detailed in Table 1 of Déry et al. (2016). The two largest rivers by volume are the Nelson River, which discharges into the western Hudson Bay, and La Grande Rivière, which discharges into the eastern James Bay. Both of these river systems are heavily fragmented with diversions and storage for hydroelectric generation (Déry et al., 2018). In total, there are 84 reservoirs within the watershed with 303.7 km3 of storage (Lehner et al., 2011). Discharge from regulated rivers accounts for up to 53% of freshwater discharge into Hudson Bay (Déry et al., 2018).

Table 1

Calibration steps and H-HYPE parameters. Process flow diagrams corresponding with each step are provided in Figures 2 and S2 to S4. DOI:

HYPE parameters Description Weighted objective function*

1: Evaporation & Runoff
            kc_corr Adjustment factor to the crop coefficients for PET 0.20 × MNS
0.20 × MCC
0.20 × MRS
0.20 × MRE
0.20 × RRE
            fc_corr Adjustment factor to the fractions of soil available for evaporation
            wp_corr Adjustment factor to the wilting point for all 3 soil layers
            fepotsnow_corr Adjustment factor to the snowpack sublimation coefficient
            srrcs_corr Adjustment factor to the recession coefficient for surface runoff
            rrcs_corr Adjustment factor to the recession coefficient for soil runoff for all 3 soil layers
            ep_corr Adjustment factor to the effective porosity for all 3 soil layers
            cmlt_corr Adjustment factor to the snowmelt (mm °C–1 day–1)
2: Evaporation & Runoff^
            kc × 5 Crop coefficient for PET 0.40 × MNS
0.40 × MRE
0.20 × RRE
            fc × 4 Fraction of soil available for evaporation
3a: Lake Discharge^^
            ilrratk × 3 Exponent for internal lake discharge curve 0.40 × MNS
0.10 × MCC
0.30 × MRS
0.10 × MRE
0.10 × RRE
            ilrratp × 3 Coefficient for internal lake discharge curve
            gicatch × 3 Fraction of runoff directed to internal lake
            olrratk × 4 Exponent for outlet lake discharge curve
            olrratp × 4 Coefficient for outlet lake discharge curve
3b: Dams & Reservoirs
            qprod Production flow (m3 s–1) Manual calibration for individual gauges
            qamp Amplitude of production flow (m3 s–1)
            qpha Phase of production flow (day)
4: Routing
            damp Fractional delay in river 0.60 × MNS
0.10 × MCC
0.20 × MRS
0.10 × MRE
            rivvel Celerity of flood in river (m s–1)

* Objective functions: MNS – mean Nash-Sutcliffe Efficiency (Nash and Sutcliffe, 1970); MKG – median subbasin KGE (Gupta et al., 2009); MCC – mean subbasin Pearson correlation coefficient (component of KGE); MRS – mean subbasin error in standard deviation (component of KGE); MDV – mean subbasin deviation of volumes (component of KGE; Gupta et al., 1999); RDV – regional relative error (data from all stations combined into a single vector).

^ Multipliers for Evaporation and Runoff parameters correspond to the number of significant landcover and soil types requiring a unique parameter calibration.

^^ Multipliers for Lake Discharge correspond to the number of lake types requiring a unique parameter calibration (Section 4.2).

The HBDB spans several climatic zones. Mean annual air temperature ranges from 4°C (Canadian Prairies and upper midwestern United States) to –12°C in Nunavut. In northern and southern regions of the basin, climate is drier (~200 mm annually) with wetter conditions for the Boreal Forest region (~800 mm annually). Generally, the HBDB is characterized by long, cold winters having significant snowpack accumulation ranging from 100 mm mean annual snow water equivalent (SWE) in the Prairies, to 400 mm SWE in northern Québec (Déry et al., 2005). Snow cover typically begins by early October (mid-November in the south) and stays until mid-June (mid-April in the south) (McKay and Gray, 1981). The HBDB has a primarily nival river regime. Further details on the HBDB watershed and physiography are found in Stadnyk et al. (2019).

3. Methods

3.1 Arctic-HYPE hydrological model

A-HYPE (Andersson et al., 2015; Gelfan et al., 2016) is a regional implementation of the HYPE hydrological model (Lindström et al., 2010), developed by the Swedish Meteorological and Hydrological Institute (SMHI). A-HYPE is set up for hydrological modeling of the pan-arctic landmass that discharges into the Arctic Ocean system, spanning North America and Eurasia. In addition to the base version of HYPE, A-HYPE includes representations of cryospheric processes including glacier accumulation and melt, lake and river ice, and a combined air temperature and radiation index snowmelt module (Gelfan et al. 2017).

3.2 H-HYPE model

The HBDB is set up within A-HYPE as 6,668 subbasins having mean (median) area of 597.8 km2 (351.0 km2), ranging from 2.0 km2 in size, to a maximum of 24,431 km2. Streamflow routing occurs between hydrologically connected subbasins. Within each subbasin, process calculations (e.g., evaporation, snowmelt, infiltration and runoff) are performed for combinations of different landcovers and soil classes. Calculations are performed at a daily time interval, forced by precipitation, mean air temperature, maximum air temperature and minimum air temperature from the HydroGFDv2 gridded reanalysis (Berg et al., 2018). The phase of precipitation depends on a combination of the minimum and maximum daily air temperature: all snowfall below 0°C, all rainfall above 1°C and mixed phase between 0°C and 1°C. Gridded HydroGFDv2 forcing at 0.5° resolution is interpolated to A-HYPE by assigning subbasin forcing to that of the nearest centroid cell. For larger subbasins containing more than one gridded forcing, an area weighted average, or Thiessen polygon-type methodology is used. Eight landcover types are used: crops, forest, open vegetation, bare, open water, glacier, wetland and urban (European Space Agency, 2014). Seven soil types are used (Nachtergaele et al., 2012): coarse, medium, fine, organic, shallow/rock, glacial and urban. Table S1 lists the data sources used in A-HYPE, including inputs for snow and ice (Krenke 2004), glaciers (WGMS, 2012), and lakes (Lehner and Döll 2014).

3.2.1 Regulated reservoirs

There are multiple regulated reservoirs within the HBDB as identified by the Canadian Dam Association criteria (Natural Resources Canada, 2008), documented within Stadnyk et al. (2019). Tefs et al. (in revision) develop a generalized daily regulation routine for HYPE that can be applied to predict specific reservoir releases for 10 of the largest reservoirs operated within the HBDB, targeting monthly operational levels. There remain, however, numerous reservoirs with human-altered regulation that do not have sufficient data to develop stage-discharge relationships for, or which reside in ungauged basins (or where hydrometric data are unavailable).

Using a modification of the existing sine wave relationship coded in HYPE (Lindström et al., 2010), we derive a more suitable and generalized reservoir release relationship for the HBDB in consultation with hydropower utilities having knowledge of the major regulation within the basin (Tefs et al., in revision; Stadnyk et al., 2019). Outflow from regulated reservoirs is predicted using a parameterized specified desired discharge (qprod), which is satisfied if there is sufficient volume in the reservoir. qprod can have a uniform value throughout the year or vary sinusoidally over the year as:


where qprodmean is the mean annual production flow (m3 s–1), qamp is the amplitude of the production flow, dayno is the day of the year, days is the number of days per year (365 or 366) and qpha sets the seasonal phase (minimum sine curve occurs at day number qpha). For reservoir levels above a maximum threshold (w0 [m]), spillway flow is assumed to occur and follows the rating curve:


where k and p are coefficients, and wlm is lake level [m]. Parameters qprod, qprodmean, qamp, qpha, k and p are manually calibrated for a number reservoirs in the HBDB. This is critical for accurate prediction of volume and timing within basins with regulated reservoirs.

3.2.2 Calibration approach

Aside from ensuring adequate process representation, it is imperative that a strong model calibration strategy is used. Recurrent topics in calibrating and validating hydrological models include the choice of calibration algorithm, objective function and split-sample periods (e.g., Refsgaard, 1997, Lindström et al., 1997, and Moriasi et al., 2007). Unaddressed in the literature is the objective selection of gauges for calibration over continental regions, such as the HBDB, where gauge density is irregular and often clustered. The relatively higher density of gauges in the southern HBDB compared to the north presents the risk of over-calibrating the model to fit hydrographs in higher-density regions. We employ a flow signature and cluster analysis within our calibration methodology to avoid overfitting more densely gauged regions, and the under-representation of sparsely gauged (northern) regions (Donnelly et al., 2016).

We perform simulations from 1961–2013, with model spin up from 1961–1970. A split calibration period is used to span both the early (relatively drier and cooler) and late (relatively wetter and warmer) periods; calibration years were 1971–1975, 1981–1985, 1991–1995 and 2001–2005, with skipped years used for validation (Section 3.2.3). Statistics were calculated on a daily basis for day in each year of calibration. Markov chain differential evolution, a genetic algorithm, was used to define uncertainty during auto-calibration (Vrugt, 2009) with 10 populations and 400 generations for each calibration step. This method selects new parameter values as a function of randomness and inherits traits from accepted previous solutions, based on probability analysis. The advantage of this approach over other differential evolutionary methods is that both a global optimum value and probability-based uncertainty estimate of that value are obtained. Calibration is performed using daily hydrometric data from the Water Survey of Canada and US Geological Survey. Of 245 hydrometric gauges, 101 gauges were used to calibrate land surface parameters and 12 are used for calibrating dams and reservoirs (shown on Figure S1). Although the hydrometric gauge on the Sylvia Grinnell River near Iqaluit on Baffin Island is not within the HDBD, its data were used for validation because it is the only gauge located near the northernmost reaches of the basin (i.e., our model domain extended further north of the HBDB extent to include this northern-most subbasin for validation purposes). Calibration steps are outlined in Table 1 and the process flow visually presented on Figure 2 (additional process diagrams are provided in Supplemental material for each operation: sensitivity analysis, guided calibration, and process improvement; Figures S2 to S4, respectively). A subset of unregulated gauges (10 in total) without large headwater lakes were used in aspects of model calibration further described in Section 4 (Table S2).

Figure 2 

H-HYPE calibration, sensitivity and processes improvement methodology. Sub-operation (yellow box) process flow diagrams are provided in Supplemental Figures S2 to S4 for Sensitivity, Guided Calibration and Process Improvement steps. DOI:

The initial parameter set was that calibrated by SMHI for the entire pan-Arctic (Gelfan et al., 2017). The first stage of calibration in this study was a regional calibration of evaporation and runoff adjustment factors (0.80–1.20) to scale the existing A-HYPE landcover and soil type dependent parameters to the regional HBDB context (step 1, Table 1). The same adjustment factor is applied to all landcover and soil types for a given parameter. Evaluation metrics comprised bias, correlation and variance errors, which were equally weighted in the Kling-Gupta Efficiency (KGE; Gupta et al., 2007) score. The Nash-Sutcliffe Efficiency score (NSE; Nash and Sutcliffe, 1970) was also evaluated. Multiple criteria during calibration were weighted according to the error components they most strongly influence, as determined by sensitivity runs (shown in Table 1, last column). Sensitivity runs were conducted separately before each calibration step listed in Table 1, with 10 populations and 200 generations for each, with the objective function an equal weighting of all error criteria (see Figure S2). For each step, the coefficient of determination (R2) between each error criterion and perturbed parameters was calculated. Evaporation parameters most affected simulation bias, runoff parameters most affected NSE and the KGE correlation and variance errors, and lake parameters most affected NSE and the KGE variance error. Details on hydrometric gauge selection for calibration are provided in Section 4.4.

Following regional adjustment factor calibration was the PET crop coefficient calibration that controls evapotranspiration loss (kc), and the fraction of soil available for evaporation (fc), across all landcover and soil types (Step 2, Table 1). Step 3 involved calibration of lake discharge from unregulated basins (Step 3a), and dams and reservoirs for regulated basins (Step 3b), which were performed in parallel. Lastly, routing parameters (damp, rivvel) were optimized in the final stage of calibration to leverage improvements in the estimation of runoff and hydrological storages from all other steps in the calibration process (Step 4, Table 1). Figure S3 provided a process flow of the guided and random calibration conducted for H-HYPE to derive the final set of parameters, and in the stepwise introduction of new processes (and the evaluation of the effectiveness of each process) to the model structure.

We employed k-means clustering to develop the prairie NCA parameterization (Section 4.1), lake clusters (Section 4.3), and for the selection of hydrometric gauges for calibration (Section 4.4). In each case, input data were normalized as standard scores (z-scores) prior to clustering. Two to 10 clusters were calculated for each application. Silhouette widths, which are a measure of how similar a point is to its own cluster compared to other clusters (Rousseeuw, 1987), were calculated for each cluster grouping to objectively select the optimal number of clusters. Clusters with highest average silhouette width and fewest negative scores were selected (Tables 2 and 3, Table S3).

Table 2

Description of five selected flow signature clusters for HBDB. DOI:

Cluster Description Number of Gauges Predominant Location

1 Moderate 86 Entire domain
2 Greatest flashiness 8 Western mountainous headwaters
3 Greatest low flow ratio; lowest high flow ratio 63 Entire domain except Nelson basin
4 Greatest flow magnitude and specific discharge; greatest variance of peak values 38 Shield region
5 Greatest high flow ratio and variance 50 Prairies

Table 3

Description of selected olake and ilake clusters. DOI:

Cluster Description Number of Subbasins Predominant Location

Olake clusters
1 Moderate relief 348 Entire domain
2 Greatest lake surface area 14 Entire domain
3 High relief 140 Canadian Shield
4 Low relief 326 Entire domain
Ilake clusters
1 Median lake count = 150 944 Low Arctic, Canadian Shield and Hudson Plains
2 Median lake count = 17 4479 Boreal Plains and Prairies
3 Median lake count = 407 150 Low Arctic and Canadian Shield

Lastly, flow signatures were used to capture specific (sometimes regional) similarities from measured hydrometric data to assist with calibration gauge selection. Eight flow signatures were calculated using daily observed discharge based on the method from Pechlivanidis and Arheimer (2015; Figure S5; Table S4) and were applied to select distinct groupings of gauges for calibration (Section 4.4, Table 2). Gauges used for calibration, relative to flow signatures and NCA sites, are shown on Figure S1.

3.2.3 Validation

After each auto-calibration step, a split sample temporal validation of the optimized parameter set was performed before starting the next auto-calibration step. This resulted in five validation steps using the same objective functions as derived for calibration. Validation years were 1976–1980, 1986–1990, 1996–2000 and 2006–2013.

4. Model enhancements and results

A-HYPE model enhancements were incrementally implemented as part of the development of the H-HYPE model, which are further described below. Model performance from calibration is presented sequentially, for each enhancement. Figures with boxplots show the median, 25th and 75th percentiles at the hinges, 95% confidence interval at the whisker extents, and outliers beyond.

4.1 Prairie non-contributing area

The new parameterization strategy leverages the existing NCA map from the PFRA (1983), which includes anthropogenically drained or altered landscapes. The NCA (shown on Figure 1) is that for a return period of two years; therefore, it is likely that the PFRA NCA is underestimated for low runoff conditions, and potentially overestimated for the highest runoff conditions. An increase in the NCA decreases the effective drainage area of a subbasin.

We identify thresholds in simulated daily runoff at which different NCAs (or effective drainage areas) are applied. The analysis uses A-HYPE (prior to H-HYPE calibration) evaluated at 10 gauges with upstream NCA but without large upstream lakes or regulation (Table S2, Figure S1). Although the Water Survey of Canada database ( indicates that these are natural rivers, it cannot be confirmed that their drainage areas are completely void of agriculture withdrawals or altered drainage. Although there is considerable variability in the relationship between simulated runoff and observed specific discharge for these 10 gauges (Figure 3), a Pearson’s correlation coefficient of 0.87 (p < 0.001) indicates that, even when neglecting the influence of complex pothole connectivity dynamics (contributing to the scatter, or variability), runoff is a suitable for parameterizing NCA within A-HYPE.

Figure 3 

A-HYPE simulated daily runoff (original parameter set) relative to measured daily specific discharge for 10 gauges used in NCA parameterization development. A black line represents the linear regression through the simulation data with associated Pearson’s r correlation. Scatter indicates the complex pothole runoff generation dynamics of the region. DOI:

The relationship between simulated runoff and the error in simulated specific discharge is examined using differences between the logarithms of simulated and measured discharge to give comparable weighting to high and low flows (Figure 4). K-means clustering (Section 3.2) is performed to identify runoff thresholds at which different NCA schemes may be applied. Two runoff thresholds (three clusters) are selected: <0.351 mm day–1, 0.351 to 1.071 mm day–1, >1.071 mm day–1; thresholds are taken as the minimum values of runoff clusters 2 and 3.

Figure 4 

Simulated daily runoff and discharge error (a function of Qsim (simulated) and Qobs (measured) discharge) clusters. Data for all 10 gauges used in NCA parameterization development are shown. DOI:

Sensitivity runs were performed, using KGE as a performance metric, to find an NCA scheme for each of the three runoff threshold groupings. Each runoff threshold has a multiplier, αi, that is applied to the published NCA in subsequent simulations:

NCAnew={α1NCApub,  runoff<0.351 mm/dayα2NCApub,  0.351 mm/dayrunoff<1.071 mm/dayα3NCApub,  runoff1.071 mm/day

where NCAnew is the updated NCA and NCApub is the published PFRA NCA. Best fit multipliers are found to be: α1 = 1.15, α2 = 0.65 and α3 = 0.25. The new NCA parameterization reduces subbasin drainage area, and therefore runoff values within each timestep. It also corrects the underestimation of PFRA NCA and effective drainage for low runoff, and overestimation for higher runoff. The new runoff threshold based NCA parameterization improves model performance across the 10 gauges (Figure 5). NCA processes and parameterization were introduced into the H-HYPE model, and the model recalibrated (Figure S4). Mean deviation of runoff volumes (Dv) improve from –13% to –8%, and median KGE improves from 0.21 to 0.28, relative to a simple NCA mask. Results confirm that NCA modifications impact median to low flow more significantly than peak flow, evidenced by similarity among NSE scores.

Figure 5 

Improvement in model error statistics by adding NCA parameterizations to the A-HYPE model: (a) Dv, (b) KGE, (c) NSE. Daily performance statistics for 10 gauges used in NCA parameterization development are shown for years 1971–2013. DOI:

4.2 Lake clusters

Lakes are represented in HYPE as either outlet lakes (olakes) or internal lakes (ilakes). Olakes are located at the outlet of subbasins. Ilakes are located conceptually within a subbasin and receive a portion of subbasin runoff. Olakes, in particular, exert a strong control on discharge characteristics. Outflow from olakes and ilakes in HYPE are parameterized using rating curves:

outflow=k(wlmlake depth)p

where k and p are coefficients, and wlm is lake level (m); lake depth is taken from Kourzeneva (2010). In H-HYPE, 84% of the subbasins contain ilakes and 12% contain olakes. Few of these lakes are gauged directly at the mouth, therefore, it is necessary to generalize lake outflow rating curve parameters in some way.

Groupings of ilakes and olakes with similar physiographic characteristics were made using k-means cluster analyses (Section 3.3.2). Lakes with similar physiography have similar discharge characteristics and can be given common calibrated rating curve parameters (k and p). Combinations of a number of physiographic parameters are considered; for olakes: olake surface area, slope (10 km buffer around olake), relief (10 km buffer around olake), elevation (10 km buffer around olake) and upstream area; for ilakes: ilake area, ilake count per subbasin, mean ilake area, standard deviation of ilake area, subbasin slope, subbasin relief.

It is best to avoid a large number of lake clusters to limit equifinality during calibration (Beven, 2001). Table 1 identifies three parameters × number of ilake clusters and two parameters × number of olake clusters that are used in the model setup and calibration. Silhouette analyses (Table S3) indicate that the best selections of lake clusters were four olake clusters based olake area and relief, and three ilake clusters based on ilake count alone (Figure 6 and Table 3). Subbasins with high ilake count and relief were most common in the Canadian Shield and low Arctic (Figure S6), as would be expected resulting from localized internal (to the subbasin) storage in isolated and perched wetlands and lakes within the Canadian Shield and poorly draining sub-Arctic terrain. Cluster regions, assigned at the subbasin level, determine how ilakes and olakes are parameterized within a given region (Figure S6a, b). Selection of the optimal number of lake clusters was made by finding the number of clusters which returned the smallest number of negative silhouette scores and the largest average silhouette width by ranking these data and taking the highest combined rank (Table S3), with distributions of characteristics for the optimal selection shown in Figure 6.

Figure 6 

Characteristics of (a) olake clusters, and (b) ilake clusters. DOI:

Model performance results from two alternate lake cluster selections were compared to those based on the above cluster formations for validation of the methodology: Alternative 1) four olake clusters based on olake area and three ilake clusters based on ilake area; and Alternative 2) four olake clusters based on olake area, slope and elevation and three ilake clusters based on ilake area, standard deviation of ilake area, elevation and relief. Lake parameters were recalculated for the alternate lake clusters (i.e., parameters from calibration steps 1 and 2 (Table 1) are retained, but step 3a is repeated separately). Results from the calibration of olake and ilake cluster parameters are presented in Section 4.5.1.

4.3 Frozen soils

Four representations of frozen soils in A-HYPE are considered: 1) none, 2) frozen soils limiting infiltration, 3) frozen soils limiting subsurface runoff, and 4) frozen soils limiting both subsurface runoff and infiltration. Soil layer temperatures are calculated following the equations presented in Pers et al. (2016). Prior to this study, there was no representation of frozen soils in A-HYPE (Option 1).

Option 2 employs the parametric expression for estimating infiltration into frozen soils from Zhao and Gray (1999). The relationship relates infiltration to total soil saturation and temperature at the beginning of snowmelt, the soil surface saturation during melt and the infiltration opportunity time estimated following Zhao and Gray (1997). Infiltration calculations are grouped into three categories:

Restricted – infiltration is completely restricted due to impermeable surface conditions such as ice lens formation. This is identified by the soil temperature falling below a threshold for lens formation (default: –10°C) and sufficient recent snowmelt (default: 5 mm).

Limited – capillary flow predominates and infiltration is primarily controlled by soil physical properties.

Unlimited – gravity flow predominates and all water infiltrates.

The parametric equation is used for the limited infiltration case:

INF=CS02.92(1SI)1.64(273.15TI273.15)0.45t00.44    TI273.15K

where INF is the potential infiltration capacity (mm), C is an empirical constant with value taken from Gray et al. (2001), S0 is the soil surface saturation (assumed equal to 1.00 during infiltration), SI = θI0 is the premelt pore saturation of the upper soil layer with θI being the volumetric soil moisture at the start of infiltration and θ0 being the porosity, TI is the pre-melt temperature of the upper soil layer (K), and t0 is the infiltration opportunity time (hours).

Option 3 employs the freezing point depression equation (e.g., Fuchs et al., 1978; Flerchinger and Saxon, 1989). As soil temperatures fall below freezing, water nearest soil particles remains in liquid form due to absorptive and capillary forces. The soil matric potential has some dependence on soil texture, expressed by the empirically derived Clapp-Hornberger parameters (Cosby et al., 1984). Combining the equation relating water suction to temperature when ice is present (e.g., Black and Tice, 1988) with the Clapp and Hornberger (1978) equation for suction as a function of liquid water content gives the upper limit of liquid water content for subfreezing temperatures:


where θl,max is the maximum volumetric liquid water content (m3 m–3), θS is the saturation volumetric water content or porosity (m3 m–3), LF is the latent heat of fusion (3.34 × 105 J kg–1), Tfreeze = 273.15 K, Tsoil is soil layer temperature (K), g is gravitational acceleration (9.81 m s–2), ψs is soil water potential at saturation (m), and b is the shape coefficient of the soil water potential-moisture curve. ψs and b are empirically derived soil texture dependent parameters with values taken from Cosby et al. (1984). Equation 6 is applied to each soil layer in A-HYPE, with only θl,max being available for runoff.

The model is calibrated separately for each of the four frozen soil representations using the same methodology described in Section 3.2 (shown on Figure S4). Five populations of 300 generations were used, with equally weighted mean KGE and mean NSE as the objective function. Only soil runoff parameters are calibrated (srrcs_corr, rrcs_corr, ep_corr, wp_corr and fc from Table 1). Performance statistics for the best parameter set from each calibration chain were combined for evaluation and all gauges in the HBDB are assessed (refer to Section 4.5.1).

4.4 Flow signatures and calibration gauge selection

For a large watershed such as the HBDB, the hydrometric gauges selected for calibration should be non-biased with respect to drainage area covered by distinct flow signatures. K-means cluster analyses (Section 3.3.2) are used to group hydrometric gauges based on similar flow signatures (Section 3.3.1, Figure S5). A grouping of five flow signature clusters was chosen to maximize mean silhouette width with the fewest negative scores. Resulting flow signatures follow expected regional streamflow patterning: greatest flashiness in higher sloping mountainous regions, greatest normalized high flow ratios due to NCAs across the prairies, and greatest flow magnitudes in wettest regions with the largest drainage areas (Table S2).

An appropriate number of area-weighted stations from each flow signature cluster are selected for calibration (Figure S1). The particular calibration stations are manually selected with the goal of minimizing overlapping upstream drainage area, maximizing spatial coverage, and maximizing the range of years covered. In practice, selection of stations for calibration is limited in most cases to those available in, and required from each, flow signature cluster. The ilake and olake clusters are well-covered by the calibration gauges: olake cluster 1 (45.2%), olake cluster 2 (35.7%), olake cluster 3 (47.9%), olake cluster 4 (20.6%), ilake cluster 1 (40.5%), ilake cluster 2 (37.5%), ilake cluster 3 (32.5%).

To verify the flow signature approach to gauge selection, model calibration is repeated for four different random selections of 101 gauges. Calibration steps 1–3a (Table 1) were repeated; model performance results are compared in section 4.5.

4.5 Calibration results

4.5.1 Frozen soils, lake clusters and flow signatures

The ilake and olake clusters were selected using the silhouette analyses to optimize model performance (NSE and KGE scores) among the alternate lake clusters (Figure S7 only NSE is shown). Median NSE and KGE are 0.04–0.07 and 0.04–0.08 better than the alternate clusters. There are negligible differences in runoff volume (Dv) scores between the different clusters because bias is more so controlled by the evaporation and runoff parameters.

For frozen soil representation, option 3 (only frozen subsurface) provides the best performance across all evaluation metrics, i.e. NSE, KGE and Dv (each of mean, median, 1st quartile and 3rd quartile of statistics; Figure 7 – only KGE shown). When compared to no frozen soils, option 3 provides slight improvements in mean evaluation scores (2.9% Dv, 0.04 KGE and 0.06 NSE improvements). Option 3 was retained for further modeling and calibration. Including frozen soil effects on infiltration (option 2) degrades model performance across all measures (18.5% Dv, 0.14 KGE and 0.23 NSE degradations compared to no frozen soils).

Figure 7 

Comparison of KGE scores for each representation of frozen soils across 245 hydrometric gauges in the HBDB for years 1971–2013. DOI:

The flow signatures approach to gauge selection for calibration provides better performance largely in validation (i.e., long term simulation) than the random selection of gauges across all measures of model performance. Median NSE, KGE and Dv are 0.03–0.11, 0.04–0.16 and 3.4%–8.0% better than the random gauge selection for validation. Lower quartiles improved by 0.08–0.17, 0.03–0.20, and 2.2%–7.4% during validation when using the flow signatures.

4.5.2 Final model performance

Calibration improves all performance measures, with similar scores for calibration and validation indicating a verifiable parameter set. Calibration improves the variance of performance across all gauges (i.e., a lower difference between the third and first quartiles). There is substantial improvement in simulated volume with median Dv error decreasing from 27.4% to 2.7% and variance reduced (all gauges, both calibration and validation time periods; Figure not shown). This is an indication of general improvement of estimation error at all gauges across the HBDB domain. There is also increased skill in the timing and variance of discharge, with median NSE improving from 0.26 to 0.43 and median KGE improving more substantially from 0.39 to 0.59 (Figure 8). This result indicates a reduction of the long-term error estimation was more substantial than the estimation of individual peak flows among gauges.

Figure 8 

Final model performance for KGE error across all hydrometric gauges associated with each model calibration step. DOI:

Calibration of evaporation and runoff parameters produces a marked improvement in all measures (Figure 8; only KGE shown), particularly in reducing volume errors. There is negligible improvement in volume error from subsequent calibration steps. Lake calibration improves the simulated variance of discharge, as reflected in improved KGE scores. Lake calibration also degraded some NSE scores, but higher scores are obtained following dam and routing calibration. Routing calibration has little effect on model error statistics, indicating the strong control of lake discharge parameters on routing performance.

Spatially, calibration improves model performance across the whole HBDB. Volume error and NSE improve across all regions of the HBDB with few degradations, with comparatively smaller improvements in NSE than KGE (KGE not shown). Some regional distinctions in statistical improvement are apparent (Figure 9). Spatial plots indicate lower Dv values in northern rivers, in Québec, and parts of the southern prairies. Conversely, there are gauges in the southwest HBDB where volume errors are high (Figure 9b). Figure 10 presents the final KGE error across the HBDB domain, used as a final assessment of model performance for the long-term predictive capability to estimate mean and variance in freshwater discharge to Hudson Bay.

Figure 9 

Spatial plots of the change in model performance statistics pre- and post-calibration: (a) Deviation of runoff volume (original parameters), (b) Deviation of runoff volume (calibrated parameters), (c) Nash-Sutcliffe Efficiency (original parameter), (d) Nash-Sutcliffe Efficiency (calibrated parameters). DOI:

Figure 10 

Spatial analysis of model performance based on the KGE metric. DOI:

5. Discussion

The results presented in Section 4.5.2 demonstrate improvement in model performance from the process-based enhancements added to the HYPE model for the HBDB. Including an NCA representation markedly reduces streamflow overestimation across the prairies (Figure 5). The developed parameterization does not, however, simulate prairie pothole fill and spill dynamics, which are the primary runoff generation control in the prairie environment (Stichling and Blackwell, 1957; Shook and Pomeroy, 2011; Shook et al., 2013). It is expected that the simulation of large-scale prairie pothole dynamics will aid in simulating the timing of discharge and peak flow performance (NSE); whereas the parameterization developed here is primarily a volume-based (storage) correction. There have been several algorithms developed over the past decade to specifically tackle runoff generation in the Prairie Pothole Region, including probability distribution models (Mekonnen et al., 2014; 2016), fill and spill DEM or depression-driven methods (Shook et al., 2013; Nasab et al., 2017), and hybrid approaches combining physically-based and stochastic (artificial neural networks and physically based models) (Mekonnen et al., 2015). But a fully realistic representation of the complexity and dynamic nature of the prairie potholes in large-scale hydrological modeling experiments still requires the use of a routing product derived from high-resolution digital terrain models, and exponentially increases computational resource requirements (Muhammad et al., 2019), which is not practical or possible for the vast HBDB. The NCA representation we develop in this study balances improved model accuracy (for low flow approximation) with efficiency and realistic input data requirements.

Model improvements realized from the lake clusters were demonstrated during calibration (Figure 8). Combined calibration of the olake and ilake cluster parameters improves the efficacy of H-HYPE to correctly simulate the variance of observed discharge, with the variance error component of KGE reflecting this result. Calibration improves the mean and median KGE variance error components to 0.96 (from 0.70) and 0.98 (from 0.77), respectively, which approach the optimal value of 1.00. This results from improvements to the amount of land-based storage, or simulated volume of runoff within the model (ilakes), resulting in a more constrained simulated discharge value downstream of olakes.

With global warming frozen soils have become an indicator of change in the natural ecohydrological system. The implementation of a frozen soil scheme in cold region land surface models or spatially distributed hydrological models is required for more realistic long-term projections of runoff impacted by changing soil moisture and soil temperature profiles (Wang et al., 2010; Qin et al., 2016). Frozen soils are distinct because of unique hydrological characteristics relative to thawed soils (e.g. lower hydraulic conductivity, higher soil water content, seasonal freeze/thaw dynamics, and redistribution of frozen/thawed water) (Cherkauer and Lettenmaier, 1999; Qin et al., 2016). In this study, the freezing point depression representation of subsurface runoff provides a small improvement in discharge simulations compared to no representation of frozen soils (Figure 7). These improvements are generally greater at higher latitudes, where H-HYPE performs poorest. There is otherwise not a clear spatial pattern to model improvement. The Zhao and Gray (1999) parameterization tested for reducing infiltration into frozen soils degrades model performance across the HBDB in this study as it produced earlier and flashier runoff, which was not commensurate with observed discharge records. It is noted, however, that regulated discharge signals (downstream of reservoirs) produce flattened hydrograph responses that do not mimic naturally flashy hydrographs typical of frozen soils (i.e., particularly in the Prairie basins). In these regions, where frozen soils and regulated reservoirs co-exist, there is likely error in the modeling process resulting from both the predicted (anthropogenic) reservoir releases and (natural) frozen soil infiltration that cannot be distinguished without more data to validate frozen soils. Pitman et al. (1999) showed that runoff simulations at large scales can be poorer when suppressing infiltration into frozen soils and suggest this is due to the presence of preferential drainage pathways such as macropores.

Our calibration methodology improves model performance, but there are identified areas for further improvement. The H-HYPE model set-up has eight landcover types and seven soil types, with common landcover and soil parameters applied to various areas across the whole HBDB. We suspect that further discretization of soil and land classes, and a more regional approach to model calibration, could improve model performance. NSE scores are lowest (<0) in the lake-heavy Churchill River Basin where a high fraction of ilake cluster 3 exists. Further improvements to lake storage-discharge curve relationships is also likely required in this region to improve model performance. Representations of more than a single ilake and olake per subbasin may be warranted.

Model performance metrics indicate the H-HYPE model performs well at estimating long-term mean and variability in streamflow (i.e., KGE scores), but still struggles to estimate the volume, timing (or both) of peak flow, evidenced by smaller improvements to the NSE scores. This is particularly true, and most prevalent (i.e., region with lowest NSE scores) across the Prairie region of the Western Hudson Bay where hydrographs are characteristically flashy. This region is notoriously difficult to simulate using hydrologic models and has been referred to as the “graveyard of hydrologic models” (Muhammad et al., 2019). Moriasi et al. (2015) suggest a hydrologic model is satisfactorily calibrated when NSE>0.45 at the watershed scale. The HBDB is an extremely large and hydrologically diverse (and complex) region, and yet we are able to achieve this standard for more than half the subbasins in the H-HYPE domain. Moreover, given the emphasis on long-term and large-scale mean and variability in freshwater discharge to Hudson Bay, KGE scores better reflect the desired performance of our model. Based on KGE performance, we determine the H-HYPE model to be satisfactorily fit-for-purpose for estimating the freshwater discharge to Hudson Bay for the BaySys project.

6. Conclusions

This study summarizes the initial development of the H-HYPE hydrological model, developed to estimate freshwater discharge across the HBDB terrestrial domain and into Hudson Bay; a largely ungauged continental scale region. This modeling was undertaken as part of the larger BaySys project to provide non-stationary, dynamic time-series of freshwater flux for the NEMO ocean model (Ridenour et al., 2019), which will then be used to generate scenarios of sea ice, biogeochemical, and ecosystem processes. A key aspect of the BaySys project is to separate the relative contributions of human and climate-driven impacts to the freshwater discharge to Hudson Bay. To that end, several key developments were needed for the model to be capable of simulating and discerning climate-driven hydrology from those imposed by anthropogenic control, or water resources management. Our publication documents those developments, and a stepwise approach to calibrating the baseline freshwater discharge scenarios for the BaySys project.

We made three landscape scale storage-based improvements to the original A-HYPE model to facilitate improved separation of natural from human-driven changes to the water cycle, including the 1) representation of prairie pothole NCAs, 2) parameterization of lakes and regulated storage-discharge relationships, and 3) representation of frozen soils in the model structure.

The NCA parameterization provided marked reduction in simulated discharge bias across the Prairie region. As it is principally a volume-based correction without representing the dynamic fill and spill runoff generation process, it is recognized that a more refined understanding of large-scale prairie runoff generation mechanisms will further improve hydrological simulations in this region. Our approach to identifying NCA thresholds is applicable in other Prairie region hydrological studies, and the identified thresholds transferable to other studies within this region.

The use of lake clusters improved the ability of H-HYPE to simulate the variability of observed discharge. The lake clustering technique, based on physiographic characteristics, allows the transfer of calibrated rating curve parameters to ungauged regions. This approach is generalizable and can be applied in other, regional or continental scale ungauged basins with numerous lakes (e.g., Viglione et al., 2013). By identifying lake similarity characteristics (e.g., Fig. 6) that dominate the lake runoff response, and calibrating that response among gauged basins or lakes in each cluster, the identified parameters can then be transferred to ungauged regions identified as having similar hydrologic behaviour. Further, more refined (region-specific) calibration of the most lake-heavy subbasins in the HBDB, such as the western Canadian Shield, are needed to further improve model performance. Since the goal of this study was to derive a Hudson Bay-wide model and calibration that resulted in a good match of total freshwater flux to the Bay, region-specific calibration was out of scope for this project, and there were insufficient data to support such an endeavour. Though the timing and volume of freshwater flux delivered to the Bay are important for local-scale circulation, the scope of the BaySys project (given the size of the domain and scale of analyses) was to analyse variability on a monthly basis, including the separation of regulation and climate-driven influences (Barber, 2014). The nature of how H-HYPE freshwater discharge is integrated into the NEMO ocean model also spatially aggregates discharge on a regional basis, meaning that discharge from several adjacent rivers is combined into one ocean polygon. Given these spatial and temporal constraints imposed on BaySys analyses, we find the H-HYPE model to be appropriate for simulating the dynamics of the Hudson Bay freshwater discharge response. A recent study by Ridenour et al. (2019) compared the H-HYPE freshwater discharge inputs, and their impact on ocean variables, to reanalysis discharge datasets and found comparable performance.

Reducing subsurface runoff due to soil freezing provided only minor improvements in simulated discharge; however, limiting infiltration into frozen soils was found to degrade discharge simulations. We suggest that the parameterization of infiltration into frozen soils at large scales should remain an open area of research, and the algorithms identified in our study tested in other cold regions at varying spatial scales. An emerging critical area of research is on the parameterization of deep, frozen soils (i.e., permafrost). Several studies have been published looking at the effects of frozen soils and permafrost under changing climates (e.g., Evans et al., 2018; Connon et al., 2014; Smith et al., 2010), but further research into the impact on surface runoff, and for methods of regional data collection (i.e., remote sensing) are warranted and could potentially improve discharge simulations in the Prairie and most northern regions of the HBDB, where the model performs most poorly. Unfortunately, at the time of this study, there were a lack of frozen soil and permafrost data to further validate such modeling efforts.

The use of flow signatures for an unbiased calibration gauge selection may also be applied to other large, sparsely gauged basins. This approach avoids overfitting parameter sets to regions where gauge density is greatest (e.g., the southern HBDB) in gauges where there is uneven gauge coverage. There are significant regions of the HBDB that remain ungauged, such as Foxe Basin (furthest north), which remain problematic and of concern. In the HBDB, there was a significant decrease in the number of gauged basins with increasing latitude (Coulibaly et al., 2013), explaining the poorer model performance in the north. It is in these areas where we rely on the inherent value of semi-distributed hydrological models, where calibrated parameters are transferred to ungauged areas with the same landcover and soil classifications (Dornes et al., 2008). Of note is the non-independence of certain flow signatures, which could introduce bias in the flow signature clustering. A future application of this method may consider self-organizing maps in place of flow signatures (e.g., Ley et al., 2011). Challenges remain in improving model performance over the northern and north-eastern parts of the basin, where a regionalized approach to model calibration is needed and an increased discretization of landcover and soil classes.

Data Accessibility Statement

All simulation output used in this research is available on Compute Canada, and made available through the BaySys data repository (open access).

Supplemental file

The supplemental file for this article can be found as follows:


Financial support for this research was provided by Manitoba Hydro and the Natural Sciences and Engineering Research Council of Canada through the BaySys project. Kristina Koenig, John Crawford, Mark Gervais, Phil Slota, Mike Vieira and Shane Wruth (Manitoba Hydro) provided helpful advice throughout this work and beneficial reviews prior to submission. Mark Gervais provided shapefiles for Figure 1. Charlotta Pers and Jörgen Rossberg (Swedish Meteorological and Hydrological Institute) assisted with HYPE code understanding and drainage basin set up, respectively. We thank those along the way, including our reviewers, who have contributed constructive comments to improve this manuscript.

Funding information

This research was funded by the Natural Sciences and Engineering Research Council (NSERC) of Canada, a national granting agency. Research was conducted as part of the BaySys collaborative research and development NSERC grant (D. Barber), which was jointly supported by Manitoba Hydro (cash contribution) and the following agencies: ArcticNet, Ouranos, Hydro Quebec, and various Canadian Universities (in-kind).

Competing interests

The authors have no competing interests to declare.

Author contributions

  • Contributed to conception and design: TAS, MM, SJD
  • Contributed to acquisition of data: SJD, MM, DG, KI, BA
  • Contributed to analysis and interpretation of data: MM, TAS, SJD, HA, AT
  • Drafted and/or revised the article: MM, HA, AT, SJD, TAS
  • Approved the submitted version for publication: TAS, SJD


  1. Andersson, JCM, Pechlivanidis, IG, Gustafsson, D, Donnelly, C and Arheimer, B. 2015. Key factors for improving large-scale hydrological model performance. European Water 49: 77–88.

  2. Barber, D. 2014. BaySys – Contributions of climate change and hydroelectric regulation to the variability and change of freshwater-marine coupling in the Hudson Bay System. Retrieved from [Accessed 2017-11-05].

  3. Berg, P, Donnelly, C and Gustafsson, D. 2018. Near real-time adjusted reanalysis forcing for hydrology. Hydrol. Earth Sys. Sci 22: 989–1000. DOI: 10.5194/hess-22-989-2018

  4. Beven, K and Freer, J. 2001. Equifinality, data assimilation, and uncertainty estimation in mechanistic modeling of complex environmental systems using the GLUE methodology. J. Hydrol 249: 11–29. DOI: 10.1016/S0022-1694(01)00421-8

  5. Black, PD and Tice, RA. 1988. Comparison of soil freezing curve and soil water curve data for windsor sandy loam. US Army Cold Regions Research and Engineering Report , 88-16. 42p. DOI: 10.21236/ADA202365

  6. Canadian Council of Ecological Areas. 2004. Ecozones of Canada. Accessed 11 Apr 2015.

  7. Cherkauer, KA and Lettenmaier, DP. 1999. Hydrologic effects of frozen soils in the upper Mississippi River basin. J. Geophys. Res 104(D16): 19599–19610. DOI: 10.1029/1999JD900337

  8. Cherkauer, KA and Lettenmaier, DP. 2003. Simulation of spatial variability in snow and frozen soil. J. Geophys. Res 101(D22): 8858. DOI: 10.1029/2003JD003575

  9. Connon, RF, Quinton, WL, Craig, JR and Hayashi, M. 2014. Changing hydrologic connectivity due to permafrost thaw in the lower Liard River valley, NWT, Canada. Hydrol. Process 28: 4163–4178. DOI: 10.1002/hyp.10206

  10. Cosby, BJ, Hornberger, GM, Clapp, RB and Ginn, TR. 1984. A statistical exploration of the relationships of soil moisture characteristics to the physical properties of soils. Wat. Resour. Res 20: 682–690. DOI: 10.1029/WR020i006p00682

  11. Coulibaly, P, Samual, J, Pietroniro, A and Harvey, D. 2013. Evaluation of Canadian National Hydrometric Network density based on WMO 2008 standards. Can. Wat. Resour. J 38: 159–167. DOI: 10.1080/07011784.2013.787181

  12. Déry, SJ, Stadnyk, TA, MacDonald, MK and Gauli-Sharma, B. 2016. Recent trends and variability in river discharge across northern Canada. Hydrol. Earth Sys. Sci 20: 4801–4818. DOI: 10.5194/hess-20-4801-2016

  13. Déry, SJ, Stadnyk, TA, MacDonald, MK, Koenig, K and Guay, C. 2018. Flow alteration impacts on Hudson Bay river discharge. Hydrol. Process. 32(24): 3576–3587. DOI: 10.1002/hyp.13285

  14. Déry, SJ, Stieglitz, M, McKenna, EC and Wood, EF. 2005. Characteristics and trends in river discharge into Hudson, James and Ungava Bays, 1964–2000. J. Climate 18: 2540–2557. DOI: 10.1175/JCLI3440.1

  15. Donnelly, C, Andersson, JCM and Arheimer, B. 2016. Using flow signatures and catchment similarities to evaluate the E-HYPE multi-basin model across Europe. Hydrol. Sci. J 61: 255–273. DOI: 10.1080/02626667.2015.1027710

  16. Dornes, PF, Tolson, BA, Davison, B, Pietroniro, A, Pomeroy, JW and Marsh, P. 2008. Regionalisation of land surface hydrological model parameters in subarctic and arctic environments. Phys. Chem. Earth 33: 1081–1089. DOI: 10.1016/j.pce.2008.07.007

  17. European Space Agency. 2014. Land Cover CCI Product user guide version 2. Accessed 26 June 2017.

  18. Evans, SG, Ge, S, Voss, CI and Molotch, NP. 2018. The role of frozen soil in groundwater discharge predictions for warming alpine watersheds. Wat. Resour. Res. 54: 1599–1615. DOI: 10.1002/2017WR022098

  19. Evenson, GR, Golden, HE, Lane, CR and D’Amico, E. 2016. An improved representation of geographically isolated wetlands in a watershed-scale hydrologic model. Hydrol. Process 30: 4168–4184. DOI: 10.1002/hyp.10930

  20. Flerchinger, GN and Saxton, KE. 1989. Simultaneous heat and water model of a freezing snow-residue-soil system I. Theory and development. Trans. ASAE 32: 0565–0571. DOI: 10.13031/2013.31040

  21. Fortin, J-P, Turcotte, R, Massicotte, S, Moussa, R, Fitzback, J and Villeneuve, J-P. 2001. Distributed watershed model compatible with remote sensing and GIS data. I: Description of model. J. Hydrol. Eng 6: 91–99. DOI: 10.1061/(ASCE)1084-0699(2001)6:2(91)

  22. Fuchs, M, Campbell, GS and Papendick, RI. 1978. An analysis of sensible and latent heat flow in a partially frozen unsaturated soil. Soil Sci. Soc. Am 42: 379–385. DOI: 10.2136/sssaj1978.03615995004200030001x

  23. Gelfan, A, Gustafsson, D, Motovilov, Y, Arheimer, B, Kalugin, A, Krylenko, I and Lavrenov, A. 2017. Climate change impact on the water regime of two great Arctic rivers: modeling and uncertainty issues. Clim. Change 141: 499–515. DOI: 10.1007/s10584-016-1710-5

  24. Gray, DM, Toth, B, Zhao, L, Pomeroy, JW and Granger, RJ. 2001. Estimating areal snowmelt infiltration into frozen soils. Hydrol. Process 15: 3095–3111. DOI: 10.1002/hyp.320

  25. Gupta, HV, Kling, H, Yilmaz, KK and Martinez, GF. 2009. Decomposition of the mean squared error and NSE performance criteria: implications for improving hydrological modeling, J. Hydrol 377: 80–91. DOI: 10.1016/j.jhydrol.2009.08.003

  26. Gupta, HV, Sorooshian, S and Yapo, PO. 1999. Status of automatic calibration for hydrologic models: Comparison with multilevel expert calibration. J. Hydrol. Eng 4: 135–143. DOI: 10.1061/(ASCE)1084-0699(1999)4:2(135)

  27. Hayashi, M. 2013. The cold vadose zone: hydrological and ecological significance of frozen-soil processes. Vadose Zone J. 12. DOI: 10.2136/vzj2013.03.0064

  28. Huang, SL, Young, C, Abdul-Aziz, OI, Dahal, D, Feng, M and Liu, SG. 2013. Simulating the water budget of a Prairie Potholes complex from LIDAR and hydrological models in North Dakota, USA. Hydrol. Sci. J 58: 1434–1444. DOI: 10.1080/02626667.2013.831419

  29. Kourzeneva, E. 2010. External data for lake parameterization in Numerical Weather Prediction and climate modeling. Boreal Env. Res 15: 165–177. DOI: 10.4296/cwrj1301062

  30. Kouwen, N. 1988. WATFLOOD: A Micro-Computer based Flood Forecasting System based on Real-Time Weather Radar. Can. Wat Resour. J 13: 62–77.

  31. Kouwen, N. 2016. WATFLOOD/CHARM Canadian Hydrological and Routing Model. Accessed 26 June 2017.

  32. Krenke, A. 2004. Edited by National Snow and Ice Data Center. 1998, updated. Former Soviet Union Hydrological Snow Surveys, 1966–1996. Boulder: National Snow and Ice Data Center. Accessed 26 June 2017.

  33. Lehner, B and Döll, P. 2004. Development and validation of a global database of lakes, reservoirs and wetlands. J. Hydrol 296: 1–22. DOI: 10.1016/j.jhydrol.2004.03.028

  34. Lehner, B, Reidy Liermann, C, Revenga, C, Vörösmarty, C, Fekete, B, Crouzet, P, Döll, P, Endejan, M, Frenken, K, Magome, J, Nilsson, C, Robertson, JC, Rodel, R, Sindorf, N and Wisser, D. 2011. Global Reservoir and Dam Database, Version 1 (GRanDv1): Dams, Revision 01. Palisades, NY: NASA Socioeconomic Data and Applications Center (SEDAC). DOI: 10.7927/H4N877QK

  35. Ley, R, Casper, MC, Hellebrand, H and Merz, R. 2011. Catchment classification by runoff behavior with self-organizing maps (SOM). Hydrol. Earth Sys. Sci 15: 2947–2962. DOI: 10.5194/hess-15-2947-2011

  36. Lindström, G. 2016. Lake water level for calibration of the S-HYPE model. Hydrol. Res 47: 672–682. DOI: 10.2166/nh.2016.019

  37. Lindström, G, Johansson, B, Persson, M, Gardelin, M and Bergström, S. 1997. Development and test of the distributed HBV-96 hydrological model. J. Hydrol 201: 272–288. DOI: 10.1016/S0022-1694(97)00041-3

  38. Lindström, G, Pers, C, Rosberg, J, Strömqvist, J and Arheimer, B. 2010. Development and testing of the HYPE (Hydrological Predictions for the Environment) water quality model for different spatial scales. Hydrol. Res 41: 295–319. DOI: 10.2166/nh.2010.007

  39. Manak, DK and Mysak, LA. 1989. On the relationship between Arctic sea-ice anomalies and fluctuations in Northern Canadian air temperature and river discharge. Atmos.–Ocean 27: 682–691. DOI: 10.1080/07055900.1989.9649361

  40. McClelland, JW, Déry, SJ, Peterson, BJ, Holmes, RM and Wood, EF. 2006. A pan-arctic evaluation of changes in river discharge during the latter half of the 20th century. Geo phys. Res. Lett 33: L06715. DOI: 10.1029/2006GL025753

  41. McKay, GA and Gray, DM. 1981. The distribution of snowcover. In: Gray, DM and Male, DH (eds.), Handbook of Snow: Principles, Processes, Management and Use . Toronto: Pergamon Press, pp. 153–190.

  42. Mekonnen, MA, Wheater, HS, Ireson, AM, Spence, C, Davison, B and Pietroniro, A. 2014. Towards an improved land surface scheme for prairie landscapes. J. Hydrol 511: 105–116. DOI: 10.1016/j.jhydrol.2014.01.020

  43. Mekonnen, BA, Nazemi, A, Mazurek, KA, Elshorbagy, A and Putz, G. 2015. Hybrid modeling approach to prairie hydrology: fusing data-driven and process-based hydrological models, Hydrol Sci. J 60(9): 1473–1489. DOI: 10.1080/02626667.2014.935778

  44. Mekonnen, BA, Mazurek, KA and Putz, G. 2016. Incorporating landscape depression heterogeneity into the Soil and Water Assessment Tool (SWAT) using a probability distribution. Hydrol. Process 30: 2373–2389. DOI: 10.1002/hyp.10800

  45. Mengistu, SG and Spence, C. 2016. Testing the ability of semidistributed hydrological model to simulate contributing area. Wat. Resour. Res 52: 4399–4415. DOI: 10.1002/2016WR018760

  46. Mlynowski, TJ, Hernández-Henríquez, MA and Déry, SJ. 2011. An evaluation of hydrometric monitoring across the Canadian pan-Arctic region, 1950–2008. Hydrol. Res 42: 479–490. DOI: 10.2166/nh.2011.105

  47. Moriasi, DN, Arnold, JG, Van Liew, MW, Bingner, RL, Harmel, RD and Veith, TL. 2007. Model evaluation guidelines for systematic quantification of accuracy in watershed simulations. Trans. Am. Soc. Agri. Bio. Eng 50: 885–900. DOI: 10.13031/2013.23153

  48. Moriasi, DN, Gitau, MW and Daggupati, P. 2015. Hydrologic and water quantity models: performance measures and evaluation criteria. Transactions of the ASABE 58(6): 1763–1785. DOI: 10.13031/trans.58.10715

  49. Muhammad, A, Evenson, GR, Stadnyk, TA, Boluwade, A, Jha, SK and Coulibaly, P. 2019. Impact of model structure on the accuracy of hydrological modeling of a Canadian Prairie watershed. J. Hydrol.: Reg. Stud 21: 40–56. DOI: 10.1016/j.ejrh.2018.11.005

  50. Nachtergaele, F, van Velthuizen, H, Verelst, L and Wiberg, D. 2012. Harmonized world soil database version 1.2. Laxenburg, Austria: FAO, Rome and IIASA.

  51. Nasab, MT, Zhang, J and Chu, X. 2017. A new depression-dominated delineation (D-cubed) method for improved watershed modeling. Hydrol. Process 31: 3364–3378. DOI: 10.1002/hyp.11261

  52. Nash, JE and Sutcliffe, JV. 1970. River flow forecasting through conceptual models part I — A discussion of principles, J. Hydrol 10: 282–290. DOI: 10.1016/0022-1694(70)90255-6

  53. Natural Resources Canada. 2008. Atlas of Canada 1,000,000 National Frameworks Data, Hydrology – Dams. Phillips, RW, Spence, C and Pomeroy, JW, 2011. Connectivity and runoff dynamics in heterogeneous basins.

  54. Pechlivanidis, IG and Arheimer, B. 2015. Large-scale hydrological modeling by using modified PUB recommendations: the India-HYPE case. Hydrol. Earth Sys. Sci 19: 4559–4579. DOI: 10.5194/hess-19-4559-2015

  55. Pers, C, Temnerud, J and Lindström, G. 2016. Modeling water, nutrients, and organic carbon in forested catchments: a HYPE application. Hydrol. Process 30: 3252–3273. DOI: 10.1002/hyp.10830

  56. Pitman, AJ, Slater, AG, Desborough, CE and Zhao, M. 1999. Uncertainty in the simulation of runoff due to the parameterization of frozen soil moisture using the Global Soil Wetness Project methodology. J. Geophys. Res.: Atmos 104: 16879–16888. DOI: 10.1029/1999JD900261

  57. Pomeroy, JW, Gray, DM, Brown, T, Hedstrom, NR, Quinton, WL, Granger, RJ and Carey, SK. 2007. The cold regions hydrological model: a platform for basing process representation and model structure on physical evidence. Hydrol. Process 21: 2650–2667. DOI: 10.1002/hyp.6787

  58. Prairie Farm Rehabilitation Administration (PFRA) Hydrology Division. 1983. The Determination of Gross and Effective Drainage Areas in the Prairie Provinces. Regina: Agriculture Canada, PFRA Engineering Branch Hydrology Report 104, 22p.

  59. Qin, Y, Lei, H, Yang, D, Gao, B, Wang, H, Cong, Z and Fan, W. 2016. Long-term change in the depth of seasonally frozen ground and its ecohydrological impacts in the Qilian Mountains, northeastern Tibetan Plateau. J. Hydrol 542: 204–221. DOI: 10.1016/j.jhydrol.2016.09.008

  60. Refsgaard, JC. 1997. Parameterisation, calibration and validation of distributed hydrological models. J. Hydrol 198: 69–97. DOI: 10.1016/S0022-1694(96)03329-X

  61. Ridenour, N, Hu, X, Jafarikhasragh, S, Landy, JC, Lukovich, JV, Stadnyk, TA, Sydor, K, Myers, PG and Barber, DG. 2019. Sensitivity of freshwater dynamics to model resolution and river discharge forcing in Hudson Bay Complex. J. Marine Sys 196: 48–64. DOI: 10.1016/j.jmarsys.2019.04.002

  62. Rousseeuw, PJ. 1987. Silhouettes: A graphical aid to the interpretation and validation of cluster analysis. J. Comp. App. Math 20: 53–65. DOI: 10.1016/0377-0427(87)90125-7

  63. Saucier, F and Dionne, J. 1998. A 3-D coupled ice-ocean model applied to Hudson Bay, Canada: The seasonal cycle and time-dependent climate response to atmospheric forcing and runoff. J. Geophys. Res 103(C12): 27689–27705. DOI: 10.1029/98JC02066

  64. Saucier, F, Senneville, S, Prinsenberg, S, Roy, F, Smith, G, Gachon, P, Caya, D and Laprise, R. 2004. Modeling the sea ice–ocean seasonal cycle in Hudson Bay and Hudson Strait, Canada. Clim. Dyn 23: 303–326. DOI: 10.1007/s00382-004-0445-6

  65. Shook, K, Pomeroy, JW, Spence, C and Boychuk, L. 2013. Storage dynamics simulations in prairie wetland hydrology models: evaluation and parameterisation. Hydrol. Process. 27: 1875–1889. DOI: 10.1002/hyp.9867

  66. Shook, KR and Pomeroy, JW. 2011. Memory effects of depressional storage in northern prairie hydrology. Hydrol. Process. 25: 3890–3898. DOI: 10.1002/hyp.8381

  67. Slater, AG, Pitman, AJ and Desborough, CE. 1998. The validation of a snow parameterization designed for use in general circulation models. International J. Climat. 18: 595–617. DOI: 10.1002/(SICI)1097-0088(199805)18:6<595::AID-JOC275>3.0.CO;2-O

  68. Smith, S, Romanovsky, V, Lewkowicz, A, Burn, C, Allard, M, Clow, G, Yoshikawa, K and Throop, J. 2010. Thermal state of permafrost in North America: a contribution to the international polar year. Permafrost Periglac. Process 21: 117–135. DOI: 10.1002/ppp.690

  69. Spence, C and Woo, MK. 2003. Hydrology of subarctic Canadian Shield: soilfilled valleys. J. Hydrol 279: 151–166. DOI: 10.1016/S0022-1694(03)00175-6

  70. Stadnyk, TA, Déry, S, MacDonald, M and Koenig, K. 2019. Theme 1: Physical Environment. Iv. The Freshwater System. In: Kuzyk, ZZ and Candlish, L (eds.), From Science to Policy in the Greater Hudson Bay: An Integrated Regional Impact Study (IRIS) of Climate Change and Modernization . ArcticNet, 113–155.

  71. Stichling, W and Blackwell, SR. 1957. Drainage areas as a hydrologic factor on the Canadian Prairies. International Union of Geodesy and Geophysics Proceedings 111, 12p.

  72. Viglione, A, Parajka, J, Rogger, M, Salinas, JL, Laaha, G, Sivapalan, M and Blöschl, G. 2013. Comparative assessment of predictions in ungauged basins – Part 3: Runoff signatures in Austria. Hydrol. Earth Sys. Sci 17: 2263–2279. DOI: 10.5194/hess-17-2263-2013

  73. Vrugt, JA, Ter Braak, CJF, Diks, CGH, Hyman, JM, Robinson, BA and Higdon, D. 2009. Accelerating Markov Chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling. Int. J. Nonlinear Sci. Num. Sim. 10: 271–288. DOI: 10.1515/IJNSNS.2009.10.3.273

  74. Wang, L, Koike, T, Yang, K, Jin, R and Li, H. 2010. Frozen soil parameterization in a distributed biosphere hydrological model. Hydrol. Earth Sys. Sci 14: 557–571. DOI: 10.5194/hess-14-557-2010

  75. WGMS. 2012. Fluctuations of Glaciers 2005–2010 (Vol. X): Zemp, M, Frey, H, Gärtner-Roer, I, Nussbaumer, SU, Hoelzle, M, Paul, F and Haeberli, W (eds.). Zurich, Switzerland: ICSU (WDS)/IUGG (IACS)/UNEP/UNESCO/WMO, World Glacier Monitoring Service. Based on database version. DOI: 10.5904/wgms-fog-2012-11

  76. Woo, MK and Winter, TC. 1993. The role of permafrost and seasonal frost in the hydrology of northern wetlands in North America. J. Hydrol 141: 5–31. DOI: 10.1016/0022-1694(93)90043-9

  77. Wood, EF, Lettenmaier, DP and Zartarian, VG. 1992. A land-surface hydrology parameterization with subgrid variability for general circulation models. J. Geophys. Res 97: 2717–2728. DOI: 10.1029/91JD01786

  78. Zhao, L and Gray, DM. 1997. A parametric expression for estimating infiltration into frozen soils. Hydrol. Process 11: 1761–1775. DOI: 10.1002/(SICI)1099-1085(19971030)11:13<1761::AID-HYP604>3.0.CO;2-O

  79. Zhao, L and Gray, DM. 1999. Estimating snowmelt infiltration into frozen soils. Hydrol. Process 13: 1827–1842. DOI: 10.1002/(SICI)1099-1085(199909)13:12/13<1827::AID-HYP896>3.0.CO;2-D

comments powered by Disqus