Abstract
Abstract
Subarctic wetlands span almost one-sixth of the Canadian wetlands and have been acknowledged as an important ecotone between arctic tundra and boreal forest. Their hydrological complexity is determined by the climatological and physiographical characteristics with regard to the spatiotemporal distribution of water resources. In this study, two semi-distributed hydrological models, SLURP (semi-distributed land use-based runoff processes) and WATFLOOD™, were employed to understand their effectiveness in modeling the hydrology of subarctic wetlands. Comparisons of their delineation approaches, formulations, parameters, and simulation results indicated that both models were capable of simulating the hydrological processes. However, differences were also observed. Besides their different segment delineation approaches, snowmelt and spring peak flows simulated by SLURP were 4–7 days earlier than those estimated by WATFLOOD because SLURP predefines snowmelt rates as variables, whereas WATFLOOD applies constants in the degree-day method. Due to the lack of considering the existence of permafrost and numerous seasonal ponds, both models tended to underestimate the spring peak flows. Evapotranspiration estimated by the Morton complementary relationship areal evapotranspiration method adopted in SLURP was lower than that calculated by WATFLOOD. Summer runoff only appeared during intense rainfall events, and its concentration was much faster in both models as compared with the observed records, which may be attributed to the variations of permafrost depth, soil water storage capacity, and seasonal pond levels. These findings will be helpful in improving the modeling quality of the two models and understanding the hydrologic features of subarctic wetlands.
Introduction
Semi-distributed hydrological models have been widely used as simplified, conceptual, and mathematical representations of the hydrologic cycle. Unlike lumped conceptual rainfall-runoff models which describe natural environment as a point without any spatial variability, or fully distributed models that always have excessive data demand in grid element scale, semi-distributed hydrological models fit between these two types of models on the basis of a trade-off between process resolution, data availability, and model accuracy. They provide a mechanism for evaluating the impacts of geographical conditions and climatic variations on water resources. The spatial heterogeneity is taken into account by employing various watershed delineation concepts, such as hydrological response unit (HRU) and grouped response unit (GRU) approaches. The HRU approach derives segments from the homogeneity of physiographic and geographic parameters. It has been widely applied with geographic information systems (GIS) and remotely sensed data to facilitate hydrological modeling (See et al., 1992; Greene and Cruise, 1995; Chen, 2007, 2008). Units with similar hydrological response can be grouped by overlying different layers of physiographic information that preserves the heterogeneity of the watershed (Flugel, 1995). Bongartz (2003) stated that the HRU approach can delineate scale-independent homogeneous modeling entities according to spatially distributed physiographical catchment properties. However, it is also sensitive toward variations in land-use change when generating associated runoff. Shaw et al. (2005) claimed that the HRU approach could be impractical for macro-scale basins due to the large amount of segmentation. Hydrological modeling of large basins usually requires distinct hydrological response areas to be grouped within segments. By subdividing a watershed into multiple grid segments, the GRU approach groups of all areas with the similar land cover such that a grid segment may contain a number of distinct GRUs, regardless of location. Hydrological response of each land cover within each grid segment is calculated and then weighted by the fractions of land cover areas (Leon et al., 2001). It allows minority land covers (e.g., marsh, bare soil) to be represented while permitting a large segment that efficiently addresses the problem of computational limitations (Shaw et al., 2005). Kite et al. (1994) applied the GRU approach to simulate the discharge from a 1,600,000 km2 watershed by delineating it into only six computational segments. Much discussion on their performance in modeling different watersheds has been continuing in the literature (Filoso et al., 2004; Thompson et al., 2004; Wu and Johnston, 2007; Schmalz et al., 2008; Hattermann et al., 2008).
Many studies have examined the hydrologic features of wetlands using the HRU and GRU approaches (Shaffer et al., 1999; Mansell et al., 2000; Chen et al., 2003; Chen et al., 2004; Price et al., 2005; St Laurent and Valeo, 2007; Hattermann et al., 2008). Nonetheless, such efforts in modeling Canadian subarctic wetland dominated regions have been limited. There are still knowledge gaps in understanding the feasibility and effectiveness in modeling the water cycle of subarctic wetlands and the interactions among hydrological processes, such as snowmelt and spring peak, permafrost and wetland water storage, precipitation intensity and runoff concentration. To better understand the hydrological processes within subarctic wetlands, in this study, the SLURP (semi-distributed land use-based runoff processes) model and the WATFLOOD model were targeted and tested in a Canadian subarctic wetland—the Deer River watershed in northern Manitoba. As two of the major Canadian hydrological models using the HRU and GRU approaches respectively, SLURP and WATFLOOD have been extensively applied in many basins ranging in size from a few hectares to million square kilometres (Haberlandt and Kite, 1998; Fassnacht et al., 1999; Su et al., 2000; Pietroniro et al., 2006; Shin and Kim, 2007; Dibike and Coulibaly, 2007; Armstrong and Martz, 2008; Chen et al., 2010; Jing et al., 2010; Li et al., 2011; Jing and Chen, 2011). However, they have not been parallelly applied in the subarctic region, and their performance in simulating wetland hydrology is still unknown. This study attempts to help fill the above gaps, understand how these two hydrological models, particularly their watershed delineation approaches, would fit the requirements of wetland simulation, and explore their strength and shortcomings in modeling subarctic wetlands.
Study Areas and Field Investigation
The study site—the Deer River watershed is located in north Hudson Bay Lowlands (HBL, 57°55′ N 94°46′ W), 70 km south of Churchill, Manitoba, Canada. It has a drainage area of 5,048 km2 with elevations ranging from 16 to 232 m, most of which is covered by tundra, shrub, deciduous, and evergreen forests. The Deer River, as one of the largest tributaries of the Churchill River, originates in the southwest and flows toward the northeast (Fig. 1). The watershed is a broad polygonized peat plateau that consists of high- and low-centered polygons. Volumetric water content of peat soil and hummock is maintained at 60% to 90%. Vegetation is predominantly lichen-heath, in which lichen coverage varies with sub-basins and ranges from 67% to 83% (Wessel and Rouse, 1994; Jing et al., 2009). Deciduous and evergreen forests dominate the headwaters of the Deer River and its downstream channels. The hummocky terrain prevails and consists of porous peat overlying a thick layer of mineral substrate. A great number of seasonal ponds stretch over, whereas continuous permafrost underlies the study area with an active depth of ∼1 m by late August. The primary soils are brunisolic static cryosol, brunisols, brunisolic turbic cryosol, and organo cryosol according to Mills et al. (1976). Mean annual precipitation at nearby Churchill is (1978–2007) 462 mm with maximum in summer (mean July=63 mm; mean August=74 mm; mean September=71 mm), and mean annual temperature is −6.5°C. Winters are long and cold with average temperature varying around −20°C from November to April inclusive. Winter runoff is extremely low and is sustained by groundwater discharge. Snowmelt starts in early or mid May and ceases by the beginning of June, which generates annual peak runoff. Soil water deficits are common in summer and fall due to the intensive evapotranspiration.

