Trade-offs between economic profitability, erosion risk mitigation and biodiversity in the management of uneven-aged Abies alba Mill. stands

Abstract

Key message

Multi-objective forest planning methods were used to assess the trade-offs between three ecosystem services: timber production, erosion protection and biodiversity. The use of trade-off analysis helps to define proper weights for the management objectives and evaluate the feasibility of obtaining economic profit from timber while controlling the erosion risk and maintaining biodiversity of Abies alba Mill.

Context

Forests provide several ecosystem services (ES), many of which are in trade-off. The assessment of relationships between ES helps to optimize forest management.

Aims

This study analyses the trade-offs between timber production, erosion protection and biodiversity in uneven-aged mountain forest of Abies alba Mill.

Methods

Multi-functional forest management was optimized using a simulation-optimization software. Timber production was measured with present value, erosion risk was modelled as a function of stand structure, and biodiversity was dealt with by setting a minimum number of 10 large trees (dbh > 60cm) per hectare as a requirement. Optimizations were conducted for slopes of 10%, 45% and 80% and cutting cycles of 10, 30 and 50 years. Discount rate 2% was used in all optimizations. Trade-offs were evaluated with production possibility boundaries.

Results

We found relevant trade-offs between timber production and erosion protection, which depended on the slope steepness and the length of the cutting cycle. Trade-offs were marginal on 10% slope (50-year cycle) and large on 80% slope (10-year cycle). Biodiversity constraint reduced both economic profits and erosion protection values.

Conclusion

In multifunctional mountain forest, defining proper weights for ecosystem services that are in trade-off is important for satisfying different management objectives in a sustainable way.

1 Introduction

Ecosystem services (ES) are benefits that ecosystems provide for human societies through provisioning, regulation, cultural and supporting services (MEA 2005). Particularly mountain forests deliver benefits in all four categories: timber, fuel, food, genetic resource, biodiversity, climatic regulation, soil and water regulation, erosion control and carbon storage, as well as cultural and recreational amenities (e.g. skiing and hiking environments, scenic beauty and sense of remoteness) (de Groot et al. 2010; EME 2011; Price et al. 2011). The focus in traditional forest management has been on economic profit and the provision of timber and fibre (Kareiva et al. 2007; Messier et al. 2015), while non-timber products and the other uses of forests such as recreation, erosion control and water regulation have been often ignored (Kräuchi et al. 2000; Rodríguez et al. 2006; Kurttila et al. 2018). During the recent years, the topic has gained importance resulting in several studies about surveying the supply and evaluating management alternatives for the simultaneous production of several forest ES (e.g. Blattert et al. 2017; Díaz-Balteiro et al. 2017; Mina et al. 2017; Kurttila et al. 2018). Many of forest ES have negative relationship e.g. timber and biodiversity (Lanfond et al. 2017), meaning there are trade-offs between them, when the use or the management for one service decreases the benefits supplied by another service (Rodríguez et al. 2005). In addition, positive relationships, synergies, exist when both services can be produced together e.g. timber and mushrooms (Palahí et al. 2009). Trade-offs and synergies between different ES and their spatial and temporal scales and reversibility cause challenges to managing forest for several ES (MEA 2005; Rodríguez et al. 2006; Bennett et al. 2009). Therefore, knowledge about environmental conditions and their influence on various ecosystem services on the smaller local scale are required to define reasonable provision targets for the forest ES (Mina et al. 2017).

Multi-objective forest planning aims at several simultaneous management objectives. However, managing forests simultaneously for multiple ES can be complex due to different relationships between the ES (Bennett et al. 2009). Methods to study the distribution of forest-based ES in the landscape, for instance identifying hotspot areas where a specific service has to be prioritized and where trade-offs and synergies may take place, are key factors to planning how and where forest should be managed for various potential ES (de Groot et al. 2010; Raudsepp-Hearne et al. 2010; Roces-Díaz et al. 2018). Demand for ES on forest lands is the basis to define strategic forest plans, as it should define our management objectives, for instance maximizing timber production, conserving habitats or maintaining the recreational values of the forest (Pukkala 2002). Several studies have assessed management alternatives and their effects using Pareto frontier methods e.g. between timber production and habitat availability (Mönkkönen et al. 2014); carbon storage, timber and cork production (Borges et al. 2014); timber production, biodiversity conservation and protection against avalanches, rockfalls and landslides (Lafond et al. 2017); and timber production and non-wood forest products (Kurttila et al. 2018). Numerical methods allow the comparison of different management alternatives while incorporating multiple objectives simultaneously, by using e.g. multi-objective utility functions, penalty functions or constraints (Pukkala 2002; Blattert et al. 2017; Lafond et al. 2017).

