Prediction of snowmelt derived streamflow in a wetland dominated prairie basin

The Cold Regions Hydrological Modelling platform (CRHM) was used to create a prairie hydrological model for Smith Creek Research Basin ( ∼400 km2), eastcentral Saskatchewan, Canada. Physically based modules were sequentially linked in CRHM to simulate snow processes, frozen soils, variable contributing area and wetland storage and runoff generation. Five “representative basins” (RBs) were defined and each was divided into seven hydrological response units (HRUs): fallow, stubble, grassland, river channel, open water, woodland, and wetland. Model parameters were estimated using field survey data, LiDAR digital elevation model (DEM), SPOT 5 satellite imageries, stream network and wetland inventory GIS data. Model simulations were conducted for 2007/2008 and 2008/2009. No calibration was performed. The model performance in predicting snowpack, soil moisture and streamflow was evaluated against field observations. Root mean square differences (RMSD) between simulation and observations ranged from 1.7 to 25.2 mm and from 4.3 to 22.4 mm for the simulated snow accumulation in 2007/2008 and 2008/2009, respectively, with higher RMSD in grassland, river channel, and open water HRUs. Spring volumetric soil moisture was reasonably predicted compared to a point observation in a grassland area, with RMSD of 0.011 and 0.009 for 2008 and 2009 simulations, respectively. The model was able to capture the timing and magnitude of peak spring basin discharge, but it underestimated the cumulative volume of basin discharge by 32% and 56% in spring 2008 and 2009, respectively. The results suggest prediction of Canadian Prairie basin snow hydrology is possible with no calibration if physically based models are used with physically meaningful model parameters that are derived from high resolution geospatial data. Correspondence to: X. Fang (xif382@mail.usask.ca)


Introduction
The prairie region of Canada (the Prairies) lies in the southern part of provinces of Alberta, Saskatchewan, and Manitoba and is a portion of the vast Prairie Pothole Region of North America (Winter, 1989).The Canadian Prairies are characterized by relatively low precipitation especially in the southwest part and are highly subject to frequent and severe droughts (Nkemdirim and Weber, 1999;Fang and Pomeroy, 2007).Annual precipitation in the prairie region of Saskatchewan ranges from 300 to 400 mm (Pomeroy et al., 2007a), approximately one third of which occurs as snowfall (Gray and Landine, 1988).The Canadian Prairies are a cold region and exhibit typical cold region hydrology typified by continuous snowcover and frozen soils throughout the winter.Great variation in hydrology exists across the prairie region of Saskatchewan, with fairly well-drained, semi-arid basins in the southwest part and with numerous wetlands and lakes development in the sub-humid north central and eastern parts (Pomeroy et al., 2007a).
Important hydrological characteristics of the prairie region of Saskatchewan are long periods of winter (usually four to five months) with occasional mid-winter melts (common in the southwest and rare in the northeast) and a snowcover modified by wind redistribution and sublimation of blowing snow (Pomeroy et al., 1993).The blowing snow process is affected by the interaction of local topography and surficial vegetation cover with regional wind flow patterns (Pomeroy et al., 1993;Fang and Pomeroy, 2009).High surface runoff derives from spring snowmelt, which is 80% or more of annual local surface runoff (Gray and Landine, 1988), and occurs as a result of frozen mineral soils at the time and a relatively rapid release of water from melting snowpacks (Gray et al., 1985).Meltwater infiltration into frozen soils can be restricted, limited, and unlimited depending on soil infiltrability (Gray et al., 1985;Zhao and Gray, 1997).Deep soils are characterized by good water-retaining capacity and high X.Fang et al.: Prediction of snowmelt derived streamflow unfrozen infiltration rates (Elliott and Efetha, 1999).Most rainfall occurs in spring and early summer from large frontal systems and the most intense rainfall in summer is associated with convective storms over small areas (Gray, 1970).During summer, most rainfall is consumed by evapotranspiration (Armstrong et al., 2008).Evapotranspiration occurs quickly via wet surfaces such as water bodies, wetted plant canopies and wet soil surfaces and relatively slowly from unsaturated surfaces such as bare soils and plant stomata (Granger and Gray, 1989).
The Canadian Prairies are characterized by numerous small wetlands as known locally as "sloughs" or "potholes"; these depressions formed from previous glaciations of the landscape.The majority of the depressional wetlands do not naturally integrate to any natural drainage system (LaBaugh et al., 1998) and are often internally drained, forming closed basins (Hayashi et al., 2003); in normal hydrological conditions these basins are termed non-contributing areas (Godwin and Martin, 1975).These wetlands occasionally connect to one another during wet conditions through the "fill and spill" mechanism (van der Kamp and Hayashi, 2009).The water balance of these wetlands is influenced by redistribution of snow by wind from adjacent upland areas, precipitation, evapotranspiration, snowmelt runoff, groundwater exchange, and antecedent status of soil and depressional storage (Fang and Pomeroy, 2008;van der Kamp and Hayashi, 2009).Depending on the water balance, these wetlands vary from being shallow and seasonal to deep and permanent.The depressional wetlands are important hydrological elements as they have large storage capacities (Hayashi et al., 2003) which can regulate peak runoff.They are also valuable habitats for migratory waterfowl (Smith et al., 1964).However, hydrology of these wetlands is very sensitive to changes in air temperature, seasonal precipitation and other climatic variability (Poiani et al., 1995;Fang and Pomeroy, 2008;van der Kamp et al., 2008).Land use alteration in surrounding upland areas can produce noticeable impacts on snowpack trapped by wetland vegetation, surface runoff to wetlands, and wetland pond level (van der Kamp et al., 2003;Fang and Pomeroy, 2008).50 to 75% of the original Prairie wetlands have been filled, levelled, and drained since European settlement (Dahl and Johnson, 1991;Gleason and Euliss, 1998), which has been implicated as a cause for downstream flooding (Rannie, 1980;Hubbard and Linder, 1986).
Substantial efforts have been made to investigate hydrological processes governing prairie wetlands in terms of surface and subsurface hydrological processes, dynamics of wetland storage, and surface runoff (Woo and Rowsell, 1993;Hayashi et al., 1998;Berthold et al., 2004;Spence, 2007;van der Kamp and Hayashi, 2009).Hydrological modelling systems have been developed to focus on predicting water balance for large scale basins with considerable wetland storage (Vining, 2002;St. Laurent and Valeo, 2007;Wang et al., 2008), whereas physically based models integrating more cold regions hydrological processes have been assembled to simulate hydrological processes for the individual closed wetland basins (Su et al., 2000;Pomeroy et al., 2007b;Fang and Pomeroy, 2008).In light of the hydrological and ecological importance of prairie wetlands, the objectives of this paper are to: (1) develop a physically based hydrological model for a Canadian Prairie basin with large wetland areas; (2) derive all model parameters using field survey data, digital elevation model (DEM), satellite imageries, stream network and wetland inventory GIS data (no calibration); (3) evaluate the model performance in simulating winter snow accumulation, estimating spring soil moisture, and predicting basin streamflow.

