RESTORATION OF A DEGRADED BOG HYDROLOGICAL REGIME USING SYSTEM DYNAMICS MODELING

In scope of biodiversity and sustainable ecosystem development the restoration of the bog ecosystem is important, because by reducing the drainage effect on the bog, the negative impact on adjacent intact or relatively intact raised bogs and other wetland hydrological regimes is lowered. To restore a degraded bog hydrological regime, it is necessary to fill up the drainage ditches and cut out part of the forest stand. While researching scientific literature the author has obtained no evidence that there is a system dynamics model developed in order to simulate the tree cutting intensity in a degraded bog after filling up the drainage ditches with the aim to speed up the restoration of hydrological regimes, thus this approach is an innovative way of restoring the hydrological regime of degraded bogs. In previous studies the author concluded that the STELLA® system dynamics model is an appropriate tool to model the hydrological regime of bog. As a result of this research there is a STELLA® system dynamics model developed which through mathematical relationships helps to better understand the bog water cycle and to determine the consequences of any intervention on the bog ecosystem, primarily the effect of tree cutting. While running this STELLA® system dynamics model by changing the leaf area index, changes in the peat layer moisture level can be observed, which allows to predict the tree cutting intensity in order to reach the desired peat layer moisture level. By changing input data, this STELLA® system dynamics model could be used in other restoration projects of degraded raised bogs. UDC Classification: 574; DOI: http://dx.doi.org/10.12955/cbup.v6.1301


Introduction
After construction of drainage ditches, thus changing the natural hydrologic regime of bogs, many bogs in Latvia have been degraded, thus also leaving a negative impact on the adjacent unaffected and slightly affected hydrological regime of raised bogs and other wetlands. After drainage work, with groundwater level decreasing, trees begin to grow faster and larger in volume, filling the previous bog area and its surroundings The large volume of tree foilage, through the interception of rainwater, prevents it from reaching the ground because it evaporates back into the atmosphere. In addition, trees favor water absorption from the soil trough their roots further enhancing the drainage effect. To restore a degraded bog, first of all, it is necessary to fill up drainage ditches, but if the bog is overgrown with trees, this action may not be enough, and it is indispensable to cut down part of the trees. Experimentation with the intensity of forest thinning in reality would take years to establish the optimal scale of intervention in an ecosystem. This is because one would need to determine the intensity of forest thinning in order to reach the desired result but by cutting down as few trees as possible as well as to increase saturation of soil with water that would contribute to the resurgence of the bog and swampy forest biotopes. This is the reason why system dynamics modeling lwas chosen to solve this problem, because it allows to predict the influence of tree cutting on the ecosystem without expensive and time consuming experiments in the real life. Leaf area index (LAI) is one of the most important parameters influencing forest stand evaporation, rainfall, and snow interception. By running this STELLA model and changing LAI it is clearly visible how the forest stand interacts with the water level of the peat layer. This article describes what the components of the system dynamics model are and how these work mathematically..

Literature review
Studies of other authors' experiences have helped to determine elements which shape the bog ecosystem and interact with each other. Based on the exploration of elements forming the bog ecosystem system, a dynamics model in the STELLA® environment has been constructed. The two main articles which helped to understand the hydrological processes of soil are "Simulation of the hydrological processes on reconstructed watersheds using system dynamics" (Elshorbagy et al., 2007) and "System dynamics approach to assess the sustainability of reclamation of disturbed watersheds" (Barbour et al., 2005). This degraded bog system dynamics model uses the Aston (1979) interception equation which is noticed as an appropriate interception model in the publication "Different views on tree interception process and its determinants", because it includes the effect of rainfall characteristics (Klamerus-Iwan, 2014). Evaporation from the forest stand is based on the leaf area index (LAI) and is governed by logic described in the publication "A generic system dynamics model for simulating and evaluating the hydrological performance of reconstructed watersheds" (Carey et al., 2009). As it is a time consuming process to calculate LAI manually, a remote sensing approach was used as described the article "Derivation and validation of Canada-wide coarseresolution leaf area index maps using high-resolution satellite imagery and ground measurements" by Brown et al. (2002). Evaporation from open water depends on water vapor pressure differences and wind speed as Penman (1956) described in his publication "Estimating evaporation". The channel flow from the bog lake occur according to Manning's (1891)  o Leaf distribution coefficient.