Taking into account multiple objectives in forest management is particularly relevant in mountain areas, which are ecologically vulnerable and sensitive to changes. They provide goods and services to both mountain and lowland communities, for instance fresh water, carbon storage, protection against gravitational hazards (Price et al. 2011; Blattert et al. 2017). During the past decades, the rural population of European mountainous areas has declined dramatically (Collantes and Pinilla 2004), resulting in the abandonment of traditional management practices both in agriculture and forestry, reduction in primary production (timber and food) and changes in forest management that have led to overstocked stands and increased erosion (Oliva and Colinas 2007; EME 2011). In mountain areas, the forest cover is essential for maintaining the hydrological balance, protecting the soil against sediment loss and protection against landslides and rockfalls (Price et al. 2011; Lanfond et al. 2017). In fact, soil protection is one of the most important regulating services of mountain forests but its economic effect is difficult to assess as it does not produce direct income but indirectly affects e.g. timber production and watershed protection (Price et al. 2011). Analyzing the effect of forest site and stand characteristics on soil erosion can help recognize forests susceptible to erosion and determine management practices that can enhance soil protection (Selkimäki et al. 2012).

In European mountain forest, Abies alba Mill. (silver fir) is the most important indigenous fir from both economic and ecological points of view: its timber is used for construction, furniture and plywood and the species is important for soil protection as well as for maintaining forest biodiversity as it provides habitats for fauna (Wolf 2003; Tinner et al. 2013). Although the species has shown some signs of decline in several mountain forests in Europe, in the Spanish Pyrenees, it has increased its dominance over the pines in mixed stands because pines cannot regenerate under the shading canopies of Abies alba (Oliva and Colinas 2007). As a shade-tolerant species, it is typically managed as uneven-aged forest. Trees can remain suppressed for decades and start growing well when the suppression is reduced (Blanco et al. 1997). Traditional selection methods for shade-tolerant species do not maintain high amounts of dead wood and old large trees thus reducing biodiversity (Angers et al. 2005; Rosenvald and Lõhmus 2008). In northern Spain, the typical management method for Abies alba has been cutting all trees over 35 cm in diameter every 10 years (Solano et al. 2007). However, during the past 30 years, socio-economic changes have taken place in the region and forest management practices have reduced while other activities such as recreation and tourism have been gaining importance. The lack of management in Abies alba stands has led to overstocked stands and rather regular structures, which compromises the multifunctionality of these mountain forests characterized by protective, biodiversity and recreational functions (Aunós and Blanco 2006; Tinner et al. 2013). In addition, the natural regeneration of Abies alba is threatened due to lack of management, as competition decreases regeneration in overstock stands. Therefore, silvicultural intervention would be needed to sustain the natural regeneration of silver fir (Oliva and Colinas 2007).

Uneven-aged forest management can potentially promote multiple ES (Bennett et al. 2009), but the trade-offs and synergies between ES should be taken into account. The assessment of trade-offs between ES that can show efficient management options is efficient when considering several ES’s simultaneously (Garcia-Gonzalo et al. 2015). The aim of this research was to assess the trade-offs between three ecosystem services (timber production, erosion protection and biodiversity) in uneven-aged mountain forest of Abies alba in Spanish Pyrenees mountain where these services are in decline due to lack of management. In addition, erosion occurrence has been high according to Spanish National Forest Inventory (DGCN 2005). First, erosion risk was modelled with a logistic function based on a classification tree of a previous study (Selkimäki et al. 2012). The obtained model was integrated into a stand simulation and optimization software. The management objectives were to (1) get revenues from timber, (2) reduce the risk of soil erosion and (3) enhance biodiversity.

2 Material and methods

2.1 Study area

The study area includes mountainous forest in the Spanish Pyrenees, where silver fir (Abies alba Mill.) forms pure stands or sometimes mixtures with Fagus sylvatica, Pinus sylvestris and Pinus uncinata. Forests dominated by silver fir are scattered through mountain altitudes ranging from 1000 to 2000 m a.s.l. where the mean annual precipitation is between 1000 and 1500 mm (Blanco et al. 1997; Gracia et al. 2004). Silver fir typically grows on northern aspects, on steep slopes and has high productivity of about 8.7 m³ ha ⁻¹a⁻¹ (Solano et al. 2007). Silver fir occupies about 12,400 ha, of which about 9500 ha (76%) are on slopes steeper than 40% (Gracia et al. 2004). Typical management method for shade-tolerant silver fir has been cutting all trees over 35 cm in diameter every 10 years (Solano et al. 2007). Abandonment of the mountain rural area has led also to abandonment of forest management. The regeneration of silver fir has been natural but the mortality rate of regeneration is often high due to lack of management and the consequent intense competition (Oliva and Colinas 2007).

2.2 Erosion risk modelling

The forest variables and soil erosion observations were obtained from the third Spanish National Forest Inventory (NFI) for the region of Catalonia (DGCN 2005). The NFI is based on a systematic sample of permanent plots distributed on a 1 × 1 km square grid and about 10-year measurement interval. For every plot, the presence of erosion has been visually evaluated up to 60 m from the centre of the plot and coded into six categories: no erosion, surface erosion, rill erosion, V-shape gully erosion, U-shape gully erosion and landslides. Surface erosion was observed in 689 NFI plots and the number of no-erosion plots was 8872. For modelling, we selected pure or almost pure Abies alba (> 80% of basal area) plots, which presented no erosion or surface erosion, resulting in 56 plots (Table 1). The slope of these plots was obtained from a digital elevation model (DEM) of 30 × 30 resolution. These data were used to fit a logistic model for the probability of erosion. The site and stand variables used as predictors are based on classification tree results of a previous study (Selkimäki et al. 2012).