Study site and field observations
The study was conducted in the Smith Creek Research Basin (SCRB), which is located between the Rural Municipalities of Churchbridge and Langenburg in the east-central Saskatchewan, Canada approximately 60 km southeast of the City of Yorkton shown in Fig. 1a.The SCRB was initially estimated to have a gross area of about 445 km 2 based on a Ducks Unlimited Canada (DUC) basin delineation shown in Fig. 1b.Agricultural cropland and pasture are the dominant land uses, with a considerable area left to natural wetlands, native grassland and deciduous woodland.Soil textures mainly consist of loam (Saskatchewan Soil Survey, 1991).The basin is characterized by low relief with elevations varying from 490 m above sea level near the basin outlet area at the south end to 548 m in the north end upland; slopes are gentle and range from 2 to 5%.The 30-year (1971The 30-year ( -2000) ) annual average air temperature at Yorkton Airport is 1.6 • C, with monthly means of −17.9 • C in January and +17.8 • C in July; the 30-year mean annual precipitation at Yorkton Airport is 450.9 mm, of which 106.4 mm occurs mostly as snow from November to April (Environment Canada, 2009).Frozen soils and wind redistribution of snow develop over the winter, and snowmelt and meltwater runoff normally occur in the early spring with the peak basin streamflow usually happening in the latter part of April.The spring snowmelt runoff is the main annual streamflow event in the basin and much of this runoff accumulates in the seasonal wetlands and roadside ditches.Many water control structures such as road culvert gates exist in the basin and are operated by local farmers to regulate the runoff in their cropland areas; the gates are closed during extremely high runoff periods (i.e. during fast snowmelts or intense rain storms) but remain open otherwise.
Instrumentation at SCRB consists of a streamflow gauge, main meteorological station, network of 10 rain gauge stations, and network of seven wetland water level transducers shown in Fig. 1b.The main meteorological station (SC-1) was set up in July 2007 and includes the measurements of air temperature, radiation (incoming short, long, outgoing short, and long-wave), relative humidity, wind speed and direction, soil moisture (0-40 cm), soil temperature (0-20 cm), snow depth, rainfall, and snowfall.Snowfall was corrected for wind-undercatch using the algorithm of MacDonald and Pomeroy (2007).These data were collected for two field seasons: 2007-2008 and 2008-2009.A stream depth gauge located at the basin outlet (05ME007) is operated by Water Survey of Canada at a site with a stable rating curve and has been used to estimate basin streamflow discharge since 1975.Field surveys of soil properties and vegetation were conducted in the fall of 2007 and 2008.Soil samples were collected from the 18 field transects located nearby the rain gauge and water level stations and were later used to determine the soil moisture and porosity.These transects were selected to represent characteristic basin land uses: summer fallow, grain stubble, grassland, woodland, wetland, and drainage channel.Vegetation height, type, and density were recorded from the same field transects.In addition, snow surveys were taken from the same field transects over the winter of 2007-2008 and 2008-2009.Each survey was comprised of 420 samples of snow depth and 102 samples of snow density; the depth and density were used to estimate the water equivalent of snowpack.

