Simulating the Effects of Thinning Events on Forest Growth and Water Services Asks for Daily Analysis of Underlying Processes

Abstract

Forest growth function and water cycle are affected by climatic conditions, making climate-sensitive models, e.g., process-based, crucial to the simulation of dynamics of forest and water interactions. A rewarded and widely applied model for forest growth analysis and management, 3PG, is a physiological process-based forest stand model that predicts growth. However, the model runs on a monthly basis and uses a simple soil-water module. Therefore, we downscale the temporal resolution to operate daily, improve the growth modifiers and add a responsive hydrological sub-model to represents the key features of a snow routine, a detailed soil-water model and a separated soil-evaporation calculation. Thereby, we aim to more precisely analyze the effects of thinning events on forest productivity and water services. The novel calibrated 3PG-Hydro model was validated in Norway spruce sites in Southern Germany and confirmed improvements in building forest processes (evapotranspiration) and predicting forest growth (biomass, diameter, volume), as well as water processes and services (water recharge). The model is more sensitive to forest management measures and variability in soil water by (1) individualization of each site’s soil, (2) simulation of percolation and runoff processes, (3) separation of transpiration and evapotranspiration to predict good evapotranspiration even if high thinning is applied, (4) calculation in daily time steps to better simulate variation and especially drought and (5) an improved soil-water modifier. The new 3PG-Hydro model can, in general, better simulate forest growth (stand volume, average diameter), as well as details of soil and water processes after thinning events. The novel developments add complexity to the model, but the additions are crucial and relevant, and the model remains an easy-to-handle forest simulation tool.

1. Introduction

Rising afforestation necessitates the adaptive management of the provision of ecosystem services from natural forest resources in response to the observed and foreseen changes in climatic conditions [1,2,3]. Forests are facing and will continue to deal with an increase in droughts and waterlogging—two of the most limiting factors to growth and production system [2]. Therefore, as manifested in the Global Forest Goals Agenda of 2030 [4], well-balanced water management plays a crucial role not only for forest growth but for sustaining water quality and water yield of the whole catchment area [2,5,6]. This makes the development and application of reliable models of forest and water interactions indispensable to the achievement of intelligent planning and decision making [5,7].

Forest growth (carbon) and water cycles are affected by climatic conditions, making climate-sensitive models, e.g., process-based, crucial to the simulation of dynamics of forest and water interactions. A rewarded and widely applied model for forest growth analysis and management is 3PG [8,9]. 3PG is a physiological process-based forest stand model that predicts growth and operates in annual or monthly time-steps [10,11]. The underlying hydrological model is a simple water balance consisting of one soil-water content (SW) ranging between a minimum and maximum, which influences plant growth by a soil-moisture-dependent modifier. Input consists of precipitation, infiltrating until the maximum SW is reached, with excess water either pooling or running off. Output of the soil is evapotranspiration by the stand. In an updated version of 3PG (3PG vsn 2.7: [11,12]), the hydrological processes were improved by including a snow routine and separating evapotranspiration in the canopy and the understory [13].

Even though relevant forest growth dynamics display in years or centuries, important hydrological processes are taking place rapidly on a smaller timescale, from hours to weeks [14,15]. This makes a temporal downscaling of 3PG necessary for a good representation of forest and water interactions. Especially when it comes to intense precipitation and runoff formation, short time steps are needed for a more accurate modelling [16].

Simulating forest growth responses to SW availability and displaying forest-water interactions (e.g., water yield) demands an accurate prediction of SW [5,17]. Therefore, the simple soil model of 3PG is not sufficient because, firstly, plants primarily depend on water of an effective root zone, whereas water processes affect soil properties [18,19]. Secondly, these processes are complex and highly texture-dependent [6,17,20,21].

Evaporation and transpiration are calculated as one mass flow in 3PG, which can cause simulated evapotranspiration and thus SW to significantly deviate from real conditions, especially when the stand is thinned or a clear-cut [22,23]. Therefore, in order to achieve a better simulation of SW prediction that is more sensitive to stand management applications, a separated calculation of soil evaporation is necessary [24,25]. The van-Genuchten-Mualem model (VGM) is the most used soil parametrization for water flow modelling in porous media and is flexible in application. Hence, it is a good approach for modelling soil-water processes [26,27,28].

3PG has been processed intensively, e.g., for spatial analysis with GIS [29], a large-scale application for analyzing ecosystem services [30], or applied to a catchment [31], as well as by enhancing the water balance with a hydrological sub-model [32]. However, a sufficient representation of forest soil layers and water processes is still in development [13].

Developing a new and more detailed hydrological sub-model requires a different hydrological approach. This approach should keep the model and its applicability simple by adding the least possible number of new parameters while improving the mentioned deficiencies. Nevertheless, this updated version requires the addition of new temporal resolution to the model so that it can use daily input data and a set of new parameters for the soil. Without daily climate data, precipitation and hydrological processes would be displayed inaccurately in their variability, and the same would be true for SW and forest growth [14,15,16]. The new soil parameters are essential to the application of the model to different soil types and the simulation of the variability of soil hydrological characteristics because soils differ significantly in terms of texture [21,33]. To simplify this, we integrated VGM parameter presets for each soil class of the original 3PG and left VGM parametrization optional.

The ultimate goal of the model development was to support the discovery of adaptive thinning regimes to cope with the climate change—especially drought events—and to incorporate forest water services in the analysis, in addition to timber production and carbon storage. Therefore, the new version of the 3PG-Hydro model was developed to investigate (1) the interaction between water cycles and forest growth and (2) how forest management may deal with the effects of climate change on these interacting processes. To check the performance of our new version, we tested 3PG-Hydro on three stands of Norway spruce (Picea abies (L.) H. Karst.) and validated it for forest-growth and soil-water-content prediction. Hereby, we focused on comparing the results of the 3PG-Hydro against those of 3PG vsn. 2.7 on both daily and monthly resolution.