Table 1 Means and ranges of site and stand characteristic in non-eroded and eroded plots of pure Abies alba

The model for the probability of erosion was as follows:

$$ {E}_{\mathrm{risk}}={\left(1+{e}^{-\left[-12.356+0.04 BA+0.301H+0.134\mathrm{dead}\%+0.13\mathrm{slope}\right]}\ \right)}^{-1} $$

(1)

where E_risk is the probability of erosion in a given stand, BA is the basal area (m² ha⁻¹), H is the mean height on the trees (m), dead% is the percentage of standing dead trees and slope is the steepness of slope in degrees of the inventory plot calculated from the DEM (later converted to percent for the optimization software). The model had good predictive power (Nagelkerke R² = 0.448). The overall percentage of correct classifications was 82.1% (with the cut value of 0.5), the Hosmer and Lemeshow test statistic was 0.462 and the area under the receiver operating characteristic curve (ROC) was 0.836. For optimizations, the probability of erosion (E_risk) was converted to “erosion protection” (EP) value calculated as 1-E_risk.

2.3 Optimized cases

We studied multi-objective forest management of uneven-aged Abies alba stands aiming at (1) profit from harvested trees measured by present value (PV), (2) reducing the risk of soil erosion (in other words, maximizing the erosion protection) and (3) maintaining a minimum number of 10 large trees per hectare for biodiversity reasons. We simulated the dynamics and management of uneven-aged Abies alba stands on three slopes: moderate 10%, steep 45% and very steep 80%. The selection of slopes aimed to present the natural distribution range of the species as well to cover high erosion risk areas on steep mountain slopes. For each slope, management was optimized with three cutting cycles: 10, 30 and 50 years. The selection of cutting cycles was done starting from the typical 10-year cycle and until longer cycles in order to reach the large diameter classes over 60 cm. The harvesting was done at the end of the cutting cycles as selective felling aiming at returning the original diameter distribution to maintain the uneven-aged structure of the stand.

The site parameters represented the average conditions of typical Abies alba stands. Elevation was assumed to be 1600 m a.s.l., aspect north and soil cambisol. These variables are used as predictors in the growth models. The growth index (a predictor of diameter increment) was assumed to be 1, which means that the growth level corresponded to the mean growth of trees that were used in growth modelling (tree measurements in the 2nd and 3rd NFI plots in Catalonia). The growth index measures site quality and is calculated as the mean ratio between measured past growths of sample trees and the predicted past growths of the same trees (Trasobares and Pukkala 2004a).

2.4 Optimization method

Multi-objective optimization methods and simulations have been used in several studies for assessing the trade-offs and synergies between forest ecosystem services (e.g. Mönkkönen et al. 2014; Lafond et al. 2017; Mina et al. 2017). Multi-criteria decision analysis is commonly used in the field of forest management planning when multiple objectives must be evaluated simultaneously dealing typically with both economic and ecological issues such as timber production versus biodiversity or non-timber products (e.g. Pukkala 2002; Díaz-Balteiro and Romero 2008). The utility weighting takes into account different objectives and makes the simulation and optimization software a flexible management tool (Blattert et al. 2017), allowing the analysis of the trade-offs between the objective variables.

We used the simulation-optimization software Rodal (Palahí and Pukkala 2003) to optimize the structure and management of stands. The software uses individual-tree models for diameter increment, tree height and survival and stand level models for ingrowth to simulate stand dynamics. We used the growth and yield models developed for Abies alba in the Spanish Pyrenees by Trasobares et al. (unpublished). These models are based on data from repeated tree measurements made in the 2nd and 3rd NFI (ICONA 1993a and 1993b; DGCN 2005) for the region of Catalonia (Appendix 1). The aim was to optimize the steady-state diameter distribution of the stand after cutting. The optimization algorithm was the direct search method of Hooke and Jeeves (1961), which is commonly used in stand level analyses (e.g. Roise 1986; Palahí and Pukkala 2003; Trasobares and Pukkala 2004b). It is an unconstrained optimization method for minimizing or maximizing an objective function. Constrained problems can be formulated by using penalty functions (Bazaraa et al. 1993). The diameter distribution was described with the Weibull function (Bare and Opalach 1988) and the optimized decision variables were the following four parameters that describe the diameter distribution of trees:

• 1. Total number of trees per hectare

• 2. Weibull parameter b

• 3. Weibull parameter c

• 4. Maximum retained diameter (cm)