Cold regions hydrological modelling platform
The Cold Regions Hydrological Modelling platform (CRHM) was used to set up a prairie hydrological model to predict water balance for the SCRB.The development of CRHM involved many decades of hydrological research in the cold, semi-arid environment of the Canadian Prairies.CRHM is a state-of-the-art, physically based hydrological model which uses a modular, object-oriented structure (Pomeroy et al., 2007b).Within CRHM, component modules represent basin descriptions, observations or physically based algorithms for calculating hydrological processes, including redistribution of snow by wind, snowmelt, infiltration, evaporation, soil moisture balance, and runoff routing.
These processes are simulated on landscape units called hydrological response units (HRU).HRUs are defined as spatial units of mass and energy balance calculation corresponding to biophysical landscape units, within which processes and states are represented by single sets of parameters, state variables, and fluxes.HRUs in the Prairies typically correspond to agricultural fields, grassland, forest woodland, and bodies of water (Fang and Pomeroy, 2008).CRHM has shown good simulations in mountain basins (Dornes et al., 2008), boreal forest and arctic basins (Pomeroy et al., 2007b), a semi-arid, well-drained prairie basin (Fang and Pomeroy, 2007), and a wetland prairie basin (Fang and Pomeroy, 2008).
A set of physically based modules was assembled in a sequential fashion to simulate the hydrological processes relevant to the SCRB (Fig. 2).The key modules include the radiation model of Garnier and Ohmura (1970), Prairie Blowing Snow Model (Pomeroy and Li, 2000), albedo model of Gray and Landine (1987), Energy-Budget Snowmelt Model (Gray and Landine, 1988), Gray's expression for snowmelt infiltration (Gray et al., 1985), Green-Ampt infiltration model (Ogden and Saghafian, 1997), Granger and Gray's (1989) www.hydrol-earth-syst-sci.net/14/1/2010/ Hydrol.Earth Syst.Sci., 14, 1-16, 2010  unsaturated surface actual evaporation model, Priestley and Taylor's (1972) evaporation expression for wetlands, and a Muskingum streamflow routing model (Chow, 1964).A new wetland module was developed by modifying a soil moisture balance model, which calculates soil moisture balance and drainage (Dornes et al., 2008) to include depressional storage and pond surface water storage.This model was modified from an original soil moisture balance routine developed by Leavesley et al. (1983).The changes are to make this algorithm more consistent with what is known about prairie water storage and drainage (Pomeroy et al., 2007a).A flowchart of this module is shown in Fig. 3.The soil moisture balance model divides the soil column into two layers; the top layer is called the recharge zone.Inputs to the soil column layers are derived from infiltration of both snowmelt and rainfall.Evaporation only occurs from the recharge zone, and water for transpiration is taken out of the entire soil column.Excess water from both soil column layers satisfies groundwater flow requirements before being discharged to subsurface flow which represents flow in macropores that occurs in cracking clay and very coarse soils.Two components, depression and wetland pond, were added to the soil moisture balance model to simulate wetland drainage.Depressional storage represents small scale (sub-HRU) transient water storage on the surface of upland agri-cultural fields, pastures and woodlands.Wetland pond storage represents water storage that dominates a HRU in wet to moderate conditions, though the pond can be permitted to dry up in drought conditions.The inputs to depressional storage are from surface runoff and overland flow after the soil column is saturated.After the depressional storage is filled, overland flow is generated via the fill-and-spill process (Spence and Hosler, 2007), in which over-topping of the depression results in runoff but minimal leakage of water from the depression to sub-surface storage is permitted before it overtops.Evaporation is permitted from depressional storage.Wetland pond storage works in a similar manner to depressional storage, except that the pond area does not have a soil column, and inputs are derived from upland surface runoff and subsurface lateral unsaturated flow fed by infiltration.

Model parameterisation
A pre-processing procedure was taken to estimate the values of model parameters.The procedure was essentially a model parameterisation based on field observations, lookup table values, and stream network and wetland inventory datasets.The parameterisation procedure also employed land use classification using satellite images and automated GIS procedure using high spatial resolution DEM.

