A generic system dynamics model for simulating and evaluating the hydrological performance of reconstructed watersheds

A generic system dynamics watershed (GSDW) model is developed and applied to five reconstructed watersheds located in the Athabasca mining basin, Alberta, Canada, and one natural watershed (boreal forest) located in Saskatchewan, Canada, to simulate various hydrological processes in reconstructed and natural watersheds. This paper uses the root mean square error (RMSE), the mean absolute relative error (MARE), and the correlation coefficient (R) as the main performance indicators, in addition to the visual comparison. For the South Bison Hills (SBH), South West Sand Storage (SWSS) and Old Aspen (OA) simulated soil moisture, the RMSE values ranges between 2.5–4.8 mm, and the MARE ranges from 7% to 18%, except for the D2cover it was 26% for the validation year. The R statistics ranges from 0.3 to 0.77 during the validation period. The error between the measured and simulated cumulative actual evapotranspiration (AET) flux for the SWSS, SBH, and the OA sites were 2%, 5%, and 8%, respectively. The developed GSDW model enables the investigation of the utility of different soil cover designs and evaluation of their performance. The model is capable of capturing the dynamics of water balance components, and may used to conduct shortand longterm predictions under different climate scenarios.


Introduction
Hydrological models have been adopted, modified, and applied to solve a wide spectrum of problems.The difficulty of modeling watershed hydrology lies primarily in that the response of the watershed system is strongly controlled by Correspondence to: A. Elshorbagy amin.elshorbagy@usask.caits spatial and temporal heterogeneity, and this heterogeneity cannot be precisely known or described.In general, the focus of watershed modeling studies has been on rainfall-runoff relations (Beven, 2001).Over the past few decades, countless number of watershed models has been developed around the world for a variety of applications.However, the main challenge remains in applying limited and imperfect knowledge of hydrological processes while providing an acceptable prediction of the real world (Schaacke, 2002).Hydrological processes such as soil moisture redistribution and evapotranspiration (ET), are intricately linked; therefore, the understanding of their mutual interaction could lead to a more accurate simulation of the processes responsible for land-atmosphere interaction (Mahmood and Hubbard, 2003).
As natural ecosystems are complex, their characteristics and dynamic properties depend on many interrelated links between climate, soil and vegetation (Rodriguez-Iturbe et al., 2001).Soil and climate control vegetation dynamics, while in turn vegetation modulates the water balance (Porporato and Rodriguez-Iturbe, 2002;Arora, 2002).Vegetation acts as the intermediate link between soil and the atmosphere via evapotranspiration, as well as affecting soil hydraulic and mechanical properties.Moreover, it affects both the surface energy budget and soil storage in the root zone (Falkenmark, 1997).Several models have been used to simulate soil-atmospheric-vegetation interaction, e.g.Soil-Vegetation-Atmosphere Transfer schemes (SVAT).SVAT model is used to simulate, energy, carbon, water fluxes, and typically assuming static vegetation (Arora, 2002).Other models used to simulate agriculture management scenarios, e.g.SWAT, are devoted to reproducing crop's growth, nutrient and pesticide practice, yet are data intensive to implement (Neitsch et al., 2002).Recently, Quevedo and Francés (2008) used a conceptual dynamic vegetation-soil model (called HORAS) for arid and semi-arid zones.The HORAS model Published by Copernicus Publications on behalf of the European Geosciences Union.
N. Keshta et al.: GSDW model for evaluating reconstructed watersheds consists of two reservoirs; the first for the interception process and the second for near-surface soil moisture.The HO-RAS model excelled in demonstrating the adaptation of vegetation to various climatic and soil conditions.However, it simulates only the upper capillary water-soil content related to where vegetation occurs unless it is coupled with a more comprehensive hydrological model (Quevedo and Francés, 2008).
Infiltration, soil moisture redistribution, and ET are the main hydrological processes affecting the behavior of natural and restored (reconstructed) watersheds in arid and semiarid regions.Hence, the proper simulation of these processes is vital to the accurate representation of the hydrology of both watersheds (Elshorbagy et al, 2007;Quevedo and Francés, 2008).Despite the importance of ET and soil moisture redistribution in defining the water balance of the arid and semiarid regions relative to the rainfall-runoff relation, there is limited literature available on the simulation of these processes.In addition, both ET and soil moisture dynamics have an important role in the ecological behavior of reconstructed watersheds following mining.
The rapid growth of the oil sands industry results in large disturbances to the natural ecosystem as soil and overburden materials are removed to provide access to mining materials.The mining process is followed by a remediation process, through which the disturbed landscape is recovered with the intent to replicate the performance of natural watersheds functions such as habitat function (hosting aquatic ecosystems), production function (e.g., biomass), and carrier function (for dissolved and suspended material), this process is also known as land reclamation (Haigh, 2000;Barbour et al., 2004).The adverse impact of disturbing the natural ecosystem can be intensified by climate change and its projected consequences.Consequently, it is crucial to simulate and predict, as accurately as possible, the hydrological behaviour of the reconstructed watersheds.Watershed models provide a vital tool that can assist in achieving this goal by simulating the hydrological behaviour of a variety of possible soil cover designs.The aim of this paper is to develop a generic system dynamics watershed model (GSDW), which provides a reliable, simple, and comprehensive tool that facilitates the assessment of the sustainability of various reconstructed/natural watersheds.The validity of the proposed model is assessed thorough its capability in reproducing the hydrological behaviour of the reconstructed and natural watersheds.The proposed model can be used to aid in decision making and contribute to the understanding of the nonlinear and complex hydrological processes of the reconstructed systems, especially in arid and semi-arid regions.In those regions, land reclamation is affected by the local climates where potential evapotranspiration is greater than the annual precipitation.Subsequently, the designed soil covers should have the ability of minimizing runoff, and retaining soil moisture for the growing season.

