Dynamic contributing areas, various fill-and-spill mechanisms and cold-region processes make the hydrological modelling of the Prairies very challenging. Several models (from simple conceptual to advanced process-based) are available, but the focus has been largely in reproducing streamflow. Few studies have assimilated soil moisture and other hydrological fluxes for improved simulation, but the emphasis has been predominately on simulating contributing areas. However, previous research has shown that the contributing areas are dynamic, and can vary from one year to the next, depending on hydro-meteorological conditions. Therefore, the areas deemed non-contributing can also occasionally contribute to streamflow. In this study, we introduce a progressive two-stage calibration strategy to constrain soil moisture in non-contributing areas. We demonstrate that constraining soil moisture in non-contributing areas can result in improved hydrological simulations and more realistic process representations. The Nash–Sutcliffe efficiency (NSE) values for simulated soil moisture in contributing areas increased by 68% at 20 cm and 25% at 50 cm soil depths during validation when non-contributing areas were constrained. This further led to increases in NSE values in streamflow simulation during calibration (6%) and validation (12%). Our findings suggest that soil moisture in non-contributing areas should be properly constrained for improved modelling of Prairie catchments.

The Prairie region covers approximately 900,000 km2 and spans from north-central Iowa in the United States to central Alberta in Canada (Daniel & Staricka 2000). Agriculture is a key land use of the region, and hence, the economic prosperity of the region is heavily dependent on water (Pomeroy et al. 2009). The general topography of the Prairies is a hummocky terrain that consists of millions of closed depressions known as ‘prairie potholes’ or simply ‘ponds’. Wetlands develop in the depressions around the ponds and provide ecological services, including habitat for biodiversity (Johnson et al. 2010), flood water storage (Ehsanzadeh et al. 2012) and reduction of sedimentation and nutrient loading (Nachshon et al. 2013).

In Canada, the Prairies are located in the provinces of Alberta, Saskatchewan and Manitoba. The Canadian Prairies have a cold and semiarid climate with seasonally frozen ground, and low-angled, undulating topography. Annual precipitation ranges from approximately 300–400 mm, and the majority of the region sees continuous snow cover and frozen soils during the 4–5 months of winter (Pomeroy et al. 2009). Non-contributing areas in the Prairies refer to closed watersheds around the potholes/ponds since under normal hydro-meteorological conditions, the water never reaches the receiving river systems (see Figure 1 for the map of non-contributing areas). Instead, the adjacent network of ponds receives the surface runoff from the non-contributing areas, which exchange water by a fill-and-spill mechanism and ultimately feed a terminal pond. Terminal ponds lose water to evaporation and, potentially, to a local ground water flow system (Daniel & Staricka 2000). However, connectivity in the non-contributing areas is dynamic, leading to variable contributions to nearby river systems, especially during extreme flood events.

Figure 1

Static map of non-contributing areas of the Prairies. The shape file for non-contributing areas was obtained from Agriculture & Agri-Food Canada (2013).

Figure 1

Static map of non-contributing areas of the Prairies. The shape file for non-contributing areas was obtained from Agriculture & Agri-Food Canada (2013).

Close modal

Mekonnen et al. (2014) provides a summary of water balances of individual potholes in the Prairies. The sources of inflows include snowfall and snow redistribution by wind from areas of sparse vegetation to potholes, snowmelt runoff from upland areas and potential overland runoff during intense rainfall events. Overland flows from upland areas in summer are rare. The outflows occur mainly in the form of evaporation and infiltration. However, in extreme wet conditions, fill-and-spill can occur connecting one pothole to another and even contribute runoff to nearby streams. The contribution of groundwater can be minimal when underlying tills have low hydraulic conductivity. The runoff volume is dynamic and can change significantly between one year to another for similar hydro-meteorological conditions due to antecedent soil moisture and water storage (Shaw et al. 2012).

Prairie wetland hydrology is extremely complex and thus presents unique challenges in simulating their hydrological responses. Don M. Gray (known as the father of Canadian hydrology) famously called Canadian Prairies as ‘a graveyard of hydrological models’ since no existing hydrological models at the time were able to satisfactorily simulate Prairie hydrology (Shook 2012). The major modelling challenge is due to non-contributing areas that occasionally contribute to stream networks. Runoff does not occur as a simple series of cascades to the basin outlet due to irregular distribution and morphology of potholes. The contributing and non-contributing areas in the basin are dynamic and can change each year. For any year, the area of the basin that would contribute flow to the total streamflow would depend upon snowmelt, rainfall, water storage in the pond, antecedent soil moisture conditions and basin topography. But, for improved water management in Prairies, it is essential to have modelling capabilities that can simulate the changes in water storage in wetland complexes over time.

Different methods and models have been proposed in the literature for modelling potholes and non-contributing areas of the Prairie landscape. Su et al. (2000) carried out some minor modifications to the semi-distributed model SLURP to simulate water-level variations in Prairie potholes. A majority of research has been on improving the Soil Water Assessment Tool (SWAT) to simulate the Prairie landscape. Yang et al. (2010) introduced a concept of ‘hydrologic equivalent wetland’ of Wang et al. (2008) to represent wetlands within the SWAT to study wetland restoration scenarios. The same model was later used by Perez-Valdivia et al. (2017) to assess the impact of wetland drainage using hypothetical scenarios of non-contributing drainage areas. Mekonnen et al. (2015) merged the SWAT with artificial neural networks (ANNs), so that contributing areas were simulated by the SWAT and non-contributing areas by ANNs. They reported that the fusion of process-based models with data-driven models improved model performance. Other modifications to the SWAT were also done by Evenson et al. (2016), Mekonnen et al. (2016) and Muhammad et al. (2018) for improved simulation of the Prairie landscape. The Cold Region Hydrological Model (CRHM) has also been used to simulate the Prairie landscape developing Prairie relevant modules, which have been reported to improve hydrological simulations (e.g. Fang et al. 2010; Cordeiro et al. 2017). Similarly, other models, such as the Pothole Complex Hydrological Model (Liu & Schwartz 2011), the Wetland Digital Elevation Ponding Model and the Pothole Cascade Model (Shook et al. 2013) and other conceptual models (e.g. Shook & Pomeroy 2011; Shaw et al. 2012), have also been applied. New advances include the application of powerful distributed numerical models such as HydroGeoSphere and Parflow (e.g. Nussbaumer 2014; Amado et al. 2018).

