A Hydrological Concept including Lateral Water Flow Compatible with the Biogeochemical Model ForSAFE

Abstract

The study presents a hydrology concept developed to include lateral water flow in the biogeochemical model ForSAFE. The hydrology concept was evaluated against data collected at Svartberget in the Vindeln Research Forest in Northern Sweden. The results show that the new concept allows simulation of a saturated and an unsaturated zone in the soil as well as water flow that reaches the stream comparable to measurements. The most relevant differences compared to streamflow measurements are that the model simulates a higher base flow in winter and lower flow peaks after snowmelt. These differences are mainly caused by the assumptions made to regulate the percolation at the bottom of the simulated soil columns. The capability for simulating lateral flows and a saturated zone in ForSAFE can greatly improve the simulation of chemical exchange in the soil and export of elements from the soil to watercourses. Such a model can help improve the understanding of how environmental changes in the forest landscape will influence chemical loads to surface waters.

1. Introduction

Forests regulate the availability and quality of water resources. They generally reduce the total runoff through canopy interception and transpiration. Forests also change the velocity of the flow by enriching the soil with organic matter and changing the structure of the soil with their root system. By reducing the runoff and altering nutrient cycles, forests also change the quality of the soil water and runoff [1,2,3]. Therefore, environmental changes affecting forest dynamics, such as climate change or management changes, also affect water resources in forest catchments and eventually the loads of chemicals in inland and coastal waters [4,5,6,7,8].

Biogeochemical models can be used to evaluate the effects of environmental changes on forest ecosystems [8,9,10,11,12]. They can include the simulation of hydrology of the system, but the level of complexity of the water flows varies significantly among models [13,14]. To succeed in simulating how nutrients are transported to rivers and surface waters, biogeochemical models need to include vertical and lateral flows in the soil.

The biogeochemical model ForSAFE is a mechanistic model designed to simulate the dynamic responses of forest ecosystems to environmental changes [11,15,16,17]. The hydrology module in ForSAFE was designed to simulate processes at the forest site level. It included the vertical flow of water in the soil profile but not lateral water flows.

The aim of this paper is to present the hydrology module developed to include lateral water flow, with the future aim of including it in the ForSAFE model. The inclusion of lateral flow will allow the simulation of water and chemistry transport from the forest ecosystem to the stream. In addition, it will simulate a saturated zone which will better represent moisture content in the soil and the processes in deeper soil layers as well as in the riparian zone. Moreover, in the current ForSAFE model, the simulation of dynamics linked to fast hydrological response is limited by the fact that processes are simulated on a monthly time step. For this reason, the new hydrology module simulates dynamics at the daily time step. The hydrology module was calibrated and evaluated against data collected at Svartberget in the Vindeln Research Forest in Northern Sweden. The scope of the test in Svartberget was a proof of the concept rather than providing a best-fit to the measurements of water storage and flows by optimizing model parameters.

2. Model Description

The improved hydrology concept including lateral flow was designed to be compatible with the current structure of the ForSAFE model.

ForSAFE is a mechanistic model of the dynamics of forest ecosystems that combines chemical, physical and physiological processes [11,15]. It integrates the three basic material and energy cycles in a single model: the biological cycle representing the processes involved in tree growth; the biochemical cycle including uptake, litter decomposition and soil nutrient dynamics; and the geochemical cycle including atmospheric deposition and weathering processes.

In ForSAFE, each modelled forest unit includes a single soil column with different layers denoting the soil horizons. At present, the soil moisture content is calculated on a monthly basis for each of the defined soil layers. The hydrology module only simulates the vertical flows of water, including evapotranspiration and percolation. To include lateral water flow, it was necessary to represent the differences of soil properties not only with depth, but also along the slope. This was achieved by arranging several soil columns next to each other to discretely represent the vertical and horizontal structure of the soil in a forest transect. In order to test and better illustrate the new concept, the hydrology module was developed in the visual modeling environment STELLA (Figure 1) [18].

Figure 1. Simplified illustration (considering only two soil columns) of the hydrology as conceptualized in STELLA. The boxes represent water stocks in elements storing water, the blue arrows represent water flows from and to water stocks, the stand-alone circles are converters containing data or equations and the red arrows are connectors indicating which factors are regulating the flows. The striped boxes represent conveyors that transfer water after a specified number of time units.

The main changes compared to the original module in ForSAFE are:

- increased time resolution from monthly to daily time steps

- the water flow from a given soil layer is constrained by the amount of water that receiving layers can accept vertically and horizontally

- the water flow is given a velocity regulated by the soil conductivity, controlling the amount of water that can move within each time step (per day)

- inclusion of water movement along a slope, i.e., surface runoff and lateral flow

- the soil hydraulic properties are assessed as a function of soil texture

The following sections describe in detail how the water flows and stocks are calculated in the new hydrology module as compared to that in ForSAFE.

2.1. Water Inputs to the Soil

Water enters the system as precipitation on a daily time step (Figure 2, Water inputs). When entering the soil, the precipitation (Prec, m∙d⁻¹ of water) goes first through the same snow routine as in the ForSAFE model. The precipitation is stored entirely in a snowpack (Snowin, m³·d⁻¹ of water equivalents) when the average daily air temperature (Tm, °C) is below −5 °C. Between −5 °C and 2 °C, it is assumed that water inputs are a mix of snowfall and rainfall and only the snowfall is stored in the snowpack. The fraction of snowfall is reduced linearly with increasing air temperature until the 2 °C point, when all the precipitation is assumed to be in the form of rainfall [19].