Location of the Deer River watershed and Chesnaye sub-basin (with monitoring stations).
Due to the limited accessibility of the Deer River watershed, a representative sub-basin in the lower reach of the Deer River, the Chesnaye sub-basin, was selected for modeling parameter investigation during 2006–2008 (Fig. 1). Its elevation slightly varies around 52 m with many seasonally connected ponds stretching over. Vegetation cover is mainly tundra and shrub with little coniferous forest along the streams. The Canada VIA railway goes through the sub-basin from north to south and, thus, provides transportation. A monitoring network of four stream gauging stations (5, 6, 7, 10) and one automated weather station (Rail Spur, 58°09′38″ N 94°08′35.4″ W) was operated during 2006–2008 (Jing, 2009). Data from the weather station were scanned by a Campbell Scientic data logger (model CR1000) and stored as hourly averages. Frost table and surface soil moisture at multiple transects (2, 4, 6 and 8 m away from the bank) of each stream gauging station were measured using steel pole and the SM200 soil moisture sensor. Streamflow was monitored at the gauging stations using HOBO® water pressure transducer and Sontek® ADV Flowtracker. Two helicopter trips were also carried out in June and October 2007 to collect information about vegetation coverage, topographic and hydrological conditions across the watershed, and, particularly, in the upper reach of the river.
The 3-year observations at Rail Spur showed a significantly proportional relationship between daily air temperature and evapotranspiration in summer; meanwhile, precipitation played an important role in elevating evapotranspiration with an average lag of 1 day. Soil layers became more saturated as getting closer to the stream channels. On account of the extraordinarily high hydraulic conductivity, locations close to streams received more infiltrated water from streams because the extent of this infiltration process was inversely proportional to distance. After the major recharge period during snowmelt, soil moisture contents kept declining throughout summer, which was mainly due to the intensive evapotranspiration. A reciprocal relationship between the permafrost table and its distance to stream channels was observed. Moreover, the declining trend of permafrost indicated that air temperature acted as the dominant factor of its fluctuation. Most of the small or moderate rainfall events during summer barely generated notable runoff due to the descending frost table, high soil porosity, and intensive evapotranspiration.
Methodology
The SLURP and WATFLOOD hydrological models
The SLURP model version 11 (Kite, 1997) is a semi-distributed, conceptually continuous hydrological model. This daily time-step model was originally developed for meso- and macroscale basins with intermediate complexity, which incorporates necessary physical processes without compromising the simplicity of calculation. It requires the target watershed to be delineated into multiple HRUs—aggregated simulation areas (ASAs) by topographic parameterization. Each ASA is further divided into subareas with different types of land cover based on vegetation, soil, and physiographical conditions. A vertical water balance module is sequentially applied to each type of land cover within each ASA to account for precipitation, canopy interception, snowmelt, infiltration, surface runoff, and groundwater outflow. Precipitation is intercepted by canopy, whereas the amount of precipitation that passes through the canopy is counted as snow or surface storage. Snow melt is calculated using either the degree-day method (Anderson, 1973) or the simplified energy budget method (Kustas et al., 1994). Snowmelt rate is interpolated (parabolically) between the rates in January and July. Infiltration is governed by the Philip formula (Philip, 1954). Three methods from Morton (1983), Spittlehouse (1989), and Granger (1995) are available for estimating evapotranspiration. Surface runoff appears when water remaining in the aerated soil layers after infiltration exceeds the depression storage capacity. Subsurface flow (interflow and groundwater flow) is simulated at a rate depending on subsurface storage and the water transfer coefficient. Manning's equation is applied to route surface and subsurface flow from each land cover to the nearest channel within each ASA. Runoff routing between ASAs is sequentially carried out by using either the hydrological storage techniques or the Muskingum-Cunge channel routing method. Accumulated flow from the outlet of each ASA is routed to the downstream ASA and finally to the outlet of the whole watershed.
WATFLOOD, as a semi-distributed and physically based hydrological model, has been widely used to forecast flood events or simulate watersheds without sacrificing the distributed features and computational efficiency (Kouwen, 2008). As differentiated from SLURP, WATFLOOD is constructed based on the concept of GRUs, which allows its application to large watersheds where similar vegetated areas within each grid segment are grouped as one land cover type for water balance calculation. WATFLOOD allows hourly time increment and simulates both vertical and horizontal water balance in the natural environment, including interception, infiltration, evapotranspiration, snow accumulation and ablation, interflow, recharge, baseflow, and overland and channel routing (Kouwen et al., 1993). Precipitation interception is computed by the Linsley's model (Linsley et al., 1949) with minor changes with regard to the interception evaporation (IET). Snowmelt is estimated by the temperature index routine that allows refreezing (Anderson, 1973) or the radiation temperature index model. In addition, snow sublimation is considered using an adjustable sublimation factor. Potential evapotranspiration (PET) is calculated using either the Priestley-Taylor (Priestly and Taylor, 1972) or the Hargreaves' equation (Hargraeves and Samani, 1982) based on input availability. Actual evapotranspiration (AET) is either presumed as the potential rate when soil moisture is at a level of saturation or reduced to a fraction of the potential rate if soil moisture is below the saturation point. Infiltration process is governed by the Philip formula (Philip, 1954) and the Darcy's law. A ground-water depletion module is applied to diminish base flow. Channel routing is accomplished using the storage routing technique, whereas base flow is calculated by a nonlinear storage-discharge function. Moreover, routing through lakes is available using a used-specified power or polynomial function.
Delineation approach
One of the fundamental differences between SLURP and WATFLOOD is the selection of criteria for watershed discretization. A key question in watershed delineation is to determine what constitutes a hydrological homogeneous area. Usually, answers include the uniformity of physiographic, geographic, and climatological conditions within the target watershed. However, the assumption of uniformity may be violated due to the attempt to maintain a manageable number of segments as the size of the watershed increases. SLURP adopts the HRU approach and delineates a watershed into multiple ASAs based on the homogeneity of physiographic and geographic parameters such that the response to climatic forces would be uniform (Fig. 2a). An HRU can also be arbitrarily defined by land cover, soil type, slope, stream gauges, or even elevation bands. Raster grids that have upstream drainage area larger than the predefined critical source area are determined as the channel network in SLURP. A modification to the channel network is further applied to eliminate any links shorter than the preset minimum source channel length values. A notable advantage of this concept is that its outputs are available in raster format and are able to be compared with satellite-based models. Each ASA is not homogenous and contains multiple land covers with independent flow routing calculations based on their mean distances to the nearest streams. SLURP is capable of recalculating downstream flow values based on transient internal system diversion. However, watershed delineation may not be accurate under some specific topographical conditions, such as board plains where elevation hardly varies. Moreover, it might become difficult when applied to macro-scale basins where a large amount of segmentation could occur in heterogeneous basins.