Recent developments have extended to the representation of dynamic contributing areas in land-surface hydrological modelling systems, allowing not only the simulation of large Prairie catchments but also coupling to atmospheric models. Mekonnen et al. (2014) introduced a runoff-generation algorithm based on probability density functions (PDMROF) to represent the spatial variation of pothole storages in the land-surface hydrological model, MESH. A pilot evaluation showed MESH-PDMROF, not only improved simulations of streamflow but also captured the dynamics of contributing areas. Mengistu & Spence (2016) further tested the ability of MESH-PDMROF to replicate contributing areas and found satisfactory results when compared with independently derived contributing areas from satellite imagery. However, the focus of earlier hydrological modelling studies have only been on replicating streamflows and in a few cases mapping of dynamic contributing areas without explicitly considering other important fluxes, despite strong consensus in the hydrology community that single-objective calibrations (using only streamflow) may not always result in the right internal responses of the catchment and may suffer from parameter equifinality (Beven & Freer 2001; Rakovec et al. 2016).

Recently, Budhathoki et al. (2020) simulated an experimental Prairie catchment in Saskatchewan using the MESH-PDMROF and showed that incorporating soil moisture information in model calibration improves not only the replication of streamflow but also other fluxes such as evapotranspiration, and soil moisture in deeper soil layers. However, their study only considered soil moisture information from a contributing area of the basin. Literature review also showed that the state variables from non-contributing areas are not generally considered in hydrological modelling studies, despite the fact that contributing areas in the Prairies are dynamic and can significantly vary from one year to another, which is a major motivation for this study. In this study, we present a progressive two-stage calibration method in which soil moisture in non-contributing areas is calibrated first against observed soil moisture. In the second stage, soil moisture and streamflow in contributing areas are calibrated while using calibrated parameters for non-contributing areas (obtained from the first stage). The modelling results show that constraining soil moisture information in non-contributing areas can result in improved hydrological simulation in a Prairie catchment.

Study basin description

The Brightwater Creek (BWC) basin lies within the sub-basin of the South Saskatchewan River basin. The digital elevation map of the study site shows that the elevation of the basin ranges from 547 to 646 m. The gross drainage area of the BWC basin is approximately 900 km2 at the gauging station 05HG002 (operated by the Water Survey of Canada), of which only 282 km2 is contributing since about 75% of the total area in BWC consists of numerous potholes that does not contribute surface runoff to the streamflow (Agriculture & Agri-Food Canada 2013), see also Figure 2. The contributing drainage area is expected to contribute flow to the main stream channel during a flood with a return period of 2 years (Godwin & Martin 1975). Between 2009 and 2014, the mean annual precipitation in the basin was about 330 mm. The streamflow is intermittent and mostly occurs during the spring snowmelt period. The landcover of the basin predominantly consists of crops with some patches of grass. The soil texture in the area ranges from loam to clay loam.

Figure 2

Brightwater Creek Basin in Saskatchewan, Canada.

Figure 2

Brightwater Creek Basin in Saskatchewan, Canada.

Close modal

The semi-distributed hydrologic land-surface scheme

The hydrological modelling in this study was carried out using Modelisation Environmentale-Surface et Hydrologie (MESH), which is a semi-distributed physically based land-surface modelling system developed by Environment and Climate Change Canada (Pietroniro et al. 2007). It uses the Grouped Response Unit (GRU) approach to combine areas of similar hydrological behaviour to reduce complexity and heterogeneity in the drainage basin. The vertical energy and water balance of soil, vegetation and snow is computed using the Canadian Land Surface Scheme (CLASS) (Verseghy 1991; Verseghy et al. 1993) or the Soil, Vegetation and Snow (SVS) Scheme (Husain et al. 2016). The model's WATROF (Soulis et al. 2000) or PDMROF (Mekonnen et al. 2014) module simulates the lateral movement of soil and surface water to the drainage system. The WATROF uses a sloped hillslope concept where the overland flow and the interflow are represented using Manning's and Richard's equations, respectively. The PDMROF is based on the concept of the probability density approach (Moore 2007) to parsimoniously represent the runoff and storage. The streamflow through the river channels is routed using the WATROUTE module. WATROUTE mainly relies on flow directions and elevation data, and routes the surface runoff and recharge to the basin outlet using a continuity equation and Manning's formula (Kouwen 2016).

MESH-PDMROF model set-up

The MESH model was set up for the BWC basin with an outlet at Brightwater Creek near Kenaston (05HG002) (see Figure 3(a)). The drainage database was prepared using the Environment and Climate Change Canada's Green Kenue software (Canadian Hydraulics Centre 2010). The digital elevation data for the model were obtained from the Geobase database at the scale of 1:50,000, and the land cover data were derived from the Commission for Environmental Cooperation (CEC) in 30 m resolution (http://cec.org/tools-and-resources/map-files/land-cover-2010-landsat-30 m). Within the MESH model, the CLASS was selected for the vertical exchange of water and energy within a grid cell, the PDMROF for lateral soil (sub-surface) and surface water movement and WATROUTE for routing streamflow in river channels.

Figure 3

The study site. (a) Brightwater Creek Basin with tipping bucket station locations. (b) GRU delineation for contributing and non-contributing areas in Green Kenue for MESH (adapted from Budhathoki et al. (2020)).

Figure 3

The study site. (a) Brightwater Creek Basin with tipping bucket station locations. (b) GRU delineation for contributing and non-contributing areas in Green Kenue for MESH (adapted from Budhathoki et al. (2020)).