_Snowin=MAX(0, MIN(_Tm-2-7,1))_{× Prec×A1,j}

where A_1,j is the area of the first layer in soil column j. Generally, in the model, the first index of a variable (i) represents the number of the layer, starting at the soil surface and the second index (j) the number of the column, starting from the column next to the stream.

Figure 2. Illustration in STELLA of the main processes regulating water flows in the model: water inputs to the soil, percolation, lateral water flow and transpiration from a soil layer. The connectors (red arrows) indicate which factors are included in the calculation of the flows.

When the average daily temperature is greater than 0 °C, the water stored in the snowpack starts melting and enters the soil. The melting rate is determined by a degree-day factor (CMELT, m·°C⁻¹·d⁻¹) [19]. The sum of rainfall and snow melt at each time step (Surf H2O, m³ water) corresponds to the maximum amount of water that can infiltrate in the upper soil layer. The actual infiltration in the first layer of each soil column j (Infilt_1,j, m³·d⁻¹ water) is limited by the level of saturation of the soil and cannot be greater than the volume of empty soil pores in the upper layer.

_Infilt1,j=MIN(SurfH2O,(_{POR1,j -Moist1,j})_×Vol1,j)

where POR_1,j (m³ water m⁻³ soil) is the volume of pores per unit of soil volume, Moist_1,j (m³ water m⁻³ soil) is the volume of water per unit of soil volume and Vol_1,j (m³ soil) is the volume of the upper layer. The total volume of a soil layer is given by the thickness (Δz_i,j, m), the length (Δx_i,j, m) and the width of the layer (Δy_i,j, m). The impediment created by the coarse fragments is considered by reducing the total volume proportionally to the coarse fraction (CF_i,j):

_Voli,j = Δ_zi,j×Δ_xi,j×Δ_{yi,j ×(1-CFi,j})

As a consequence, the presence of coarse fragments such as boulders or stones reduces the amount of pores in the soil and the water volume that can be retained in each soil layer.

The water that cannot infiltrate because the soil is saturated becomes surface runoff, and is directed towards the next column downhill or, in the case of the last column downstream, to the stream.

2.2. Water outputs from the soil

The water can leave each soil layer through three processes:

1. Percolation

2. Lateral flow

3. Transpiration

The percolation and lateral flows are, respectively, the vertical and horizontal water movements between soil layers, while transpiration is the water extracted by plants through their roots. The flows in the model are not regulated by a gradient of potential as determined by the laws of physical hydrology. The amount of water flowing in the soil at each time step is dependent on the hydraulic properties of the soil layers, as explained in the following sections.

The hydraulic properties of the soil used to assess the water outputs are:

1. the porosity, the fraction of pores in the soil volume (POR, m³ water m⁻³ soil)

2. the field capacity, the water content held in the soil when free drainage by gravity has stopped (FC, m³ water m⁻³ soil)

3. the permanent wilting point, the soil moisture content at which plants cannot extract more water from the soil (WP, m³ water m⁻³ soil)

4. the saturated hydraulic conductivity, the rate at which water moves through the pores in the saturated soil (Ksat, m·d⁻¹)

5. the unsaturated hydraulic conductivity, the rate at which water moves through the pores in the unsaturated soil (Kh, m·d⁻¹)

6. the slope of the soil moisture characteristic, which is the relation between water tension and volumetric water content (λ)

The soil hydraulic properties are calculated for each soil horizon from the texture and the density of the soil, based on pedo-transfer functions [20,21,22].

All the hydraulic properties are representative for the volume occupied by the fine earth, i.e., soil particles smaller than 2 mm in diameter and pores.

2.2.1. Percolation

The downward vertical flow is known as percolation and it concerns the volume of soil water above field capacity, i.e., the water that can be moved by the force of gravity.

As in ForSAFE, the maximum volume of water that can percolate from a soil layer (Perc_i,j, m³·d⁻¹) is the water content above field capacity [19] (Figure 2, Percolation). In the new hydrology module, the percolation from layer i to layer i + 1 in soil column j is given by:

_Perci,j={MIN((_{Moisti,j -FCi,j})×MIN(_Voli,jdt_{, Khi,j ×Avi,j})_{, EPi+1,j}), _{Moisti,j >FCi,j}0,_Moisti,j≤_FCi,j

where Moist_i,j (m³ water m⁻³ soil) is the soil moisture content, FC_i,j (m³ water m⁻³ soil) the field capacity, Vol_i,j (m³) the volume of the percolating layer, Kh_i,j (m∙d⁻¹) the unsaturated hydraulic conductivity and Av_i,j (m²) the cross-sectional area in the direction of the vertical flow. EP_i+1,j is the capacity of the underlying layer given by the volume of empty pores at each time step (dt).

_EPi+1,j = MAX((_{PORi+1,j -Moisti+1,j})×_Voli+1,j, 0dt)

The unsaturated hydraulic conductivity (Kh, m·d⁻¹) is calculated according to the function by Saxton and Rawls [21].

K_hi,j=(Ksa_ti,j×^{(Mois_ti,jPO_Ri,j)3+2λ})

where λ is the inverse of the slope of logarithmic tension-moisture curve:

λ=ln(F_Ci,j)− ln(W_Pi,j)ln(1500) − ln(33)

The percolation at the bottom of the modelled soil column can be divided into a fraction reaching the stream and a fraction that does not contribute to the streamflow (Figure 1). In addition, a transit time can be assigned to the bottom water flow, which can be used to better represent the time required to reach the stream.

2.2.2. Lateral Water Flow

In this model, the water flows both vertically and laterally downhill to the stream. When water above field capacity cannot percolate downwards because the layer below is saturated, it moves laterally to the next soil column downslope (Figure 2, Lateral Flow). That is, in the saturated zone and in the layers above the saturated zone, the flow is mainly horizontal. The vertical flow can also be constrained by a low permeability at the bottom of the soil column. This is simulated by reducing the conductivity that regulates the percolation from the deepest layer. The way flows are regulated in the model does not account for conditions when the flow is directed first laterally than vertically. Following Darcy’s law, this would be the case when differences of potential are greater laterally than vertically, such as when the soil downhill is much drier or has a much finer texture than the deeper soil layers.

The lateral flow (Lat_i,j , m³·d⁻¹) from a layer i in column j to a layer i in column j−1 (downstream) is given by:

_Lati,j={MIN((_{Moisti,j -FCi,j})×MIN(_Voli,jdt_{, Ksati,j}×_Ahi,j×Slp)_{, EPi,j-1}), _{Moisti,j >FCi,j}0,_Moisti,j≤_FCi,j

where Ksat_i,j (m∙d⁻¹) is the saturated hydraulic conductivity, Ah_i,j (m²) the cross-sectional area in the direction of the horizontal flow and Slp the slope (m∙m⁻¹). EP_i,j−1 is the capacity of the receiving layer given by the volume of empty pores in the next column downslope at each time step:

_{EPi,j-1 =MAX(0,(PORi,j-1 -Moisti,j-1 )×Voli,j-1})

2.2.3. Transpiration

The assessment of the transpiration (T_i,j, m³·d⁻¹ of water) from a soil layer i in soil column j is given by (adapted from [19]):

_Ti,j=PET×R_Fi,j×(MIN(MAX(0, Mois_ti,j−W_Pi,jL_Pi,j−W_Pi,j), 1))

where PET is the potential evapotranspiration, RF_i,j the root fraction, LP_i,j is the limit for potential evapotranspiration and WP_i,j the permanent wilting point of the soil layer (Figure 2, Transpiration).

The PET in ForSAFE is driven by tree photosynthesis, which in turn is dependent on nutrient availability. Since the tree growth processes and nutrient cycling are not modelled in this study, we adopted a simplified approach to estimate PET by calculating it as a function of temperature with the HBV model [23]. As in ForSAFE, the potential plant transpiration is constrained by water availability in the soil. In addition, the transpiration equals the PET above a critical point identified by the limit of evapotranspiration, LP_i,j which ranges between 50%–80% of the field capacity of the layer [24] The water availability is calculated as a function of the relative water content in the root-zone [19,24]. Moreover, the water that can be extracted from each layer is dependent on the amount of roots in that layer. In this study, the root fraction is assessed according to the data reported by Rosengren and Stjernquist [25] for different tree species. It is also assumed that the tree roots extend to soil layers up to 1 m depth.

Capillary rise due to surface tension is not considered in the model and this could limit the water available for trees when dry conditions persist for longer periods.

3. Model Test

3.1. Site Description

The STELLA model was tested on the S-transect in the Vindeln Research Forest (64°14’N 19°46’E) in northern Sweden [26] by simulating water flows and soil moisture on a daily basis.

The S-transect is located at 250 m a.s.l. in the Krycklan catchment which covers a surface of about 68 km² (Figure 3). The transect is aligned parallel to lateral flow paths towards the Västrabäcken stream. Measurements of soil moisture, groundwater level and soil water chemistry have been collected from 1997 [26,27]. The mean air temperature for the period 1981–2012 was 1.8°C, the average annual precipitation 631 mm and the average runoff 308 mm [26]. About half of the runoff occurs during a few weeks after the snowmelt. The vegetation is dominated by a mature stand of Norway spruce close to the stream and Scots pine upslope. The bedrock is gneiss and it is overlaid by glacial till, which is about 10–20 m deep in the lower part of the Svartberget catchment and 30–40 m in the middle part of the catchment. Nearest to the bedrock there is dense basal till that is overlaid by a less dense melt-out till of about 1–3 m depth [28]. The soil in the S-transect is an organic rich Histosol within 10 m of the Västrabäcken stream, with well-developed Podzols further upstream.

Figure 3. Location of the S-transect in Sweden and in the Krycklan catchment.

3.2. Transect Model

The soil transect was modelled as series of soil columns (C) representing the vertical and horizontal structure of the soil from the stream to the water divide.

The transect was designed as seven soil columns of different length (C1 to C7) and 1.5 m depth to represent discretely the change of soil type from the riparian zone to the mineral soil upstream (Figure 4). Data on grain size distribution at different points along the transect were used to define the texture of the seven soil columns, each divided in seven layers ([29], Table S1).

Figure 4. Representation of the S-transect in the hydrology model when the water divide is assumed at 80 m from the stream. C: modelled soil column. S: point of measurement. The number of the sampling points represents the distance from the stream (4, 12 and 22 m).

Measurements of daily precipitation, average temperature and grain size distribution were used as inputs to the model. Model parameters and results were tested against hydrological measurements taken at three different stations along the transect which are located at 4, 12 and 22 m from the stream (S4, S12, S22) (Figure 4).

3.3. Parameterization of Soil Hydraulic Properties

Measurements of soil hydraulic properties are not commonly available. Therefore, the model is structured to calculate them based on soil texture—an input required also by the ForSAFE model- and the pedo-transfer functions [20]. In this study, the soil hydraulic properties were assessed for the different layers of all soil columns based on texture data collected along the entire transect (Table S2).

The functions were tested against measurements available along the transect. Measurements of hydraulic properties were available only at the three sampling point S4, S12 and S22 at selected depths. Porosity and the soil moisture characteristic (pF curves) were measured on cylinder samples of the soil using the pressure plate method, except for the 150 m tension for which a pressure chamber was used [27]. The saturated hydraulic conductivity was measured with a double-ring infiltrometer [30]. Texture data from the same or comparable sampling points were used to estimate the hydraulic properties with the pedo-transfer functions.

The calculated porosity is generally in good agreement with the measured data (Figure 5, POR). Some discrepancies are observed in organic soil layers where the porosity is high due to the high organic matter content [31,32,33]. A possible reason for this difference is that the volumetric water content can vary substantially within short distances at the interface between organic and mineral soil, due to the large difference in physical properties between those two types of soil. When the estimated porosity in a layer is higher than the actual, water content at saturation will be overestimated by the model. This could have the effect of overestimating the water stored in the soil, thus underestimating the streamflow. The opposite will happen when the calculated porosity is underestimated.

Figure 5. Testing of estimated hydraulic properties against measurements at the sampling points S4, S12, S22 of the S-transect. POR: porosity; WP: permanent wilting point; Ksat: saturated hydraulic conductivity. The dotted lines and their slope are an indication of the discrepancy between estimates and measured data. When the slope is > 1 the measured data are higher than the estimates and when <1 measurements are lower than estimates.

The estimated permanent wilting point is also comparable to measurements (Figure 5, WP). When different, the water volume retained at the wilting point is a low fraction of the soil volume. Therefore, the difference in terms of water content is small and it should have a minor impact on the estimate of water storage and fluxes.

Greater discrepancies are found between the calculated saturated hydraulic conductivity and measurements (Figure 5, Ksat). In upper layers, the calculated Ksat is generally higher than the empirical data and lower in deeper layers. Therefore, in agreement with other studies [34,35], the results indicate that the estimated Ksat decreases more rapidly with depth and is higher in upper layers than measurements. These discrepancies will have an effect on the velocity of the flow: when the estimated Ksat is lower than the actual, water will reach the stream later than in reality. Both calculated values and measurements indicate that layers with high organic matter (>10%) have a lower saturated hydraulic conductivity.

3.4. Sensitivity Analysis

A sensitivity analysis was performed to evaluate a set of model parameters that could influence the model outputs. These were also the parameters that were considered for the model calibration. A one-factor-at-the-time method (OFAT) was used to test the model sensitivity to parameter perturbation [36]. Two objective functions were used in the analysis: Bias and Root Mean Square Error (RMSE) (Table S3). The functions evaluate the discrepancy between model results when parameters are varied within a certain range relative to a baseline value (Table 1). The ranges were defined based on available information from the site or literature, as explained in the following paragraphs.

[ Table omitted. See PDF. ]

A Lidar-based digital elevation model was used to define the hillslope up to the water divide (WD). With the 1 m DEM, the transect extended 80 m to the water divide, while using a 5 m DEM realigned the transect slightly and pushing the water divide to 140 m [37]. The two different estimates of the water divide compared to a mean value of 110 m were used to test the effect of changes of hillslope length on the model results.

The model simulates the hillslope hydrology up to 1.5 m depth. Without further assumptions the percolation at the bottom of the column is only regulated by the calculated Ksat of the bottom layer (Kbott). As a consequence, the lower permeability in soil deeper than 1.5 m would not be considered in the simulation. The effect of changing the permeability of the bottom layer was tested in the sensitivity analysis by varying the Kbott within a range between zero (e.g., impermeable rock) and a value equal to the Ksat of the bottom layer (no further impediment at the bottom of the column).

A further unknown factor is the type of connection existing between the water percolating from the bottom of the columns and the river (FrStream). When evaluating the model, the simulated streamflow is compared to a measured value taken at a certain point along the river. Ideally all the water from the simulated area discharges at the point of measurement (FrStream = 1), but part of the deep flow can also discharge in the area downstream from the point of measurement [38]. Since the transect is aligned parallel to lateral flow paths, most of the bottom percolation should contribute to the stream flow measured at the dam downstream. In the sensitivity analysis the fraction to the stream was varied between 0.8 and 1. As a baseline scenario, we considered a linear decrease of the fraction to the stream from 1 in the column next to the river to 0.8 at the water divide.

The change of the degree-day factor was used to test the sensitivity of the model to the simulation of snow melt (Cmelt). In ForSAFE a value of 1.5 mm °C⁻¹·d⁻¹ was used in several studies in Sweden [39]. A variation of ±0.5 mm °C⁻¹·d⁻¹ was included in the sensitivity analysis.

The effect of transpiration was evaluated by varying the limit for potential evapotranspiration, LP, that is calculated as a fraction of the FC, ranging from 0.5 to 0.8 of FC [24,40].

Other factors that can affect model results are the estimated hydraulic properties, which determine the amount of water stored and flowing in the soil. In the sensitivity analysis we tested the effects of varying the porosity (POR) and the saturated conductivity (Ksat) of ±10% of the calculated value. According to the pedo-transfer functions the change of porosity affects also the field capacity and the wilting point (on average ±12% and ±14%, respectively).

Model results are also affected by the uncertainty of the input data, such as texture, PET and climate data, but these effects were not included in the sensitivity analysis.

We evaluated the model sensitivity by analyzing the modelled water storage and streamflow (Table 2 and Figure 6). We considered three different components of the water storage: the snowpack, the saturated and the unsaturated storages in the soil. We also analyzed the changes in daily streamflow (Q) and the high and low peak values by comparing the 90% and 10% percentile values of daily streamflow.

Figure 6. Change of water storage and streamflow components when parameters are varied compared to a baseline. The change is expressed in terms of Bias and RMSE.

[ Table omitted. See PDF. ]

A reduction of the conductivity regulating the bottom percolation has the most significant impact on both the simulated water storage and the streamflow. By reducing the Kbott, it is possible to increase the saturated storage in the soil, i.e., simulate a higher groundwater level. In addition, a reduced permeability at the bottom significantly changes the streamflow, which is reduced on average on a daily basis, but increased at high flow.

To a minor extent, the water storage is also affected by the calculated porosity: with higher porosity the soil is on average more unsaturated.

Other relevant factors for the simulation of streamflow are:

- the fraction of bottom percolation to the stream (FrStream): when increased the overall streamflow increases, as indicated by the Bias;

- an increase of snow melting factor (Cmelt) affects significantly the distribution of the flow during the year (increase of RMSE), but has a minor effect on the overall water balance.

3.5. Calibration

As a following step, the model results obtained when applying baseline parameter values were compared against measured data describing soil moisture and streamflow. These included dynamic soil moisture measurements with TDR and the streamflow data at the Västrabäcken stream water site (site 2, [26]) from the period 2011–2014. Soil TDR were available at 5, 12 and 22 m from the stream which are comparable to columns C2, C3 and C4 in the model. The TDR data of the second and fourth layer were available for all the three measurement points and therefore compared to the modelled results (Table 3). The discrepancies between modelled results and measurements were compared in terms of Bias, RMSE and normalized mean error (NME) (Table S3). Almost all the parameters are significantly underestimated. By adopting baseline parameter values, the soil is most of time unsaturated and the streamflow on a daily basis and at high peak flows lower than measurements. The only parameter that is greatly overestimated is the 10% percentile of the streamflow, i.e., the low flow values.

[ Table omitted. See PDF. ]

This comparison suggests that a high conductivity at the bottom of the soil columns compromises the simulation of a saturated zone and releases a constant base flow which is much higher than reality.

Therefore, in the calibration phase, the Ksat of the bottom layer was reduced by two orders of magnitude. A further decrease of Ksat was necessary in column C4 (15–25 m from the stream) where the calculated Ksat was at least twice as large as in the other columns. As a result, a total decrease of 2 × 10⁻² was applied n C4 and the Ksat at the bottom of the modelled transect was reduced from 0.33 m·d⁻¹ to 0.003 m·d⁻¹ on average. The adjustment was necessary to avoid an excessive loss of water through bottom percolation that would compromise the simulation of deep saturated layers, i.e., of the groundwater level. We justify this choice by assuming that a very low conductivity at the bottom of the column would include the lower permeability in soil deeper than 1.5 m that is not considered in the simulation.

The results from the calibration are evaluated in the following sections by analyzing the components of the water storage and of the water flows. The results on the dynamic soil moisture and runoff are compared against the TDR and streamflow measurements on a daily basis.

4. Model Results and Discussion

4.1. Simulated Water Storage and Fluxes

The runoff (R) in the entire simulation period is 36% of the precipitation (P) (Table 4). A smaller share of the runoff (11%) is simulated as deep percolation and does not reach the stream, while most of it contributes to the stream flow. Transpiration plays a significant role in the water balance by using about 60% of the total precipitation. Over the entire simulation period, the water storage in the soil increases by 4%. The increase of storage mostly happens in the first year of simulation, because it is assumed that all layers have an initial water content equal to field capacity. Therefore, a period of stabilization of the model is required to reach the saturation level in deeper soil layers.

[ Table omitted. See PDF. ]

Figure 7 illustrates the magnitude and the changes of the water storage over time. The storage is divided into storage in the saturated and unsaturated zones and in the snow pack. The period of stabilization required to reach the saturation level in deeper soil layers is clearly shown by the absence of a saturated zone at the beginning of the simulation. The stabilization ends some months after the first snowmelt, when the soil starts receiving water inputs. In the following years (2012–2014), the model simulates a seasonal change of the storage which reflects the accumulation of water in the snow pack, the increase of soil water storage after snow melt, both in the unsaturated and saturated zone, and the decrease of soil moisture in summer, due to plant transpiration.

Figure 7. Water storage in the saturated and unsaturated zone in the soil and in the snow pack over the simulation period.

4.2. Soil Moisture

Figure 8 compares the modelled soil moisture to measured TDR data at different depths at the three sampling points S4, S12, S22. The figure shows the simulation of soil moisture in the period 2012–2014, i.e., excluding the period of stabilization.

Figure 8. Modelled soil moisture as compared to the TDR data at increasing distance from the stream. The three modelled soil columns up to 25 m from the stream (C2, C3, C4) are compared to the measurements at the sampling points at 4, 12 and 22 m (S4, S12, S22). For each soil profile, simulations and measurements are compared at different soil depths. Modelled layers of similar depths are compared in the same row. Measurements were not available for layers that are missing. Differences between model and measurements are evaluated in terms of NME (%) and Bias (m3·m−3). Dashed grey lines represent the water content at field capacity.

The model simulates permanent saturated layers at depths similar to measurements and levels of soil moisture which are comparable to the TDR data.

However, some of the simulated moisture dynamics differ significantly from the measurements. The TDR measurements in the upper layers show a decrease of moisture in winter that reaches values below field capacity (e.g., in layers up to 20 cm depth). In the same period, the model predicts water contents at the lowest close to field capacity, since transpiration does not occur in winter. The decrease of the TDR values is mainly caused by soil and water frost in the upper soil layers, which is not shown in modelled results, since frozen water is not simulated by the model.

Other discrepancies between model and measurements can be attributed to different soil hydraulic properties at the points of measurement and in the modelled soil layer. Due to its structure, the model assumes that all soil properties are homogeneous at any given point within a soil layer, as compared to measurements, which represent a specific point in depth and in distance from the stream. When modelled and measured saturation levels differ, the calculated porosity for the entire layer is most likely different from the porosity at the point of measurement. Moreover, the parameterization of the soil hydraulic conductivity at the bottom of the soil columns can result in soil permeability different from the actual one in the deep soil. The difference can produce some of the discrepancies between the fluctuations of modelled soil moisture and of the TDR data in deeper layers.

4.3. Streamflow

The comparison of the modelled and measured streamflow shows that the model tends to overestimate the annual water flow at the stream (Figure 9). In 2013–2014, the cumulative modelled streamflow over time is about 10% higher than measurements.

Figure 9. Modelled as compared to measured streamflow. Differences between model and measurements are evaluated in terms of NME (%) and Bias (mm). The graph on the right compares the cumulative streamflow in 2013-2014 when measurements were available.

The overestimation of the annual streamflow is mainly caused by the modelled base flow which is higher than the measured one: about 0.60 mm·d⁻¹ against the observed 0.25 mm·d⁻¹ in winter. As highlighted by the sensitivity analysis, the magnitude of the modelled base flow is influenced by the assumptions made to model the bottom percolation: a reduced Ksat at the bottom of the soil columns have the effect of reducing the overall streamflow and the 10% percentile value. In addition, the model does not simulate soil and water frost, as shown by the mismatch between measured and modelled soil moisture in winter. The unfrozen soil water could contribute to the winter flow and thereby overestimate the streamflow in winter.

The results presented in Figure 9 also highlight that the model tends to underestimate the peak flows in spring. Also in this case, the high flow values are mostly affected by the Ksat at the bottom of the columns: a reduction of Ksat increases the 90% percentile value of the streamflow. Other factors that could contribute to the seasonal differences between modelled and measured streamflow are a snow routine that does not fully capture the processes of snow accumulation and melting, the parameterization of plant transpiration and the lack of preferential flow paths in the model.

As shown in the sensitivity analysis, the parameterization of snow melting mainly changes the distribution of the flow and affects the peak flows: a faster snow melt anticipates and increases the peak flows by concentrating the release of water from the snow pack in a shorter time. Similarly, it can be expected that by modifying the parameters regulating snow accumulation, the yearly distribution of the streamflow would change. An increase of precipitation stored in the snowpack could reduce the amount of water discharged in winter and increase the streamflow in spring.

Transpiration affects the flow during the vegetation period. A lower transpiration could reduce the uptake in spring and contribute to simulate higher peak flows. However, the overall effect of the simulation of the streamflow would be less significant than the factors previously discussed.

The lack of preferential paths contributes to the underestimation of the conductivity of the soil at high precipitation events causing a water flow more distributed over time and thereby a possible underestimation of the discharge at high flow. The difference between modelled and real soil conductivity could have an effect on the residence time of water, i.e., the time spent by water in the soil, which is relevant to predict the amount of chemicals transported to the stream. By calculating the mean residence time (T_R, days) as the ratio between the simulated soil water storage and the streamflow over the period 2012–2014, we found that the mean residence time is related to the streamflow (Sf, mm·d⁻¹) according to the following exponential function: T_R = 496.72 × Sf^−0.899 (R² = 0.93); the median T_R is 846 days.

5. Conclusions

As compared to the previous hydrology module in ForSAFE, the new hydrology concept allows for the simulation of a saturated and an unsaturated zone in the soil, as well as a water flow that reaches the stream in a way that is consistent with measurements. It simulates with a good approximation the saturated zone and the level of soil moisture at different depths. In addition, it captures part of the flow peaks observed in the stream and it simulates a total annual streamflow comparable to the measured one.

The most relevant differences compared to measurements that will affect the chemical transport to the stream are that the model simulates a higher base flow in winter and lower flow peaks after snowmelt. The discrepancies are mainly caused by the assumptions made to regulate the percolation at the bottom of the soil columns. To simulate a saturated zone in the soil and a seasonal change of streamflow it is necessary to reduce the conductivity of the deepest soil layer to include to the lower permeability of soil below the simulated depth. The differences relative to measurements observed in this study will mainly have an effect on the timing with which water reaches the stream. The residence time of water in the soil affects the biogeochemical processes in the soil and eventually the amount of chemicals that are transported to the river.

The discrete structure of the model also limits its capability to represent water storage and transport which in reality are continuous. However, it makes the hydrology module suitable for coupling to the ForSAFE model. The simulation of lateral flows and a saturated zone in ForSAFE will couple several processes in the forest ecosystem with water flows to the steam. This will greatly improve the simulation of chemical exchange in the soil and chemical transport from the soil to watercourses.

Future work including chemical transport will help clarify the impact of water residence times and help identify possible improvements in the modeling of hydrology in ForSAFE.

Supplementary Material

Table S1: Description of the soil columns in the S-Transect model; Table S2: Calculated soil hydraulic properties for the modelled soil columns of the S-transect; Table S3: Objective functions.

Acknowledgments

The authors wish to thank for the financial support granted by the projects “ForWater” and “Managing Multiple Stressors in the Baltic Sea” funded by FORMAS. We also thank Nino Amvrosiadi for the valuable discussions that contributed to the model development and Martin Erlandsson for sharing the data on soil texture.

Author Contributions

The model was developed by Giuliana Zanchi and Salim Belyazid. Lin Yu and Harald Grip contributed to a further refinement of the model. Kevin Bishop and Stephan Köhler provided the experimental data used to test the model. All authors contributed to data analysis, paper writing and proofreading.

Conflicts of Interest

The authors declare no conflict of interest.

References and Notes

References

1. Chang, M. Forest Hydrology: An Introduction to Water and Forests, 2nd ed.; CRC Press: Boca Raton, FL, USA, 2006.

2. Neary, D.G.; Ice, G.G.; Jackson, C.R. Linkages between forest soils and water quality and quantity. For. Ecol. Manag. 2009, 258, 2269–2281.

3. FAO. Forest and Water-International Momentum and Action; FAO: Rome, Italy, 2013; p. 86.

4. Muller, F.; Chang, K.-C.; Lee, C.-L.; Chapman, S. Effects of temperature, rainfall and conifer felling practices on the surface water chemistry of northern peatlands. Biogeochemistry 2015, 126, 343–362.

5. Hellsten, S.; Stadmark, J.; Pihl Karlsson, G.; Karlsson, P.E.; Akselsson, C. Increased concentrations of nitrate in forest soil water after windthrow in southern sweden. For. Ecol. Manag. 2015, 356, 234–242.

6. Taboada-Castro, M.M.; Rodríguez-Blanco, M.L.; Diéguez, A.; Palleiro, L.; Oropeza-Mota, J.L.; Taboada-Castro, M.T. Effects of changing land use from agriculture to forest on stream water quality in a small basin in nw spain. Commun. Soil Sci. Plant Anal. 2015, 46, 353–361.

7. Schelker, J.; Öhman, K.; Löfgren, S.; Laudon, H. Scaling of increased dissolved organic carbon inputs by forest clear-cutting—what arrives downstream? J. Hydrol. 2014, 508, 299–306.

8. Xiao, W.; Chun, Z.; Qing, J. Impacts of climate change on forest ecosystems in northeast china. Adv. Clim. Chang. Res. 2013, 4, 230–241.

9. Aber, J.; Neilson, R.P.; McNulty, S.; Lenihan, J.M.; Bachelet, D.; Drapek, R.J. Forest processes and global environmental change: Predicting the effects of individual and multiple stressors we review the effects of several rapidly changing environmental drivers on ecosystem function, discuss interactions among them, and summarize predicted changes in productivity, carbon storage, and water balance. Bioscience 2001, 51, 735–751.

10. Campbell, J.L.; Rustad, L.E.; Boyer, E.W.; Christopher, S.F.; Driscoll, C.T.; Fernandez, I.J.; Groffman, P.M.; Houle, D.; Kiekbusch, J.; Magill, A.H.; et al. Consequences of climate change for biogeochemical cycling in forests of northeastern north americathis article is one of a selection of papers from ne forests 2100: A synthesis of climate change impacts on forests of the northeastern us and eastern canada. Can. J. For. Res. 2009, 39, 264–284.

11. Belyazid, S.; Westling, O.; Sverdrup, H. Modelling changes in forest soil chemistry at 16 Swedish coniferous forest sites following deposition reduction. Environ. Pollut. 2006, 144, 596–609.

12. Gaudio, N.; Belyazid, S.; Gendre, X.; Mansat, A.; Nicolas, M.; Rizzetto, S.; Sverdrup, H.; Probst, A. Combined effect of atmospheric nitrogen deposition and climate change on temperate forest soil biogeochemistry: A modeling approach. Ecol. Model. 2015, 306, 24–34.

13. Homann, P.S.; McKane, R.B.; Sollins, P. Belowground processes in forest-ecosystem biogeochemical simulation models. For. Ecol. Manag. 2000, 138, 3–18.

14. Maxwell, R.M.; Putti, M.; Meyerhoff, S.; Delfs, J.-O.; Ferguson, I.M.; Ivanov, V.; Kim, J.; Kolditz, O.; Kollet, S.J.; Kumar, M.; et al. Surface-subsurface model intercomparison: A first set of benchmark results to diagnose integrated hydrology and feedbacks. Water Resour. Res. 2014.

15. Wallman, P.; Svensson, M.G.E.; Sverdrup, H.; Belyazid, S. Forsafe—An integrated process-oriented forest model for long-term sustainability assessments. For. Ecol. Manag. 2005, 207, 19–36.

16. Yu, L.; Belyazid, S.; Akselsson, C.; van der Heijden, G.; Zanchi, G. Storm disturbances in a swedish forest—A case study comparing monitoring and modelling. Ecol. Model. 2016, 320, 102–113.

17. Zanchi, G.; Belyazid, S.; Akselsson, C.; Yu, L. Modelling the effects of management intensification on multiple forest services: A swedish case study. Ecol. Model. 2014, 284, 48–59.

18. Richmond, B. An Introduction to System Thinking with Stella; Isee Systems Inc: Lebanon, NH, USA, 2004.

19. Lindström, G.; Gardelin, M. Chapter 3.1: Hydrological Modelling-Model Structure. In Modelling Groundwater Response to Acidification; Swedish Meteorological and Hydrological Institute (SMHI): Norrköping, Sweden, 1992; pp. 33–36.

20. Balland, V.; Pollacco, J.A.P.; Arp, P.A. Modeling soil hydraulic properties for a wide range of soil conditions. Ecol. Model. 2008, 219, 300–316.

21. Saxton, K.E.; Rawls, W.J. Soil water characteristic estimates by texture and organic matter for hydrologic solutions. Soil Sci. Soc. Am. J. 2006, 70, 1569–1578.

22. Reid, I. The choice of assumed particle density in soil tests. Area 1973, 5, 10–12.

23. Bergström, S. The HBV Model: Its Structure and Applications; Swedish Meteorological and Hydrological Institute (SMHI): Nörrköping, Sweden, 1992; p. 35.

24. Dingman, S.L. Physical hydrology, 2nd ed.; Waveland Press, Inc.: Long Grove, IL, USA, 2008; p. 646.

25. Rosengren, U.; Stjernquist, I. Gå på djupet!: Om rotdjup och rotproduktion i olika skogstyper; SUFOR: Alnarp, Sweden, 2004; p. 55.

26. Laudon, H.; Taberman, I.; Ågren, A.; Futter, M.; Ottosson-Löfvenius, M.; Bishop, K. The krycklan catchment study—a flagship infrastructure for hydrology, biogeochemistry, and climate research in the boreal landscape. Water Resour. Res. 2013, 49, 7154–7158.

27. Nyberg, L.; Stähli, M.; Mellander, P.-E.; Bishop, K.H. Soil frost effects on soil water and runoff dynamics along a boreal forest transect: 1. Field investigations. Hydrol. Process. 2001, 15, 909–926.

28. Lindqvist, G.; Nilsson, L.; Gonzalez, G. Depth of Till Overburden and Bedrock Fractures on the Svartberget Catchment as Determined by Different Geophisical Methods; University of Luleå: Luleå, Sweden, 1989; p. 16.

29. Erlandsson, M.; Uppsala University, Uppsala, Sweden. Personal communication, 2015.

30. Bishop, K.H. Episodic Increases in Stream Acidity, Catchment Flow Pathways and Hydrograph Separation; University of Cambridge: Cambridge, UK, 1991.

31. Yi, S.; Manies, K.; Harden, J.; McGuire, A.D. Characteristics of organic soil in black spruce forests: Implications for the application of land surface and ecosystem models in cold regions. Geophys. Res. Lett. 2009.

32. Mitsch, W.J.; Gosselink, J.G. Wetlands, 5th ed.; Wiley: Hoboken, NJ, USA, 2015.

33. Kellner, E. Wetland—Different Types, Their Properties and Functions; Department of Earth Sciences/Hydrology: Uppsala, Sweden, 2003; p. 62.

34. Bishop, K.; Seibert, J.; Nyberg, L.; Rodhe, A. Water storage in a till catchment. Ii: Implications of transmissivity feedback for flow paths and turnover times. Hydrol. Process. 2011, 25, 3950–3959.

35. Rodhe, A. On the generation of stream runoff in till soils. Nord. Hydrol. 1988, 20, 1–8.

36. Cuo, L.; Giambelluca, T.W.; Ziegler, A.D. Lumped parameter sensitivity analysis of a distributed hydrological model within tropical and temperate catchments. Hydrol. Process. 2011, 25, 2405–2421.

37. Amvrosiadi, N.; Bishop, K.; Grabs, T.; Beven, K.; Seibert, J. Water storage dynamics in a till hillslope: The foundation for modeling flows and residence times. Hydrol. Process. 2015, submitted.

38. Grip, H.; Rodhe, A. Vattnets väg från regn till back; Hallgren & Fallgren: Uppsala, Sweden, 2000.

39. Bengtsson, L. Snowmelt estimated from energy budget studies. Nord. Hydrol. 1976, 7, 3–18.

40. Grip, H.; Halldin, S.; Lindroth, A. Water use by intensively cultivated willow using estimated stomatal parameter values. Hydrol. Process. 1989, 3, 51–63.