Basin physiographic parameters
For modelling large basins such as SCRB, CRHM has a new "representative basin" (RB) feature, in which a set of physically based modules are assembled with a number of HRUs to represent a sub-basin.The RB can be repeated as necessary in a basin, with each sub-basin having the same modules but differing parameter sets as needed.Streamflow output from a number of RBs is then routed along the main stream through lakes, wetlands and channel.The SCRB was divided into five sub-basins that are represented by five RBs (Fig. 4).An automated basin delineation technique, "TOPAZ" (Garbrecht andMartz, 1993, 1997) was used to extract the sub-basin (Fig. 5a).A 1-m LiDAR DEM was resampled to 50-m to provide a more computational efficient input for the TOPAZ program.TOPAZ channel and sub-basin segments were generated, and the sub-basin segments were aggregated to five sub-basins illustrated in Fig. 5a.Within each RB, seven hydrological responses units (HRUs) were derived from the supervised land use classification based on two SPOT 5 10-m multispectral images that were acquired on 5 July 2007 and 1 October 2008 (Fig. 5b).The summer image was used mainly for separating vegetation and nonvegetation features, while the fall image was used to separate cropland and natural vegetation.Areas for fallow, stubble, grassland, open water, woodland, and wetland HRUs were determined from SPOT 5 land use classification; areas for river channel HRU was estimated from Ducks Unlimited Canada (DUC) drainage network GIS data.The average elevation for HRU at different sub-basins was determined from DEM and HRU classification.The latitude for the basin is the geographic centre of SCRB and was measured from GPS.The average ground slope of HRU was approximated from the reported slope values in Saskatchewan Soil Survey (1991).

Blowing snow and frozen soil parameters
Blowing snow fetch distance is the upwind distance without disruption to the flow of snow.A computer program "FetchR" (Lapen and Martz, 1993) was used to estimate the fetch for the large exposed areas (i.e.fallow, stubble and grassland HRUs) from the DEM and vegetation classification, resulting in fetches of 1000 m, 1000 m, and 500 m respectively.For river channel, open water, woodland and wetland HRUs, a 300 m fetch length was assigned.The vegetation height, stalk density and stalk diameter were calcuwww.hydrol-earth-syst-sci.net/14/1/2010/ Hydrol.Earth Syst.Sci., 14, 1-16, 2010 lated based on vegetation survey measurements.The distribution factor parameterizes the proportional allocation of blowing snow transport from aerodynamically smoother (or windier) HRU to aerodynamically rougher (or calmer) ones and was decided according to observed prairie landscape aerodynamic sequencing to favour deposition in wetland and river channel HRUs (Fang and Pomeroy, 2009).A frozen soil infiltration parameter, initial fall soil saturation, was determined from the soil porosity and volumetric fall soil moisture.The soil porosity was estimated from soil texture, which is predominately loam in the basin.Volumetric fall soil moisture was approximated from gravimetric measurement of soil survey samples.

Wetland and soil module parameters
For the soil column, the maximum water holding capacity was determined from multiplying the rooting zone depth by soil porosity; the initial value of available water in the soil column was estimated by multiplying the maximum water holding capacity by volumetric fall soil moisture content.The soil recharge layer is the shallow top layer of the soil column, approximately 60 mm; the initial value of available water in the soil recharge layer was determined by the product of the maximum water holding capacity and volumetric fall soil moisture content.It should be noted that the model treats river channel, open water, and wetland HRUs as having no soil column, and sustaining permanent surface ponding.Subsurface and groundwater drainage factors control the rate of flow in the subsurface and groundwater domains; these rates are slow in the prairie environment (Hayashi et al., 1998) and were estimated from the saturated hydraulic conductivity based on soil texture.
An automated procedure involving LiDAR DEM and various ArcGIS tools was used to extract initial depth, area and volume of surface depression which were in turn input into a depth-area-volume relationship, yielding final depth, area and volume of surface depression.The basin LiDAR DEM was resampled from its original 1-m spatial resolution to 10-m, and a "fill pits" ArcGIS algorithm by Martz and Garbrecht (1998) was used to created a depressionless DEM from the 10-m LiDAR DEM; both were used as inputs in the ArcGIS 3D spatial analyst "cut/fill"."Cut/fill" detects changes in the area and volume of a surface between two times due to addition or removal of material.If a surface is characterized as "cut" from erosion, it is categorized into "net loss", and if a surface is identified as "fill" from deposition, then it is regarded as "net gain", and "unchanged" is another category if there is no change on the area and volume of a surface.Using both the original DEM and the depressionless Hydrol.Earth Syst.Sci., 14, 1-16, 2010 www.hydrol-earth-syst-sci.net/14/1/2010/ X. Fang et al.: Prediction of snowmelt derived streamflow 7 DEM in the "cut/fill" created a virtual surface during two periods and generated only one category "net gain".The DUC sub-basin wetland GIS inventory and the basin "cut/fill" surface depressions were input in the ArcGIS "intersect" tool, producing area and volume of sub-basin "cut/fill" surface depressions in the wetland area.The sub-basin supervised land use classification and the basin "cut/fill" surface depressions were together input in the ArcGIS "intersect" tool, generating the area and volume of sub-basin "cut/fill" depressions for each land use, which were then filtered by DUC sub-basin wetland GIS inventory to "Erase" the wetland portion.The final results were area and volume of sub-basin "cut/fill" surface depressions in the upland area.The volume of "cut/fill" surface depressions (V 3-D cut/fill [m 3 ]) results from the product of depth (d 3-D cut/fill [m]) and area (A 3-D cut/fill [m 2 ]), thus the depth of "cut/fill" surface depressions was calculated based on Eq. ( 1): Then, a simplified depth-area-volume relationship (Brooks and Hayashi, 2002) was used to calculate the maximum surface depression volume (V max [m 3 ]) according to Eq. ( 2): where A max [m 2 ] and d max [m] are the maximum surface area and depth of depressions, respectively, and p[−] is the shape coefficient of depressions.Rearranging the Eq. ( 2), the maximum surface depression storage sd max [mm] was estimated based on Eq. (3): where d max is estimated from the depth of "cut/fill" surface depressions d 3-D cut/fill calculated by Eq. ( 1).d 3-D cut/fill was assumed to be the maximum for the depressions in the upland area, but was adjusted for the depressions in the wetland area due to the inability of the LiDAR signal to penetrate water stored in the permanent wetland.The average fall depth from the monitored wetlands shown in Fig. 1b was added to get d max in the wetland area.A 3-D cut/fill was assumed to be the maximum.The shape coefficient p varied with area of each wetland; for the wetland smaller than 10 000 m 2 , p=1.72 was used, the average value estimated from Smith Creek wetland volume analysis (Minke et al., 2010).Values of p, 3.3 and 6, as discussed by Hayashi and van der Kamp (2000) were used for medium size wetlands (i.e. 10 000 m 2 ≤A max ≤ 100 000 m 2 ) and large size wetlands (i.e.A max >100 000 m 2 ), respectively.The maximum surface depression storage in the wetland and upland areas was determined from average value of individual "cut/fill" surface depression storage in these areas using Eq. ( 3).The maximum storage of river channel HRU was estimated from the DUC drainage networks GIS data assuming that the channel has parabolic cross-section.For the river channel, open water and wetland HRUs, the initial surface depression storage was approximated by the product of the maximum storage and the average percentage of fall storage capacity of the monitored wetlands.The initial surface depression storage for the upland area was set as zero due to its ephemeral nature of storage and typical dry antecedent condition in the fall.The estimated values of both initial and maximum surface depression storage are shown in Table 1.