Contrastingly, WATFLOOD employs the GRU approach, which divides the predefined grid segments into a number of distinct GRUs based on their topographic and climatological similarity (Fig. 2b). This approach is ideal to couple to remotely sensed data and GIS. Each segment has identical digital elevation model grids and various land covers with determined ratios, because different types of land cover can have a wide range of response characteristics. The GRU approach calculates the response for each land cover within each segment using standard distributed modeling techniques and, subsequently, weights the response by land cover area fractions. The results from each land cover are then summed to determine the response of each gird segment and the basin. The GRU approach is significantly based on the percentage of each land cover within each grid segment rather than their location. In addition, it is less data-intensive and more efficient than the HRU approach over large spatial scales. These features allow the size of grid units to be determined by a size where routing time within the units are much shorter than the overall basin routing time. Another prominent advantage of the GRU approach is that it can incorporate necessary hydrological features without compromising the simplicity of computation and introducing uncertainties caused by watershed delineation. Additionally, increasing the size of grid segments can remarkably reduce computation and input parameterization effort. However, an inherent weakness of this concept is that heterogeneity may be lost, and only land covers differentiation within one grid could be derived. Moreover, determination of flow direction may also be inaccurate if grid size is too large.
Snowmelt
Snowmelt is defined as surface runoff produced from melting snow in spring and is usually an important fraction of the annual water cycle in subarctic wetlands. The rate of snowmelt is determined by the amount of energy available to the snowpack and is usually dominated by net radiation. Cold snowpack has negative energy balance, whereas warming results in the isothermal condition and additional energy starts the melt. The most commonly used snowmelt-runoff model is the degree-day method or temperature index routine, (Anderson, 1973) which has been successfully verified worldwide over a number of basins and climates. Due to the lack of net radiation data, it was selected for both models. Although this method has inadequate accuracy in simulating small-scale snowmelt rates, it does produce acceptable results at basin scale (Pohl et al., 2005). SLURP uses an increasing parabolic interpolation from values specified in January and July to determine the degree-day factor (Table 1). The later the date is, the greater melt rate. WATFLOOD models snow-free and snow-covered area separately. Although the temperature index routine may not be effective for hourly time increments as the melting process is mainly controlled by radiation (Rango and Martinec, 1995), recent studies reported good modeling results using this routine (Donald, 1992; Seglenieks, 1994; Hamlin, 1996). The degree-day factors of each land cover are calibrated as constants in WATFLOOD in reference to Kouwen (2008) and Stadnyk (2008) (i.e., 0.13–0.17 mm/day/°C). This difference, if under the same ambient conditions (e.g., date, time, temperature, and radiation) in the melting season (April-June), may result in more rapid snow depletion in SLURP.
Variables are defined in the Notation.
SLURP, semi-distributed land use-based runoff processes; HRU, hydrological response unit; GRU, grouped response unit.
Evapotranspiration
Evapotranspiration is defined as the sum of vaporization from both vegetations and open water bodies, and plant transpiration. It is affected by a number of factors, including vegetation growth stage and maturity, percentage of soil cover, soil moisture, solar radiation, humidity, air temperature, and wind speed. The Penman equation has been widely used to estimate PET, which specifies the maximum evapotranspiration capacity if there were abundant water available (Penman, 1948). However, PET is a theoretical concept and is unlikely to be achieved in most scenarios. Many hydrological models calculate AET from reducing PET, because some water will be lost due to infiltration and surface runoff.
In this study, the Morton complementary relationship areal evapotranspiration (CRAE) model was selected to calculate AET in SLURP (Morton, 1983) due to the lack of sunshine hours and soil heat flux data. As a modification of the Penman equation, it computes PET by replacing the wind function with a vapor transfer coefficient to solve the energy balance and aerodynamic equations. Wet-environment evaporation is derived from an empirical estimation with regard to saturation vapor pressure, temperature curve, and net radiation. AET is computed from a complementary relationship as shown in Table 1. This method has wide appeal, because it does not require many parameters that may affect evapotranspiration and recognizes that AET is often less than PET. However, soil moisture feedback is not included and AET needs to be set to zero if surface and soil water storage depletes. In addition, this method is not land cover dependent.
In the absence of hourly radiation data, the Hargreaves equation was applied in WATFLOOD to empirically estimate PET from daily air temperature (Table 1). Solar input (Ra) is simplified as a function of latitude and Julian day, whereas Ct is regressed from measurements of relative humidity made in the United States (Kouwen, 2008). It is a lumped empirical estimation that may not accurately reflect the real conditions depending on the location of the study area. Nonetheless, this equation consistently produces accurate evapotranspiration approximations in an increasing number of documented studies in spite of not considering dewpoint temperature, vapor pressure, and wind conditions (Hargraeves and Samani, 1982; Mohan, 1991). PET is then further reduced to AET by applying soil moisture (UZSI) and temperature (FPET2), vegetation (FTALL), and peat layer effect (ETP) coefficients when IET is equal to zero. If IET is greater than zero, AET is further adjusted based on the difference between values of PET and IET. This reduction approach appears to be promising, because it addresses the feedback of soil moisture and land cover effects. In spite of this, more specific coefficients may raise the challenges of data input and quality, because they have to be either measured or calibrated.
Modeling input
A 3-arc-second (90 m resolution) digital elevation model of the Deer River watershed was obtained from the National Map Seamless Server of U.S. Geological Survey (USGS, 2008). The study area was delineated into 206 ASAs in SLURP and 54 × 48 square grids in WATFLOOD. Land cover datasets (1 km resolution) were obtained from the Systeme Probatoire d'Observation dela Tarre earth observation satellite system (SPOT IMAGE and VITO, 2008) and reclassified into six land cover classes, including water, impervious, marsh, shrub, coniferous trees, and deciduous trees (Table 2) based on some previous studies in Canadian subarctic region (Stow et al., 2004; St Laurent and Valeo, 2007). Historical streamflow data (1978–1997) were obtained from Water Survey Canada and further compared with modeling outputs at the D. River N. Belcher station (58°0′48″ N 94°11′53″ W, ID: 06FD002, Fig. 1). Meteorological data (1978–1997) at the Churchill-A Climate Station (ID 5060600) were provided by Environment Canada.
Sensitivity analysis
To understand which parameters contribute most to the output, sensitivity analysis was conducted for SLURP and WATFLOOD by adjusting all the model parameters (one-factor-at-a-time) and evaluating their impacts on modeling results. The base values of these parameters were determined in reference to the field investigation in the Chesnaye sub-basin, manuals of both models, and other researchers' work in the HBL where the study area is located (Kite, 1997; Su et al., 2000; Metcalfe and Buttle, 2001; Kite, 2002; Thorne, 2004; Woo and Thorne, 2006). Based on the assumption of independent interrelationship, each parameter was individually adjusted by ±5%, ±15%, and ±30% while keeping all other ones at their initial values (adjustment of each parameter was made to all land covers at one time), and evaluated by the fluctuation of the 10-year (1978–1987) logarithmic Nash and Sutcliffe efficiency (NSEIn) and the correlation coefficient (R) to derive a limited number of influential ones. The 10-year NSEIn and R were calculated by the following equations:
where N is the number of days; Qi is the daily observed flow in 10 years (m3/s); Qm is the daily modeled flow in 10 years (m3/s); and Qaverage is the 10-year mean observed flow (m3/s). The NSEln describes model performance in the case of discharge simulations during low-flow periods. The closer it is to 1, the better the performance of the model. The R statistic describes the degree of colinearity between observed and modeled flow. It tends to the unity the more simulated flows have the similar dynamics as observations. The most sensitive parameters to be calibrated in SLURP are initial contents of snow store and slow store, maximum infiltration rate, retention constants for fast store and slow store, maximum capacity for fast store and slow store, precipitation factor, rain/snow division temperature, and snowmelt rates in January and July. Contrastingly, the most influential parameters in WATFLOOD include lower zone drainage function exponent, base temperature for snowmelt, porosity of the wetland or channel bank, reduction in soil evaporation due to tall vegetation, crude snow sublimation factor, melt factor, and manning's roughness coefficient.
Modeling validation
To minimize the difference between observed and simulated runoff rates and compensate for errors in measurement, a 10-year (1978–1987) calibration was conducted for the influential parameters in SLURP and WATFLOOD using their built-in optimization modules and manual adjustment. Both the NSEln criterion and the correlation coefficient (R) were used as statistical measures of the goodness of fit of the calibration results at the D. River N. Belcher station. Tables 3 and 4 summarize the final values of model parameters obtained after the calibration of SLURP and WATFLOOD, respectively. Modeling verification was performed for the next 10-year period (1988–1997) to ensure both models accurately reflect the hydrological features of the study area. Tables 5 and 6 show the models' outputs and performance measures using parameters declared in Tables 3 and 4.
AET, actual evapotranspiration.
Results and Discussion
The first year of a simulation period is an initialization period due to the uncertainty of the initial values for snow cover and water storages. This explains the low NSEIn value (8%) of WATFLOOD in the calibration year 1978 with regard to the lack of initial water storage settings. On the other hand, SLURP utilizes and optimizes specific parameters for initial water storage, which, therefore, enables a good fit between simulated and observed hydrographs in the start year.
SLURP and WATFLOOD allowed fair reproduction of runoff generation dynamics during both high- and low-flow periods. They were able to capture the subarctic wetland character of the discharge regime of the Deer River, governed by snow accumulation and low flow in winter, by typical snowmelt runoff peak in spring, and by intensive evapotranspiration as well as soil water deficit in summer. Application of SLURP to the study area in 1988 resulted in an NSEIn of 72% and an R of 0.92, whereas WATFLOOD produced an NSEIn of 48% and an R of 0.56. The year 1988 was selected as an example, because it was a median year from a meteorological perspective (Table 5). Of the deficiency in NSEln, the majority was attributed to simulation during spring snowmelt and summer, whereas the minority was contributed during fall and winter. Overall, SLURP had better performance than WATFLOOD, which might be attributed to the differences in snowmelt and evapotranspiration, as well as model complexities.
Snowmelt runoff
In subarctic wetlands, snowmelt runoff usually commences in May with an increase in runoff during the melting period corresponding to the saturation of organic soil layers and a sharp decline after the snow disappearance. Fluctuations of melt rate exist due to the varying rhythm of melt responding to the seasonal trend. Saturation begins 12–18 days after the start of snowmelt infiltration, indicating most of the snowmelt peaks occur in May or June and lead to pronounced flow maxima in streams. Van der Linden and Woo (2003) applied SLURP to a typical subarctic catchment, the Liard basin in northwest Canada, and reported spring peak flows in early June. Similar trends were observed by other studies conducted in subarctic wetlands (Woo et al., 2000; Carey and Woo, 2001).
Figure 3a shows the Julian days of observed and simulated snowmelt peaks at the D. River N. Belcher station. It can be concluded that most of the simulated peaks were 2–11 days later than observed records. This can probably be explained by the existence of permafrost that blocked the percolation of water and, therefore, increased observed runoffs. In spring, when the permafrost layer is shallow, meltwater could infiltrate the surface porous soil; however, infiltration is retarded by permafrost and a perched saturated zone is formed toward the ground surface to produce continuous lateral runoff (Carey and Woo, 2001). These effects are not considered in either SLURP or WATFLOOD, and meltwater is routed to interflow and groundwater instead of being discharged as instant runoff. Therefore, the modeling results implied a number of days' delay to convey meltwater from snow cover to surface runoff. Another interesting finding was that peaks modeled by SLURP were ∼8 days earlier than those modeled by WATFLOOD. Snowmelt base temperatures were calibrated in both models, and their values were close to each other. Instead of having preset melt rates as in WATFLOOD, SLURP employs parabolically increasing melt factors that might intensify the melt process and result in earlier peaks.