2. Materials and Methods

2.1. Upgrade Overview

We implemented 3PG-Hydro in 3PG vsn. 2.7 using the model described by Sands (2004 & 2010) and Meyer et al. (2018) [11,12,13]. The changes consist of downscaling the temporal resolution to operate daily, improving the growth modifiers and adding a hydrological sub-model. The hydrological sub-model represents the core upgrade of 3PG-Hydro with the key features of a snow routine, a detailed soil-water model and a separated soil-evaporation calculation (Figure 1).

2.2. Temporal Resolution

The downscaling from monthly to daily processing was achieved by using daily input data (temperature, precipitation, solar radiation and optional vapor pressure deficit (VPD)) and transforming the time-dependent equations, which affected every main routine of 3PG. The equations, from growth modifiers (temperature, VPD, soil water, frost, age and physiological modifier) to calculation of leaf area index (LAI), canopy conductance, net primary production (NPP), gross primary production (GPP), evapotranspiration and water balance till stand growth and dying (stem, foliage, roots), were downscaled temporally to daily calculations. Besides the original options for monthly and annual output, outputs are now available in daily resolution.

2.3. Upgraded Modifiers

Apart from downscaling the growth modifiers, the frost modifier, as well as the soil-water modifier, was upgraded. The improved frost modifier labels a day as a frost day if daily mean temperature drops below 0 °C, rather than having numbers of frost days per month as an input. This makes the application easier by demanding one less input variable. The original equation of the soil-water modifier, f_θ [9], was upgraded by making it only dependent on the usable field capacity of the effective root zone (ER) [18,34], making f_θ more sensitive to changes in the ER moisture:

$$f_{\theta }=\frac{1}{{\left(1+\frac{1-Minimum\left(\frac{{\theta }_{ER}-{\theta }_{wp}}{{\theta }_{fc}-{\theta }_{wp}},1\right)}{c_{sw}}\right)}^{p_{sw}}}$$

(1)

where θ_ER is the volumetric water content of the ER (m³ × m⁻³), θ_wp is the wilting point (m³ × m⁻³) and θ_fc is the volumetric wilting point of the ER (m³ × m⁻³), with all three being new variables of 3PG-Hydro (explained in Section 2.4.2). c_sw is the soil-water constant (−), and p_sw is the soil water power (−), both original parameters of 3PG.

2.4. Hydrological Sub-Model

2.4.1. Snow Routine

An accurate snow routine had been missing in 3PG vsn. 2.7 and was therefore implemented. To keep it simple, the routine consists of a snow balance, changing precipitation to snow if mean temperature drops below 0 °C and melting it vice versa. Snow gets separately accumulated as intercepted snow by canopy and understory, as well as throughfall, which is calculated in the similar way as rain [9]. The melting process is calculated by the day-degree equation (Appendix A.1.1) [35,36], which needs a new parameter, the snowmelt factor, c_melt, because of its spatial variability [37,38]. This consists of melting and direct evaporation of accumulated intercepted snow, as well as melting and transferring of accumulated throughfall snow in the water model.

2.4.2. Soil-Water Model

The soil-water model represents the core of 3PG-Hydro by simulating the processes for forest growth, as well as soil-water interactions, on a more detailed scale. To achieve this, the single-box soil model of 3PG vsn. 2.7 was enlarged to two boxes, representing the effective root zone (ER), where the majority of primary roots and therefore water uptake occurs, and the deep root zone (DR), which is mainly affected by soil-water processes but can also be affected by deep-root water uptake (Figure 1). The depth of the ER depends on species, age and soil type but generally varies between 0.5 and 2 m [18,34,39,40,41]. The depth of the DR is defined as the water table. The processes affecting the ER water content, SW_ER (mm), and the DR water content, SW_DR (mm), are presented by two water balances. In all following equations, mm is defined as liter × m⁻²:

$$S{W}_{ER}\left(t\right)=S{W}_{ER}\left(t-1\right)+I\left(t\right)-ET\left(t\right)-PC\left(t\right)+C{R}_{ER}\left(t\right),t\ge 1$$

(2)

$$S{W}_{ER}S{W}_{DR}\left(t\right)=S{W}_{DR}\left(t-1\right)-Trans{p}_{DR}\left(t\right)+PC\left(t\right)-C{R}_{ER}\left(t\right)-DP\left(t\right)+C{R}_{DR}\left(t\right),t\ge 1$$

(3)

where t is time (day), I is infiltration (mm × day⁻¹), ET is evapotranspiration (mm × day⁻¹), PC is percolation from ER to DR (mm × day⁻¹), CR_ER is capillary rise from DR to ER (mm × day⁻¹), Transp_DR is deep-root transpiration (mm × day⁻¹), DP is deep percolation as outflow of the DR (mm × day⁻¹) and CR_DR is capillary rise from water table to DR (mm × day⁻¹). SW_ER(t = 0) and SW_DR(t = 0) are initial water contents of each layer, respectively. Infiltration is calculated as the daily sum of throughfall, snowmelt, pooled soil water and irrigation infiltrating in the ER to a threshold, either being the saturated hydraulic conductivity, K_sat, of 1 m² area of the ER (mm × day⁻¹) or until SW_ER reaches the maximum soil water content of the ER, SW_ERsat (mm), with exceeding water becoming runoff (Appendix A.1.2). K_sat is a new optional parameter with a preset for each 3PG soil type: sand, sandy loam, clay loam and clay (

Table A1. VGM preset soil parameters. Values from Jabro (1992), Gootman et al. (2020), Hodnett & Tomasella, (2002) and Nemes et al. (2001) were considered and further modified.

Soil Type	θsat (m3 × m−3)	θres (m3 × m−3)	n (−)	α (m−1)	Ksat (m × day−1)
sand	0.38	0.02	1.55	4	3.50
sandy loam	0.4	0.08	1.35	3.5	1
clay loam	0.44	0.1	1.25	2.8	0.4
clay	0.5	0.12	1.1	2.4	0.1

) [42,43]. SW_ERsat is calculated in the initialization of 3PG-Hydro:

$$S{W}_{ERsat}=\frac{V_{ER}*\left(1-f_{skER}\right)*{\theta }_{sat}}{A}*1000$$

(4)

where V_ER is the soil volume of 1 m² area and ER height (m³), A is the area of 1 m², f_skER is the fraction of skeleton in the ER (−), θ_sat is the volumetric saturation water content (m³ × m⁻³) and the 1000 conversion factor from m³ × m⁻² to mm. SW_DRsat is calculated similarly but with V_DR being the soil volume of 1 m² area and DR height (m³) and f_skDR being the fraction of skeleton in the DR (−) (Appendix A.1.3). The fraction of skeleton in soil is important because of its influence on absolute soil water content [2], a making it a mandatory parameter for both ER and DR to be set before simulations. θ_sat, along with θ_r, the volumetric residual water content (m³ × m⁻³), the Van Genuchten n (−) and Van Genuchten α (m⁻¹) are the initialization parameters for the soil characteristics based on van-Genuchten-Mualem equations (VGM) [27,28]. All parameters are part of the preset for each 3PG soil type (Table A1) [44,45]. However, they can optionally be changed. With these parameters, the volumetric field capacity, θ_fc, and volumetric wilting point, θ_wp, respectively, (m³ × m⁻³) are calculated by VGM in the initialization of 3PG-Hydro (Appendix A.1.4) and are essential for further simulation of soil and water processes.

Soil properties change with depth due to compaction and pore-size shrinkage [46,47,48]. Therefore, soil properties should display a difference between ER and DR. To simulate this change by depth but also keep the parametrization simple, we implemented it in an internal calculation, where θ_fc is decreased by 20% [45,49,50] and K_sat is decreased by a factor, depending on soil texture (Appendix A.1.5) [45,51,52].

The percolation, PC (mm × day⁻¹), and capillary rise, CR_ER (mm × day⁻¹), per day and area of 1 m² are modeled as a single vertical stationary flow based on the Darcy-Buckingham-Law and Richards equation:

$$PC=\{\begin{array}{cc}\hfill -K_u\left({\Theta }_{ER}\right)*\left(\frac{\partial h}{z_{ER}}+1\right)*1000,& {\theta }_{ER}<{\theta }_{fc} and \frac{\partial h}{z_{ER}}>-1 \hfill \\ \hfill 0,& otherwise\hfill \end{array}$$

(5)

$$C{R}_{ER}=\{\begin{array}{cc}\hfill -K_u\left({\Theta }_{ER}\right)*\left(\frac{\partial h}{z_{ER}}+1\right)*1000,& {\theta }_{ER}<{\theta }_{fc} and {\theta }_{DR}>{\theta }_{fcDR} and \frac{\partial h}{z_{ER}}<-1\hfill \\ \hfill 0,& otherwise \hfill \end{array}$$

(6)

where K_u(Θ_ER) is the unsaturated hydraulic conductivity (m × day⁻¹) as a function of the ER water content, Θ_ER (m³ × m⁻³), ∂h is the pressure-head gradient of the pressure head of the ER, h_ER, and DR is h_DR, all respectively (m). z_ER is the height of the ER (m), and the term +1 represents the gravitational potential and 1000 the conversion factor from m³ × m⁻² to mm. Percolation is subjected to the condition of the ER volumetric water content, θ_ER being higher than volumetric field capacity, θ_fc [19,20]. The matric potential below the gravitational potential is expressed by the term ∂h × z_ER⁻¹ > −1. The conditions for capillary rise are θ_ER being below θ_fc, the DR volumetric water content, θ_DR, above DR volumetric field capacity, θ_fcDR, and the matric potential being higher than the gravitational potential, expressed by the term ∂h × z_ER⁻¹ < −1. The deep percolation out of the DR DP_DR and the capillary rise from the water table to the DR CR_DR are calculated in a similar way and under similar conditions but with K_uDR(Θ_DR) being the unsaturated hydraulic conductivity of the DR (m × day⁻¹) as a function of the DR water content Θ_DR (m³ × m⁻³), ∂h being the pressure head of the DR h_DR and the water table (at saturation) h₀ = 0, all respectively (m), and z_DR being the height of the DR (m) (Appendix A.1.6). The needed variables for the conditions and equations of PC, CR_ER, DP and CR_DR are calculated step-by-step as θ_ER and θ_DR, Θ_ER and Θ_DR (Appendix A.1.7), as well as h_ER and h_DR, by transposing the VGM water-content equation to the pressure head as a function of water content, h_ER(Θ_ER) (m) (h_DR(Θ_DR), Appendix A.1.8):

$$h_{ER}\left({\Theta }_{ER}\right)=\frac{1}{\alpha }*{\left({\Theta }_{ER}^{\frac{-n}{\left(n-1\right)}}-1\right)}^{\frac{1}{n}}$$

(7)

where α is the Van-Genuchten α (m⁻¹), n is the Van-Genuchten n (−) and Θ_ER is the water content of the ER (m³ × m⁻³). The values of h_ER(Θ_ER) and h_DR(Θ_DR) are positive for simplification [27] rather than being negative. This was considered in further application. At last, the values of K_u(Θ_ER) and K_uDR(Θ_DR), respectively (m × day⁻¹), are calculated by VGM for each time step, decreasing with lower water content and depending on soil texture (Appendix A.1.10).

As mentioned earlier, in the DR, there is the possibility of water uptake by deep roots if SW_ER cannot cover the transpiration demand of the canopy. This shift/supply by the DR is simulated by a fraction of the required transpiration water, Transp_DR (mm × day⁻¹):

$$Trans{p}_{DR}=Transp*(f_{dr}*Maximum\left(\frac{W{R}_t}{W{R}_{amt}},1\right))$$

(8)

where Transp is transpiration of the canopy (mm × day⁻¹), f_dr is maximum fraction deep root transpiration (−), WR_t is weight of root dry biomass per tree (kg × tree⁻¹) and WR_amt is the average mature tree-root dry biomass (kg × tree⁻¹). f_dr can cover up to 20% of total transpiration, depending on species mass of roots in deeper soil [53]. The ratio between WR_t and WR_amt is used to apply an age dependence to the deep-root uptake, reducing it for younger trees. f_dt and WR_amt are new mandatory species parameters, but both can be derived from the literature [53,54].

2.4.3. Soil Evaporation

To enhance the prediction of forest management effects on evapotranspiration and use the opportunity given by a more accurate water content of the upper soil, evapotranspiration was partitioned into transpiration and soil evaporation. Transpiration is calculated by the Penman-Monteith-equation (P-M), as in the original 3PG [9], but using canopy-intercepted radiation instead of total radiation as input. The soil evaporation, E_soil (mm × day⁻¹), is calculated with a simple adaptation of the 3PG original P-M:

$$\begin{gathered} E_{soil}=h*\frac{g_S*\left(s*{\phi }_{na}+g_{aS}*{\rho }_a*c_{pa}*D\right)}{\lambda *(\left(\gamma +s\right)*g_S+\gamma *g_{aS}} \\ with{\phi }_{na}=a+b*RADsoil \end{gathered}$$

(9)

where g_s is the soil conductance (m × s⁻¹) substituted for canopy conductance from the original equation [55], g_as is the soil aerodynamic conductance (m × s⁻¹) substituted for the aerodynamic boundary-layer conductance from the original equation and RADsoil (W × m⁻²) as part of the equation for soil absorbed radiation, φ_na, is the residual radiation, being the solar radiation subtracted by intercepted radiation. The other variables are similar to those in the original 3PG P-M and can be read in Landsberg & Sands [9]. g_s is calculated similarly to canopy conductance in the original 3PG:

$$g_S=max{g}_S*exp\left(\frac{{\theta }_{ER}}{{\theta }_{sat}}\right)$$

(10)

where maxg_s is the maximum value of g_s (m × s⁻¹), being a new mandatory parameter to be set before simulations varying between 10⁻³ and 10⁻⁵ m × s⁻¹, depending on soil type [56] and able to be adjusted if soil evaporation is over- or underestimated. maxg_s decreases exponentially by the saturation of the soil, expressed as the ratio between the volumetric soil water content of the ER (m³ × m⁻³), θ_ER, and the volumetric saturation water content (m³ × m⁻³), θ_sat [57]. g_as is the inverse of aerodynamic resistance and can be estimated for wind speeds in forests to around 0.02 m × s⁻¹ [58,59]; however, it can be changed optionally. Because not only soil conductance but also total soil evaporation, E_soil, are directly linked to soil moisture and get hyperbolically reduced due to drying [58,60,61,62], E_soil was subdivided into two conditions: a wet soil and a drying soil:

$$E_{soil}=\{\begin{array}{cc}\hfill E_{soil},& days without Rain\le 1\hfill \\ \hfill E_{soil}*0.01*{\left(\left(dayswithoutrain-1\right)-10\right)}^2,& days without Rain>1\hfill \end{array}$$

(11)

where the conditions of days without rain ≤ 1 (day) expresses the state of wet soil, being a day with rain prolonging to one day without rain, where E_soil is not affected, and days without rain > 1 (day) expresses the state of a drying soil, starting on the second day without rain and reducing E_soil hyperbolically by a reduction equation, where the soil is dry after 10 days without rain and E_soil becomes 0.

2.5. Model Application

We applied 3PG-Hydro on three stands of Picea abies in southwestern Germany at the locations called by their region: Conventwald, Heidelberg and Ochsenhausen. Conventwald stand was 76 years, Heidelberg was 84 years and Ochsenhausen was 77 years old at the start of each simulation. All sites were planted with a unique provenances of Norway spruce from the Black Forest area. The stands and their site characteristics were monitored over a period of 20, 16 and 14 years. The management of the sites before the monitored periods is unknown. Thinning was applied at all three sites with low intensities (

Table A2. Three thinning regimes applied at the sites. Stems are the stems before thinning took place. The percentage is the proportion of thinned stems to stems before thinning.

	Age	Stems	Thinned Stems	Method
Conventwald	80	528	60 (11.4%)	Below
	85	468	52 (11.1%)	Below
Heidelberg	89	380	20 (5.3%)	Below
	91	360	15 (4.2%)	Below
Ochsenhausen	77	496	4 (0.8%)	Middle
	79	492	8 (1.6%)	Middle
	80	484	4 (0.8%)	Middle

). Model input for each stand was species-specific parameters of Picea abies, site soil characteristics (calibrated for Conventwald, preset of sandy loam for Heidelberg and clay loam for Ochsenhausen (Table A1)), climate data (daily precipitation, daily mean, minimum and maximum temperature, and daily solar radiation) and the thinning regimes. We analyzed the output of forest growth (stand volume and mean diameter) and forest processes (transpiration and evaporation), as well as soil processes (SWC, percolation, deep percolation and runoff). Stand and site data, climate data for Conventwald, as well as soil moisture data at Conventwald (measured over a period of 3-years), were provided by the regional forest research institute (FVA-BW). Climate data for the Heidelberg and Ochsenhausen sites were derived from the open-data platform of the German weather service (DWD). The species-specific parameters for Picea abies were adjusted based on the values found in the literature [34,63,64,65].

2.5.1. Performance Evaluation

For performance evaluation of 3PG-Hydro, we used the residual mean square error RMSE and the Nash-Sutcliff-Efficiency NSE. RMSE is a common statistical method for measuring the difference between simulated and observed data, with a value of 0 being no error (Appendix A.1.10). For the validation of soil water content, we adjusted the RMSE to evaluate how the simulation underestimates the first quartile of total observed data, RMSE_Q1, overestimates the third quartile of total observed data, RMSE_Q3, as well as performs in between the first and third quartile, RMSE_Q:

$$RMS{E}_{Q1}=\sqrt{\frac{1}{n}{\displaystyle \sum }_{i=1}^n\{\begin{array}{cc}\hfill 0,& Y_i^{sim}\ge Y_i^{obsQ1}\hfill \\ \hfill {(Y_i^{obsQ1}-Y_i^{sim})}^2,& Y_i^{sim}<Y_i^{obsQ1}\hfill \end{array}}$$

(12)

$$RMS{E}_{Q3}=\sqrt{\frac{1}{n}{\displaystyle \sum }_{i=1}^n\{\begin{array}{cc}\hfill 0,& Y_i^{obsQ3}\ge Y_i^{sim}\hfill \\ \hfill {(Y_i^{obsQ3}-Y_i^{sim})}^2,& Y_i^{sim}>Y_i^{obsQ3}\hfill \end{array}}$$

(13)

$$RMS{E}_Q=\sqrt{\frac{1}{n}{\displaystyle \sum }_{i=1}^n\{\begin{array}{cc}\hfill 0,& Y_i^{obsQ3}\ge Y_i^{sim}\ge Y_i^{obsQ1}\hfill \\ \hfill {(Y_i^{obsQ1}-Y_i^{sim})}^2,& Y_i^{sim}<Y_i^{obsQ1}\hfill \\ \hfill {(Y_i^{obsQ3}-Y_i^{sim})}^2,& Y_i^{sim}>Y_i^{obsQ3}\hfill \end{array}}$$

(14)

where n is the number of observations, Y_i^sim is the simulated value, Y_i^obsQ1 is the value of the first quartile of observed data and Y_i^obsQ3 is the value of the third quartile of observed data. The NSE is a test for model performance, comparing residuals of simulated and observed data to variance of observed data [66]:

$$NSE=1-\frac{{{\displaystyle \sum }}_{i=1}^n{(Y_i^{obs}-Y_i^{sim})}^2}{{{\displaystyle \sum }}_{i=1}^n{(Y_i^{obs}-{\overline{Y}}^{obs})}^2}$$

(15)

where n is the number of observations, Y_i^obs is the observed data, Y_i^sim is the simulated value and Ῡ^obs is the mean of total observed data. Even though the equation of NSE is similar to that of the coefficient of determination, R², the difference is that NSE evaluates the performance of a model simulation instead of a statistical model, and therefore, NSE can be below 0, signifying the simulation is worse than the mean of total observed data. Values above 0 mean better performance than the mean, with 1 being perfect fit. For the validation of soil water content, we also adjusted NSE to quantify how well the simulation predicts the observed data by not underestimating low values (first quartile), NSE_Q1, not overestimating high values (third quartile), NSE_Q3, and performing in between NSE_Q:

$$NS{E}_{Q1}=1-\frac{{{\displaystyle \sum }}_{i=1}^n\{\begin{array}{cc}\hfill 0,& Y_i^{sim}>Y_i^{obsQ1}\hfill \\ \hfill {(Y_i^{obsQ1}-Y_i^{sim})}^2,& else\hfill \end{array}}{{{\displaystyle \sum }}_{i=1}^n{\left(Y_i^{obsQ1}-{\overline{Y}}^{obsQ1}\right)}^2}$$

(16)

$$NS{E}_{Q3}=1-\frac{{{\displaystyle \sum }}_{i=1}^n\{\begin{array}{cc}\hfill 0,& Y_i^{sim}<Y_i^{obsQ3}\hfill \\ \hfill {\left(Y_i^{sim}-Y_i^{obsQ3}\right)}^2,& else\hfill \end{array}}{{{\displaystyle \sum }}_{i=1}^n{\left(Y_i^{obsQ3}-{\overline{Y}}^{obsQ3}\right)}^2}$$

(17)

$$NS{E}_Q=\frac{NS{E}_{Q1}+NS{E}_{Q3}}{2}$$

(18)

where n is number of observations, Y_i^sim is the simulated value, Y_i^obsQ1 is the value of the first quartile of observed data, Ῡ^obsQ1 is the mean of total the first quartile of observed data, Y_i^obsQ3 is the value of the third quartile of observed data and Ῡ^obsQ3 is the mean of total third quartile of observed data.

2.5.2. Calibration

The soil parameters θ_sat, θ_r, n, α and K_sat were calibrated for prediction error of simulated volumetric soil water content red to observed data over a period of 2 years (number of observations N = 730). Calibration was done by the Monte-Carlo method [67], where 5000 parameter sets were randomly selected of a range resembling that of values for soils from sand to clay loam, evaluating the prediction on RMSE and RMSE_Q to get the minimal divergence of simulated and observed data. Furthermore, a sensitivity analysis was done by the Hornberger-Spear-Young method (HSY) for the 1000 parameter sets with least RMSE and RMSE_Q with a behavioral rate of 15% (150 best outcomes) [68]. Sensitivity was tested with a two-sample Kolmogorov-Smirnov Test (KS-Test), defining sensitivity and slight sensitivity as significance p ≤ 0.01 and p ≤ 0.05, respectively, and critical value D_0.01 = 0.15 and D_0.05 = 0.12, respectively [69].

2.5.3. Validation

Forest growth was simulated with 3PG-Hydro, as well as with 3PG vsn. 2.7, on daily and monthly calculation resolution. We validated each simulation for all stands over the monitored period to see not only how the hydrological sub-model performs but also how the downscaling of temporal resolution affected forest-growth prediction. We simulated soil water content of the ER (depth = 0.6 m) with 3PG-Hydro daily, with the calibrated soil parameters as input. For comparison purposes, we simulated the soil water content of the ER with 3PG vsn. 2.7 on daily resolution. We set the 3PG vsn. 2.7 soil parameters (minimum and maximum available soil water) to the calibrated values of θ_r and θ_sat, respectively. We compared the results to measured data at 0.3 m and 0.6 m depth over a 3-year period at Conventwald (N = 1096). To evaluate the effect of temporal downscaling on soil-water processes, we compared percolation prediction by 3PG-Hydro daily to prediction by 3PG-Hydro monthly. Evapotranspiration (ET) simulation was not validated due to non-availability of forest ET data. However, to get an overview of the ET performance, it was simulated with 3PG-Hydro daily and 3PG vsn. 2.7 daily for all stands over the monitored periods and visually evaluated. Furthermore, we simulated ET of a clear-cut (0 trees × hectar⁻¹) area for each stand and visually compared the results to data of actual ET of grass obtained from nearby areas to check how the separated soil evaporation predicted ET without trees. Furthermore, we did a water balance for all stands over the monitored periods.

3. Results

3.1. Calibration

Depending on evaluation with RMSE or RMSE_Q, the calibration provided two parameter sets (

Table 1. Calibration parameter sets. θ_sat (m³ × m⁻³), θ_r (m³ × m⁻³), n (−), α (m⁻¹), K_sat (m × day⁻¹), RMSE & RMSE_Q (m³ × m⁻³).

	θsat	θr	n	α	Ksat	RMSE	RMSEQ
On RMSE	0.38	0.143	1.737	0.647	2.632	0.026	0.015
On RMSEQ	0.38	0.149	1.817	0.663	1.931	0.027	0.014

):

RMSE and RMSE_Q, respectively (m³ × m⁻³), showed almost the same deviation. Apparently, θ_sat, θ_r, n and α display a minor difference, while K_sat shows a greater difference. For further evaluation, we decided to use the parameters from the calibration on RMSE_Q because our aim was to achieve good performance within water-content range from the first to third quartile. A sensitivity was not visible by examining the outcome of the 5000 parameter sets (Figure A2) but by analysis with HSY (Figure A3). The KS-Test of the parameters obtained with calibration on RMSE led to sensitivity of θ_sat, θ_r, α and K_sat and no sensitivity of n. Calibration on RMSE_Q led to sensitivity of θ_sat and θ_r as well as slight sensitivity of α and K_sat and no sensitivity of n (Figure A3).

3.2. Validation

3.2.1. Forest Growth

The evaluation of stand-volume (SV) simulation (Figure 2) indicates a better performance of 3PG-Hydro compared to 3PG vsn. 2.7, with almost perfect prediction at Heidelberg and Ochsenhausen. Comparing the prediction of daily versus monthly calculations, the deviation from 3PG-Hydro is of minor significance at Conventwald and Heidelberg, but calculations are better at Ochsenhausen for 3PG vsn. 2.7. The better performance of the daily resolution is evident. For mean diameter at breast height (DBH), a different picture emerges, showing very good and similar performance of all models at Conventwald, better performance of 3PG-Hydro and especially 3PG-Hydro on monthly resolution at Heidelberg and at Ochsenhausen. The worst performance was achieved by 3PG-Hydro daily calculations, while the best and almost perfect prediction what achieved by 3PG-Hydro monthly calculations. However, the RMSE of 3PG-Hydro daily at Ochsenhausen is only 2%, compared to 0.8% of 3PG-Hydro monthly. Generally, DBH simulation by 3PG-Hydro is better on monthly resolution, while 3PG vsn. 2.7 performs better on daily resolution.

3.2.2. Thinning

Thinning was applied with very low intensities at Heidelberg, Ochsenhausen and Coventwald (Table A2). Figure 3 displays the residuals (observed data—simulated data) of SV and DBH for the years after a thinning event until the next thinning event. This shows a poorer performance of 3PG-Hydro for low thinning at the Conventwald site (Figure 3). However, the available data for this site are very few (N = 2) and should be interpreted very carefully. For the sites with more data and very low thinning, the performance of 3PG-Hydro was clearly better. 3PG-Hydro underestimates SV and DBH growth after thinning, while 3PG vsn 2.7 overestimates SV and DBH growth after thinning. This is important when it comes to correctly planning thinning to achieve certain yields.

3.2.3. Soil Water Content

Visual reflection of the SW prediction (Figure 4) clearly shows the better predictions of 3PG-Hydro, following the behavior of the observed data and staying within the range of the first (Q1) and third quartile (Q3), with a few exceptions. In contrast, 3PG vsn. 2.7 highly overestimates the observed data and is generally out of range. These observations reflect in the performance evaluation (

Table 2. Performance evaluation of 3PG-Hydro and 3PG vsn. 2.7. Observations, N = 1096. RMSE respectively (m³ × m⁻³).

	RMSE	RMSEQ1	RMSEQ3	RMSEQ	NSE	NSEQ1	NSEQ3	NSEQ
3PG-Hydro	0.026	0.01	0.007	0.012	−0.14	0.83	0.93	0.88
3PG vsn. 2.7	0.09	0.008	0.067	0.068	−11.8	0.88	−5.2	−2.16

), presenting a total RMSE on the mean observed data for 3PG-Hydro of 10% compared to 36% for 3PG vsn 2.7. Within the range of Q1 and Q3, the RMSE_Q of 3PG-Hydro is almost six times lower than that of 3PG vsn 2.7. Validated on the mean observed data (NSE), 3PG-Hydro performs worse than the mean of all observed data, but this indicates no failure because the performance within the range is very good (NSE_Q = 0.88). In comparison, 3PG vsn. 2.7 not only performs poorly on mean data (NSE) but also within the range (NSE_Q), with NSE_Q3 clearly indicating the overestimation of high values.

3.3. Evaluation

3.3.1. Percolation

The advantage of daily resolution is better depiction of percolation after high and low precipitation events, whereas monthly resolution takes the mean and thus misses the dynamic variability. This is especially the case for high daily precipitations (Figure A4). At Conventwald, the importance of daily resolution is visible due to the strong reaction to higher rainfall. At Heidelberg and Ochsenhausen, the difference appears minor but is still evident for very significant precipitation events. (Figure A4). The discrepancies at the sites are due to the different soil characteristics and the overall precipitation dynamics.

3.3.2. Evapotranspiration

The mean monthly evapotranspiration (ET) simulated on daily resolution by 3PG-Hydro and 3PG vsn 2.7 show very similar behavior between the separated and the original one-mass-flow calculation at all three sites, with only small deviations (Figure A4). In contrast, the clear-cut simulations clearly show the better performance and the advantage of a separation. 3PG-Hydro estimates the evaporation close to the value of actual evaporation of grass in terms of periodicity and quantity, whereas 3PG vsn. 2.7 simulates values close to 0 (Figure A5).

3.3.3. Water Balance

The water balance shows all input, output and storage contents (mm) summed up over the whole observed periods of 20, 16 and 14 years, respectively for each site (

Table 3. Water balances over the whole periods at all three stands. All values in (mm), (%) indicating the share of input or output, respectively. (*) mean SW_ER and SW_DR are not part of the balance (storage) but displayed for interpretation purposes.

	Input
	Precipitation		Soil Water
	Rain	Snow	Initial Water ER	Initial Water DR	Deep Capillary Rise
Conventwald Heidelberg Ochsenhausen	18,806 (95.7%)	302 (1.5%)	66 (0.3%)	484 (2.5%)	0 (0%)
	11,565 (93.8%)	241 (2%)	100 (0.8%)	420 (3.4%)	0 (0%)
	13,062 (90.6%)	776 (5.4%)	93.5 (0.6%)	490 (3.4%)	0 (0%)
	Output
	Evapotranspiration			Soil processes
	Transpiration	Intercept. evap.	Soil evap.	Runoff	Deep percolation
Conventwald Heidelberg Ochsenhausen	11,446 (59.8%)	3089 (16.1%)	2622 (13.7%)	11 (0.1%)	1970 (10.3%)
	6529 (55.7%)	1603 (13.7%)	2442 (20.8%)	0 (0%)	1140 (9.7%)
	9624 (70.3%)	2260 (16.5%)	1613 (11.8%)	0 (0%)	188 (1.4%)
	Storage
	Final SWER	Final SWDR	Snow storage	Mean SWER *	Mean SWDR *
Conventwald Heidelberg Ochsenhausen	66	454	0	61	458
	106	489	17	73	429
	116	621	1.5	92	607

). Input is dominated by rain, with low contribution of snow, depending on, site and not affected by deep capillary rise. Output is clearly dominated by evapotranspiration, with shares of around 90% at Conventwald and Heidelberg, as well as 98.6% at Ochsenhausen, with remaining portion being deep percolation, while runoff plays no role. The end-of-simulation storage depicts water being stored in snow beyond simulation time, as well as a soil water increase of around 25% in the ER (SW_ER) and DR (SW_DR) at Ochsenhausen, whereas changes at the other stands were minor. Comparison to the mean SW_ER and SW_DR of the whole period to check if these changes were only of temporary showed that, indeed, the change of SW_ER was temporary but that SW_DR increase is robust. The comparison also pointed out a SW_ER decrease of 27% at Heidelberg.

4. Discussion

4.1. The Model

Environmental models combining accurate forest-growth prediction and good SW depiction are scarce but of high importance when it comes to simulation of the effects of future climate change scenarios on forest growth. Our results reflect that the previous version of 3PG did not fulfill the combination of both accurate forest-growth prediction and good SW depiction. However, the updated 3PG-Hydro model achieved a reliable performance, coupling forest growth and soil-water-processes by temporal downscaling, as well as by adding a hydrological sub-model with upgraded modifiers. Daily resolution predicted forest growth better, except for smaller deviations on DBH prediction due to the small magnitude and therefore higher sensitivity of the values. 3PG-Hydro simulated very low thinning better than 3PG vsn. 2.7 because of its daily resolution and more accurate soil water prediction, therefore achieving a higher sensitivity. However, the results are limited due to the low density of data used to analyze the effects of thinning events on forest growth. Higher thinning regimes could not be tested because of the lack of data for validation. However, due to the good separation of evapotranspiration (Figure A5), we assume a better performance of 3PG-Hydro. Further, the particular advantage of a daily simulation is the depiction of variability in precipitation, especially extreme rainfall, and consequently, SW, both being necessary for climate change simulations, where droughts and intense rainfall events are expected to increase [70,71] (Trenberth, 2011; Teuling, 2018). This advantage is reflected in the reaction of percolation to intense rainfall (Figure A4). To further illustrate this point, we created a very extreme scenario. We applied 3PG-Hydro daily to simulate stand volume with precipitation data, where every month had 29 days without rain, followed by one day with 150 mm of rain. We applied 3PG-Hydro monthly with the same precipitation data (interpreted as 150 mm × month⁻¹, daily mean of 5 mm × day⁻¹). We found that the monthly calculation diverges extremely from the daily calculations as a consequence of misrepresenting precipitation details and thus the soil-water processes (Figure 5).

Furthermore, heavy rainfalls lead to runoff formation, which is illustrated in Figure 6, where runoff was predicted by 3PG-Hydro (daily and monthly resolution) with each of the four soil-type presets (Table A1). We used a precipitation scenario of 29 days without rain and 2 days of extreme rainfall (150 mm, 80 mm) as daily input, whereas monthly resolution treated the data as 230 mm × month⁻¹ (daily mean of 7.4 mm × day⁻¹). No runoff was predicted with 3PG-Hydro on monthly resolution due to precipitation being very low. The simulation with daily resolution led to runoff with each soil type in response to the event of 150 mm × day⁻¹. Further, the dependence of runoff magnitude on soil type is evident, indicating the importance of soil texture implementation in the model, which is also reflected in the sensitivity analysis, with four of the five soil parameters being sensitive or slightly sensitive (Section 3.1).

In general, the hydrological sub-model predicted forest growth better than 3PG vsn. 2.7 (Figure 2) and can be verified with a very satisfying outcome, which is also due to its qualified SW prediction at Conventwald (Figure 4). The performance in between the range of observed data (NSE_Q = 0.88) is good, especially because, after all, 3PG is a forest-growth model and a complex coupling of many processes, therefore posing difficulty in combining it with other sub-models. Furthermore, the deviations of SW prediction are marginal because only low overestimations of SW were observed, which is extremely important when it comes to modelling water resources in times of drought. On the other hand, 3PG vsn. 2.7 constantly overestimated SW. Certainly, the maximum soil water parameter of 3PG vsn 2.7 can be set lower, improving the overestimation of SW, but this fitted value would diverge from real conditions (θ_sat). The hydrological processes of runoff, percolation (PC), deep percolation (DP) and capillary rise (CR_ER, CR_DR) can only be vaguely analyzed and interpreted. First, their complexity is very simplified by the sub-model, and second, there is a lack of observed data to compare. Runoff did almost not appear at any site (Table 3), mainly for two reasons: one is the high infiltration (K_sat) and storage capacity of the simulated soils, and the other is that 3PG-Hydro only simulates infiltration based on daily precipitation values, but runoff generation is highly dependent on hourly intensities [16]. This leads to a general underestimation of runoff by 3PG-Hydro. PC, being directly influenced by infiltration and displaying the variability of precipitation and soil texture, appeared in high rates at Conventwald (Table 3). Therefore, the accuracy of SW_ER simulation can be assessed as good. At Heidelberg, PC also played a major role (Table 3) due to the specific soil characteristics of relatively high hydraulic conductivities (sandy loam, Table A1). At Ochsenhausen, it only played a minor role (Table 3) due to the high storage capacity and low hydraulic conductivity of the set soil (clay loam, Table A1). The outflow to the water table is classified as deep percolation (DP) and not groundwater recharge due to it being only one-dimensional (vertical) and therefore not representing important lateral-flow processes, resulting in no accurate recharge. However, DP values can be used to assess forest management (e.g., thinning, rejuvenation, species) influences on water yield. DP showed a similar pattern as PC due to the specific soil characteristics and was of higher importance at Conventwald and Heidelberg but played no role at Ochsenhausen (Table 3). Capillary rise only becomes important in soils under very dry conditions, which was not the case at the test sites (Table 3).

The separation of ET was successfully carried out because of only small differences between the ET simulations of 3PG-Hydro and 3PG vsn. 2.7 (Figure A5). This implies that separation did not lead to over- or underestimation of ET. The great advantage of a separated ET calculation is the improvement in prediction of forest management effects (e.g., clear-cut, Figure A5) on forest and water processes. The generally high shares of ET on output (Table 3) are a peculiarity of 3PG due to ET being calculated on potential VPD, with all intercepted rain gets fully evaporated [9]. For further improvement of 3PG-Hydro, a surface model could be integrated to simulate lateral soil water flow, runoff could be enhanced by adding hourly precipitation intensities for heavy rainfall events and evapotranspiration calculation could be further improved. However, it is important to consider the trade-off between model precision and relevance and carefully add complexity to the relatively easy-to-handle 3PG and 3PG-Hydro models.

4.2. Advantages for Forest Management Application

The motivation for develop ing 3PG-Hydro was to enhance the simulation of forest management measures on both growth and water cycle, especially in the view of a changing climate and drought. 3PG-Hydro can effectively help to investigate how forest management measures (e.g., thinning) (1) can improve individual tree growth by limiting decrease in SW in times of drought and (2) increase the water yield (deep percolation and runoff) for ecological, as well as economical (commercial water), reasons. This is achieved by (1) individualization of each site’s soil characteristics with a preset or by specific parametrization, (2) simulating percolation and runoff processes, (3) separating transpiration and evapotranspiration to predict good evapotranspiration, even if high thinning is applied, (4) calculating in daily time steps to simulate variation and especially drought better and (5) improving the soil-water modifier. This results in a growth simulation more sensitive to forest management measures and variability in soil water. This is a crucial base for evaluating forest water service and studies of its cost-effectiveness [72].

5. Conclusions

3PG-Hydro achieved its set purpose of coupling higher temporal resolved forest-growth prediction with a reliable soil-water model. In its entity, 3PG-Hydro performed very well for all validated and evaluated upgrades compared to 3PG vsn. 2.7. Moreover, the novel development (daily temporal resolution, detailed soil-water processes) aligns with the specific research questions and the crucial need to analyze forest-growth interactions with water and soil processes. 3PG-Hydro can, in general, better simulate forest growth (stand volume, average diameter), as well as details of soil and water processes after thinning events. However, the marginal benefit is not very high regarding intensive data needs and increased model complexity. Overall, 3PG-Hydro offers a good balance of input (low data intensity) and output (relevant decision parameters) crucial in analysis of forest water interactions under changing climate conditions. This is essential to the study of forest drought and analysis of thinning effects on forest growth, water and carbon values.