Routing parameters
For the routing sequence within RBs (Fig. 6a), runoff in the upland area of fallow, stubble, and grassland is routed to the upland woodland, and then is routed to wetland, open water, and river channel.Runoff from the wetland is accumulated in the open water, which connects to the river channel.This routing sequence represents a characteristic surface runoff flow pattern on the flat prairie landscape and describes the runoff sequence observed throughout Smith Creek basin during spring snowmelt runoff.Similarly, Su et al. (2000) had to modify the original routing sequence (i.e.uplands routing to channels and then out of basin) in the SLURP model to a sequence that was appropriate for a prairie basin (i.e.uplands routing to wetlands then to channels and finally out of basin).The routing sequence between RBs (Fig. 6b) was determined from the direction of flow between the sub-basins of Smith Creek along the main channel of Smith Creek as controlled by topography.
Muskingum routing module (Chow, 1964) was used for both routing within and between RBs.For the routing within RBs, the routing length is the distance from each HRU to the main channel; for the routing between RBs, the routing length is the main channel length in each sub-basin, and both types of routing length were estimated from DUC drainage networks GIS data.Manning's equation (Chow, 1959) was used to calculate the average flow velocity; the parameters used in the equation include hydraulic radius, longitudinal friction slope, and Manning's roughness coefficient.Hydraulic radius was determined from flow depth based on the channel shape.Longitudinal friction slope was calculated from the average change in elevation over a routing length using the DEM and DUC drainage networks GIS data.Manning's roughness coefficient was estimated based on the channel's condition.From the average flow velocity and routing length, the storage constant was estimated.The dimensionless weighting factor controls the level of attenuation, ranging from 0 (maximum attenuation) to 0.5 (no attenuation), and can be approximated by a variety of methods (Wu et al., 1985;Kshirsagar et al., 1995).However, due to lack of information for the approximation, medium value, 0.25 was assumed for the basin.
www.hydrol-earth-syst-sci.net/14/1/2010/ Hydrol.Earth Syst.Sci., 14, 1-16, 2010  A weighted routing distribution parameter is used to partition the amount of runoff between HRUs and the values were determined from a modified Hack's law length-area relationship (Granger et al., 2002).The parameter is multiplied times the outflow from each HRU to distribute this outflow as inflow to the downstream HRU.For each non-river channel HRU, the land use polygons from the supervised classification were used to extract total polygon area and the longest linear length within the polygon.The extracted area and longest length were graphed on a log-log plot to generate the modified Hack's law length-area relationship: where L (km) is Hack's law length for each HRU and A (km 2 ) is total area for each HRU.For the river channel HRU, the original Hack's law length-area relationship (Hack, 1957) was used: where n is number of samples, X o , and X s are the observed and simulated values, respectively.The RMSD is a weighted measure of the difference between observation and simulation and has the same units as the observed and simulated values.The MB indicates the ability of model to reproduce the water balance; a positive value or a negative value of MB implies model overprediction or underprediction, respectively.

