Assessment of the Sensitivity of Streamflow Simulations to Changes in Patch Resolution Using GIS Based Hydro-Ecologic Model

Eight different patch configurations were investigated to analyze the effect of patch characterization/formation in streamflow simulation, using the Regional Hydro-Ecologic Simulation Systems (RHESSys) model. It is investigated for eight different patch configurations of a subcatchment of the Turkey Lakes Watershed, Ontario. The model’s hydrological parameters are calibrated for each of these patch configurations and the performance of the simulations is evaluated. Results indicate that both the nature of the flow simulation and the calibrated parameter values are sensitive to patch configuration. The best simulation results were obtained for the patch configuration with the highest spatial variation of climate, stream network and hillslope conditions across the subcatchment. Different patch configurations also lead to markedly different calibrations of the model’s hydrological parameters (54.26 < k < 119.13; and 1.02 < m < 2.28), which has implications for the physical interpretation and transferability of the calibrated parameter values.


Introduction
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 [1].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 [1] [2].
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 [3]- [5].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 [1].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 [6] [7] and TOPKAPI [8].
This grid parameterization approach, however, is considered as inadequate by some authors.[9] 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 [10].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 [11]; the Representing Element Area concept [12]; hydrologically similar units or "patches" [13]; the Hydrological Response Units concept of [14]; and basin-scale model equations (e.g.[15]).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 [1] [5] [16]- [21].[17] 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) [22], 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 [22].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.

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 [23] within a distributed watershed context to represent the spatial distribution of hydrologic and ecophysiological processes and characteristics of a catchment [22].RHESSys is composed of three process submodels: a meteorological model (based on MT-CLIM; [24]), an eco-physiological model (based on BIOME-BGC; [25]), and two distributed hydrologic models.TOPMODEL [26] 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 [27] 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 [28].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 [22].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 [22].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 GISgenerated 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.

Study Area Description
This research was conducted on a subcatchment of the Turkey Lake Watershed (TLW).The TLW has an area of 10.5 km 2 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 [30].[31] 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 [32].
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 [33].Details on the physical and biological characteristics of the watershed are given by [34].
TLW has been used as an experimental basin since 1980 to study anthropogenic perturbations in Canadian Shield ecosystems in central Ontario [30].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 m 2 .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.

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 [34].Topographic attributes, including both primary (slope, aspect, specific contributing area) and secondary (wetness index) attributes, were derived from established terrain analysis techniques [26] [36]- [38].
Two sets of stream networks were derived from the DEM, using specific contributing area thresholds of 2000 m 2 and 5000 m 2 .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 m 2 and 5000 m 2 .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 [39].A detailed description of the physiological parameters for the vegetation classes is provided by [29].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 [22].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).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; [40]).For each patch configuration, a Monte-Carlo calibration was applied, in which 75 different combinations of m and k were simulated.

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.02and 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   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.

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

Figure 2 .
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) [29].

Figure 3 .
Figure 3. Location of subcatchment C38 within Ontario and the turkey lake watershed.

Figure 5 .
Figure 5. Maps of the eight patch characterization scenarios used in this study project.

Figure 6 .
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.

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

Table 1 .
Range of values of input files for patch creation.

Table 2 .
Tablethatshows the best m and k combinations and the resulting total simulated flow simulated for the eight patch scenario projects.