The simulation-optimization system includes two software components: optimization algorithm where the Hooke and Jeeves algorithm is used to find the optimal values for the decision variables and simulation software where the stand development is simulated and the values of the objective variables calculated (more detailed description in Appendix 2). First, the optimization software starts from a candidate diameter distribution of the trees which was defined by the Weibull parameters (500 trees ha⁻¹, Weibull b equal to 15, Weibull c equal to 1, maximum retained dbh 50 cm). After that, the stand is let to grow for one cutting cycle (10, 30 or 50 years) after which a cutting is simulated aiming at returning the original diameter distribution. The erosion protection value (EP) was calculated at the beginning and every 10 years until the end of the cutting cycles. Utility function was used to maximize timber production and erosion protection. Penalty functions were added to the objective function to guarantee the sustainability of diameter distribution and to force a sufficient number of large trees:

$$ \mathrm{OF}=\mathrm{Utility}-{P}_{\mathrm{d}}-{P}_{\mathrm{lt}} $$

(2)

where OF is the objective function, P_d is the penalty due to a difference between the initial diameter distribution and the post-cutting distribution after one cutting cycle. P_lt is the penalty due to insufficient number of large trees.

The utility function was as follows:

$$ \mathrm{Utility}={w}_1{u}_1\left(\mathrm{PV}\right)+{w}_2{u}_2\left(\mathrm{EP}\right) $$

(3)

with

$$ {u}_{\mathrm{i}}=\frac{\left({q}_{\mathrm{i}}-{\min}_{\mathrm{i}}\right)}{\left({\max}_{\mathrm{i}}-{\min}_{\mathrm{i}}\right)} $$

(4)

where w₁ and w₂ are the weights and u₁ and u₂ are the sub-utility functions of the two objective variables; u_i is the sub-utility function, q_i is the quantity, min_i is the minimum value and max_i is the maximum value of the objective variable. The sum of the weights was equal to 1, and the sub-utility functions scale also the objective variables to range 0–1. We used five combinations of weights for PV and EP, representing the preferences assigned to each management objective. These combinations were 1/0, 0.75/0.25, 0.5/0.5, 0.25/0.75 and 0/1. The production possibility boundaries for PV and EP were calculated by changing the weights to obtain different points of the curves. The solutions are therefore Pareto-optimal, meaning that it is not possible to increase one criterion without decreasing the other (Branke et al. 2008).

2.5 Economic parameters

Economic parameters included timber prices (Eq. 5), harvesting (Eq. 6) and transportation (Eq. 7) costs and the discount rate, which were used to calculate the PV (Eq. 8). In all optimizations, the discount rate was 2%. It has been proposed by Díaz Balteiro and Prieto Rodríguez (1999) as it is very close to the rate of the return of public depth in Spain. Data for timber price as well as harvesting and transportation costs were obtained from the Regional Forest Service. The road side timber price was calculated from dbh according to Solano et al. (2007) (Eq. 5). The harvesting (Eq. 6) and transportation (Eq. 7) cost functions have been presented in earlier studies (Solano et al. 2007; González-Olabarria and Pukkala 2011; Pascual et al. 2018). In this study, distance to road was assumed to be 500 m in all cases to represent the average forwarding distance. Additionally, the entry cost for cutting was set at 100 € ha⁻¹. Variables affecting the costs were tree diameter, slope and distance to road. Timber price and the effect of slope on costs are visualized in Fig. 1.

Fig. 1

Timber price and harvesting (H_cost) and transportation cost (T_cost) as a function of diameter for three different slopes (10, 45 and 80%)

The roadside timber price and the harvesting and transportation costs were calculated using the following functions:

$$ {T}_{\mathrm{price}}=\frac{{\mathrm{dbh}}^2}{1.69+0.52\mathrm{dbh}+0.019{\mathrm{dbh}}^2}\kern0.5em $$

(5)

$$ {H}_{\mathrm{cost}}={e}^{3.406-0.568\ln \left(\mathrm{dbh}\right)+0.01\mathrm{slope}} $$

(6)

$$ {T}_{\mathrm{cost}}={e}^{4.396-0.110\mathrm{dbh}+0.012\mathrm{slope}+0.001\mathrm{distance}} $$

(7)

where T_price is the price (€ m⁻³), H_cost is the harvesting cost (€ m⁻³), T_cost is the transportation cost (€ m⁻³), dbh is the diameter of the harvested tree measured at 1.3 m (cm), slope is the slope of the stand (%) and distance is the distance to nearest road (m).

The PV was calculated as:

$$ {\mathrm{PV}}_{\mathrm{t}}=\frac{{\mathrm{VG}}_{\mathrm{t}}}{\left({\left(1+i\right)}^t-1\right)} $$

(8)

where PV_t is the present value for a t-year cutting cycle (€ ha⁻¹), VG_t is the net income obtained from trees harvested every t years (€ ha⁻¹), i is the discount rate and t is the length of the cutting cycle.

The aim was to keep the stand in a steady-state structure. If the stand structure differed at the end of cutting cycle (after simulating the cutting) from the original structure at the beginning of the simulation period, a penalty was subtracted from the objective function (Eq. 2) (Bazaraa et al. 1993). The penalty was calculated as follows:

$$ {P}_{\mathrm{d}}=0.01\sum \limits_{i=1}^4{\left|{N}_{{\mathrm{end}}_{\mathrm{i}}}-{N}_{{\mathrm{i}\mathrm{ni}}_{\mathrm{i}}}\right|}^2 $$

(9)

where N_endi, N_inii are the numbers of trees per hectare of the ith diameter class. The four diameter classes used in the penalty function were < 20 cm, 20–34.99 cm, 35–59.99 cm and ≥ 60 cm.

The biodiversity constraint was also implemented as a penalty function in order to guarantee a sufficient number of large trees. The presence of large trees has been suggested as an indicator of biodiversity in the forests of Catalonia (Camprodon 2001) and a threshold of ten trees per hectare has been suggested to indicate of abundance of habitat trees (Bütler et al. 2013). Therefore, the example of biodiversity constraint was defined as minimum number of 10 trees ha⁻¹ with dbh > 60 cm. The penalty was set so that every missing tree reduced the utility function by 20%. The number of large trees was checked at the beginning and after every 10-year interval until the end of the cutting cycle. Penalty was calculated always when the number of large tree was less than the required minimum.

3 Results

3.1 PV and EP values in different scenarios

According to the results (Table 2), the highest PVs values were obtained when profit from timber was the only objective, decreasing gradually as more importance was given to erosion protection. The highest PV value was obtained for 10% slope. PV decreased as the steepness of the slope increased due to increasing harvesting and transportation costs. In addition, the PV values decreased as the cutting cycles were longer, partly due to the discounting effect. The lowest PV values were obtained when all weight was on erosion protection. In these cases, the PV values were clearly negative (results not shown) due to the small average size of the optimal stand structure and, as a consequence, high harvesting and transportation costs per cubic meter and low roadside value of harvested trees. When the biodiversity constraint was included, leading to the maintenance of many large trees in the stand, both erosion protection and PV values were lower compared with solutions obtained without the biodiversity constraint.

Table 2 Volumes (V) and number of trees (N) before and after harvesting and the net incomes, PV and EP by slope, cutting cycle (CC) and the weights of management objectives

In all the scenarios, the erosion protection values improved when more weight was given on EP. Erosion protection values were the highest on 10% slope and were getting smaller as the steepness of the slope increased. Ten-year cutting cycles had the largest EP values on all three slopes, whereas longer cutting cycles showed slightly smaller EP values. The scenarios including the biodiversity constraint resulted in lower EP values compared with unconstrained scenarios, especially on the steepest 80% slope.

3.2 Trade-offs between erosion protection and timber production

The production possibility boundaries represent efficient (Pareto-optimal) production of PV and EP with and without the biodiversity constraint. The production possibility boundaries revealed that there are increasing trade-offs between EP and PV when slope gets steeper and when the cutting cycle gets shorter. In all three slopes, the curve has a concave shape, suggesting increasing rate of transformation (Fig. 2). This means there are losses in timber production in exchange for better erosion protection values. Very concave curves mean that for some combinations in the left part of the curve, trade-offs are indeed small. For a 50-year cutting cycle for example enhancing the erosion protection from 0.4 to 0.7 results in only a small reduction of PV. However, further increase in erosion protection from 0.7 to 0.8 leads to large reductions in PV. Therefore, the trade-off varies according to the starting situation of the weight combination for the objectives. Combining all the curves into the same graph (Fig. 3) shows that on 80% slope, the curves are very concave, suggesting that a joint production of the two ES would almost certainly be optimal. The trade-off curve shows that on slope 10% there are only minimal improvements in erosion protection at the expense of losses in money. In those conditions, prioritizing only one ES might also be the most efficient management choice depending on the importance of PV and EP. Adding the biodiversity constraint did not alter the main trends in trade-offs. PV and EP values were, in general, lower with biodiversity constraint. However, the weight combination of 0.25 PV and 0.75 EP resulted in most cases slightly higher profit from timber but still lower erosion protection compared with no constraint scenarios. The production possibility curves for scenarios including the biodiversity constraint demonstrated that the biodiversity objective resulted in lower erosion protection values for all weight combinations of PV and EP.

Fig. 2

Production possibility boundaries for 10-, 30- and 50-year cutting cycles on a 10%, b 45% and c 80% slopes representing the trade-offs between erosion protection and present value (PV) without and with biodiversity constraint (continuous presence of at least 10 trees per ha with dbh > 60 cm). Note that the x-axis is cut in diagrams a and b

Fig. 3

Production possibility boundaries a without and b with biodiversity constraint for different slopes and cutting cycles

3.3 Optimal diameter distributions

