Reconciling forest profitability and biodiversity conservation under disturbance risk: the role of forest management and salvage logging

We employed data from a field-based research composed by 109 plots in the Southern Black Forest- Germany, covering a forest connectivity (<50%, 50%–75% and >75%) and a forest structure (<5 habitat trees/ha, 5–15 habitat trees/ha and >15 habitat trees/ha) gradient. The inventory recorded tree identity and diameter at breast height (DBH) of all trees, the height of a representative 7% of the trees and the volume of standing and lying deadwood. This dataset was used initialize the forest growth model. The richness of forest birds, orders of flying insects, tree microhabitats and bat activity were assessed in the same plots. Birds seen or heard in the plots were recorded in three surveys during the spring of 2017. Bats were recorded via ultrasonic sound recorders and flying insects were collected with modified window traps in the same period (see Knuff et al 2019 for details). Finally, TreMs were evaluated and recorded based on a detailed catalogue for TreMs classification by Larrieu et al 2018 (for details see supplementary 1 (available online at http://stacks.iop.org/ERL/15/0940a3/mmedia) and Storch et al 2020).

The forest simulations were carried out using the semi-empirical growth model Sibyla (Fabrika 2005), in a 50-year simulation period (2016–2065). Sibyla is a distance-dependent individual tree model that simulates forest stand dynamics at tree, stand and landscape level. The model includes the description of growth, mortality and detailed management options. Forest growth is dependent on multiple climatic variables (atmospheric CO₂- and N₂O- concentration, average temperature, temperature range, precipitation) that inform a potential growth modifier. This is subsequently combined with a tree competition index and a tree vitality index to yield the actual increment.

We simulated four thinning intensities, applying thinning from above (see supplementary 2), namely the business-as-usual (BAU) intensity (calculated as the volume removal in percentage of the standing stock), increase or decrease in thinning intensity and no thinning. These thinning schemes were tailored to the main species in the study region, according to the German National Forest Inventory data (https://bwi.info). Each thinning intensity, referred to as management regimes in the remainder of the manuscript, was driven by three climate change scenarios given by the HadGEM2-ES global climate model, forced by the representative concentration pathways 4.5, 6.0 and 8.5 (https://www.isimip.org).

We retrieved from the forest simulations (for each thinning intensity and climate change scenario) the list of trees in the plots at each time step of the simulation period. These information included tree growth, biomass, the list of trees remaining in the plot, trees removed by management, dead trees due to natural mortality and dead trees due to windstorms. Moreover, the list of trees also included information regarding the species identity, DBH, height, total volume and volume per wood assortment class.

We included deadwood dynamics in the simulations via a post-processing routine, considering deadwood input from natural mortality and windstorms. Both sources of deadwood were assigned to the lying deadwood and snag pools. For deadwood arising from natural mortality, we differentiated between standing and lying deadwood, according to the proportions retrieved from the German National Forest Inventory in the past 10 years (https://bwi.info). For deadwood created by windstorms, a critical slenderness coefficient was applied (height/DBH ratio) as a decision criterion for stem breakage or uprooting. Uprooted trees were allocated to the lying deadwood pool, whereas broken trees were allocated to the snag pool (standing portion of the stem) and lying deadwood pool (broken portion of the stem). Finally, we computed at each period the actual lying and standing deadwood using the decay models developed by Meyer et al (2009) (details in supplementary 2).

We used the wind damage model proposed by Albrecht et al (2012) to estimate the probability of total stand damage (equation (1)), using the species and dominant height as predictors. We set the other predictors in the original model to its mean value and adding the respective constant to the intercept of the model. Subsequently, we randomly assigned total damage events (equal to 75% of the standing stock) across the simulation plots based on the calculated probability. We calculated the mortality of partially damage stands using equation (2), for each plot and species. Finally, the wind damage probability was included in the forest growth simulations using a Monte Carlo approach endogenous to the growth model Sibyla (considering ten repetitions), for each climate scenario and thinning intensity.

$$\begin{equation}TotalDamage = logit^{-1} ({b_0} + {b_1}Species + {b_2}hdom)\end{equation} \tag{ 1 }$$

$$\begin{equation}Damage = logit^{-1} ({b_0} + {b_1}Species + {b_2}hdom + {b_3}TQ)\end{equation} \tag{ 2 }$$

Where: $TotalDamage$: probability of total damage; $Damage$: amount of damage (in %) for partially damaged stands; ${b_i}$: predictor's coefficients retrieved from Albrecht et al (2012); $Species$: tree species group; $hdom$: dominant height (mean height of 100 largest trees in a stand); $TQ$: thinning quotient of past 10 years.

We fitted empirical models for each biodiversity indicator against stand predictors affected by forest management and windstorm events, as well as environmental proxies (ranges are given in supplementary 1). We used the data collected in the field surveys to regress species richness of birds, bats and insect orders, as well as the richness of TreMs, against the stand height, basal area, density, deadwood volume, number of snags, share of conifers in the stand (in terms of the standing volume), altitude and mean temperature during the vegetation season. We applied GLMMs, with the plots as random effects, due to potential spatial autocorrelation among them. All indicators of stand's structural complexity were assumed to potentially influence all biodiversity indicators, either through the direct availability of resource and roosting opportunity or by influencing related forest taxa.

To calibrate the models, we applied Bayesian inference, employing the R package brms (Bürkner 2017). We used a Poisson likelihood to estimate the richness of all indicators (equation (3)) and uniform priors in the interval [−10, 10] for all predictors, which were normalized to avoid scaling issues. We ran three chains for each model with 20 000 iterations and discarding the initial 10 000. Finally, environmental parameters were included in the analysis via the average temperature in the vegetation season (May to September) and the Altitude. The former is important to manifest the effects of climate change on species richness and the latter to express the variation in forest types and management systems applied in the study region. To predict the biodiversity indicators, we tracked the stand structural parameters (used as predictors) at each time step from the forest simulations described in section 2.2. Moreover, we assumed no change in species composition, i.e. no immigration or emigration beyond the study site, as well as no change in the tree species composition in forest management (which would affect TreMs).

$$\begin{equation} Bird\_Richness,\\Bat\_Richness, = \exp ({b_0} + {b_1}Height + {b_2}BasalArea\\ Insect\_Richness, \,\,+ {b_3}Deadwood + {b_4}Density + {b_5}SnagCount\\ TreM\_Richness \,\,\,+ {b_6}ConiferShare + {b_7}Altitude + {b_8}AvgTemp) \end{equation} \tag{ 3 }$$

Where: ${b_i}$: model's coefficients; $TreM\_Richness$: tree microhabitat richness (N/ha); $Bird\_Richness$: bird richness(N/ha); $Bat\_Richness$: approximated bat richness(N/ha); $Insect\_Richness$: insect order richness(N/ha); $Height$: mean plot height (m); $BasalArea$: plot basal area (m²/ha); $Deadwood$: lying deadwood volume (m³/ha); $Density$: plot density (N/ha); $SnagCount$: number of snags (N/ha); $ConiferShare$: volume share of conifers (%); $Altitude$: plot altitude (m a.s.l.); $AvgTemp$: mean temperature in the vegetation season in 2017 and 2018 (°C).

Harmonizing wood production with the provision of other ecosystem services remains an important challenge in the management of forest resources (Piponiot et al 2019). To address this issue and account for potential different preferences related to forest profitability and biodiversity conservation goals of multiple agents or stakeholders, we developed alternative preference scenarios (table 1). Here we adopted a normative approach to the forest planning problem and aimed to derive alternative management responses that optimally balance economic and conservation objectives. These solutions may serve as a benchmark for decision-makers, which can adapt the framework proposed to local contexts and individual preferences.

To find these optimal responses in each preference scenario, we employed robust optimization models that simultaneously take into account economic and ecological indicators of forest management under windstorm risk, climate change and economic parameters uncertainty (interest rate and wood price). To this end, we applied a similar approach to Augustynczik et al (2018) and used a Bernstein approximation scheme to construct Bernstein uncertainty sets (Ben-Tal et al 2009). This is a more general form of ellipsoidal uncertainty sets, with the advantage of offering a same level of protection with reduced conservativeness. Thereby, we developed robust counterparts of the planning problems, finding solutions that guaranteed the performance of the objectives in worst-case scenarios, in the face of parameter fluctuations (details in supplementary 3).

For the economic analysis, we computed the NPV of forest management in the 50-years simulation period (details in supplementary 4). We considered the standing stock at the beginning of the simulation period as the initial investment and the returns were given by discounted thinning revenues and the value of the standing stock at the end of the simulation period. Based on the economic output of forest management and the corresponding impact on the biodiversity indicators considered here, we defined the optimal management portfolio considering four preference scenarios (table 1).

Uncertainty related to climate change (3 RCPs), windstorm risk (30 realizations), economic parameters (Monte Carlo draw of interest rate and wood prices) and on the estimated species richness (via a posterior sampling) were included in the optimization models for policies B to D. For the uncertainty in species and TreM richness, we sampled parameter values for all predictors (equation (3)) from their posterior distribution considering one standard deviation from the mean values. For the interest rate uncertainty, we sampled uniformly from the interval [0.005, 0.01] and for wood price we applied a variation of up to plus and minus 20% of the current price for each assortment class and species. Furthermore, we created Compromise Programming models for scenarios C and D, to balance the various biodiversity indicators (details in supplementary 3 and supplementary 5).

The richness of forest birds was most responsive to environmental parameters, namely the average temperature during vegetation season and altitude. Conversely, richness of bats and TreMs were more responsive to stand structure and thus to forest management interventions (figure 1). Finally, moderate effects of both environment and stand structure were observed for insect order richness, with altitude, basal area and stand height having the strongest effects on this indicator.

Zoom In Zoom Out Reset image size

The increase in average temperature during the vegetation season positively affected the richness of bats and insect orders, while an increase in this parameter caused a reduction in the richness of bird species and TreMs. A reduced richness of birds and insects was recorded in higher altitudes, whereas bats showed the opposite trend. The relationships between plot altitude and species richness might be mediated by the different species composition and climatic parameters, e.g. an increased share of conifers and lower stand height.

Tree density and plot height were important parameters for the evaluation of bat richness, which improved in sparse stands. Conifer-dominated stands displayed lower bat richness, whereas lying deadwood volume positively affected this taxon. TreMs responded positively to stands with higher basal area and lower height, accompanied by a negative impact of stand density, i.e. TreM di versity was increased in stands with larger trees. Conversely, bird richness was weakly affected by management parameters.

We highlight that the direction of the biodiversity indicators' responses to the structural parameters of forest stands showed a considerable degree of uncertainty. The posterior distributions for several structural characteristics exhibited positive probability density over positive and negative values. Additionally, the explained variance of the fitted models ranged from 0.19 for birds to 0.48 for TreMs (for details on the models' diagnostics see supplementary 1).

All biodiversity indicators analyzed in this study showed a positive response to the occurrence of wind disturbances (figure 2). The richness of forest birds increased up to 4% due to the creation of deadwood and reduction in basal area, compared to the baseline case (no disturbances). The occurrence of disturbances had also a positive impact on the predicted richness of insects and bats (5% and 7% increase in richness at the end of the planning horizon, respectively), due to the reduction in density, basal area and increase in deadwood amounts. The results also indicate that these patterns may be transient for some biodiversity indicators, due to the decay of deadwood over time. Increased deadwood volume and a reduction in the basal area benefited insect order richness, increasing this biodiversity indicator in damaged stands. These patterns were less pronounced in managed plots, since windstorm events led to lower amounts of deadwood and snag creation in these scenarios.

Zoom In Zoom Out Reset image size

Management had important implications to windstorm damage and economic outcomes. Increased thinning intensity led, in general, to lower economic damage (details in supplementary 5). Salvage logging buffered the impacts of disturbances, bringing the biodiversity indicators closer to the baseline scenario (0% change in richness) (figure 2). For bat richness, however, the positive influence of disturbances were largely maintained as basal area and density, important parameters for this indicator (see figure 1), were unaffected by salvage logging. Moreover, revenues generated by this management intervention enabled a marginal recovery of economic value from windfelled trees (108 to 2803 EUR ha⁻¹, depending on the management scenario).

The application of different preference scenarios resulted in a major divergence in the optimal portfolios obtained. In general, scenarios aiming at forest profitability (A and B) predominantly applied increased and BAU thinning intensities, whereas the opposite pattern occurred in biodiversity-oriented scenarios (C and D), corresponding to a 33%, 36%, 30% and 21% removal of the standing stock in each period for preference scenarios A to D, respectively. Moreover, a strong management diversification was promoted in scenarios including risk, especially in the balanced solution (C) (figure 3).

Zoom In Zoom Out Reset image size

Apart from the impacts on forest biodiversity, the application of different management portfolios led to a variation in forest profitability under windstorm risk. In the face of uncertain conditions, preference scenario (B) (risk mitigation) outperformed the baseline preference (A), with higher mean NPV (2000 and 1680 EUR ha⁻¹ EUR/ha, respectively) and VaR (1642 and 985 EUR ha⁻¹, respectively). When windstorm damage was taken into account, a reduction in standing stock occurred during the simulation period for scenarios focused on profitability, especially during the initial periods of the planning horizon (right panels in figure 3). This was accompanied by an increase in the harvesting volume, reducing the exposure of these plots to windstorm damage.

The biodiversity indicators were also sensitive to the optimal portfolios applied, i.e. the optimal solutions of preference scenarios (A) to (D) (figure 4). Bird richness decreased up to 8%, depending on the preference scenario. Increases in bat richness ranged from 17%–20% and up to 5% for insects. TreM richness reduced up to 6% in scenarios (A) and (B). As expected, the balanced (C) and biodiversity maximization preference scenarios (D) showed the best outcomes for all indicators. Additionally, the risk mitigation scenario (B) displayed slightly better biodiversity indicator values than the baseline.

Zoom In Zoom Out Reset image size

The biodiversity indicators investigated in this study displayed opposing responses to several stand structural parameters and environmental drivers, with forest bats and insects benefitting from higher temperatures, whereas bird richness was negatively correlated with this parameter. Similarly, disturbance events benefitted indicators with positive response to decreased basal area and increased deadwood availability (e.g. birds and bats). Thom et al (2017) also reported diverging responses from forest taxa to climate change impacts on a temperate forest landscape in Austria. A positive influence of climate change to the richness of ground vegetation, Hemiptera, Hymenoptera, Mollusca, Symphyta and Syrphidae was observed. Negative responses, in terms of species richness, were detected for Araneae, Carabidae and saproxylic beetles. In agreement with this study, the authors found a positive influence of disturbance events on the richness of all taxa considered.

Our results show mixed effects of climate change on the indicators studied. Rebelo et al (2010) also report increasing species richness of bats in temperate European forests in the period 2050–2060 using species distribution models. However, a reduction towards the end of the century in more extreme climate scenarios can be expected. Birds may also benefit from increased temperatures, given their ecological adaptability (Stephens et al 2016). Bale et al (2002) and Robinet and Roques (2010) indicate that increased temperatures in temperate regions may benefit some insect groups, resulting in more generations within a year, due to faster development times. Positive responses of bat and TreM richness can be expected after windstorm events, due to the creation of gaps (increasing foraging opportunity for forest bats) and microhabitats related to tree damage (e.g. broken braches and crowns) (Loeb and O'Keefe 2011; Thorn et al 2016).

Salvage logging marginally counterbalanced the economic losses caused by windstorms. Apart from the recovery of economic value, salvage logging has the goal of mitigation bark beetle outbreaks after wind disturbances (Leverkus et al 2018). Such outbreaks are likely to become more frequent under climate change, due to the positive impacts of increased summer temperatures on the development rates of insects (e.g. Seidl and Rammer 2017).

The application of salvage logging remains controversial, however, due to insufficient evidence on the effectivity of this management action to control bark beetle outbreaks and the negative effects on biodiversity (Dobor, et al., 2020). Salvage logging dampened the positive impact of disturbances on deadwood-dependent biodiversity indicators in our analysis. Lindenmayer and Noss (2006) identify this trade-off between economic goals and biodiversity generated by salvage logging and indicate that this practice leads to alterations in the structural complexity and removal of legacies, negatively affecting the taxa associated with these structures. Thus, salvage logging must be conducted in a sensible manner, optimally planning the allocation and intensity of wood removals.

The inclusion of windstorm risk and economic uncertainty in the optimization model required a strong management diversification (in terms of thinning intensities applied). Such patterns are ubiquitous in the literature, since the risk related to wood prices and disturbances can be thereby diluted (e.g. Paul et al 2019, Knoke et al 2020a). We considered solely the modification of thinning intensity to cope with these risks in the optimization model. Forests in the study area are majorly managed under a continuous-cover system, and in the time period considered in this study, thinning intensity is a determinant factor for forest profitability and risk exposure. Nevertheless, for other management systems and longer planning horizons, further options are available to support forest resilience. For example, the promotion of mixed stands and modification of rotation lengths may reduce vulnerability and exposure to biotic and abiotic risks (Díaz-Yáñez, et al 2019; Neuner et al 2015, Knoke et al 2020a). Such options may be seamlessly implemented in the framework proposed here, e.g. by creating additional management strategies as combinations of thinning schemes, rotation lengths and species composition (e.g. Augustynczik et al 2020).

Typically, the consideration of windstorm risk in forest planning problems demands an increase in the share of early harvestings and reduction of the standing stock of forest stands (e.g. Amacher, et al 2009; Augustynczik et al 2018). This modification of stand structure may cascade to the availability of habitats for different forest taxa, raising concerns on the impact of such management interventions on forest biodiversity. The results obtained in this study indicate that the risk-mitigation scenario led to poorer outcomes, in terms of species richness, compared to the biodiversity maximization or balanced scenarios. Nevertheless, it still performed better than the baseline management for the mature continuous-cover systems studied here.

We used the NPV as an economic decision criterion in the optimization models and considered a range of different preference scenarios to describe agents' priorities in the management of forest resources (scenarios (A) to (D)). Hence, our results may provide benchmark solutions for decision-makers, which can adapt these to local contexts. It is important to consider, however, that the NPV maximization assumes risk neutrality or a marginal project in the agents' portfolio, and different streams of revenues may generate different utility levels (see Knoke et al 2020b). Furthermore, the optimization model developed here may deviate from individual preferences and attitudes towards conservation and profitability goals. Although it is impractical to derive a single solution that is optimal for all preferences and stakeholders simultaneously, such standard solutions may shed light on the optimal response under various preferences combinations and provide general management guidelines. We highlight that the development of richer decision models may provide additional insights on optimal forest management strategies, e.g. considering utility maximization models tailored to local preferences, welfare maximization models and the inclusion of decision-making biases in the decision process (e.g. Dhami 2016, Augustynczik et al 2020, Knoke et al 2020b).

Furthermore, additional robustness metrics may be applied to balance different management goals and biodiversity indicators, and warrant attention in future studies. For example, Knoke et al (2020b) proposes a framework for multiple criteria robust optimization for forest planning, using ellipsoidal uncertainty sets and minimizing deviations of parameters of interest in relation to their reference scenarios. Similarly, the framework of Robust Decision-Making is well suited to address problems of this nature, allowing to search and select robust solutions that meet multiple criteria across a wide range of scenarios or minimizing the regret related to alternative management solutions (e.g. Kasprzyk et al 2013).

We addressed here the impacts of environmental drivers and disturbances on indicators related to multiple forest taxa and TreMs. Although these biodiversity indicators span over a variety of habitat requirements, this is not an exhaustive list of forest taxa and the inclusion of additional indicators is warranted. Moreover, the responses of the different indicators analyzed here was based on (log)linear empirical relationships to stand structure and environmental parameters. Some of the parameters included in our model might diverge from this assumption and display an optimal response (e.g. to temperature). Nevertheless, given a reasonable distance from the optimal condition, linear approximations are able to adequately capture the indicators' responses. We also did not include potential interactions among different indicators. In particular, the abundance of insects may define resource availability for several bat and bird species, and an increase in its population due to increased disturbance activity or warmer climate may enhance the responses of other taxa.

We highlight that the balanced and biodiversity maximization scenarios ((C) and (D), respectively) applied a same weight to the different biodiversity indicators and could, therefore, trade the loss of one indicator by a gain in the other, if preferences for indicators change. Different biodiversity indicators may have different sensitivity towards management interventions, and some respond more to climate drivers than to stand structure. Thus, additional requirements aiming to limit losses related to specific indicators may be included in our framework via the addition of constraints to limit this behavior.