Improved modelling of a Prairie catchment using a progressive two-stage calibration strategy with in situ soil moisture and streamflow data

The Prairie region covers approximately 900,000 km² 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.

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.

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 km² at the gauging station 05HG002 (operated by the Water Survey of Canada), of which only 282 km² 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.

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).

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.

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.

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.

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.

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.

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).

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.

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 m³/m³ at 20 cm soil depth and 0.03 m³/m³ at 50 cm soil depth. In the model validation period, NSE values of 0.80, bias of +5% and RMSE of 0.03 m³/m³ 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 m³/m³ were achieved. Interestingly, the model performed better during validation than calibration for both 20 and 50 cm soil depths (see also Figure 4).

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 m³/m³ 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 m³/m³ 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 m³/m³ at 20 and 50 cm soil depths, respectively.

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.

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 m³/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.

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.

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).

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.