The simulation results also demonstrated that the amount of peak flows was 30% less than the historical records (Fig. 3b). Permafrost represents an over-winter surface storage of groundwater. When temperature increases and snow starts to melt, water is also released as stream runoff from the thaw of ice contents within the organic soil layers (Woo and Thorne, 2006). The numerous seasonal ponds may also contribute to melt peaks, because they have a vast water equivalent storage capacity that can be discharged when air temperature rises. In addition, from the hydrological point of view, permafrost and a discontinuity in the soil properties restrict water percolation and, therefore, excludes the effects of soil water storage. However, meltwater could be intercepted as surface and groundwater storage in both models, resulting in alleviated spring peaks compared with observed records.
Evapotranspiration
The average annual AET estimated by SLURP and WATFLOOD during the modeling period were 202.1 and 334.6 mm, respectively. Rouse (1998) estimated the water balance of a hummocky sedge fen in northern HBL and reported an annual AET rate of 269 mm. Figure 4 shows the temporal trends in observed annual runoff at the D. River N. Belcher station and cumulative AET modeled by the Morton CRAE method in SLURP, the Hargreaves equation WATFLOOD, and Rouse (1998). It can be concluded that SLURP may underestimate AET, whereas WATFLOOD may overestimate it. Summertime AET estimated by the Hargreaves equation was dramatically higher than that from the Morton CRAE method. Take August 1995 as an example, the monthly cumulative AET were 128 and 51 mm for WATFLOOD and SLURP, respectively. This difference may be attributed to a number of possible reasons. First, inputs for PET computations were not the same. SLURP uses a modified version of the Penman equation, whereas WATFLOOD only estimates by temperature, location on the earth, and the Julian day. Second, PET calculated by both models may be subject to inaccuracy due to the lack of applications and experiences in subarctic wetlands, such as the lack of soil wilting points, soil water contents, net radian flux, water depth in the canopy, and saturated vapor pressure. Third, the reduction rates of PET in both models were different. WATFLOOD appears to have a more realistic concept, because it employs soil moisture feedback and land-cover effects. However, these extra factors may cause arbitrary uncertainties during the calibration process.