Close modal
In this study, the PDMROF module was used because it better represents the variable nature of contributing and non-contributing areas of the Prairies (Mengistu & Spence 2016). The PDMROF was specifically designed by Mekonnen et al. (2014) to simulate the complex hydrological behaviour of the Prairies, by employing the concept of the probability density approach of Moore (2007) to parsimoniously represent their runoff and storage processes. The basic equation is given as follows:
where F(C) is the distribution function of storage capacity, C indicates the actual storage capacity, is the maximum storage capacity and B represents a shape factor parameter. The two parameters, and B, control the degree of spatial variability of the storage capacity across the basin. The unit of C is ‘m’ and B is a unitless parameter. A low value of ‘B’ indicates a highly disconnected pothole distribution, whereas a higher one suggests a well-connected pothole distribution (Mekonnen et al. 2014).

The vegetation parameters for the model were obtained from the CLASS manual (Verseghy 2009) and various literature sources (e.g. Pan et al. 2017; Rokaya et al. 2019), and the soil texture data were obtained from the Agriculture and Agri-Food Canada's Real-Time In-Situ Soil Monitoring for Agriculture Network metadata (http://wcag.fieldvision.ca/) for different soil layers. The model was set up at a grid resolution of 0.05° (approx. 4 km), resulting in 66 grid cells. The grid size was selected taking into account the spatial resolution of other input data and computational requirements. Some input data such as land cover data were available at a finer spatial resolution (e.g. 30 m), but other input data such as meteorological forcing data had coarser spatial resolutions (e.g. 15 and 10 km), thus from a basin heterogeneity and computational perspective, 4 km resolution was deemed suitable.

Two GRUs, contributing and non-contributing, were created with similar crop vegetation on both GRUs based on landuse data (see Figure 3(b)). The soil layer in the model was conceptualized by six layers at the depths of 10, 30, 70, 130, 300 and 400 cm from the surface (a total depth of 400 cm), in order to get the average soil moisture data at 5, 20, 50, 100, 215 and 350 cm. The observed soil moisture data were available at 5, 20 and 50 cm soil depths; therefore, the above-mentioned soil layer discretization facilitated the comparison of simulated soil moisture at 20 and 50 cm soil depths with the field observations. Although measured soil moisture data at 5 cm were also available, it was not used in this study since Yassin et al. (2017) and Budhathoki et al. (2020) found that at 5 cm, soil moisture was dominated by streamflow and therefore, not sensitive to soil moisture calibration.

Meteorological forcing data