Modeling of reconstructed watersheds
In general, the literature of reconstructed watersheds emphasizes geotechnical perspectives, and Julta (2006) noted that most publications targeted modeling an individual component of the hydrological cycle.For example, HELP (Hydrological Evaluation of Landfill Performance) is a widely applied landfill water budget model (Yalcin and Demirer, 2002).Berger et al. (1996) used HELP to simulate the water balance of a landfill cover system, where they achieved good lateral drainage simulation, yet failed to model the linear leakage of comprehensive soil liners.The model only applies to simple covers, considers only grass as the vegetation type, and has a poor performance in estimating the long-term hydrologic processes (Berger, 2000).
The root zone water quality model (RZWQM) was used to simulate the volumetric soil water content of the reconstructed slopes of the South West Sand Storage (SWSS) in northern Alberta (Mapfumo et al., 2006).This model is used extensively in many agricultural studies; however, it tends to overestimate the saturated hydraulic conductivity and accordingly underestimates the surface soil moisture contents, especially in wet conditions.Furthermore, the soilatmosphere model (SoilCover) (Geoanalysis Ltd., 2000) was used, by Shurniak (2003), to predict moisture movement in a variety of reconstructed soil cover systems.Shurniak (2003) recommended that the overall cover thickness to be more than 0.6 m to improve plant survival.
Finally, Elshorbagy et al. (2005) developed a site-specific system dynamics watershed (SDW) model to simulate hydrological processes using a daily time-step in reconstructed watershed in northern Alberta, Canada.This model was extended by Elshorbagy et al. (2007) to simulate three prototype watersheds, however, the model remained site-specific.Elshorbagy and Barbour (2007) presented a probabilistic approach, using the SDW model, to assess the long-term hydrologic performance of three inclined reconstructed watersheds.The validated SDW model was used along with the available meteorological historical data to generate continuous simulated records of the daily depth-averaged soil moisture content.These records were used to estimate the maximum annual moisture deficits as indicators of the hydrologic performance of the considered watershed.The probabilistic approach was used to quantify the predictive uncertainty of the SDW model.However, the efficiency of this approach depends on the reduction of the predictive uncertainty of the model, which can be mitigated through a generic model that can simulate various reconstructed and natural watersheds under potential and uncertain changes of the prevailing climatic conditions.

GSDW model development and formulation
The proposed GSDW model is a lumped conceptual model capable of simulating various components of watershed Hydrol.Earth Syst.Sci., 13,[865][866][867][868][869][870][871][872][873][874][875][876][877][878][879][880][881]2009 www.hydrol-earth-syst-sci.net/13/865/2009/ hydrology.This model is an upgrade/generalization of the existing site-specific SDW model, which was developed by Elshorbagy et al. (2005Elshorbagy et al. ( , 2007)).The main drawback of the previous SDW model was that it did not account for a canopy interception.Moreover, the SDW model has no flexibility in choosing the number of soil layers (three layers only), layer thickness, and topographic inclination.The GSDW model uses sets of meteorological, vegetation, and hydrological data to evaluate different hydrological processes on a daily basis.The model is entitled "generic" in the sense that it can be implemented on a wide spectrum of watersheds, soil cover alternatives, and topographic conditions in semiarid regions for the purpose of assessing the performance of reconstructed watersheds.The system dynamics simulation environment (STELLA), (HPS, 2001), was used for modeling the watershed as a dynamic system in a user-friendly environment.
The system dynamics (SD) approach is based on the understanding of the complex relationships existing among the different elements within the considered system.In general, the SD approach can be defined as: "a theory of system structure and a set of tools for representing complex systems and analyzing their dynamic behavior" (Forrester, 1980a,b).Ford (1999) defined the SD approach as a method of analyzing problems in which time is an important factor.The main issue in using SD is to understand the system and its boundaries by identify its key building blocks, and the proper representation of the physical processes through relatively accurate mathematical relationships.SD models have the potential of implementing a combination of empirical formulations and physically based concepts and also allows for building on a tentative knowledge of the relation between two parameters by incorporating a qualitative relationship between those parameters.(Elshorbagy et al., 2007).The proposed GSDW model will have the ability to simulate relevant hydrological processes, e.g., canopy interception, evapotranspiration, surface runoff, lateral interflow, infiltration, and soil moisture redistribution in unsaturated/saturated layers, based on the surface energy and water balances.Particular attention is given to the parameterization, which is kept as simple as possible and reliant on widely available data.A schematic diagram of the major processes modelled by the proposed GSDW model is shown in Fig. 1.
Figure 1 shows the simple daily water balance of the GSDW model, which consists of three storage components: (1) canopy storage, (2) surface storage, and (3) soil storage.The system dynamics hypothesis of the developed GSDW model is represented in a causal-loop diagram (Fig. 2).The feedback loops illustrate the mutual interaction among the different factors affecting watershed hydrological processes.Negative and positive signs denote the type of relationship between corresponding variables.Figure 2   lowing outlines modifications and improvements suggested to the SDW model, whose details can be found in Elshorbagy et al. (2005Elshorbagy et al. ( , 2007)).

Canopy storage
Interception losses range from 10-40% of gross precipitation for different vegetation types (Dingman, 2002).Therefore, consideration of interception losses will improve AET predictability of the developed model.One challenge posed by the incorporation of vegetation processes in hydrological models is the incorporation of new parameterizations.parameters representing the aforementioned physical characteristics, such as the canopy storage capacity (S c ), trunk storage capacity (S t ), trunk evaporation as a fraction of the total evaporation (ε), and the canopy drip as a fraction of the drainage (1−ρ).These parameters are based on detailed information of the vegetation structure.The leaf area index (LAI) is used as an indicator of the canopy interception, which can be deduced based on the ratio of the canopy shaded areas to the bare areas.The evaporation rate from the canopy (E c ) is computed as the sum of both the trunk and the canopy evaporation.The Penman equation is used to compute the rate of evaporation (E p ) of the intercepted water.Equations (1), (2), and (3) are the mathematical representation of the evaporation rates from different canopy components: where E cc is evaporation from leaves (mm), C c (t) is the actual amount of water stored on the canopy leaves in mm, E ct is the trunk evaporation, C t (t) is the actual amount of water stored on the trunk in mm, and F is the fraction of area covered by the forest canopy.The main practical drawback of the Valente et al. (1997) model lies in its extensive data requirements.
The van Dijk and Bruijnzeel (2001) model; based on a modification of Gash et al. (1995) interception model, retains some of the simplicity of the empirical approaches.It is based mainly upon the LAI, and the canopy storage (S c ).The model assumptions are: (i) the relative evaporation rate Ē/ R can be expressed as a function of LAI, (ii) the canopy storage capacity (S c ) is linearly related to LAI, and (iii) the LAI was treated as static value during the growing season for each case study.This approach is represented in the GSDW model by the following: where S L denotes the specific leaf storage (the depth of water retained by the leaf per unit LAI).The S L -values, as suggested by Pitman (1989), range between (0.4-5.88), c is the canopy cover fraction, k is the extinction coefficient and it depends on the leaf inclination angle and distribution, the kvalues ranges between 0.2 and 0.8, and E p is the Penman potential evapotranspiration.The GSDW model provides the user with the flexibility to use either one of the previous two approaches based on available data, in addition to a third selection, where the canopy interception is not incorporated due to the lack of information regarding the canopy coverage, or the absence of canopy in the case of newly reconstructed watersheds.