Temporal trends in observed annual hydrograph at the D. River N. Belcher station and cumulative actual evapotranspiration modeled by SLURP, WATFLOOD, and Rouse (1998).
Rainfall-runoff response
Both historical records and model results showed that most light and moderate rainfall events in summer (July–September) were not able to generate notable runoff due to the combined effects from interception, depression and subsurface storage, permafrost, and evapotranspiration. Many lichens, mosses, small shrubs, and conifers flourish in the summer. They have considerable capacities to intercept a great amount of precipitation and this leads to non-negligible water loss for the drainage basin caused by evaporation. Depression storage is mainly referred to the numerous seasonally connected ponds that behave as runoff buffers due to their fluctuating water storage capacities. Highly porous soil and descending frost table allows more water to penetrate into deep soil layers in summer. These distinguished features of subarctic wetlands not only affect snowmelt runoff but also alleviate peak streamflows and extend runoff concentration after rainfall events. Precipitation is retained in soil layers and ponds for a longer period and are more likely to be released as evapotranspiration or groundwater flow rather than surface runoff. As the most significant natural means of cycling water back into the atmosphere, evapotranspiration also tends to be intensified due to higher air temperature and longer daylight period in the summer and, therefore, further reduces surface runoff. These combined factors resulted in the fact that only heavy or continuous rainfall events were able to generate countable runoff.
A high-intensity rainfall event brought large amount of water to the wetland in a short period. After the surface soil layer was saturated, excess water generated flashy runoff, resulting in quicker runoff response. To quantitatively assess the rainfall-runoff response in subarctic wetlands, 44 summer heavy storms (daily amount is more than or equal to 20 mm) that occurred during 1978–1997 were analyzed. Observed streamflow records at the D. River N. Belcher station indicated an average time lag of 3–6 days between the peaks of rainfall and runoff in summer. The modeling results showed that this delay was only 1–3 days, which could be attributed to some unique features of subarctic wetlands (Fig. 5). The variation of permafrost depth, soil water storage capacity, and seasonal pond level can prolong the runoff concentration. Lack of consideration of these buffering effects in both models resulted in those more rapid runoff response. The mean values and standard deviations of rainfall-runoff response in SLURP (1.44 and 0.74 days) and WATFLOOD (1.84 and 0.72 days) indicated that WATFLOOD has relatively longer runoff concentration time, which may be attributed to its depression storage subroutine.

