Received 12 December 2015; accepted 3 April 2016; published 6 April 2016
Fundamental physical processes in hydrology such as rainfall, evaporation, transpiration, infiltration, surface and ground water flow and connectivity can be precisely described by physical laws at small spatial and temporal scales  . However, one of the unresolved key challenges in hydrology is how the laws can be applied at scales by far larger and temporally and spatially more heterogeneous than the small scale they have been proved valid   .
Over the last few decades, considerable studies have been made in using large scale numerical hydrologic models developed to understand environmental processes and predict environmental changes at various spatial and temporal scales  -  . Most of these studies have employed GIS-coupled physically based distributed hydrologic models. In these models, the spatial variability is taken into account by dividing the study area into several smaller elements or regular mesh pixels over which point scale laws or measurement data at point or local scales are upscaled to predict water dynamics at larger scales  . There is a general assumption that change in scales can be accommodated by using effective parameters that should be determined during calibration of these models. Examples are MIKE-SHE   and TOPKAPI  .
This grid parameterization approach, however, is considered as inadequate by some authors.  describes the use of effective parameter values as inadequate to the scaling problem. In addition, effective parameter based distributed models are heavily over parameterized and computationally not efficient, as they result in the application of a set of model parameters to each grid cell where the small scale physical laws are employed  . Hence, various researchers have proposed the concept of using homogeneous modeling entities, which are aggregated area of similar hydrologic behavior on the basis of hydrology, land use, soil, vegetation, slope and so on. Examples include the simple scaling and multiscaling framework  ; the Representing Element Area concept  ; hydrologically similar units or “patches”  ; the Hydrological Response Units concept of  ; and basin-scale model equations (e.g.  ). The hydrologically similar units approach as compared to the grid approach is believed to be more efficient computationally as specific set of models parameters is applicable to each type of entities rather than individual grids.
Despite these various attempts, a number of recent reviews have identified that scaling issues are at the heart of most hydrologic problems which are not yet resolved    -  .  noted that many of the better known scaling procedures in hydrologic modeling did neglect the important components of the complex hydrologic processes that we observe in the field. For example, the amount of information contained in hydrologically similar units like patches is highly dependent on the various assumptions considered during their formation. The way that these modeling entities have been formed or characterized may have impact on simulation results of the hydrologic regime.
This study, therefore, aims to investigate how the characterization of hydrologically similar units affects simulation output in hydrological models. The simulations are performed using the Regional Hydro-Ecologic Simulation Systems (RHESSys)  , in which “patches” represent the finest scale of the spatial hierarchy, i.e. the basic modeling entities on which energy and water components are simulated. These patches are formed by the intersection of each coarser level of the hierarchy and other attributes of the watershed  . It is important to understand the sensitivity of the model to patch characterization before applying the model for the purpose of investigating the effect of potential management practices in catchment’s hydrologic regimes, both water quantity and quality, investigated at the outlet of a catchment. Particularly, this study investigates the impact of patch characterization on simulated flow regime, and on the calibration of RHESSys’s main hydrological parameters: saturated hydraulic conductivity, k, and, decay of saturated hydraulic conductivity with depth, m.
2. Description of the Hydrological Model
RHESSys is a hydro-ecological, GIS-based model that is designed to simulate fluxes of water (evaporation, transpiration, stream flow, and soil water), nutrients (net photosynthesis, plant growth), and carbon in the landscape. The model uses the Soil-Plant-Atmosphere Continuum concept  within a distributed watershed context to represent the spatial distribution of hydrologic and ecophysiological processes and characteristics of a catchment  . RHESSys is composed of three process submodels: a meteorological model (based on MT- CLIM;  ), an eco-physiological model (based on BIOME-BGC;  ), and two distributed hydrologic models. TOPMODEL  is a “quasi”-distributed hydrologic model that distributes soil water according to the statistical distribution of the topographically defined soil wetness index of the patches. An adapted version of DHVSM  is a distributed hydrologic model, where saturated subsurface through flow and overland flow are explicitly routed between contiguous patches. The explicit routing, whilst computationally not as efficient as TOPMODEL, produces hydrologic patterns that more closely matched empirical data  . Therefore, for this study, the explicit hydrologic routing approach is used.
The RHESSys model is structured as a spatially nested hierarchical representation of the landscape (Figure 1), in which different hydrological, climatological, and ecosystem processes are partitioned. Different processes are modeled at each of these levels within the hierarchy. At the base of the hierarchy is the drainage network. Within the drainage network, levels in the hierarchy, range from climate zones (representing orographic effects on precipitation), hillslopes (representing horizontal water, carbon, and nitrogen transfers), and patches (representing vertical water, C, and N cycling). A detailed description of the parameters and processes at each level of the hierarchy is provided in  . The patches, which represent the finest scale of the spatial hierarchy, are formed by the intersection of each coarser level of the hierarchy and other attributes of the watershed  . Therefore the patches are homogeneous units in terms of climate, hillslope, soils, riparian areas, vegetation, and land use.
The RHESSys model requires a large number data inputs, including both spatial data, in the form of GIS- generated maps, and nonspatial data, in the form of tables (Figure 2). The maps, which are usually GIS generated, represent spatially variable properties, such as topography, vegetation, soil, and land use. The tables include temporally variable climate data (daily precipitation, maximum and minimum temperatures) and default values for physiological characteristics of the vegetation and physical characteristics of the soils.
Figure 1. The spatial hierarchy within RHESSys  .
Figure 2. Required inputs for the RHESSys model. Inputs to the model include spatial data (maps) generated within a GIS and non-spatial data (tables)  .
3. Study Area Description
This research was conducted on a subcatchment of the Turkey Lake Watershed (TLW). The TLW has an area of 10.5 km2 and located within sugar maple dominated forest of the Algoma Highlands of Central Ontario (47˚80'30''N, 84˚82'50''W; Figure 3), on the Precambrian shield approximately 60 km north of Sault Ste. Marie  .  noted that the climate within the TLW is continental and is strongly influenced by the proximity of Lake Superior, with mean annual precipitation and temperature of 1200 mm and 5.0˚C respectively, for the period 1981 to 2006. A snow pack typically persists from late-November, early December through to late-March, early April. Peak stream discharge occurs during snowmelt and again in October to November during autumn storms. Almost 35% of the average precipitation is in the form of snowfall and spring snowmelt accounts for 30% - 60% of annual runoff  .
Vegetation in TLW is relatively undisturbed and is predominantly mature to overmature tolerant hardwood forest (approximately 90% sugar maple, Acer saccharum Marsh.) that underwent a light selective harvesting (“high grading”) in the 1950s  . Details on the physical and biological characteristics of the watershed are given by  .
TLW has been used as an experimental basin since 1980 to study anthropogenic perturbations in Canadian Shield ecosystems in central Ontario  . It is divided into a number of subcatchments each having hydrometric stations at their outlet. Catchment (C38), which is used herein, is one of these subcatchments (Figure 3), with an area of 64,625 m2. Most of the area within the subcatchment has an elevation of nearly 405 m above sea level (Figure 4(a)). In some places the elevation can go as high as 450m that makes the total relief within the watershed to be approximately 45 m. The soil (Figure 4(b)) is predominantly silty loam (71%), followed by sandy loam (15%), and wetland area (14%). Outside of the wetland, which is located in the centre of the catchment, land use is entirely sugar maple forest (Figure 4(c)). The subcatchment outlet is identified as a weir point (47˚2'52''N, 84˚24'29''W) where continuous discharge measurements are available.
Figure 3. Location of subcatchment C38 within Ontario and the turkey lake watershed.
(a) (b)(c) (d)
Figure 4. Spatial properties of subcatchment C38. (a) DEM; (b) Soil distribution; (c) Land use and (d) Leaf-area index.
4. Model Setup
The stream flow at the outlet of C38 was modelled for eight patch configuration scenarios. The simulation results were compared to see the effect of patch resolution on the flow regime at the outlet of the catchment.
Spatial data were derived from terrain analysis of a 5m resolution Digital Elevation Model (DEM) obtained from the Ontario Base Map of the Ontario Ministry of Natural Resources (OMNR), and preprocessed using the Terrain Analysis System  . Topographic attributes, including both primary (slope, aspect, specific contributing area) and secondary (wetness index) attributes, were derived from established terrain analysis techniques   -  .
Two sets of stream networks were derived from the DEM, using specific contributing area thresholds of 2000 m2 and 5000 m2. The stream networks were buffered to consider the heterogeneity in topography of riparian area near the streams during the creation of patches. Similarly, two sets of hillslopes were derived from the DEM, also with thresholds of 2000 m2 and 5000 m2. Finally, two sets of climate zones were created. Since there is no statistical correlation denoting orographic/elevation effect on precipitation, climate zones were determined arbitrarily. As such, two climate zones were created, one having one climate zone and the other having two climate zones. The latter was created by reclassifying the processed DEM into two zones: areas with elevation from 405 m to 427.4 m classified in the first zone and areas with elevation between 427.4 m and 450 m classified in the second zone.
A digital vegetation map was adapted from the Forest Resource Inventory (FRI) surveyed by OMNS in 1995. Vegetation classes were amalgamated to include conifer, hardwood, poplar, and black spruce. Canopy Leaf Area Index (LAI) was derived from LANDSAT Thematic Mapper data using the methods of  . A detailed description of the physiological parameters for the vegetation classes is provided by  . Soil attributes were derived from data obtained from Forest Landscape Productivity Survey (FLaPS) maps created and updated by the OMNR. In RHESSys soils are classified based on their textural classes. The model contains default values for the hydraulic characteristics of most common soil textural classes.
Patches, which represent a unique homogeneous modeling units, are created by overlaying spatial layers of the climate zone, hill slope, vegetation, soils and buffered stream networks  . Since there is only one soil map, and one vegetation map, eight different patch configurations (PC1 to PC8) can be obtained by overlying these with different combinations of the climatic zones, hill slopes and stream networks (Table 1; Figure 5). The number of patches in each configuration depends on the overlay and varies between 93 and 168 (Table 1).
Figure 5. Maps of the eight patch characterization scenarios used in this study project.
Table 1. Range of values of input files for patch creation.
Temporal climate data (daily precipitation, daily maximum temperature, and daily minimum temperature) for a 5 year calibration period (June 1986 to July 1991) were obtained from a station close to the southeastern corner of TLW (Figure 3) where daily meteorological data around the watershed have continuously been measured since 1980.
For each of the eight patch configurations, the model was calibrated for optimal flow regime simulation. All calibrations were based on continuous flow measurement at the outlet of subcatchment, covering the entire simulation period (June 1986 to July 1991). Before calibration, the catchment hydrologic and nutrient condition was stabilized by spinning up of the model since the year 1600. The calibrations optimize the soil parameters that define soil transmissivity: the saturated hydraulic conductivity, k, and the decay of saturated hydraulic conductivity with depth, m. Model performance is measured with the Nash-Sutcliffe efficiency coefficient, E, which indicates the agreement between the observed and simulated discharges (E = 0 represents poor agreement and E = 1 represents complete agreement;  ). For each patch configuration, a Monte-Carlo calibration was applied, in which 75 different combinations of m and k were simulated.
5. Results and Discussion
The Nash-Sutcliffe efficiencies, resulting from the Monte-Carlo simulations for each of the eight patch configurations, are shown in Figure 6. For each of the patch configurations there appears to be a systematic trend in the relation between E and k, although the nature of that relation varies with the different patch configurations. On the other hand, there is no discernable trend in the relation between E and m. Thus, the hydraulic conductivity, k, has a stronger influence over E, than the decay in hydraulic conductivity with depth, m. Most patch configurations result in acceptable simulation performances (E > 0.5), for most combinations of k and m. Exceptions to this are PC2 (−1.4 < E < 0), PC6 (−1.0 < E < 0.2), which indicates that the combination of low stream network threshold and coarse hillslope thresholds is not suitable. In PC8, with low thresholds for both stream network and hillslope, most k and m combinations resulted in good simulation performance (E > 0.6). Optimal values for k and m vary for each patch category (Table 2), indicating that the calibration of these parameters is sensitive to the way the patches are defined. In general, the patch configurations with two climate zones (i.e. PC5, PC6, PC7 and PC8) have slightly better performance than the corresponding patch configuration with only one climate zone (i.e. PC1, PC2, PC3 and PC4, respectively). However, there is no generalized pattern in the performance of the patch configurations with different thresholds for hillslope and stream network definition, other than that the two configurations with a high thresholds for hillslope definition and a low threshold for stream network definition (PC2 and PC6) both perform poorly (E = −0.02 and E = 0.20, respectively). The total simulated flow and the relative difference to the total observed flow during the simulation period also reflect these observations (Table 2).
These quantitative results are a first indication that model simulation is sensitive to the way patches are defined. For each patch configuration, the optimal m and k values (Table 2) were selected for further analysis and flow simulation at the outlet of the subcatchment. Qualitative comparison of the observed and simulated daily
Figure 6. Maps showing distribution of the Nash Sutcliffe efficiency values that correspond to selected 75 k (a) and m (b) values generated during model calibration.
Table 2. Table that shows the best m and k combinations and the resulting total simulated flow simulated for the eight patch scenario projects.
specific discharge (Figure 7) shows that the model captures the observed discharge variations reasonably well, again with the exception of PC2 and PC6. However, the model seems to be overly responsive in 1986 and 1988, simulating higher peak flows than were observed during the summer periods. The model also under predicts the flow at the beginning of the simulation period, especially in scenarios with river network threshold (PC1, PC3, PC5 and PC7). Additionally, the spring runoff in the year 1988 was consistently underpredicted, which may be due to a delay in the melting of snow cover. Overall, however, the timing of flow events is well captured in all cases, with the exception of PC2 and PC6. Qualitative comparison of the yearly simulated and observed runoff, based on water years (from June to May in the following year), shows that the total annual flow is consistently underpredicted in most patch configurations, except PC2, PC6 and PC8 (Figure 8). Configurations PC2 and PC6 resulted in overprediction of the total yearly runoff for all the simulation years. Again, it is clear that patch scenario 8 has the best simulation result, with the least difference between the simulated and the observed annual average runoff values.
Figure 7. Daily time series curves of the observed and simulated flows at the outlet of the c38 watershed for the eight patch classification scenarios.
Figure 8. Comparison of the yearly total and average observed and simulated flows at the outlet of the c38 watershed for the eight patch classification scenarios.
Eight patch configuration scenarios for subcatchment 38 of the Turkey Lake Watershed were calibrated for hydrological flow simulation with the RHESSys model. Most patch configurations could, after calibration, be made to give reasonable estimation of observed flow patterns for the five-year water simulation period, with the exception of the patches defined by a high threshold for hillslope definition and a low threshold for stream network definition (PC2 and PC6). The best simulation result was obtained for patch configuration PC8, which was formed by combining layers of more classified climate zones and layers of stream network and hillslope derived using smaller threshold value. Similarly, the smaller threshold area specified during stream network delineation gave a relatively detailed stream network that in turn influenced the detail of hillslopes generated within the subcatchment.
Based on these results, it can be concluded that the RHESSys model flow simulation is sensitive to patch classification. Patches formed from finely aggregated patch forming inputs layers give better results as compared to patches formed from coarsely aggregated inputs. From the point of view of computation efficiency and memory, it may be very interesting to find the threshold level at which patch forming inputs should be aggregated so that flow simulation is not affected under this threshold.
The most notable result, however, is that the model’s k and m parameters are dependent on the patch configuration and vary notably between the different patch configurations (54.26 < k < 119.13; and 1.02 < m < 2.28). This dependency poses problems for the physical interpretation of these parameter values, and for their transferability to neighboring catchments.
Patch configuration in RHESSys modeling applications is often arbitrary or undocumented. This study, however, has shown that the result and the accuracy of flow simulations could be affected by the assumptions made during patch formation process.
The authors acknowledge Catchment Research Facility (CRF) of Western University, Canadian Forest Services of Natural Resources Canada, and National Water Research Institute and Meteorological Services of Environment Canada for providing the TLW datasets used in this study. We also acknowledge Dr. Marco Van De Wiel for his assistance with editing the manuscript.
 Li, L., Ngongondo, C.S., Xu, C.-Y. and Gong, L. (2013) Comparison of the Global TRMM and WFD Precipitation Datasets in Driving a Large-Scale Hydrological Model in Southern Africa. Hydrology Research, 44, 770.
 Zhang, X., Drake, N.A., and Wainwright, J. (2013) Spatial Modelling and Scaling Issues. In: Wainwright, J. and Mulligan, M. Eds., Environmental Modelling, John Wiley & Sons, Ltd., Chichester, UK, 69-90.
 Abbott, M.B., Bathurst, J.C., Cunge, J.A., O’Connell, P.E. and Rasmussen, J. (1986) An Introduction to the European Hydrological System—Systeme Hydrologique Europeen, “SHE”, 1: History and Philosophy of a Physically-Based, Distributed Modelling System. Journal of Hydrology, 87, 45-59.
 Abbott, M.B., Bathurst, J.C., Cunge, J.A., O’Connell, P.E. and Rasmussen, J. (1986) An Introduction to the European Hydrological System—Systeme Hydrologique Europeen, “SHE”, 2: Structure of a Physically-Based, Distributed Modelling System. Journal of Hydrology, 87, 61-77.
 Ciarapica, L. and Todini, E. (2002) TOPKAPI: A Model for the Representation of the Rainfall-Runoff Process at Different Scales. Hydrological Processes, 16, 207-229.
 Becker, A. (1992) Criteria for a Hydrologically Sound Structuring of Large Scale Land Surface Process Models. In: O’Kane, J.P., Ed., Advances in Theoretical Hydrology, Elsevier, Amsterdam, 97-111.
 Flügel, W.-A. (1995) Delineating Hydrological Response Units by Geographical Information System Analyses for Regional Hydrological Modelling Using PRMS/MMS in the Drainage Basin of the River Bröl, Germany. Hydrological Processes, 9, 423-436.
 Kavvas, M.L., Chen, Z.Q., Tan, L., Soong, S.T., Terakawa, A., Yoshitani, J., et al. (1998) A Regional-Scale Land Surface Parameterization Based on Areally-Averaged Hydrological Conservation Equations. Hydrological Sciences Journal, 43, 611-631.
 Arrigo, J.A.S. and Salvucci, G.D. (2005) Investigation Hydrologic Scaling: Observed Effects of Heterogeneity and Nonlocal Processes across Hillslope, Watershed, and Regional Scales. Water Resources Research, 41, W11417.
 Cerdan, O., Le Bissonnais, Y., Govers, G., Lecomte, V., van Oost, K., Couturier, A., et al. (2004) Scale Effect on Runoff from Experimental Plots to Catchments in Agricultural Areas in Normandy. Journal of Hydrology, 299, 4-14.
 Corwin, D.L., Hopmans, J. and de Rooij, G.H. (2006) From Field- to Landscape-Scale Vadose Zone Processes: Scale Issues, Modeling, and Monitoring. Vadose Zone Journal, 5, 129-139.
 Schulze, R. (2000) Transcending Scales of Space and Time in Impact Studies of Climate and Climate Change on Agrohydrological Responses. Agriculture, Ecosystems & Environment, 82, 185-212.
 Tague, C.L. and Band, L.E. (2004) RHESSys: Regional Hydro-Ecologic Simulation System—An Object-Oriented Approach to Spatially Distributed Modeling of Carbon, Water, and Nutrient Cycling. Earth Interactions, 8, 1-42.
 Running, S.W., Nemani, R.R. and Hungerford, R.D. (1987) Extrapolation of Synoptic Meteorological Data in Mountainous Terrain and Its Use for Simulating Forest Evapotranspiration and Photosynthesis. Canadian Journal of Forest Research, 17, 472-483.
 Running, S.W. and Coughlan, J.C. (1988) A General Model of Forest Ecosystem Processes for Regional Applications I. Hydrologic Balance, Canopy Gas Exchange and Primary Production Processes. Ecological Modelling, 42, 125-154.
 Beven, K.J. and Kirkby, M.J. (1979) A Physically Based, Variable Contributing Area Model of Basin Hydrology. Un modèle à base physique de zone d’appel variable de l’hydrologie du bassin versant. Hydrological Sciences Bulletin, 24, 43-69.
 Tague, C.L. and Band, L.E. (2001) Evaluating Explicit and Implicit Routing for Watershed Hydro-Ecological Models of Forest Hydrology at the Small Catchment Scale. Hydrological Processes, 15, 1415-1439.
 Sanford, S.E. (2005) Scale-Dependence of Natural Variability of Stream Flow Parameters in a Forested Drainage Basin on the Boreal Shield. M.S. Thesis, University of Western Ontario, London, Ontario.
 Monteith, S.S., Buttle, J.M., Hazlett, P.W., Beall, F.D., Semkin, R.G. and Jeffries, D.S. (2006) Paired-Basin Comparison of Hydrological Response in Harvested and Undisturbed Hardwood Forests during Snowmelt in Central Ontario: I. Streamflow, Groundwater and Flowpath Behaviour. Hydrological Processes, 20, 1095-1116.
 Murray, C.D. and Buttle, J.M. (2005) Infiltration and Soil Water Mixing on Forested and Harvested Slopes during Spring Snowmelt, Turkey Lakes Watershed, Central Ontario. Journal of Hydrology, 306, 1-20.
 Laporte, M.F., Duchesne, L.C. and Morrison, I.K. (2003) Effect of Clearcutting, Selection Cutting, Shelterwood Cutting and Microsites on Soil Surface CO2 Efflux in a Tolerant Hardwood Ecosystem of Northern Ontario. Forest Ecology and Management, 174, 565-575.
 Jeffries, D.S. and Semkin, R.G. (1982) Basin Description and Information Pertinent to Mass Balance Studies of the Turkey Lakes Watershed. Rep. 82-01, Turkey Lakes Watershed, Algoma, Ontario, Canada.
 Quinn, P., Beven, K., Chevallier, P. and Planchon, O. (1991) The Prediction of Hillslope Flow Paths for Distributed Hydrological Modelling Using Digital Terrain Models. Hydrological Processes, 5, 59-79.
 Nemani, R., Pierce, L., Running, S. and Band, L. (1993) Forest Ecosystem Processes at the Watershed Scale: Sensitivity to Remotely-Sensed Leaf Area Index Estimates. International Journal of Remote Sensing, 14, 2519-2534.
 Nash, J.E. and Sutcliffe, J.V. (1970) River Flow Forecasting through Conceptual Models Part I—A Discussion of Principles. Journal of Hydrology, 10, 282-290.