Winter snowpack prediction and comparison
The simulations of snow accumulation (snow water equivalent or SWE) during the February-April of 2008 and 2009 were evaluated against observations.For the 2008 simulation period, three comparisons during the pre-melt period: 7 February, 28 February, and 20 March and four comparisons during the melt period: 11-14 April were conducted.Three comparisons during the pre-melt period: 5 February, 3 March, and 20 March and four comparisons during the melt period: 3-9 April were carried out for the 2009 simulation period.Table 2 shows RMSD for SWE simulations in all five subbasins.For the 2008 simulation period, RMSD ranged from 1.7 to 7.9 mm for fallow, stubble, open water, and woodland HRUs, indicating generally good performance; larger RMSD were found for grassland, river channel, and wetland HRUs, ranging from 7.1 to 25.2 mm.For the 2009 simulation period, RMSD for fallow, stubble, grassland, and woodland HRUs ranged from 4.3 to 8.6 mm, while greater RMSD ranging from 7.8 to 22.4 mm were for river channel, open water, and wetland HRUs.For both simulation periods, there was a wide range of differences but generally small when comparing the RMSD of the same HRU in different sub-basins.The woodland and wetland HRUs respectively had the smallest and largest range.The woodland HRU had RMSD values from 2.7 to 3.1 mm and from 8.3 to 8.4 mm in the 2008 and 2009 simulation periods, respectively.RMSD for the wetland HRU ranged from 7.1 to 25.2 mm and from 8.4 to 22.4 mm in the 2008 and 2009 simulation periods, respectively.In addition, the river channel HRU was found to have the second largest range of RMSD, from 10.0 to 17.4 mm and from 9.5 to 17.9 mm in the 2008 and 2009 simulation periods.Interestingly, the wetland and river channel HRUs with the largest and second largest range of RMSD also obtained high RMSD values, implying model's insufficiency in simulating the SWE for both HRUs.Nevertheless, the model simulated the general sequence of wind redistribution of snow fairly well; snow was relocated from fallow and stubble fields with the pre-melt SWE ranging 30 to 75 mm (Figs.7a, b, 8a  and b) to river channels and wetlands having the pre-melt SWE from 70 to 220 mm (Figs.7d, g, 8d and g).

Spring soil moisture prediction and comparison
After a 12.6 mm rainfall occurred on 22 March 2009, ice layer formation in the cropland, grassland, and shrubby wetland areas was noticed.The snowmelt infiltration into soils was restricted with the ice layer forming above soils, and the initial fall moisture status of soil matrix was no longer valid in this case.To cope with this, the initial fall soil saturation of 2008 for fallow, stubble, and grassland HRUs was adjusted to 80% from their original measured values, and the www.hydrol-earth-syst-sci.net/14/1/2010/ Hydrol.Earth Syst.Sci., 14, 1-16, 2010  wetland HRU was set to the restricted case where no infiltration is permitted.With this adjustment, the predicted volumetric spring soil moisture from 14 April to 8 May in both 2008 and 2009 was tested against the observations from the main weather station (Fig. 9).Earlier observations cannot be used because of partially frozen soil.The simulated values were somewhat higher than observed in the 2008 simulation and somewhat lower than observed in the 2009 simulation period (Fig. 9).RMSD was 0.011 and 0.009 for the simulations of spring soil moisture in the 2008 and 2009 simulation periods, respectively.This indicates on average, that the difference between the observed and simulated volumetric soil moisture was between 1.1% and 0.9%.

Spring streamflow prediction and comparison
Spring streamflow was simulated for both 2008 and 2009, and the predicted daily mean basin discharge was compared to the observations for both simulation periods (Fig. 10).
For the 2008 simulation period, the simulation showed good timing for estimating the peak daily discharge (Fig. 10a); the peak daily discharge was two days late compared to the observed one.The observed peak daily discharge was 4.65 m 3 s −1 , which is very comparable to the predicted value (4.68 m 3 s −1 ) shown in Table 2. On average, relatively small differences between the observed daily discharge and the simulation were found; Table 2 shows that RMSD was 0.12 m 3 s −1 for the model simulation.Furthermore, the simulation predicted 27 days of spring streamflow, which is three days shorter than the observed streamflow duration.MB listed in Table 2 for the simulation was −0.32, suggesting that the model underestimated the cumulative basin discharge volume by 32%.
For the 2009 simulation period, the model predicted the same timing for peak daily discharge as the observation (Fig. 10b).The simulated peak discharge was 6.29 m 3 s −1 , which is quite similar to the observed value (i.e.6.22 m 3 s −1 ; Table 2).RMSD was 0.33 m 3 s −1 for the simulated spring basin daily mean discharge, indicating that on average, the difference between the observation and simulation was slightly higher for the 2009 simulation period compared to the 2008 simulation period.The simulated duration of spring streamflow was 20 days shorter than the observed one.For the cumulative basin spring discharge, the simulated volume was lower by 56% when compared to the observation.