Rainfall-runoff response of high intensity rainfall events in observed records, SLURP, and WATFLOOD at the D. River N. Belcher station (squares show mean values; bottom and top of the boxes describe the 25th and 75th percentiles; circles represent the minimum and maximum values; error bars indicate standard deviations; 36 rainfall events in total from 1978 to 1997 where cumulative daily precipitation exceeds 20 mm).
Conclusions
This research targeted the hydrological processes in subarctic wetlands to understand the water cycle and its unique hydrological features. Based on historical data support and extensive field investigation during 2006–2008 in the Deer River Watershed, Manitoba, two hydrological models, SLURP and WATFLOOD, were applied to simulate wetland hydrology in the typical subarctic region of the HBL. They were further compared from the aspects of delineation approaches, formulations, and parameters, as well as simulation performance, to evaluate their feasibility and effectiveness in modeling subarctic wetlands. Despite the limited applicability to large-scale basins, the HRU watershed delineation approach used in SLURP can recalculate downstream flow values based on transient internal system diversion with good computational efficiency. The GRU approach embedded in WATFLOOD is ideal to couple to remotely sensed data and GIS for large scale basins; however, the heterogeneity may be lost and only land covers differentiation within one grid could be derived. The degree-day snowmelt method is embedded in both models except that SLURP uses a date-varying snowmelt rate that may cause snow depletion to be accelerated. The modeling results indicated that the Morton CRAE method in SLURP had lower evapotranspiration outputs than those from WATFLOOD. Runoff of rain water during summer only occurred after occasional intense precipitation events. It also noted that the lack of considering the existence of permafrost and ponds as well as their interactions in both models appeared to be one of the main reasons for rainfall-runoff modeling inaccuracy. Modifications to both models, such as snowmelt algorithm, permafrost layer, and wetland routing, are recommended for both models to handle these uncertainties and further improve their performance in simulating subarctic wetlands. The findings of this study will be helpful in understanding the hydrologic features of subarctic wetlands and supporting wetland conservation, particularly under climatic change context.
Footnotes
Acknowledgments
Special thanks go to Manitoba Hydro, ArcticNet, NSERC, NFSC, and Churchill Northern Studies Center (CNSC) for funding this work. The advice and help from Dr. Kenneth Snelgrove at Memorial University of Newfoundland and Dr. Tim Papakyriakou at the University of Manitoba are also highly appreciated.
Author Disclosure Statement
The authors declare that no conflicting financial interests exist.