Surface water storage
The change in the surface water storage (SW) can be expressed by: where P (mm/day) represents precipitation, in either the form of snow or rainfall, f L1 is the infiltration rate to the top soil layer (mm/day), and O F represents overland flow in mm/day.

Soil storage
The developed GSDW model is designed to facilitate the consideration of multiple layers of soil cover, as opposed to pre-set number of covers in the SDW model.This expands the applicability of the model to simulate a wide variety of alternatives, in addition to, enhancing soil moisture predictably.Therefore, the vertical movement of the soil moisture between any two subsequent layers is described by considering layer (i) as a control volume in the water balance.For example, the change of the moisture storage in the ith layer depends on downward movement of water from the upper (ith-1) layer, evapotranspiration, interflow from the ith layer, and the downward water movement to the underlying (ith+1) layer.Therefore, the change of moisture storage in the ith layer can be expressed as follows: where S i is ith layer storage in mm,f i is the downward water movement rate of the ith layer, f i+1 is the downward water movement rate to the underlying layer in mm/day, I i is the interflow rate for of the ith layer in mm/day, and ET i is the evapotranspiration rate from ith layer in mm/day.Voinov et al. (2004) suggested that infiltration rate of the top soil layer is equal to rainfall intensity before soil saturation is reached.In fact, some studies suggested that the rate of infiltrated water from a typical cover system is correlated to the degree of saturation of the soil, soil moisture retention characteristic, and climatic factors (e.g.rainfall) (Milczarek et al., 2000;Milczarek et al., 2003).On the other hand, the Green-Ampt equation governs the vertical movement of water during the saturation stage under the condition of soil temperature being greater than zero (unfrozen soil).The infiltration capacity (rate) based on total infiltration volume is expressed by the Green-Ampt equation in the case of a thawed saturated soil (Dingman, 2002): where K si is the saturated hydraulic conductivity of the ith layer in mm/day, θ si is the saturated moisture content of the ith layer (%), θ ii is the initial moisture content of the ith layer (%), ψ i is the suction pressure head at the wetting front in the ith layer in mm, and F i is the cumulative volume of infiltration in the ith layer in mm.There exist methods for quantifying infiltration into frozen soils, however such methods are highly data-intensive (Elshorbagy et al., 2007).Some studies suggest that the frozen soil layer does not impede infiltration (Iwata et al., 2008).An empirical approach for snowmelt infiltration was suggested by Li and Simonovic (2002) and has been validated by Julta (2006).This approach is based mainly upon the idea that infiltration rates in frozen soils are influenced by temperature and temperature accumulation.The infiltrated water will gain its dynamics based upon the temperature index, where soil will refreeze if the temperature drops below zero for a certain number of days.The active temperature accumulation will be lost and will start again from zero (Li and Simonovic, 2002).Consequently, infiltration into frozen soil is computed by multiplying infiltration rate of the ith layer, f i , by an empirical coefficient, C ti .Reference is made to Elshorbagy et al. (2007) for further details.
The movement of water between any subsequent layers is limited if the upper layer moisture content is less than the residual moisture content.Moisture movement will start when the upper layer moisture is greater than the residual moisture content, and it starts contributing to the lower layer moisture until it reaches saturation.Once the lower layer reaches saturation, the maximum rate at which water can be absorbed by the lower layer will correspond to the minimum value of the saturated hydraulic conductivity of this layer and if (soil temperature of (i) greater than 0 then (Min (drainable water from layer (i−1), follow Eq. 9) else (follow Eq. 10) else (no movement of water)) else (infiltration in frozen soil) where, Eq. ( 10) is an empirical equation, which can be written as follows (Elshorbagy et al., 2007): where I ci is the coefficient of the ith layer infiltration, which is determined during calibration of the model, and t is the solution time interval.Equation ( 10) suggests that the moisture redistribution between any subsequent layers is strongly dependent on the moisture contents of both layers.In addition, no downward moisture movement is allowed if the suction of the upper layer is greater than that of the lower layer.

Evapotranspiration module
In the model, the potential evapotranspiration is computed using the Penman equation derived in Mays (2005), while an empirical formula is used for the actual evapotranspiration calculation, based on the simulated soil moisture index, and the air temperature.To calculate the actual evapotranspiration from any soil layer; an empirical formulation used by (Julta, 2006;Elshorbagy et al., 2007) takes into consideration the available moisture, air and soil temperatures.
where c p is the evapotranspiration coefficient (mm/ • C/day), S ms(i) is the effective moisture saturation in layer (i) (dimensionless), λ is the exponential coefficient that expresses the impact of water saturation on evapotranspiration, S (i) , S N(i) , and S (i)rs are the water storage, the nominal water storage (mm), and minimum storage that can be attained (residual moisture) in layer (i), respectively.Sankarasubramanian and Vogel (2002) noted that the classical actual evapotranspiration relationships perform poorly for basins with low soil moisture storage capacity.However, the previous empirical formulas provided a better simulation than the conventional Penman equation (Elshorbagy et al., 2007).These parametric empirical equations provide better estimates of the actual evapotranspiration due to their dependence on the soil moisture content of the considered layer.The GSDW model adjoins both the estimated E c from the canopy storage and from different soil layers (net AET) to model the total actual evapotranspiration (AET).

Interflow component
The interflow component is restricted to the incidence of sloping layers.The model formulations account for the angle of inclination, which affects the interflow rate.Interflow (I i , mm/day) from the ith layer is estimated using the following modified empirical formula, used by Elshorbagy et al., 2007, as follows: where C Slope is the slope coefficient that depends on the slope value.D i is the depth of the ith layer in mm, and C i is the interflow coefficient, which is also a calibration parameter.Interflow generation is restricted to two conditions; the temperature of both layers, (i) and (i−1), are above zero and the ith layer is saturated.If the previous conditions are fulfilled, then interflow is computed by multiplying the interflow coefficient by the rate of available water in the ith layer above saturation level and the slope coefficient.However, if the temperature of layer (i−1) is below 0 • C and the ith layer storage is between saturation and field capacity, then interflow is computed by multiplying the interflow coefficient by the drainable water in layer (i) and the slope coefficient.
Another alternative for computing the interflow component was incorporated to the model structure based on twodimensional kinematic storage model for subsurface flow along a steep hillslope (Sloan and Moore, 1984).This module is adopted to calculate subsurface flow in a variety of hydrological models, e.g.SWAT, based on the mass continuity equation (Neitsch et al., 2002).Sloan and Moore (1984) prescribed high hydraulic conductivities in surface layers and an impermeable or semi-permeable layer at a shallow depth.The saturated thickness (H o ) normal to the bottom of the slope is expressed as: where φ d is the drainable porosity of the soil layer (dimensionless) and is equal to the difference between the total porosity of the soil layer (dimensionless) and the porosity when the layer is at field capacity (dimensionless), Lenght Hill is the hillslope length in mm, and drainable water from the layer is the stored water above field capacity.The drainable water from the layer is calculated using the following: The interflow at the hill slope outlet, I i in (mm/day) for layer i, is: where v Lat is the velocity of flow at the outlet.
The GSDW model allows user to select between Sloan and Moore (1984) technique and the modified Elshorbagy et al. (2007) empirical formula.

Overland flow component
Overland flow (O F ) is estimated using a modified empirically based equation from the SDW model.Since the model structure uses reservoir-based mechanisms to simulate the different hydrological processes, water in excess of infiltration capacity (saturation condition) of the first layer is directed as overland flow in the summer, and it is computed as: where f L1 layer (1) infiltration rate, O F is the overland flow in mm/day.Overland flow generation is also dependant upon the air/soil temperature and the gradient of the soil cover.
4 Case studies

Reconstructed watersheds
The GSDW model is used to simulate the hydrological performance of various watersheds to validate its capability in capturing the dynamics of the various water balance components in different sites.Verifying the ability of the GSDW model to simulate both reconstructed and natural watersheds confirms the utility of the proposed model to conduct shortand long-term predictions under different climatic conditions.The reconstructed watershed study areas are located in north of Fort McMurray (57 • 39 N and 111 • 13 W), northern Alberta, Canada.The oil sands industry has developed a system for stabilizing the surface of the reconstructed soil covers that enables re-vegetation.A few reconstructed watersheds formed of various soil covers (various soil types, layering, and depths) are selected for this study: the first site includes three inclined prototype soil covers (D1, D2, and D3-Covers).The three covers were constructed with a thickness of 0.5 m, 0.35 m, and 1.0 m compromised of 0.2 m, 0.15 m, and 0.2 m of peat/mineral mix overlying 0.3 m, 0.2 m, and 0.8 m thickness of glacial till, respectively, overlying saline sodic shale.The purpose of these experimental covers is to evaluate the performance of different alternatives in terms of moisture holding capacity and sustaining the vegetation.
The three covers have a slope of 5 H:1 V with an area of 1 ha each.These covers were constructed in 1999 and seeded with barley nurse crop (Hordeum Jubatum), and tree seedlings of white spurce (Picea glauca) and aspen (Populus tremuloides) (Boese, 2003).The D-covers were used by Elshorbagy et al. (2007) to develop the site-specific SDW model.The second site is the Hill top site; a flat reconstructed cover system located on adjacent to the D-covers.(Parasuraman et al., 2007).
An intensive hydrological and meteorological measurement program is carried out on these experimental sites to monitor the evolution of the reconstructed watersheds.The hydrological variables include the matric soil suction, volumetric moisture content (measured bi-daily using TDR and soil suction sensors at different depths), and soil temperature of different soil layers, measured on hourly basis for the corresponding soil moisture measurement depths.Additional monitored variables include runoff and interflow.Measurements of the latent heat fluxes are made with the eddy covariance technique (EC) and reported in 30-min interval (Carey, 2008).A weather station is used to provide hourly meteorological measurements of air temperature (AT), precipitation (P), net radiation (NR), water vapour gradients and other meteorological variables.More details on the field instrumentation and monitoring program can be found in Boese (2003), Julta (2006) and Carey (2008).Based on the climate data from an Environment Canada meteorological station at Fort McMurray, a 30-year period (1971McMurray, a 30-year period ( -2000) ) the mean annual temperature is 0.7 • C, and the mean annual precipitation is 455.5 mm.The soil properties are as follows: the saturated hydraulic conductivities are 408 (cm/day), 50.4 (cm/day), and 0.72 (cm/day), and the porosity values are 0.5, 0.54, 0.25 for peat, till and shale, respectively (Elshorbagy et al., 2005 (Cuenca et al., 1997).The forest canopy is dominated by trembling aspen (Populus tremuloides) with an average height of 21 m and about 2 m high hazelnut (Corylus cornuta) understory interspersed with alder (Balland et al., 2006).AT and P data are collected at 30-min intervals.Based on the data from an Environment Canada meteorological station nearby Waskesiu Lake (53.92 • N, 106.07 • W), the mean annual precipitation was 467 mm.Thermocouple sensors recorded soil temperature every 30min at 0.02, 0.05, 0.10, 0.20, 0.50 and 1.00 m below the moss layer.CS615 soil moisture sensors (TDR) were used to measure the volumetric moisture content of the soil at 0.08, 0.23, 0.45, and 1.05 m below the ground surface.Net radiation was measured using a Kipp and Zonen CNR-1 net radiometer above the canopy.Measurements of the latent heat fluxes were made with the eddy covariance technique (EC) and reported with 30-min interval.LAI was measured near the flux tower using a plant canopy analyzer (PCA) (model LAI-2000).For the purpose of integrating the extracted data into the GSDW model, all the data were aggregated to a daily time-step.Additional information for the OA site can be obtained from Cuenca et al. (1997).
The values of total precipitation for the D-covers, and the SBH site are 341.2mm and 294.3 mm for years 2005 and 2006, respectively.For the SWSS site, the value of total precipitation is 285.9 mm, and 366.3 mm in the years 2005-2006, respectively.Finally, the corresponding values for the OA site were 479 mm, and 483.7 mm in the years 1999-2000, respectively.
The main objective for this study was to develop a generic system dynamics model to simulate different hydrologic processes in disturbed and undisturbed watersheds, and various reconstructed sites were simulated to evaluate the model performance over a variety of different reclamation strategies.The model was applied to inclined watersheds, as in the case of the D-covers and the SWSS sites, to horizontal terrain, as in case of the SBH and OA sites, and on a set of different soil layers with a variety of thicknesses and soil stratifications.

Results and analysis
The goodness-of-fit between the measured and simulated datasets are generally quantified using multiple performance indicators, providing different aspects of comparison.This paper uses the root mean square error (RMSE), the mean absolute relative error (MARE), the percent of error measurement in the peak (PEP), and the correlation coefficient (R) as the main performance indicators, in addition to visual comparison.Both RMSE and MARE are overall error measures, where RMSE is a real value metric and MARE a relative value metric.The RMSE is biased towards high values, while the MARE is less sensitive to high values as it does not square the error magnitude (Dawson et al., 2007).Due to these limitations, R is used as a complementary error measure that quantifies the overall agreement between the observed and predicted values.PEP focus on the peaks values, yet not to the overall agreement between the two datasets.It comprises the difference between the highest values of observed and simulated datasets, made relative to the magnitude of the highest values in the observed dataset and expressed as a percentage and for a perfect model the PEP value is zero.
As mentioned, the traditional realm of hydrologic models has been the prediction of rainfall-runoff process.These predictions are used in flood warning systems, navigation, water quality management and many water resource applications.However, land reclamation concentrates on replicating the performance of natural watersheds in terms of supporting vegetation growth.Consequently, it is crucial to simulate hydrological processes that directly impact the ecological function of the watershed.As a result, the calibration of the GSDW model was performed based on two hydrological processes connected with the ecological function of the reconstructed watersheds: soil moisture and actual evapotranspiration (AET).Calibration was performed by setting individual parameter values and executing a series of simulations.This process was repeated (trial and error) until no further improvement in the values of the error measures and the visual match between simulated and observed AET could be attained.Table 1 lists the calibration parameters used on the different study sites together with their corresponding values.The calibration of the GSDW model indicates sensitivity to the lambda coefficients (λ n), a main factor in the Hydrol.Earth Syst.Sci., 13,[865][866][867][868][869][870][871][872][873][874][875][876][877][878][879][880][881]2009 www.hydrol-earth-syst-sci.net/13/865/2009/ AET equations, and the infiltration coefficients (Icn), which directly affect the moisture distribution in each layer.
The terrestrial ecology community over the last two decades has developed models to simulate various hydrological processes.Operational applications are demanding, as they request both efficiency and robustness.Therefore, there is always a debate of the modeling approach that should be selected; model choices must be justified through extensive testing and for robustness considerations, where simple models must be preferred over more complex models when they are of equal efficiency.The next section compares the performance of the GSDW model with the Soil and Water Assessment Tool (SWAT), a widely used hydrological model.

SWAT model
SWAT is a basin-scale, continuous-time model created in the early 1990s that operates on a daily time step and is designed to predict the impact of management on water, sediment, and agricultural chemical yields in ungauged watersheds (Arnold and Fohrer, 2005).The model is physically based and capable of continuous simulation over long time periods.SWAT components include weather, hydrology, soil temperature and properties, plant growth, nutrients, pesticides, bacteria, pathogens, and land management.In SWAT, a watershed is divided into multiple sub-watersheds, which are then further subdivided into hydrologic response units (HRUs) that consist of homogeneous land use, management, and soil characteristics.Climatic inputs used in SWAT include daily precipitation, maximum and minimum temperature, solar radiation data, relative humidity, and wind speed data.The overall hydrologic balance is simulated for each HRU, including canopy interception of precipitation, partitioning of precipitation, snowmelt water, redistribution of water within the soil profile, evapotranspiration, lateral subsurface flow from the soil profile, and return flow from shallow aquifers.More-  over, bypass flow can be simulated, as described by Arnold et al. (2005), for soils characterized by cracking.SWAT has proven to be an effective tool for assessing water resource and nonpoint-source pollution problems for a wide range of scales and environmental conditions across the globe (Arnold et al., 1998;Arnold and Fohrer, 2005)  plant uptake compensation factor (EPCO) and soil evaporation coefficient (ESCO) were the most sensitive parameters that affected the simulation results.In the current study, both watersheds are considered one hydrologic response unit (HRU) that consist of homogeneous land use, management, and soil characteristics.Table 2 presents SWAT model performance indicators with regards to its ability to simulate the soil moisture content for D2 and D3 covers.In general the GSDW model results, presented in Table 2, rival the SWAT model in the performance, particularly for the top layer.However, SWAT performance in the subsequent layer was slightly better.Figure 3 shows the soil moisture dynamics for the D3 cover.SWAT was capable of simulating the snowmelt period better than the GSDW model, which may be attributed to the fact that SWAT was capable of simulating bypass flow.In contrast, the GSDW was able to reproduce the soil moisture dynamics during the growing season in the upper layer (which is the most important layer for vegetation) superior to SWAT.Cumulative AET-values were overestimated using SWAT by 17% and 15% for D2, and D3, respectively.

Simulation results of the GSDW model
As the GSDW model is an advance on the pre-existing SDW model, preliminary runs were made to validate the model performance.First, the SDW model was recalibrated for 2005 and validated on cover D3 for 2006.The SDW model MARE were 16% and 6% for validation year compared to the GSDW model values of 8% and 6% for peat and till, respectively.The RMSE for the SDW model were 13.5 and 17.3 mm, whereas, 3.1, 4.3 mm for the GSDW model, for peat and till, respectively.The simulated cumulative AET were 253 and 283 mm for the validation year, for the SDW and the GSDW models respectively, compared to 276 mm of measured AET.The results presented in Table 2 show that the model performance with regards to the D3-cover was superior to the preliminary run for the site-specific SDW model built for the D-covers by Elshorbagy et al. (2005).A second attempt was made to recalibrate the GSDW model to simulate the moisture dynamics of D1-cover and compare the results to the findings by Elshorbagy et al. (2007).Therefore, 2001 and 2002 were chosen as calibration and validation years, respectively.The MARE values of the GSDW model were 7% and 8.2%, and the RMSE values were 4.3 mm and 7.0 mm for peat and till layers, respectively.The previous values were compared to MARE of 9% and 11%, and RMSE of 7.0 and 9.0 mm for peat and till layers, respectively, for the SDW model.The improvement in model performance was mainly attributed to implementing the canopy interception module to the GSDW model.Table 2 lists the GSDW model performance indicators with regard to its ability to simulate the soil moisture content of the study areas.For the SBH, SWSS and OA sites, the model provided satisfactory results, the RMSE values ranged between 2.5-4.8mm, which indicates that the average error was not more than ±5 mm away from the mean soil moisture value.
The moisture dynamics of the thinnest soil cover (D2) had a flashy response compared to the other two D-covers.The R statistic indicated that the GSDW model captured the general trend of the soil moisture in particular for the surface peat layer where R ranges from 0.30 to 0.77 in the validation phase.The subsurface till layer of the D2-cover showed a relatively low correlation of 0.10, whereas the SWSS site showed negative correlation coefficient, which is attributed to the high spatial variability of the soil moisture measurements in the reconstructed watershed, as well as the effect of the depth-averaging for the observed values of soil moisture.For the overall water balance, these errors are very small in terms of water depth (mm).
The simulated soil moisture dynamics of the surface and subsurface soil layers are shown in Figs. 4, 5, and 6 for the SWSS, SBH, and the natural old aspen watersheds, respectively.For the three figures, in the winter period, there was no significant dynamics for the moisture content because soil during this period was frozen and the model behaved accordingly.As cited by Boese (2003), the sensors used for the measurement of soil moisture at the study sites were not operating reliably during the frozen conditions; also, the AET values should be neglected during the winter season.Therefore, the evaluation of the soil moisture behaviour in the winter season is not significant and may mislead the analysis of the model results.Therefore, only the values of the growing season are considered.As the air temperature reaches active threshold value, snow starts melting and liquid water infiltrates into the soil layers.A sudden increase of moisture content ensues once the surface layer is thawed and corresponds to the amount of snow that is accumulated when the temperature was below zero.After the snowmelt period, soil moisture in the surface layer fluctuates due to the variation of rainfall intensity and evapotranspiration.This period lasts until the temperature falls down below the active air temperature and the soil starts refreezing again in the fall.Figures 4, 5, and 6 show consistencies of soil moisture patterns related to corresponding rainfall events.The surface layer storage component was  4 9  7 3  9 7 1 2 1 1 4 5 1 6 9 1 9 3 2 1 7 2 4 1 2 6 5 2 8 9 3 1 3 3 3 7 3 6   responsive to rainfall events, whereas the responses of the subsurface layers were not as rapid.In general, the results indicated that the simulated soil moisture patterns were similar to the observed patterns.
The PEP measurement was applied to the results of the model and presented in Table 2.The negative sign indicates an over-estimate of the peaks in the simulated runs whereas the under-estimate produces a positive sign.There are several recorded increases in the observed soil moisture at the SWSS site when the soil was frozen (Fig. 4).This could be attributed to an error in the measurement or contribution of preferential flow to the soil pores.Figures 5 and 6 show increases in the simulated soil moisture storage in the 2nd layer during the summer period for both SBH (year 2006) and OA  4 9 7 3  9 7 1 2 1 1 4 5 1 6 9 1 9 3 2 1 7 2 4 1 2 6 5 2 8 9 3 1 3 3 3 7 3 6   (year 1999).This sudden increase is correlated with rainfall events of 41 mm and 60 mm, respectively.During immense rainfall events, the increase in the soil moisture storage is due to the restriction on the lateral subsurface flow movement in flat landscapes.Figure 7 presents the cumulative AET over the growing season period for the validation years measured using the EC versus the simulated AET values.The graph presents an overall agreement of the observed and the simulated cumulative AET for the three sites, in both magnitude and trend.The reasonable match between measured and simulated AET values provide another indication of the ability of the GSDW model to capture the dynamics of  the hydrological processes in the reconstructed and the natural sites.The model slightly overestimated the cumulative AET fluxes in the natural OA sites, where the measured AET values using the EC method was 338 mm and the simulated was 365 mm.The GSDW model under-estimated cumulative AET fluxes for SWSS site, where values were 319 mm and 312 mm, for observed and simulated, respectively.For the SBH site, the GSDW model underestimated cumulative AET values, with a measured cumulative AET of 276 mm and a simulated AET flux of 261 mm.The error between the measured and simulated cumulative AET values for the SWSS, SBH, and the OA sites were 2%, 5%, and 8%, respectively.Figure 8 shows observed and simulated daily AET at SWSS, SBH, and OA, respectively, where the model simulated daily AET flux relatively good.Mainly, reconstructed covers, e.g. the inclined D-covers, were designed to act as a sponge, regarding store and release abilities.The top layer was selected having a high hydraulic conductivity to increase infiltrated surface water in to the soil  matrix and minimizing runoff generation, and redistribute the trapped liquid water to the subsequent layer to reduce evaporation.The subsequent layer is relatively fine textured with clay particles dominating the soil, the clay content reduce the hydraulic conductivity and at same time allows for high porosity to allow for storing more moisture for the growing season.Figure 9 shows the observed and simulated overland flow for the inclined D1 cover for 2001, where the model triggered some runoff during snowmelt period.
The annual water balance components, during validation periods, were simulated by the GSDW model and are summarized in Table 4.The intercepted precipitation from implementing the canopy module ranged from 4% to 15% of the total precipitation, runoff varied from 0.5% to 5% of the total precipitation, and ET ranged from 78% to 98% of the total precipitation.The difference between precipitation as input, AET, interception, and runoff as outputs is the percolated water to subsequent layers.

Discussion
Many hydrological models are able to represent hydrological processes at the watershed scale, but the majority with high parameter requirements.Therefore, there is a need for www.hydrol-earth-syst-sci.net/13/865/2009/ Hydrol.Earth Syst.Sci., 13, 865-881, 2009 a tool that facilitates the simulation of the response of various soil cover designs and to evaluate their performance relying on widely available data.During the preliminary stages of the watershed development, soil moisture plays the key role in vegetation growth, especially in the root zone layers (Kilmartin, 2000).The GSDW model is able to simulate the soil moisture response with a good accuracy of less than 5 mm, on average, from observed values.Previous efforts of simulating similar sites resulted in agreement between the observed and simulated values in trend only, not in magnitude.For example, Balland et al. (2006) modelled the snow pack, soil temperature, and soil moisture in the OA site, and achieved agreement between the measured and simulated snow pack and soil temperature, yet for soil moisture, the agreement was in trend not in magnitude.They attributed this difference in model performance to local conditions and sensor surroundings, in addition to different uncertainties associated with the modelling procedure itself.
The preliminary runs showed that the GSDW model surpasses the previous SDW model in performance, which may be attributed to the introduction of the canopy interception module in the model structure.The GSDW and SWAT preliminary runs showed a comparable performance.SWAT has the ability of subdividing a watershed into hydrologic response units (HRUs), which gives the model credibility in simulating large watersheds.In the current study, both reconstructed watersheds were considered as a single homogeneous hydrologic response unit, and the simulation was conducted using a set of lumped data.The trial and error procedure, which was implemented during the model calibration, in conjunction with the previous assumptions, reduced SWAT abilities in simulating soil moisture dynamics.The aim of comparing SWAT and GSDW models was to supply the modeller with an idea of which model to use based on a set of lumped set of data.The use of complex physically based models, such as SWAT with lumped inputs, is neither effective nor economic.The GSDW model has a simpler structure, with less number of calibrated parameters, and relies on widely available meteorological data.
Evapotranspiration is a key process for water resources management, particularly in arid regions.AET depends on a large variety of factors: vegetation, soil type, topography, and the meteorological conditions.The rate of evapotranspiration is largely controlled by available energy and available soil moisture.Soil and climatic variables control vegetation dy-  namics, while in return vegetation modulates the water balance by acting as an intermediate link between soil and the atmosphere.Soil moisture and vegetation affect the thermal inertia and shortwave albedo of the surface.Furthermore, numerous studies support the assertion that in arid and semi arid regions, soil moisture flux is the key variable in soilvegetation-atmosphere continuum (Rodriguez-Iturbe, 2000;Rodriguez-Iturbe et al., 2001;Porporato et al., 2001).Weeks and Wilson (2006) conducted a study to predict the surface water balance for constructed soil covers on waste disposals.The study showed that the slope direction and angle can have a significant effect on the net radiation received by the slope, and hence affect evaporation predictions.As a result of the assumptions associated with the model structure, the daily evapotranspiration is simply correlated with moisture availability without relying on other influencing variables, e.g., characteristics of the surrounding environment and type and condition of vegetation, which in turn affects the daily AET simulations.The EC method, which is used as a direct measurement of AET, has an accuracy range from ±15 to ±20% for hourly evapotranspiration measurements and up to ±8 to ±10% for longer periods (Eichinger et al., 2003;Strangeways, 2003).However, the GSDW model managed to simulate the cumulative annual AET to a very reasonable accuracy, less than ±8% of the annual measured AET value, with minor overestimation and underestimation periods.
Generally, three sources of uncertainties in the modeling process can be distinguished: errors in input variables, model and parameterization, and algorithms of process description.To visualize one type of uncertainty associated with errors in the input variables, Gee and Hillel (1988) pointed out that precision in precipitation is seldom less than ±5%.The spatial and temporal variability of soil physical properties within the same site adds an extra level of uncertainty to the measured data.The required characteristic of the soil physical properties for the GSDW model, such as the saturated hydraulic conductivity and the pore size distribution, are subject to high degree of spatial and temporal variability (particularly in reconstructed soil covers).The saturated hydraulic conductivity in the reconstructed D-covers, e.g., increased by 400% from 2000 to 2001, as tabulated by Elshorbagy et al. (2007).All models are limited to the extent that the parameterizations of physical processes are only approximations of the true soil physics and vegetation physiological action.Therefore, a long monitoring period for reconstructed watersheds is essential, and is suggested by Rick (1995) to be seven years.This monitoring period will allow tracking of different changes and evolution encountered in reconstructed watersheds.
There were no daily observed values for the canopy losses to compare the model results, however, the adopted approaches for calculating the interception loss were validated in several previous studies.Moreover, testing other components of the water balance implicitly validated the interception component.The GSDW model does not account for macropores.Nevertheless, flow in macropores may play an important role for soil water fluxes.The soil moisture predictions during snowmelt in some case studies are not well represented, which could be attributed to the deficiency in representing the flow of water through macropores.During the snowmelt period and while the soil is still frozen, melted snow tends to bypass into soil layers through the macropores, giving a sudden increase in soil moisture.Further improvement to the GSDW model could be achieved by incorporating macropore flow, enhancing soil moisture simulation and consequently other hydrological processes.The model shows sensitivity to AET and infiltration coefficients-related calibration parameters, which confirm that AET process and soil moisture content, play the dominant role in simulating the hydrological performance of the watersheds.

Conclusions
This study presents a generic system dynamics watershed (GSDW) model.GSDW is a simple, reliable, and comprehensive tool that facilitates hydrological simulation, and thus the assessment of the sustainability of various reconstructed watersheds.It is a lumped conceptual model capable of simulating various components of watershed hydrology.The model uses sets of meteorological, vegetation, and hydrological data to evaluate different hydrological processes on a daily basis.The validity of the proposed model was assessed by evaluating the predicted soil moisture redistribution, and actual evapotranspiration in different sites reasonably well in trend and in magnitude.The simulation results showed that the model performance with regard to the D-covers was quite comparable to the findings of Elshorbagy et al. (2005Elshorbagy et al. ( , 2007)), with considerable improvements in soil moisture simulation.For the three case studies, the model provided good results, based on the three selected performance measures.Spatial and temporal variability of the soil moisture measurements and the depth-averaging procedure affected the values of the performance measures, and consequently the overall assessment of the GSDW model.However, the simulated soil moisture showed consistent behaviour with the different rainfall events, which verifies the validity of model results.
As expected, GSDW model results indicated the sensitivity of the top layer to rainfall events and other meteorological conditions, compared to the trimmed effect on the response of the subsurface soil layers.Generally, the GSDW model was capable of capturing the dynamics of the various water balance components in both reconstructed and natural watersheds.The canopy interception module, which is the main upgrade from the SDW model, allows the GSDW model to simulate the future performance of the reconstructed watersheds and long term climate scenarios.Furthermore, it allows users to compare different vegetation alternatives for future reclaimed covers.The developed GSDW model provides a vital tool, which enables the investigation of the utility of different soil cover alternative designs, hypothetical covers, and evaluation of their performance.Also, it facilitates further probabilistic analysis and scenario analysis, which provides the mining industry with a comprehensive decision support tool.

Fig. 4
Fig. 4 Simulated and observed moisture dynamics in the D3 cover for validation year 2006 using both SWAT and GSDW models; (a) peat layer; (b) till layer.

Fig. 3 .
Fig. 3. Simulated and observed moisture dynamics in the D3 cover for validation year 2006 using both SWAT and GSDW models; (a) peat layer; (b) till layer.

Interception Precipitation Surface Water Storage Overland Flow Interflow Infiltration ET Channel Flow 1 st Layer Storage Interflow ET 2 nd Layer Storage Moisture Redistribution n Layers Base Flow ET n th Layer Storage M. Red. Canopy Storage Snowfall Rainfall Evaporation Throughfall Snow Storage Snowmelt Fig. 1. A
schematic diagram of the GSDW model structure (Parasuraman et al., 2007)ted in 2001 of 0.2 m of peat/mineral mix overlying 0.8 m of till.Both the D-covers and Hill top site are located on South Bison Hill, a reclamation landform which is approximately 2 km 2 in area, rising 60 m above the surrounding landscape and has a large relatively flat top several hundred meters in diameter.The major plant species on the top of the SBH are foxtail barley (Hordeum Jubatum), and minor species include fireweed (Epilobium angustifolium)(Parasuraman et al., 2007).The third site is the South West Sand Storage (SWSS), which was constructed of 0.2-0.4m of Till/secondary cover material over 1.0 m of tailings sands.
It is currently the largest operational tailings dam in the world, approximately 40 m high and a 20 H:1 V side slope ratio.The vegetation varies with groundcover including horsetail (Equisetum arvense), fireweed (Epilobium angustifolia), and white and yellow clover (Melilotus alba, Melilotus officinalis).Tree species include Siberian larch (Larix siberica), hybrid poplar (Populus sp.hybrid), trembling aspen (Populus tremuloides), white spruce (Picea glauca) and willow (Salix sp.) The GSDW model is used to simulate the hydrological performance of a natural watershed to validate its capability in capturing the dynamics of the various water balance components in undisturbed systems.The old aspen (OA) forest site, part of the former Boreal Atmosphere Exchange Study (BOREAS), is considered in this study as it is climatically similar to the Fort McMurray region.The OA site is, roughly 1580 km 2 , located near the south end of Prince Albert National Park, Saskatchewan (53.629 • N, 106.198 • W).Because of the relatively flat topography and the homogeneity of vegetation, as well as the semi-arid climate, the site is modeled as a lumped unit or column of soil.The field instrumentation of the OA site has been providing continuous measurements since 1997 as part of the Boreal Ecosystem Research and Monitoring Sites (BERMS) program (http://berms.ccrp.ec.gc.ca).The soil is a well drained loam to clay loam.The top 0.1 m layer is an organic layer (leaf litter, plus fermentation layer); 0.07-0.3m of till mixed with sand and clay, a 0.45 m layer derived from gravely and clay enriched till, overlying, a mixture of sandy clay loam of 0.40 m.The soil properties are as follows: the saturated hydraulic conductivities are 25 (cm/day), 5.76 (cm/day), and 4.8 (cm/day), and the porosity values are 0.51, 0.45, 0.46 for A, B and C horizons, respectively ). www.hydrol-earth-syst-sci.net/13/865/2009/ Hydrol.Earth Syst.Sci., 13, 865-881, 2009 4.2 Natural watersheds

Table 1 .
Calibration parameter values for the developed GSDW model. is an exponent describing the influence of TI on soil defrosting; (b) c 1 evapotranspiration constant (mm/day • C) from the (1st) layer; and (c) λ i is an exponential coefficient, used to calculate the AET, modified from the SDW to be temperature dependant.

Table 2 .
Performance statistics of the GSDW and the SWAT models regarding soil moisture.

Table 3 .
Performance statistics of the GSDW model regarding daily AET flux.

Table 4 .
Annual water balance components for the validation years.