Discussion
The Cold Regions Hydrological Model platform (CRHM) was used to simulate the streamflow generated from snowmelt runoff for a large wetland dominated prairie basin, Smith Creek Research Basin (∼400 km 2 ).Compared to other modelling efforts using CRHM for small prairie basins (Fang andPomeroy, 2007, 2008;Pomeroy et al., 2007b), this study is the first attempt to making prediction for such a large prairie basin using CRHM.The model showed encouraging simulations of various components of the Prairie water balance.The predictions of winter snow accumulation were very similar and compared quite well with most of the distributed field observations.The simulations were able to effectively describe the prairie blowing snow sequence (Fang and Pomeroy, 2009) and relocate snow from "source" areas (e.g.fallow and stubble fields) and deposit to "sink" or "drift" areas (e.g.tall vegetated wetland area and deeply incised channels).This is a vital process in governing the water balance of prairie basins as the majority of water in the wetlands and prairie river channels has been shown previously to be the result of the redistribution of snow by wind (Fang andPomeroy, 2008, 2009) and subsequent snowmelt runoff  (Gray and Landine, 1988;Pomeroy et al., 2007a).However, large differences in the snow accumulation between the simulation and observation at the start of winter season existed for the river channel HRU (Figs. 7d and 8d).The predicted snow accumulation in some HRUs (e.g.river channel and wetland HRUs in both simulation periods, grassland HRU in the 2008 simulation period, and open water HRU in the 2009 simulation period) had relatively large discrepancies compared to the observations (Table 2).This is attributed to the HRUs setup in CRHM for each sub-basin or RBs; the determination of the seven HRUs was based upon the supervised land use classification using SPOT 5 satellite images and DUC drainage networks.The derived HRUs were the simplest way to represent land use for the prairie basin and essentially stratified the basin into land use units.This type of HRU setup strategy is similar to the prairie basin stratification technique discussed by Steppuhn and Dyck (1974).The strengths of this HRU setup strategy are that it is easy to set up model parameters and reduce computational time, but this strategy can also cause the model to lose accuracy for estimating snowpack when compared to observed snowpack information from the actual land cover.It is certainly a challenge facing by CRHM to balance the complexity of HRU setup with model simulation accuracy, and further research is needed to resolve this.
Soil moisture prediction was quite adequate for most agricultural management purposes.The trend of predicted spring soil moisture generally matched the observations except for the period after 29 April 2008.The exact reason for that is not known.In addition, the simulated spring soil moisture in 2009 consistently lower compared to the observation.This is because the simple snowmelt infiltration expression (Gray et al., 1985) and Green-Ampt infiltration expression (Ogden and Saghafian, 1997) used in the model for estimating frozen soil moisture status cannot handle the formation of ice layers from rainfall in early spring 2009.A more sophisticated snowmelt infiltration expression (Zhao and Gray, 1997) capable of dealing with ice layer formation needs to be incorporated in the future model simulations.
The model was able to predict the timing and magnitude of the peak basin streamflow discharge derived from snowmelt (Fig. 10); this is quite encouraging considering no calibration was involved.However, there are inadequacies in the model simulations or spring basin streamflow hydrographs.One is the recession limb of the simulated hydrographs.There was a small peak simulated by the model around 6 May 2008 (Fig. 10a) after the simulated hydrograph levelled off for about a week; this peak was not observed at the hydrometric station located in the basin's outlet.This inadequacy is related to the HRUs setup in each sub-basin.Although the model incorporated a weighted hydrologic routing strategy to distribute surface runoff from several upland HRUs (e.HRUs, the semi-distributed nature of model structure is a potentially major source of error causing the simulated results.Each sub-basin had seven HRUs that were generated from the supervised land use classification using SPOT 5 satellite images and DUC drainage networks, and representing the subbasin with only seven HRUs might be a bit oversimplified, particularly for the wetlands.Wetlands in this basin are more diverse than a single classification can portray; for instance, they have different conditions (i.e.newly drained, established drained, and intact), which produce a wide range of storage volumes and result in differences in storing surface runoff.Aggregating these different types of wetland to a single wetland HRU in the current model setup is likely the cause of the second peak discharge shown in Fig. 10a.This suggests that to model a prairie basin with substantial wetland drainage development, more types of wetland representation are needed in CRHM and some type of wetland routing sequence should be incorporated.In addition, the recession limb of the simulated hydrograph in spring 2009 had generally good agreement with the observed hydrograph (Fig. 10b), but the cumulative volume of spring streamflow was underestimated by 56% with about 20 days shorter duration; the underestimation of the total basin streamflow volume was also the case for the 2008 simulation.The estimated surface depression storage may be a potential error responsible for the underestimation in total basin streamflow.All depression storage was presumed to be retained by the wetland, however many of these depressions have been artificially drained by small channels that would not be evident using the DEM analysis procedures used here to determine depression storage.
There is no way to assess the simulated surface depression storage relative to the observed because wetland levels were not monitored during the snowmelt runoff period.However, drainage of some depressions would increase the discharge volume and the recession limb of the hydrograph.The next phase of modelling for this type of watershed should contain a "drained wetland" HRU that permits discharge from depression storage.Apportionment of depression storage between drained and undrained wetlands will require further topographic analysis that is beyond the scope of this initial study.
This study demonstrated a model parameterization procedure utilizing high spatial resolution LiDAR DEM, SPOT 5 satellite images, various geospatial data such as stream network and wetland inventory GIS data.The purpose was to involve automated basin parameters delineation techniques and simplified wetland depth-area-volume calculation in order to eliminate the need for parameter calibration.Through this procedure, basin physiographic parameters such as basin area and elevation and important hydrological process parameter such as blowing snow fetch distance, wetland surface depression storage, and surface runoff and channel flow routing parameters were derived successfully.Using these parameters, the water balance for a prairie basin dominated by wetlands was reasonably simulated.This modelling pro-cedure with no calibration emphasised the use of physically based models for modelling basin hydrology and can be applied to ungauged prairie basins if sufficient meteorology, basin land use, and physiography data are available.