When the management objective was not constrained by the biodiversity objectives, the optimal diameter distributions were towards smaller trees the more weight was given to erosion protection. Higher weights for timber production resulted in wider diameter distributions (Fig. 4). Scenarios including the biodiversity constraint had much wider optimal diameter distributions and less trees removed per diameter class than in cases where PV and EP were maximized without the biodiversity constraint. Logically, removed volumes (Table 2) increased in all the scenarios when the cutting cycles were longer. Both the number and the volume of harvested trees were smaller when the biodiversity constraint was included.

Fig. 4

The optimal diameter distribution for 45% slope with 50-year cutting cycle without constraint (left) and with biodiversity constraint (right)

4 Discussion

This study presents a stand level assessment of three ecosystem services in Abies alba mountain forest, namely economic profit from timber production, erosion protection and biodiversity, enhanced as structural diversity of the stand by maintaining a minimum of 10 large trees per hectare. We assess the trade-offs between timber production and erosion protection and the inclusion of biodiversity constraint on different steepness of slopes and cutting cycles. The results suggest that in the case of multifunctional mountain forest, defining proper weights for different ecosystem services is important to satisfy management objectives in a sustainable way and still provide livelihood for rural people in mountainous areas.

Our results show how the trade-offs between profit from timber and erosion protection changed according to the steepness of the slope and the length of the cutting cycle (Fig. 3). As the slope increased, the profit from timber reduced due to increasing harvesting and transportation costs. Harvesting on steep slopes is challenging economically and for environmental and safety reasons (Visser and Stampfer 2015). Including the biodiversity constraint reduced economic profitability (PV) in most cases since part of the large trees remained in the stand. Erosion protection values were also lower when biodiversity constraint was included, as the denser stand with large trees had a higher risk of erosion especially on the steepest slopes.

The results for 10% slope suggest that the forest owner could decide between enhancing biodiversity or optimizing timber production since the erosion risk is always small, meaning that timber production and biodiversity can be pursued simultaneously (Fig. 2a). Giving much weight to erosion protection results in great loss in profit, but if biodiversity is pursued together with timber production, only small losses in profit can be expected. When the terrain becomes steeper, the production possibility curves turn more concave suggesting increasing rate of transformation. For example, with a slope of 45%, choosing weights 0.75/0.25 or 0.5/0.5 for PV/EP would result in better erosion protection than when giving all weight to PV, but without having a strong impact on timber revenue (Fig. 2b). On the steepest 80% slopes, setting a biodiversity objective has a strong negative effect on erosion protection, and it should be carefully considered which objective to focus on (Fig. 2c). Timber production and biodiversity on the steepest 80% slopes lead to low level of erosion protection. If the objectives are timber profit and erosion protection, weight 0.25 for PV and 0.75 for EP would result in good erosion protection and still produce some revenues from timber. However, if high structural diversity is also important, the trade-off curves are closer to origin and the erosion protection values remain low even with full weight on EP. This means that, on the steepest slopes, pursuing the three ES at the same time may not be recommendable.

Based on our results, uneven-aged forest management with younger stand compositions resulted in better erosion protection values. However, when all weight was on erosion protection, the stand turned rather even-aged with trees only in the smallest diameter classes. This kind of stand structure produced a negative profit due to the shape of the price and cost functions resulting in unprofitable harvesting. Especially for small trees with dbh < 20 cm, the timber price is low and the harvesting and transportation costs are high. The ingrowth of trees passing the 7.5 cm diameter limit in our simulations is dependent on the basal area of the stands. When larger trees are missing, the models overestimate the number of new trees. In scenarios without biodiversity constraint and all weight on erosion protection, the optimal diameter distribution had many small fast-growing trees and almost no large trees. Therefore, the growth models, which were developed for uneven-aged stands, did not work well in these scenarios predicting too large ingrowth and volume increment.

In the Spanish Pyrenees, around 90% of the Abies alba stands had been managed in the past but currently the lack of forest management is decreasing natural regeneration due to overstocked stands (Oliva and Colinas 2007). The management in these forests should focus on removing part of the valuable larger trees (Solano et al. 2007) for which the harvesting and transportation costs are reasonable. This would promote regeneration, provide profit for forest owners and would enhance the erosion protection by allowing the development of the shrub layer and tree saplings. This would require careful planning of how harvesting and transportation should be implemented as they could damage the vegetation and soil surface, which could increase erosion (Edeso et al. 1999). On the other hand, large trees were used as indicators for biodiversity as they provide ecological services for different species (Bütler et al. 2013; Jones et al. 2018) and are important for forest recreation and culturally (Blasco et al. 2009; Tabbush 2010). There are several alternative indicators of biodiversity such as the amount of deadwood or species richness (Gao et al. 2015). Since our study concerns pure Abies alba forest, only the amount of deadwood could improve our assessment of structural diversity, as the required models and modification on the management systems allow the selection of individual forest legacies, to properly apply retention forestry recommendations (Gustafsson et al. 2012, 2019).