Interception
Rainfall interception loss is the proportion of gross rainfall that is intercepted, stored and subsequently evaporated from all parts of vegetation during or following rainfall (Cui & Jia, 2014). Degraded bogs usually are overgrown with trees, which is why rain interception is significant. In this research Aston's (1979) rain interception equation is used (mm/d): where Irain is tree canopy interception (mm/d), Cp represents the function of canopy water coverage, Smax is maximum water capacity on the canopy (mm), k corresponds to the tree crown coefficient and gives the fraction of rainfall that falls on the canopy, P is precipitation amount (Klamerus-Iwan, 2014). = 0.935 + 0.498 − 0.00575 Smax is calculated based on the Hoyningen-Huene (1981) equation (Bulcock & Hewitt, 2010). Cp is determined based on the Gash et al. (1995) equation. k depends on leaf distribution and inclination angle, it values between 0.2 and 0.8 (Dijk et al., 2001). e is the mathematical constant 2.71828. This is an appropriate equation for the calculation of rain interception because it takes into account rainfall characteristics and it is possible to change the tree crown coefficient if using it in different degraded bog cases.
In Latvia there is snow in the winter time, therefore this model has the factor of snow taken into account. Snow interception occurs till the tree crown snow capacity is reached and is calculated based on the Hedstrom and Pomeroy equation, which is represented mathematically as: where Isnow is the intercepted water (mm/d) equivalent to the snow precipitation, PS is diurnal snowfall (mm/d), f is the snow interception coefficient . Maximum snow interception capacity B is expressed as: = (5) where LAI is the leaf area index of the bog forest stand and m is maximum snow coverage of the tree crowns (mm) which has to be determined empirically. The leaf area ratio Lr is a step function of temperature: IF Ta>-1 THEN (4) ELSE IF -1<=Ta<-3 THEN 1.5*Ta+5.5 ELSE (1) Interception efficiency decreases right after air temperature drops below -1 o C and is approximately constant for temperatures less than -3 o C, leading to a discontinuous relationship (Andreadis et al., 2009). The snow which is not intercepted falls on the ground and forms the snow cover. Mathematical changes in the snow cover are represented as: where St is the snow throughfall, ES is the snow sublimation and M is the daily snowmelt. Evaporation from the forest stand In this system dynamics model, evaporation from the forest stand is based on the leaf area index (LAI). = (7) = (8) SC is the canopy storage capacity and it is a subject of evaporation. SL denotes specific leaf storage (the depth of water retained by the leaf per unit LAI). EC is evaporation from the forest stand, Cp represents the function of the canopy water coverage (3) (Carey et al., 2009). To estimate the potential evapotranspiration PET (mm/d), Türc's (1961) equation is used: where Ta is the air temperature ( o C), and Rs is the total daily solar radiation (cal cm -2 d -1 ) (Minďáš et al., 2006). Leaf area index (LAI) The Leaf area index (LAI) quantifies the amount of leaf area bearing by a tree or the whole stand, normalized by the unit of crown projected or the whole ground area (Campbell & Norman, 1989). As plants use a fraction of the photosynthetically active radiation, it is possible to use remote sensing data to calculate the LAI for large areas in a short period of time. LAI for different tree species varies. An bogs, the main tree species are coniferous, it is appropriate to use the AVHRR and VEGETATION LAI algorithms as follows: = −16.32729 + 0.58909 × − 0.00754 × 2 + 4.57542 × 10 −5 × 3 − 1.30376 × 10 −7 × 4 + 1.40002810 −10 × 5 where SR is the simple ratio, which is calculated using corrected and normalized red (reflectance in red band), NIR (reflectance in near infrared band), BC is background SR value for coniferous trees, and D is the day of the year counting from January 1, when the remote sensing data were obtained (Brown, et al., 2002). Surface water In the winter time, surface is frequently covered by snow (Csnow). Snow can either melt and compose surface water, or sublimate in the atmosphere.
To calculate the sublimation of snow, Kuzmin's (1961) empirical formula is used. The daily sublimation is represented mathematically as: where ES is snow sublimation (mm/d), v is the wind velocity at 10 m height per second, xS is the saturated vapor pressure over the snow (mb), and x is the air vapor pressure (mb) at 2 m height (Gelfan et al., 2004). x is possible to be determined through using air temperature and a vapor pressure curve. Snow which is not sublimated melts and forms the overland water. The snowmelt is expressed by Anderson's (1976) degree-day factor: where M is the daily snowmelt (mm/d), Ta is the air temperature (°C), Tb is the base melt temperature (°C) and ρs is the snow density (kg/m 3 ). Water accumulated in the surface water storage is released to the peat layer as overland flow (Elshorbagy et al., 2007). The snow density is calculated using the Pomeroy et al. (1990) where ρs is the mean snow density (kg/m 3 ) corresponding to mean snow depth d (cm) (Essery, et al., 1998). Mathematically, the diurnal change in the surface water storage is represented as follows: where SSW is the surface water storage (mm), Rt is the rain troughfall (mm/d), M is the daily snowmelt, fP represents the infiltration rate in the peat layer (mm/d), OF is the overland flow (mm/d) and LF corresponds to fraction of water compiling the lake water. If in the bog there is a marsh lake, the water supply to the lake is calculated by multiplying the surface water by LAP which represents the lake percentage of the total bog area: In this case the peat infiltration rate must be multiplied by P AP representing the percentage of the bog area which is covered by soil: Peat layer Infiltration (fp) rate is dependent on the condition of the land surface, land vegetation cover, surface soil characteristics, storm characteristics, and surface soil temperature. The main factors influencing the water saturation of the peat layer are infiltration from the surface, water movement to the underlying till layer, peat layer evaporation, and water demand of the forest stand. The daily changes of the peat layer are represented as: where SP is the peat layer (mm), fT is the rate of downward movement of water into the till layer (mm/d), AETP is the actual evaporation from the peat layer (mm/d) and Ec is amount of the water used by the forest stand to provide for their biological processes. Before saturation, water infiltrates directly into the peat layer with the same intensity as rainfall. If the peat layer becomes saturated or the rainfall intensity exceeds saturated hydraulic conductivity, then infiltration is calculated by the Green-Ampt (Ampt & Green, 1911) equation (20). It estimates infiltration capacity fP based on cumulative infiltration volume FP (Barbour et al., 2005). Infiltration in the peat layer is represented as follows: where KsP is the saturated hydraulic conductivity of the peat layer (mm/d), MP is the initial moisture deficit (%), θsP is the saturated moisture content of the peat layer (porosity) (%), θiP is the initial moisture content of the peat layer before infiltration began (%), ψP is the wetting front suction head (capillary head in the peat layer (mm), FP is the cumulative volume of infiltration in the peat layer (mm) (Elshorbagy et al., 2007), and SmP is depth of the water already infiltrated in the peat layer (mm) (Feng, et al., 2011). SnP is the maximum water storage (mm) in the peat layer, and SrP is the minimum storage, in other words, residual moisture content that can be reached in the peat layer (Barbour et al., 2005).
To simulate infiltration into the frozen soil, the infiltration is multiplied by Li and Simonovic (2002) coefficient CtP which take into account soil defrosting and refreezing: where TImax (°C) is a maximum TI point at which the peat layer surface is fully defrosted, ci (dimensionless) is an exponent for describing the influence of TI and soil defrosting, N (days) is the number of continuous days with temperature below negative point, Nn is a maximum N after which TI will be lost and the peat layer surface will refreeze again, No is a logical variable to identify the day in which temperature is higher or lower than the active temperature (Li & Simonovic, 2002). This model accounts for the accumulation of positive temperature as well as the accumulation of the number of days in which the temperature falls below zero. Parameters ci, and TImax are estimated during the calibration process (Elshorbagy et al., 2007). To calculate the number of continuous days with temperatures below the negative point (N) there must be a stock built where an input is NO and in output is the negative temperature day counter: To calculate the point at which the peat layer surface is fully defrosted (TI) there must be a stock built where an input is air temperature (Ta) and in output is the positive temperature counter:

IF T a >0 AND N<N n THEN 0 ELSE SUM T a
The evaporation from the peat layer (ETP) is estimated using Li & Simonovic (2002) empirical formulations, which is represented as follows: where ETP is the peat layer evaporation (mm/d), cp is the evaporation constant (mm/d°C) and is determined during calibration, SmP (23) is the effective moisture saturation of the peat layer (dimensionless), λ is an exponential coefficient estimated while calibrating. Till layer Moisture movement in the till layer depend on the moisture content in the overlying peat layer, soil suction in the peat and till layers, evaporation from the till layer, and percolation to the underlying shale layer. The daily changes of the till layer are represented as: where fT is the rate of downward water movement into the till layer (mm/d), InT is the till layer interflow rate (mm/d), fS is the shale percolation rate (mm/d), ant ETT is the evaporation rate (mm/d). Percolation to the lower layer is dependent on the water saturation within the surface and subsurface soil layers (Li & Simonovic, 2002). There is no movement of water from the peat layer to the till layer if the moisture content in the peat layer is less than the wilting point, and if the field capacity of the peat layer is greater than the field capacity of till layer. Water will start to percolate when the moisture content in the peat layer is greater than the residual moisture content (Elshorbagy, 2007). In the case of a fully saturated peat layer, water will percolate into the till layer following Green-Ampt's (Ampt & Green, 1911) equation (30), and (21), (22), (23), replacing subscript P (for peat) with T (for till).
In the case of a saturated till layer, the maximum rate at which water can be absorbed by the till layer will correspond to the saturated hydraulic conductivity (Elshorbagy et al., 2007). Water movement in the till layer is estimated as follows: where D P is the depth of the peat layer, fT is the rate of downward water movement into the till layer (mm/d), Δt is the iteration time interval, and IcT is the coefficient of downward water movement into the till layer, KsT is the saturated hydraulic conductivity of the till layer (mm/d), ψT is the pressure head in the till layer, MT is the initial moisture content of the till layer (%), and FT is the cumulative volume of the infiltration into the till layer. When the peat layer is frozen, there is no water movement from the peat layer into the till layer, therefore the CtP coefficient is used. Evaporation from the till layer (ETT) occur only when the peat layer temperature (Ts) is greater than 0°C and the moisture level in the till layer is higher than wilting point and is calculated as: where ETT is evaporation rate (mm/d) from the till layer, cT is the evaporation constant (mm/d°C) from the till layer and is estimated during the calibration process. The effective moisture saturation in the till layer SmT is calculated using equation (23), replacing subscript P (for peat) with T (for till). The interflow (subsurface runoff) of the till layer is one of the ways how water leaves the bog ecosystem. The interflow occurs only when the till layer is saturated, and it is estimated as follows: where ST is the till layer storage (mm), θsT is the porosity of the till layer, DT is the depth of the till layer, and CiT is the interflow coefficient (Barbour et al., 2005), which is determined during the calibration. Shale layer Moisture change in the shale layer depends on percolation from the till layer and interflow (Barbour et al., 2005). The diurnal change in the shale layer (mm/d) is represented as: where fS is the infiltration rate of the shale layer (mm/d), and InS is the shale layer interflow rate (mm/d).
The till layer will release water only when the moisture content in the till layer is greater than the wilting point moisture content of the till layer (Elshorbagy et al., 2007). In case of a fully saturated till layer, water will percolate into the shale layer following the Green-Ampt (Ampt & Green, 1911) equation. Water movement from the till layer to the shale layer is estimated as follows: where I cS is the coefficient of shale percolation and is determined during calibration, and θiS is the moisture content in the shale layer (%) (Elshorbagy, 2007). The interflow of the shale layer occurs when the shale layer is saturated and is calculated according to the following logic: where θsS is porosity of the shale layer (%), DS is the depth of the shale layer (mm), and CiS is the interflow coefficient which is determined during the model calibration.

Lake evaporation
The lake evaporation is governed by the Penman equation which is expressed as: = 0.47(0.5 + 0.01 )( − ) (38) where EL is the evaporation rate of the lake (mm/d), v is the wind velocity at 2 m height above water surface (mi/d), xS is the saturated vapor pressure above the water surface (mb), and x is the air vapor pressure at 2 m height (Penman, 1956). As the lake water temperature is not always available, it is possible to use the relative humidity (RH) formula RH=x/xW*100 to calculate the saturated vapor pressure as follows: = ⁄ * 100 (39) Channel flow In some occasions the water from the bog lake also leaves the ecosystem through channel flow not in addition to evaporation. In this case channel flow is governed by Manning formula (40) (Manning, 1891). The channel flow occurs until the lake water reaches a low water level that the lake water does not reach the bottom of the channel and the water outflow cannot continue. The channel flow is estimated as follows: IF SL≤SLmin THEN 0 ELSE (40) where Q is the channel flow (m 3/ d), n is the Manning`s roughness coefficient, A is the channel cross section area of flow (m 2 ), R is the channel hydraulic radius (m), and S is the channel slope (m/m) (Froehlich & Jobson, 1988). To convert m 3 to mm the Manning formula (40) must become a STELLA® graphical function.
Overland flow Ponded water can occur when rainfall cannot infiltrate into the soil fast enough. Ponded water is routed downhill as surface runoff (Butts & Graham, 2005). Water flow on the ground surface is calculated as: IF SP/DP> θsP THEN (41)*PAP ELSE 0 where OF is overland flow (mm/d) and CSlope is the slope coefficient which depends on the slope value.

Conclusions and discussions
While running this STELLA® system dynamics model by changing the leaf area index, changes in the peat layer moisture level can be observed. If the preferred peat layer moisture level is known, the leaf area index (LAI) must be reduced until this level is reached. The next step is to develop a methodology to convert LAI to a unit of trees per one bog ground area which would allow us to determine the right amount of trees to be cut down. By changing input data, this STELLA® system dynamics model could be used in other restoration projects of degraded raised bogs that have similar weather conditions and soil properties.