Conclusions
The Canadian Prairie pothole region is characterized by numerous post-glacial surface depressions.These surface depressions form wetlands which are important factor in controlling the water balance in prairie basins.The ability of wetlands to trap blowing snow in winter and store runoff water is a crucial feature of the hydrology, and this poses a substantial challenge to hydrological modelling.A new wetland module dealing with wetland water storage was created in the Cold Regions Hydrological Model platform (CRHM) to predict the spring water balance in Smith Creek Research Basin.Results show that the model was capable of simulating wind redistribution of snow and snowmelt, updating frozen soil moisture content, and predicting spring basin streamflow.The model presumed that all wetlands were undrained and retained their depression storage below some spill threshold.This assumption likely resulted in underestimation of discharge volumes and hydrograph recession at the basin scale.Further modelling in this region should involve HRUs that describe artificially drained wetlands.Nevertheless, this study proposed an innovative process to derive model parameters using field survey data, high spatial resolution LiDAR DEM, SPOT 5 satellite images, stream network and wetland inventory GIS data.This model parameterization process can be useful for modelling ungauged basins if high resolution information on basin characteristics is available.
Fig. 1.(a) Extent of the semi-arid glaciated northern prairie wetland region (grey shaded area) in Canada and the United States (Winter, 1989) and the location of Smith Creek Research Basin (SCRB), and (b) extent of the SC and field observation stations of rainfall (SCR), water level (LR), hydrometeorology (SC) and streamflow (SG).

Fig. 3 .
Fig. 3. Flowchart of a wetland module of soil moisture balance calculation with wetland or depression storage and fill-and-spill.

Fig. 4 .
Fig. 4. CRHM modelling structure.Five Sub-basins are simulated by modelling structure "Representative Basin" (RB); same seven hydrological response units (HRUs) exist in each RB.Modelling structure of Muskingum routing connects all five RBs.

Fig. 5 .
Fig. 5. Basin pre-processing procedures.(a) Smith Creek sub-basin generation from the TOPAZ basin delineation program using 50-m resampled LiDAR DEM and (b) Smith Creek HRU generation from Ducks Unlimited Canada drainage networks and supervised land use classification using SPOT 5 10-m multispectral images.
Figures 7 and 8 show the comparisons of the observed SWE and the simulated SWE for fallow, stubble, grassland, river channel, open water, woodland and wetland HRUs in sub-basin 1.For the 2008 and 2009 simulation periods, the model tended to overpredict the SWE on a fairly consistent basis.The predicted SWE generally matched the observations for most HRUs; except for fallow, stubble, grassland and open water HRUs during the melt period of 2008 and the fallow, stubble, and open water HRUs during the melt period of 2009.Figures 7d and 8d demonstrate large difference between the simulated and observed SWE for the river channel HRU on 7 February 2008 and 5 February 2009, after which the predicted SWE was generally in good agreement with the observations.

Fig. 9 .
Fig. 9. Comparisons of the observed and simulated volumetric spring soil moisture from the main weather station in the Smith Creek Research Basin.(a) 2008 simulation period and (b) 2009 simulation period.

Fig. 10 .
Fig. 10.Comparisons of the observed and simulated spring daily mean discharge in the Smith Creek Research Basin.(a) 2008 simulation period and (b) 2009 simulation period.

Table 1 .
Parameters of surface depression storage for the wetland module.The initial values inside parentheses are for the fall of 2008 and the initial values outside parentheses are for the fall of 2007.

Table 2 .
Evaluation of snowpack simulations with the root mean square difference (RMSD, mm SWE).

Table 3 .
Evaluation of simulating spring basin discharge with root mean square difference (RMSD, m 3 s −1 ), model bias (MB), peak discharge (m 3 s −1 ), and duration of discharge (day) in 2008 and 2009 simulation periods.Sim and Obs are simulation, and observation, respectively.