Optimizations including biodiversity constraint resulted, as expected, in lower economic profitability as several previous studies on ecosystem services have shown (e.g. Lafond et al. 2017; Mina et al. 2017). Including the biodiversity constraint, that encouraged to leave a number of large trees unharvested, resulted also in lower erosion protection values. This tendency can be explained as from higher canopies raindrops can reach erosive speed as well as raindrops tend to be larger if compared with free rainfall (Brandt 1987; Calder 2001). Additionally, the stand basal area and the density of understorey vegetation have shown strong negative correlation (Coll et al. 2011), implying that leaving large Abies alba trees on the stand will limit the amount of understorey vegetation (Blanco et al. 1997; Oliva et al. 2009), which plays an important role in the protection of soil against erosive forces (e.g. Hartanto et al. 2003; Razafindrabe et al. 2010). The role of understory vegetation as an additional protective layer against raindrop impact and the surface runoff (Calder 2001; Hartanto et al. 2003) was visible in the region, where the mean shrub cover was 19% in non-eroded and 12% in eroded Abies alba plots, according to the 3rd NFI data.

The assessment of erosion protection as one ecosystem service in the studies has often been based on potential soil loss calculations with the USLE (García-Nieto et al. 2013) or RUSLE models (Anaya-Romero et al. 2016), or on the general assumption that dense forest on steep slopes > 30% supplies this service (Roces-Díaz et al. 2018). The role of forest providing this service on steep slopes is certainly important although the risk of soil erosion is high even under forest cover due to orographic and climatic reasons. Typically, mountain forests receive the highest precipitation of both water and snow making these areas risky for erosion, landslides, snow avalanches and rockfalls (Kräuchi et al. 2000; Price et al. 2011; Mina et al. 2017). Therefore, erosion protection should be assessed on the small scale, which takes into account the forest type and stand structure, variables which have been recognized to play an important role against gravitational hazards (snow avalanches and rockfalls) in the mountain forests in central Europe (Rammer et al. 2015; Blattert et al. 2017) and have shown to influence surface erosion in other Mediterranean areas (Selkimäki et al. 2012).

Managing forest for multiple purposes needs studies about the trade-offs of ES to help choosing appropriate management options. Our results suggest that economic profits can be produced together with erosion protection by choosing proper weights for the objectives depending on the location. Similar results are suggested by Rammer et al. (2015), who studied the effect of different forest management options on rockfall protection efficiency. They found out that in business as usual management (shelterwood cutting), high economic returns were related to lower protection values whereas maximizing both objectives simultaneously (protection forest management) resulted in both high income and good protection level. This tendency was also identified in the case of fire risk in Catalonia, where the traditional management options aiming at maximizing income or reducing fire hazard were found to be non-efficient when compared with forest plans that considered a combination of both objectives (González-Olabarria and Pukkala 2011). Management in mountain forest should be regionally adapted to site and stand conditions due to high variability in environmental conditions (Mina et al. 2017). Managing mountain forests and enhancing their capabilities to produce ecosystem services are important for rural people and for mitigating the effects of climate change (MEA 2005; Mina et al. 2017). When deciding the management strategy in our case of Abies alba, it should be analyzed whether the production of all three ecosystem services simultaneously is feasible, or would it be better to modify management objectives according to site and stand conditions. For example, on very steep slopes, more weight could be given to erosion protection and reduce the expectations from timber and biodiversity. Since most Abies alba forests are located on steep slopes, Fig. 5 shows as an overall suggestion, where the three different ES should be pursued. However, timber harvesting should be carefully planned as harvesting operations and forest road construction can significantly increase sediment loss (Edeso et al. 1999; Croke et al. 2001; Madej et al. 2006). On the other hand, selective cuttings can enhance regeneration and ingrowth, which decrease the risk of surface erosion. In mountain forest, stand level models have been proven suitable (Rasche et al. 2011) and further research could evaluate erosion protection considering other forest types and tree species.

Fig. 5

Illustration of the distribution ranges of Abies alba on different slope steepness (based on NFI plots) in the Spanish Pyrenees and suggested management objectives for three ecosystem services (timber, erosion protection and biodiversity) based on the results (font size and order defines the priorities suggested)

In strategic planning, the trade-offs between different ecosystem services should be evaluated together with management priorities, restrictions and recommendations. Before implementing studies at regional or even landscape level (Raudsepp-Hearne et al. 2010), it is necessary to understand which type of management is the most suitable to attain management objectives that are complex and usually in trade-off.

5 Conclusion

In multifunctional mountain forest, defining proper weights for forest ecosystem services that are in trade-off is important for satisfying different management objectives in a sustainable way. The results of this study show that adequate forest management can enhance erosion protection while still providing revenues from timber. The importance of careful planning is emphasized, since forest structure is not the only variable affecting the soil erosion risk, but also harvesting operations and the construction of forest roads. Stand level analyses on the feasibility of obtaining different ecosystem services should be conducted to find meaningful combinations of management objectives.

Appendices

Appendix 1

Trade-offs between economic profitability, erosion risk mitigation and biodiversity in the management of uneven-aged Abies alba stand.

Growth and yield models of Abies alba:

The individual-tree diameter growth model:

$$ id10=\exp \Big(1.6681+0.6131\times \ln (d)-0.0209\times d-\ln (G)\times 0.3349-0.0494\times \frac{\mathrm{BAL}}{\ln \left(d+1\right)}-0.0022\times {\mathrm{ELE}}^2+0.3743\times \ln (GI)\times 1.1771 $$

Height model:

$$ h=\frac{\left(50.9537-0.8180\times \mathrm{ELE}-\frac{1.1734}{\mathrm{GI}}\right)}{\left(1+\frac{29.6656}{d}+\frac{130.23}{d^2}\right)} $$

Ingrowth models:

$$ \mathrm{ING}=247.455-50.149\times \ln \left({G}_{\mathrm{abies}}\right) $$

$$ \mathrm{Din}=9.583-0.13\times \ln (G) $$

The density-dependent survival model:

$$ P\left(\mathrm{survive}\right)=\frac{1}{1+\exp \Big(-\left(-4.703+\frac{37.635}{d}+0.35\mathrm{ELE}-0.03G+0.121h+0.885\ln \left(\mathrm{GI}\right)\right)} $$

Id10 is the future diameter growth (cm in 10 years), d is the diameter at breast height (cm), G is the stand basal area (m² ha⁻¹), BAL is the total basal area of trees larger than the subject tree (m² ha⁻¹), ELE is the stand elevation (100 m), GI is the stand growth index, h is the tree height (m), ING is the ingrowth (number of trees ha⁻¹) during a 10-year growth period, G_abies is the basal area of Abies alba (m² ha⁻¹) and DIN is the mean dbh of ingrowth trees (cm) at the end of a 10-year growth period.

Appendix 2

1.1 Optimization method

The study optimized the post-cutting diameter distribution of the stand. The distribution was described with the two-parameter Weibull function, total number of trees per hectare (N) and the largest diameter retained in cutting (D_max.

The probability density function of the two-parameter Weibull distribution is as follows:

$$ f(d)=\left(\frac{c}{b}\right){\left(\frac{d}{b}\right)}^{c-1}\exp \left[-{\left(\frac{d}{b}\right)}^c\right] $$

where d is the diameter, b is a location parameter and c is a shape parameter. Parameters b and c were used to calculate the frequencies of different diameter classes (using class midpoint diameter). Then, the frequencies of trees greater than D_max were set to zero. The frequencies of the remaining were scaled so that their sum was equal to N.

The growth of a stand having the above-described post-cutting diameter distribution was simulated to the end of the cutting cycle i.e. 10, 30 or 50 years. After the simulation, partial cutting was simulated so that the numbers of trees in different diameter classes were returned to the same values as in the beginning of the simulation. Incomes and costs were calculated from the diameters and numbers of removed trees. The net present value was calculated by assuming that the same cutting cycle is repeated to infinity.

The method of Hooke and Jeeves (HJ) was used to find the optimal values for N, D_max and the two Weibull parameters (b and c). Combination of these four variables is called as solution vectors, denoted as x. The HJ method needs an initial solution vector x₁. The search begins with exploratory search in the directions of coordinate axes (one of the four optimized variables is altered at a time). After completing one round of exploratory searches, the algorithm goes to pattern search if the exploratory search is successful. The pattern search alters the values of more than one variables simultaneously. If exploratory search cannot improve the solution, the step size used in exploratory search is halved and the search is repeated. The search is stopped once the step size becomes smaller than a predefined stopping criterion.

Let d₁,…,d_m be the coordinate directions and let y₁ = x₁ and k = j = 1.

The algorithm works as follows (Bazaraa et al. 1993):

Exploratory search

• Step 1: If f(y_j + Dd_j) > f (y_j) (success), let y_j + 1 = y_j + Dd_j and go to Step 2. Otherwise if f (y_j–Dd_j) > f (y_j) (success), let y_j + 1 = y_j–Dd_j and go to Step 2. Otherwise let y_j + 1 = y_j and go to Step 2 (no improvement found in direction j).

• Step 2: If j < m, replace j by j + 1 and repeat Step 1 (go to next decision variable). Otherwise go to Pattern search if f (y_m + 1) > f (x_k) (at least one successful change detected in the directions of coordinate axes). If f (y_m + 1) ≤ f (x_k) go to Step size reduction.

Pattern search

• Step 3: Let x_k + 1 = y_m + 1 and let y₁ = x_k + 1 + (x_k + 1–x_k). Replace k by k + 1, let j = 1, and go to Step 1.

Step size reduction

• Step 4: If D ≤ ε, stop. Otherwise replace D by D/2. Let y₁ = x_k, x_k + 1 = x_k, j = 1, replace k by k + 1, and repeat Step 1.

The method has three parameters: initial steps size (D), stopping criterion (ε) and α. The initial step size was different for different decision variables; it was 0.1 times the allowed range of the variable. The stopping criterion was equal to 0.01 times the initial steps size. The search was stopped once the step size was smaller than the stopping criterion for every decision variable. Parameter α was taken as 1.0.