The meteorological inputs for MESH include incoming shortwave radiation, incoming longwave radiation, precipitation, temperature, barometric pressure, specific humidity and wind speed. All of these inputs, except precipitation data, were retrieved from the Global Environmental Multi-scale (GEM) Numerical Weather Prediction (NWP) model developed by Environment and Climate Change Canada (Côté et al. 1998), which has a spatial resolution of 15 km and a temporal resolution of 1 h. The precipitation data included a hybrid dataset. The rainfall data from May to September were obtained from the tipping bucket stations NW07 and NE36 (https://www.frdr.ca/repo/handle/doi:10.20383/101.0116, also see Figure 2 for the locations of the stations) whereas, for the rest of the months, precipitation data were acquired from the Canadian Precipitation Analysis (CaPA) (Mahfouf et al. 2007), which is run by the Canadian Meteorological Centre (CMC). The CaPA covers all of North America on a 10 km grid and is available at 6 hourly and 24 hourly accumulation intervals (Fortin et al. 2015), and has been found to be a suitable forcing data for hydrological applications, particularly in the mid- to high-latitude regions with sparse meteorological observation networks. The hybrid precipitation dataset was prepared because Budhathoki et al. (2020) compared the rainfall data from CaPA and tipping bucket, and found that CaPA underestimated both mean and variability of precipitation in summer months compared with tipping bucket at the BWC. Note that tipping bucket data are not collected during the winter period. The blended precipitation product helped to reduce input data uncertainty for summer months and represent precipitation intensity and distribution as accurately as possible.

Soil moisture and streamflow data

The soil moisture data were retrieved from the Brightwater Creek Monitoring Network (https://www.frdr.ca/repo/handle/doi:10.20383/101.0116) for both the contributing area (station NE36) and the non-contributing area (station NW07). Summer data were extracted for 20 and 50 cm soil depths for both stations. Winter soil moisture data are not available at the site, since the instruments do not give reliable readings during frozen conditions as reported by Tetlock et al. (2019). Similarly, the data for spring for the study stations (i.e. NE36 and NW07) are not published online and were, therefore, acquired from the National Hydrology Research Centre of Environment and Climate Change Canada, which manages the portion of the monitoring network, including the stations used in this study. Although few additional monitoring stations were available, only stations NE36 and NW07 were selected due to issues related to data quality and computational demand. The stations used in this study have a complete quality-controlled dataset for the study period. They are also representative of the study area in terms of topography, climatic conditions and location since the entire catchment lies in the same ecozone (i.e. Prairie Ecozone) and Ecoregion (i.e. Moist Mixed Grassland). Furthermore, additional soil moisture stations would have increased the number of objective functions, making it difficult to identify the optimum solutions. The daily streamflow data were retrieved from the Water Survey of Canada through the Environment Canada's Hydat database for the basin outlet at Brightwater Creek near Kenaston (05HG002). The data were extracted for the period of 2011–2018.

Model calibration

A progressive two-stage calibration strategy was adopted for model calibration. Firstly, soil moisture at the 20 and 50 cm soil depths of the non-contributing grid cell (station NW07) was calibrated. The calibration was performed using the Optimization Software Toolkit for Research Involving Computational Heuristics (OSTRICH) (Matott 2005). For this, a multi-objective aggregate optimization using Asynchronous Parallel Dynamically Dimensioned Search (PDDS) (Tolson et al. 2014) was chosen with a weightage of 50:50% for the 20 and 50 cm soil moisture depths, respectively. The non-contributing GRU does not regularly contribute to the streamflow and hence the streamflow parameters were not calibrated in this step. Only the parameters that have been identified as sensitive to soil moisture calibration in previous studies (e.g. Haghnegahdar et al. 2017; Yassin et al. 2017) were calibrated at this stage and are presented in Table 1.

Table 1

List of calibrated parameters and their ranges for the non-contributing area

ParametersUnitDescriptionRange
DRN – Drainage index 0–1.0 
KSAT m/s Saturated surface soil conductivity 1.28 × 10−6–3.45 × 10−5 
LAMX – Annual maximum leaf-area index of the crop 4.0–6.0 
SAND-L1 Sand in soil layer 1 30–60 
SAND-L2 Sand in soil layer 2 30–60 
SAND-L3 Sand in soil layer 3 30–60 
CLAY-L1 Clay in soil layer 1 10–40 
CLAY-L2 Clay in soil layer 2 10–40 
CLAY-3 Clay in soil layer 3 10–40 
ParametersUnitDescriptionRange
DRN – Drainage index 0–1.0 
KSAT m/s Saturated surface soil conductivity 1.28 × 10−6–3.45 × 10−5 
LAMX – Annual maximum leaf-area index of the crop 4.0–6.0 
SAND-L1 Sand in soil layer 1 30–60 
SAND-L2 Sand in soil layer 2 30–60 
SAND-L3 Sand in soil layer 3 30–60 
CLAY-L1 Clay in soil layer 1 10–40 
CLAY-L2 Clay in soil layer 2 10–40 
CLAY-3 Clay in soil layer 3 10–40 

Note that if the sum of sand and clay in each soil layer does not equal to 100%, then the remaining percentage is assumed to be silt.

In the second step, a multi-objective ensemble optimization was performed with the streamflow in the contributing GRU and soil moisture at 20 and 50 cm depths in the contributing area grid cell (NE36). The algorithm used for the multi-objective optimization in OSTRICH is Parallel Pareto Archieved Dynamically Dimensioned Search (ParaPADDS) (Tolson et al. 2014). Parameters and their ranges identified by previous sensitivity analyses to be important and sensitive to Prairie catchments were calibrated. The list of parameters and their calibration ranges are provided in Table 2. For further details on MESH parameters' sensitivities and their implication in hydrological modelling, readers are referred to Haghnegahdar et al. (2017), Yassin et al. (2017) and Budhathoki et al. (2020).

Table 2

List of calibrated parameters and their ranges for the contributing area

ParametersUnitDescriptionRange
R2N1 – Channel Manning's 0.001–2.0 
R1N1 – Overbank Manning's 0.001–2.0 
CMAX Maximum surface storage capacity 0.01–10.0 
– Shape factor for Pareto distribution function 0.01–30 
DRN – Drainage index 0–1.0 
KSAT m/s Saturated surface soil conductivity 1.28 × 10−6–3.45 × 10−5 
LAMX – Annual maximum leaf-area index of the crop 4.0–6.0 
SAND-L1 Sand in soil layer 1 30–60 
SAND-L2 Sand in soil layer 2 30–60 
SAND-L3 Sand in soil layer 3 30–60 
CLAY-L1 Clay in soil layer 1 10–40 
CLAY-L2 Clay in soil layer 2 10–40 
CLAY-3 Clay in soil layer 3 10–40 
ParametersUnitDescriptionRange
R2N1 – Channel Manning's 0.001–2.0 
R1N1 – Overbank Manning's 0.001–2.0 
CMAX Maximum surface storage capacity 0.01–10.0 
– Shape factor for Pareto distribution function 0.01–30 
DRN – Drainage index 0–1.0 
KSAT m/s Saturated surface soil conductivity 1.28 × 10−6–3.45 × 10−5 
LAMX – Annual maximum leaf-area index of the crop 4.0–6.0 
SAND-L1 Sand in soil layer 1 30–60 
SAND-L2 Sand in soil layer 2 30–60 
SAND-L3 Sand in soil layer 3 30–60 
CLAY-L1 Clay in soil layer 1 10–40 
CLAY-L2 Clay in soil layer 2 10–40 
CLAY-3 Clay in soil layer 3 10–40 

Note that if the sum of sand and clay in each soil layer does not equal to 100%, then the remaining percentage is assumed to be silt.

Both PDDS and ParaPADDS are the improved parallel version of the dynamically dimensioned search (DDS) (Tolson & Shoemaker 2007), which minimizes the objective functions using an evolutionary search algorithm. The degree of agreement between simulated and observed outputs was measured using a Nash–Sutcliffe efficiency (NSE) coefficient (Nash & Sutcliffe 1970) as an objective function. The value of NSE can range between − and 1, where the efficiency closer to 1 corresponds to perfect coincidence between modelled and observed values. The calibration period included 2012–2014 and validation years included 2015–2018. The model spin-up period was 2011–2012. The model runs at 30-min time-steps, but streamflows were generated at a daily time step for comparison and analysis.

Constraining soil moisture in non-contributing areas

In the first step, soil moisture in non-contributing areas was calibrated. The calibration resulted in NSE values of 0.66 at 20 cm and 0.65 at 50 cm soil depths. The biases were +7% and −14% in 20 and 50 cm soil depths, respectively. Similarly, RMSE values were 0.04 m3/m3 at 20 cm soil depth and 0.03 m3/m3 at 50 cm soil depth. In the model validation period, NSE values of 0.80, bias of +5% and RMSE of 0.03 m3/m3 were obtained for 20 cm soil depth, whereas at 50 cm soil depth, the NSE values of 0.78, bias of −8% and RMSE of 0.029 m3/m3 were achieved. Interestingly, the model performed better during validation than calibration for both 20 and 50 cm soil depths (see also Figure 4).

Figure 4

Calibration and validation of soil moisture in non-contributing areas.

Figure 4

Calibration and validation of soil moisture in non-contributing areas.

Close modal

Simulation of soil moisture in contributing areas

In the second step, soil moisture in contributing areas along with streamflow was calibrated using calibrated parameters for non-contributing areas. Note that multi-objective calibration using ParaPADDS does not give one best solution but dominated and non-dominated (optimal) solutions with trade-offs between the objective functions. The Pareto runs (not shown here) yielded 99 optimal solutions. The ensemble simulated soil moisture (from 99 simulations) and the median of ensemble and observed soil moisture at 20 and 50 cm soil depths are presented in Figure 5. For the calibration period, the NSE of 0.51, bias of −4% and RMSE of 0.042 m3/m3 were obtained at 20 cm soil depth, whereas, at 50 cm soil depth, the NSE of 0.45, bias of −8% and RMSE of 0.05 m3/m3 were achieved when median-simulated soil moisture was compared against observed soil moisture. The validation period yielded in the NSE of 0.42 and 0.81, bias of +9% and −1%, and RMSE of 0.05 and 0.02 m3/m3 at 20 and 50 cm soil depths, respectively.

Figure 5

Calibration and validation of soil moisture in contributing areas.

Figure 5

Calibration and validation of soil moisture in contributing areas.

Close modal

Constraining the soil moisture of non-contributing areas improved the simulation of soil moisture in contributing areas, particularly during the validation period. When soil moisture of non-contributing areas was not constrained, mean NSE values of 0.44 and 0.42 during calibration and 0.25 and 0.65 during validation were obtained for 20 and 50 cm soil depths, respectively (see Budhathoki et al. (2020) for further details). Therefore, an increase in NSE values by 68% at 20 cm and 25% at 50 cm was achieved during validation after constraining soil moisture of the non-contributing area. However, the calibration period only showed some improvement in NSE values, i.e. 15% at 20 cm and 4% at 50 cm soil depths compared to the validation period.

Simulation of streamflows

The streamflow calibration resulted in NSE values ranging from 0.71 to 0.79 for the calibration period for the 99 optimal parameter sets, whereas for validation years, the values ranged from 0.48 to 0.61 (see also Figure 6). The mean bias of the ensemble was 6% during the calibration but increased to 30% during the validation period when compared with observed streamflow. Similarly, average (of ensemble) RMSE during calibration and validation were 1.0 and 1.15 m3/s, respectively. When soil moisture information in non-contributing areas were not constrained, NSE values ranged between 0.68 and 0.75 during the calibration period, and from 0.45 to 0.52 during the validation period. Therefore, constraining soil moisture in non-contributing areas yielded an increase in the NSE value by 6% on an average for the calibration period, whereas the validation period showed an increase of 12% (on an average) in the NSE value.

Figure 6

Calibration and validation of streamflow.

Figure 6

Calibration and validation of streamflow.

Close modal

Although there are not any reported studies (to the best of authors' knowledge) in the literature that have constrained soil moisture in non-contributing areas of the Prairies, other studies have constrained soil moisture in general and achieved similar improvements in streamflow simulations. Yassin et al. (2017) incorporated GRACE TWS in their model calibration and reported a slight improvement in the streamflow simulation during the validation period. Improved streamflow and soil moisture simulations were obtained by Rajib et al. (2016) when root zone soil moisture information was included in the model calibration. Similar findings were reported by Leroux et al. (2016) who assimilated remotely sensed soil moisture data in the model calibration. However, it is important to note that constraining fluxes and state variables, although improving model fidelity, does not always result in improved streamflow simulations. Note that previous studies have shown both improvements (e.g. Werth et al. 2009; Chen et al. 2017) and degradation (e.g. Livneh & Lettenmaier 2012; Rakovec et al. 2016) in their streamflow simulations when multi-objective calibration was adopted since precipitation is properly partitioned among different state variables and fluxes in multi-objective calibration.

Representation of a dynamic change in the contributing area

Although we have used a static map to separate contributing and non-contributing areas in the basin, the streamflow response of the grid cells in the basin is dynamic over time and depends on hydro-meteorological conditions. The model simulations for the first three hydrological years, i.e. 2012, 2013 and 2014 with the varying degree of hydrological conditions, are presented in Figure 7. The figure depicts the fractional contribution of each grid cell in the basin when soil moisture in the non-contributing area was either constrained or not constrained. The flow hydrograph on the left panel shows that 2012 and 2013 are relatively low flow years compared to 2014. Therefore, runoff from non-contributing areas in these low flow years is not expected since only in exceptionally wet years non-contributing areas contribute to streamflows. However, when soil moisture in non-contributing areas (see Figure 3 for contributing and non-contributing areas within the BWC basin) was not constrained, the results show that the grid cells in the non-contributing areas were still contributing for 2012 and 2013, which is an erroneous representation. But when soil moisture was constrained in non-contributing areas, the grid cells in non-contributing areas did not contribute to streamflow in low flow years such as 2012 and 2013. For 2014, which was comparatively wetter than 2012 and 2013, both approaches showed some contribution of non-contributing areas (although their ranges are different).

Figure 7

Fractional contribution of grid cells to total streamflow (a) when soil moisture in the non-contributing area was not constrained and (b) when soil moisture in the non-contributing area was constrained.

Figure 7

Fractional contribution of grid cells to total streamflow (a) when soil moisture in the non-contributing area was not constrained and (b) when soil moisture in the non-contributing area was constrained.

Close modal

Cold-region processes, variable contributing and non-contributing areas, and various fill-and-spill mechanisms pose significant challenges in modelling Prairie catchments. Different models and methods have been proposed in the literature for improved modelling of the Prairie landscapes, ranging from simple conceptual models to more advanced process-based models. However, the focus of earlier modelling studies has been largely on replicating the streamflow response. Recent studies have incorporated other observed state variables and flux data, such as soil moisture, but the attention has been only on contributing areas of the Prairies. However, previous studies have demonstrated that the contributing and non-contributing areas in the basin are dynamic and can affect the generation of total streamflow. In this study, we presented a progressive two-stage model calibration strategy where soil moistures in non-contributing areas is constrained first, and then streamflow and soil moisture in contributing areas are subsequently calibrated.

The proposed approach of constraining the non-contributing area improved the hydrological model performance. An increase in NSE values by 6% and 12% were obtained for streamflow simulations during calibration and validation periods, respectively. Larger improvements were seen on soil moisture simulations, particularly during the validation period. NSE values increased by 68% at 20 cm soil depth and by 25% at 50 cm soil depth, although the calibration period showed only a nominal increase in performance at both 20 and 50 cm soil depths. Thus, for both streamflow and soil moisture simulations, significant improvements were only obtained for the validation period. This is due to the fact that during calibration, the presence of a large number of parameters in physically based models provides a large degree of freedom; therefore, comparable results can be obtained for the periods for which the models are trained for, even without constraining other state variables. In the earlier study of the BWC, Budhathoki et al. (2020) showed that, even without using soil moisture information to constrain the model, the streamflow hydrograph can be reasonably replicated. Hence, the real test of the predictive capability of a model is not on its simulation of the calibration period (to which it is trained), but rather for the validation period, as shown in this study.

We also demonstrated that the presented methodology is able to represent the dynamic change in contributing areas in Prairies. The fractional contribution of grid cells from non-contributing and contributing areas was more realistic when soil moisture in the non-contributing area was constrained. While our results are in line with what is expected (based on current knowledge), the validation of the results remains a critical limitation in the absence of alternative data. Therefore, the model results from fractional grid contributions in different years should be taken more as proof of concept rather than viewing them as explicit values for a particular location and time. Fortunately, with new advances in science, new data sources such as remote sensing are becoming available, which could potentially provide an opportunity to validate these results. In St. Denis Creek watershed in western Canada, Mengistu & Spence (2016) were able to validate the MESH-PDMROF-generated total contributing area using the satellite imagery. Further research is recommended in evaluating and validating model generation contribution areas using remote sensing and other data products for improved understanding and simulation of Prairie watersheds.

In the simulation of soil moisture for the non-contributing area, the model appears to slightly overestimate at 20 cm soil depth and slightly underestimate at 50 cm soil depth, even though equal calibration weightages were provided to both soil layers. One important bias to note in the soil moisture simulation for both contributing and non-contributing areas is the timing of the error. Larger discrepancies are observed either at the beginning of the spring period or at the end of the summer period. These two periods are important times in the year where phase changes occur from winter to spring or from summer to fall. The biases could have stemmed partly from the inadequacy of the model to correctly simulate the energy and water balances during the phase change and partly due to errors in the observed data. Simulating soil moisture and streamflow generated during soil freeze and thaw periods have been a critical challenge in cold-region hydrology (Woo et al. 2004). Similarly, observed data are also susceptible to a range of uncertainties. In the absence of additional data for comparison, it was difficult to quantitatively assign the simulated biases to model structure errors or observed data uncertainty.

It is also important to recognize that hydrological modelling is inherently susceptible to different sources of uncertainties. We expounded on precipitation uncertainty in the ‘Meteorological forcing data’ section and our approach of using a blended precipitation product to minimize the errors since CaPA was underestimating both mean and variability of rainfall during summer months. Tipping buckets do not measure winter precipitation, and, due to the lack of alternative measured data, we were not able to evaluate if CaPA satisfactorily captures winter precipitation. Some of the errors in simulated streamflows such as underestimation of the 2015 peak event could be related to meteorological data uncertainty. Besides forcing data uncertainty, there are also errors associated with model parameters, model structure and observation data. Different combinations of parameters can give similar results, leading to equifinality (Beven & Freer 2001); therefore, parametric uncertainty poses another key challenge in hydrological modelling. To address this issue, we adopted an ensemble modelling approach (using 99 optimal parameter sets) instead of a deterministic simulation based on one optimal parameter set. Similarly, the choice of calibration period also affects the generated optimal parameters and hence, simulated streamflows. It is possible to end up with different optimal parameter sets, depending on the time period to which the model is calibrated. We calibrated the model for the 2011–2014 period and validated it for the timeframe 2015–2018. The year 2011 was considered a spin-up period to minimize the effects of initial conditions on simulated results. Our logic was to use half of the study period for calibration and other half to assess the predictive ability of the calibrated model through temporal validation. We also ensured that both relatively dry (2012, 2013, 2016 and 2017) and wet (2011, 2014, 2015 and 2018) conditions are present during both calibration and validation periods. We also acknowledge that practically, it is not possible to avoid all possible errors associated with hydrological modelling, but we used a uniform set of forcing data, model structure, parameter range and observation data between two approaches (i.e. using progressive two-stage calibration vs. calibrating both contributing and non-contributing areas at once), so that the differences observed are due to methodology, not because of forcing data or other model-related errors.

This study presents a progressive two-stage model calibration strategy where soil moisture in non-contributing areas is constrained first, and then streamflow and soil moisture in contributing areas are subsequently calibrated. Our results show that constraining soil moisture in non-contributing areas leads to improvement in hydrological process simulations (particularly during validation) and, therefore, greater confidence in hydrological modelling. The NSE values in simulated soil moisture of contributing areas increased by 68% at 20 cm soil depth and by 25% at 50 cm soil depth during validation when soil moisture in the non-contributing area was constrained. Similarly, an average increase of 6% in calibration and 12% in validation in NSE values were obtained for streamflow simulation. More importantly, the representation of dynamic variability in contributing areas was more realistic with the presented approach. The presented methodology is model and scale-independent; therefore, it can be used with other models and replicated to other small Prairie catchments or upscaled to larger Prairie watersheds. Since, the economic prosperity of the Prairie region is heavily dependent on water, the findings of this study is expected to contribute to improved water management in the region.

This research was financially supported by the Global Water Futures program. We thank Daniel Princz from Environment and Climate Change Canada for providing meteorological forcing data for model calibration and validation. We are also grateful to Erica Tetlock of Environment and Climate Change Canada and Warren Helgason and Andrew Ireson from the University of Saskatchewan for the soil moisture data. The manuscript benefited greatly from the constructive comments and feedback from three anonymous reviewers.

Agriculture and Agri-Food Canada
2013
Gross and Effective Drainage Areas for Hydrometric Gauging Stations of the AAFC Watersheds Project – 2013
.
Amado
A. A.
Politano
M.
Schilling
K.
Weber
L.
2018
Investigating hydrologic connectivity of a drained prairie pothole region wetland complex using a fully integrated, physically-based model
.
Wetlands
38
(
2
),
233
245
.
doi:10.1007/s13157-016-0800-5
.
Beven
K.
Freer
J.
2001
Equifinality, data assimilation, and uncertainty estimation in mechanistic modelling of complex environmental systems using the GLUE methodology
.
Journal of Hydrology
249
(
1–4
),
11
29
.
http://dx.doi.org/10.1016/S0022-1694(01)00421-8
.
Budhathoki
S.
Rokaya
P.
Lindenschmidt
K. E.
Davison
B.
2020
A multi-objective calibration approach using in-situ soil moisture data for improved hydrological simulation of the Prairies
.
Hydrological Science Journal
65
(
4
),
638
649
.
Canadian Hydraulics Centre
2010
Green Kenue Reference Manual
.
National Research Council
,
Ottawa, Ontario, Canada
,
340 pp
.
Cordeiro
M. R.
Wilson
H. F.
Vanrobaeys
J.
Pomeroy
J. W.
Fang
X.
2017
Simulating cold-region hydrology in an intensively drained agricultural watershed in Manitoba, Canada, using the cold regions hydrological model
.
Hydrology and Earth System Sciences
21
(
7
),
3483
3506
.
Côté
J.
Gravel
S.
Méthot
A.
Patoine
A.
Roch
M.
Staniforth
A.
1998
The operational CMC–MRB global environmental multiscale (GEM) model. Part I: design considerations and formulation
.
Monthly Weather Review
126
(
6
),
1373
1395
.
Daniel
J. A.
Staricka
J. A.
2000
Frozen soil impact on ground water-surface water interaction
.
JAWRA Journal of the American Water Resources Association
36
(
1
),
151
160
.
Ehsanzadeh
E.
Spence
C.
van der Kamp
G.
McConkey
B.
2012
On the behaviour of dynamic contributing areas and flood frequency curves in North American Prairie watersheds
.
Journal of Hydrology
414
,
364
373
.
Evenson
G. R.
Golden
H. E.
Lane
C. R.
D'amico
E.
2016
An improved representation of geographically isolated wetlands in a watershed-scale hydrologic model
.
Hydrological Processes
30
(
22
),
4168
4184
.
Fang
X.
Pomeroy
J.
Westbrook
C.
Guo
X.
Minke
A.
Brown
T.
2010
Prediction of snowmelt derived streamflow in a wetland dominated prairie basin
.
Hydrology and Earth System Sciences
14
(
6
),
991
1006
.
Fortin
V.
Roy
G.
Donaldson
N.
Mahidjiba
A.
2015
Assimilation of radar quantitative precipitation estimations in the Canadian Precipitation Analysis (CaPA)
.
Journal of Hydrology
531
,
296
307
.
Godwin
R.
Martin
F.
1975
Calculation of gross and effective drainage areas for the Prairie Provinces
. In
Proceedings of the Canadian Hydrology Symposium
.
Husain
S. Z.
Alavi
N.
Bélair
S.
Carrera
M.
Zhang
S.
Fortin
V.
Gauthier
N.
2016
The multibudget Soil, Vegetation, and Snow (SVS) scheme for land surface parameterization: offline warm season evaluation
.
Journal of Hydrometeorology
17
(
8
),
2293
2313
.
Johnson
W. C.
Werner
B.
Guntenspergen
G. R.
Voldseth
R. A.
Millett
B.
Naugle
D. E.
Olawsky
C.
2010
Prairie wetland complexes as landscape functional units in a changing climate
.
BioScience
60
(
2
),
128
140
.
Kouwen
N.
2016
WATFLOOD/WATROUTE Hydrological Model Routing and Flood Forecasting System, User's Manual
.
University of Waterloo
,
Waterloo, ON
.
Leroux
D. J.
Pellarin
T.
Vischel
T.
Cohard
J.-M.
Gascon
T.
Gibon
F.
Seguis
L.
2016
Assimilation of SMOS soil moisture into a distributed hydrological model and impacts on the water cycle variables over the Ouémé catchment in Benin
.
Hydrology and Earth System Sciences
20
(
7
),
2827
2840
.
Liu
G.
Schwartz
F. W.
2011
An integrated observational and model-based analysis of the hydrologic response of prairie pothole systems to variability in climate
.
Water Resources Research
47
(
2
),
W02504
.
Livneh
B.
Lettenmaier
D. P.
2012
Multi-criteria parameter estimation for the Unified Land Model
.
Hydrology and Earth System Sciences
16
(
8
),
3029
3048
.
doi:10.5194/hess-16-3029-2012
.
Mahfouf
J. F.
Brasnett
B.
Gagnon
S.
2007
A Canadian precipitation analysis (CaPA) project: description and preliminary results
.
Atmosphere-Ocean
45
(
1
),
1
17
.
Matott
L. S.
2005
OSTRICH: An Optimization Software Tool: Documentation and Users Guide
.
University at Buffalo
,
Buffalo, NY
.
Mekonnen
M.
Wheater
H.
Ireson
A.
Spence
C.
Davison
B.
Pietroniro
A.
2014
Towards an improved land surface scheme for prairie landscapes
.
Journal of Hydrology
511
,
105
116
.
Mekonnen
B. A.
Nazemi
A.
Mazurek
K. A.
Elshorbagy
A.
Putz
G.
2015
Hybrid modelling approach to prairie hydrology: fusing data-driven and process-based hydrological models
.
Hydrological Sciences Journal
60
(
9
),
1473
1489
.
Mekonnen
B. A.
Mazurek
K. A.
Putz
G.
2016
Incorporating landscape depression heterogeneity into the Soil and Water Assessment Tool (SWAT) using a probability distribution
.
Hydrological Processes
30
(
13
),
2373
2389
.
Mengistu
S.
Spence
C.
2016
Testing the ability of a semidistributed hydrological model to simulate contributing area
.
Water Resources Research
52
(
6
),
4399
4415
.
Moore
R. J.
2007
The PDM rainfall-runoff model
.
Hydrology and Earth System Sciences
11
(
1
),
483
499
.
doi:10.5194/hess-11-483-2007
.
Muhammad
A.
Evenson
G.
Stadnyk
T.
Boluwade
A.
Jha
S.
Coulibaly
P.
2018
Assessing the importance of potholes in the Canadian Prairie Region under future climate change scenarios
.
Water
10
(
11
),
1657
.
Nachshon
U.
Ireson
A.
van der Kamp
G.
Wheater
H.
2013
Sulfate salt dynamics in the glaciated plains of North America
.
Journal of Hydrology
499
,
188
199
.
Nash
J.
Sutcliffe
J.
1970
River flow forecasting through conceptual models part I – a discussion of principles
.
Journal of Hydrology
10
(
3
),
282
290
.
Nussbaumer
R.
2014
A Physically Based Model for Groundwater-Pond Interactions in the Canadian Prairies
.
MSc, Imperial College London
,
London
,
UK
.
Pan
X.
Helgason
W.
Ireson
A.
Wheater
H.
2017
Field-scale water balance closure in seasonally frozen conditions
.
Hydrology and Earth System Sciences
21
(
11
),
5401
5413
.
Pietroniro
A.
Fortin
V.
Kouwen
N.
Neal
C.
Turcotte
R.
Davison
B.
Evora
N.
2007
Development of the MESH modelling system for hydrological ensemble forecasting of the Laurentian Great Lakes at the regional scale
.
Hydrology and Earth System Sciences Discussions
11
(
4
),
1279
1294
.
Pomeroy
J. W.
Fang
X.
Williams
B.
2009
Impacts of climate change on Saskatchewan's water resources
.
Centre for Hydrology, University of Saskatchewan
,
Saskatoon, Canada
, pp.
1
46
.
Rajib
M. A.
Merwade
V.
Yu
Z.
2016
Multi-objective calibration of a hydrologic model using spatially distributed remotely sensed/in-situ soil moisture
.
Journal of Hydrology
536
,
192
207
.
Rakovec
O.
Kumar
R.
Attinger
S.
Samaniego
L.
2016
Improving the realism of hydrologic model functioning through multivariate parameter estimation
.
Water Resources Research
52
(
10
),
7779
7792
.
Rokaya
P.
Morales-Marín
L.
Bonsal
B.
Wheater
H.
Lindenschmidt
K.-E.
2019
Climatic effects on ice phenology and ice-jam flooding of the Athabasca River in western Canada
.
Hydrological Sciences Journal
64
(
11
),
1265
1278
.
doi:10.1080/02626667.2019.1638927
.
Shaw
D. A.
Vanderkamp
G.
Conly
F. M.
Pietroniro
A.
Martz
L.
2012
The fill–spill hydrology of prairie wetland complexes during drought and deluge
.
Hydrological Processes
26
(
20
),
3147
3156
.
doi:10.1002/hyp.8390
.
Shook
K.
2012
Complexity of Prairie Hydrology: Distinctiveness of Prairie Hydrology
.
Centre for Hydrology, University of Saskatchewan, Saskatoon, Canada. Available from: https://wiki.usask.ca/download/attachments/516948256/01_ComplexityOfPrairieHydrology_KShook.pdf (accessed on 15 February 2020)
.
Shook
K. R.
Pomeroy
J. W.
2011
Memory effects of depressional storage in Northern Prairie hydrology
.
Hydrological Processes
25
(
25
),
3890
3898
.
Shook
K.
Pomeroy
J. W.
Spence
C.
Boychuk
L.
2013
Storage dynamics simulations in prairie wetland hydrology models: evaluation and parameterization
.
Hydrological Processes
27
(
13
),
1875
1889
.
Soulis
E.
Snelgrove
K.
Kouwen
N.
Seglenieks
F.
Verseghy
D.
2000
Towards closing the vertical water balance in Canadian atmospheric models: coupling of the land surface scheme CLASS with the distributed hydrological model WATFLOOD
.
Atmosphere-Ocean
38
(
1
),
251
269
.
Su
M.
Stolte
W.
Van der Kamp
G.
2000
Modelling Canadian prairie wetland hydrology using a semi-distributed streamflow model
.
Hydrological Processes
14
(
14
),
2405
2422
.
Tolson
B. A.
Shoemaker
C. A.
2007
Dynamically dimensioned search algorithm for computationally efficient watershed model calibration
.
Water Resources Research
43
(
1
),
W01413
.
Tolson
B. A.
Sharma
V.
Swayne
D. A.
2014
Parallel implementations of the Dynamically Dimensioned Search (DDS) algorithm
. In:
International Symposium on Environmental Software Systems, 22–25 May, Prague
, pp.
573
582
.
Verseghy
D. L.
1991
Class – a Canadian land surface scheme for GCMS. I. Soil model
.
International Journal of Climatology
11
(
2
),
111
133
.
doi:10.1002/joc.3370110202
.
Verseghy
D.
2009
CLASS – The Canadian Land Surface Scheme (Version 3.4), Technical Documentation (Version 1.1). Climate Research Division, Science and Technology Branch, Environment Canada, 180
.
Verseghy
D. L.
McFarlane
N. A.
Lazare
M.
1993
Class – a Canadian land surface scheme for GCMS, II. Vegetation model and coupled runs
.
International Journal of Climatology
13
(
4
),
347
370
.
doi:10.1002/joc.3370130402
.
Wang
X.
Yang
W.
Melesse
A. M.
2008
Using hydrologic equivalent wetland concept within SWAT to estimate streamflow in watersheds with numerous wetlands
.
Transactions of the ASABE
51
(
1
),
55
72
.
Werth
S.
Güntner
A.
Petrovic
S.
Schmidt
R.
2009
Integration of GRACE mass variations into a global hydrological model
.
Earth and Planetary Science Letters
277
(
1
),
166
173
.
https://doi.org/10.1016/j.epsl.2008.10.021
.
Woo
M. K.
Arain
M.
Mollinga
M.
Yi
S.
2004
A two-directional freeze and thaw algorithm for hydrologic and land surface modelling
.
Geophysical Research Letters
31
(
12
),
L12501
.
Yang
W.
Wang
X.
Liu
Y.
Gabor
S.
Boychuk
L.
Badiou
P.
2010
Simulated environmental effects of wetland restoration scenarios in a typical Canadian prairie watershed
.
Wetlands Ecology and Management
18
(
3
),
269
279
.
Yassin
F.
Razavi
S.
Wheater
H.
Sapriza-Azuri
G.
Davison
B.
Pietroniro
A.
2017
Enhanced identification of a hydrologic model using streamflow and satellite water storage data: a multicriteria sensitivity analysis and optimization approach
.
Hydrological Processes
31
(
19
),
3320
3333
.
doi:10.1002/hyp.11267
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).