Forest management for optimizing soil protection:a landscape-level approach

Ana Raquel Rodrigues,Susete Marques,Brigite Botequim,Marco Marto and José G.Borges

Abstract

Keywords: Ecosystem services, Erosion, Optimization, Silvicultural practices

Background

Soil erosion has been identified as the most common cause of land degradation, driving worldwide actions to stop and reverse it (Lal 2014; FAO 2017; Panagos and Katsoyiannis 2019). It can be defined as the net longterm balance of sediments detachment and transport from its original location, thus involving soil structure destruction, nutrient loss and decrease in water storage capacity (FAO 2019). Since soil functions are intrinsically related to ecosystem sustainability, soil degradation may transversely affect all ecosystem services, including:support of primary production, nutrient cycling and geochemical processes; regulation of climate and water cycle; the provision of food, water, fibre, fuel and raw earth materials; and the preservation of aesthetic, spiritual and heritage values (FAO and ITPS 2015).

Soil erosion by water impacts and runoff are primarily influenced by erosivity, which is defined by both rainfall amount and intensity (Panagos et al. 2015a). Consequently, water erosion is particularly susceptible to increase under climate change scenarios that indicate shifts in the amount, and temporal and spatial distribution of precipitation (Bredemeier 2011; IPCC 2014a).

Forest ecosystems are commonly associated with positive impacts in the ecosystem water-related regulatory services, including soil protection from erosion (Durán Zuazo and Rodríguez Pleguezuelo 2008; Keenan and Van Dijk 2010; Bredemeier 2011). Yet, this crucial relation is often disregarded in forest management plans.Forest management planning is currently facing substantial challenges worldwide, mainly concerning the requirement for a range of often-conflicting purposes(Cosofret and Bouriaud 2019; Baskent 2020). Therefore,there is a growing need for the development and application of management tools that can efficiently depict possible trade-offs and synergies between different ecosystems services (ES).

Identifying and understanding the best course of action to manage forest ecosystems resources using mathematical optimization tools is relatively recent (early 1960s), but evolved significantly to the present day along with decision sciences and computer technologies (Vacik and Lexer 2014; Kaya et al. 2016). Some studies have considered soil protection as a crucial forest ecosystem service for stand-level management planning (e.g., Kele? and Ba?kent 2011; Rodrigues et al. 2020; Selkim?ki et al.2020), but few mathematical models have been specifically developed for optimization of water-related management problems at landscape-level (e.g. Dumbrovsky and Korsu 2012; Fotakis et al. 2012). Yet, while most landscape-level forest management optimization studies and applications focus on wood production or economic performance maximization, possible trade-offs with other relevant ecosystem services are frequently overlooked(Baskent and Kücüker 2010; Kaya et al. 2016).

The present work describes the design and application of a linear programming (LP) resources capability model(RCM), aiming at solving a typical forest landscape management problem: minimizing soil erosion while ensuring timber provision. A 14,765-ha study area was considered, with several alternative management prescriptions, under two climate change scenarios and over a 90-year time span, divided in 10-year periods. We hypothesized that, given an appropriate set of alternative silvicultural practices, landscape-level forest management decisions can be conducted to reduce soil loss risk without hampering wood provision goals. Timber provision and soil protection are acknowledged as conflicting ES, so trade-offs between them are analysed. By evaluating results of a practical application, for a real case study area, the simulated forest management models are compared and discussed, regarding their performance towards soil conservation and timber provision.

Methods

Study area and forest management programs

The Vale do Sousa, located in the north-western region of mainland Portugal (Lat: 41.1343, Lon: -8.2951), fairly represents the country’ forest context, namely in terms of species composition and ownership structure (Novais and Canadas 2010; ICNF 2013). It extends over 14,837 ha, 14,765 of which correspond to forest land, divided in 1373 stands that are managed by ca. 376 forest owners.

Most of the study area soils are developed over schist and granite bedrocks, being mainly classified as Umbric Leptosols and Leptic Regosols (IUSS Working Group WRB 2015). The climate is typical Mediterranean, with mild and dry summers and wet winters. Local mean annual rainfall is 1194 mm and mean annual air temperature is 13.8°C (1981–2010; Palma 2015).

Current forest management comprise mixtures of maritime pine (Pinus pinaster Ait.) and eucalypt (Eucalyptus spp.) with different proportions, pure chestnut stands (Castanea sativa Mill.) and pure eucalypt plantations. Recent wildfires and stakeholders concerns on the long-term sustainable provision of a variety of forest ecosystem services,raised the opportunity to consider alternative management models (Agestam et al. 2018;Marques et al. 2020a; Marques et al. 2020b), including:lower density maritime pine plantations, pure pedunculate oak (Quercus robur L.) and pure cork oak stands (Q.suber L.).

Several species-specific silvicultural treatments were considered in each forest management model, for a 90-year planning horizon. All stands follow an even-aged regime that includes harvests, except for cork oak, to reflect regional silvicultural standards (Table 1).

The complete set of prescriptions was matched to the landscape management units, according to their biophysical aptitude (e.g., soil type and depth, altitude) and location (policy regulations). When a conversion to a new management model was a possibility, it was scheduled to occur at the year of the first clearcut harvest of the current species. In result, most of the area was considered suitable for maritime pine plantations, both pure and mixed (13,616 ha, i.e., 92% of the landscape); chestnut, pedunculate oak and cork oak alternative management models were considered for nearly 46% of the landscape (ca. 6786 ha); and pure and mixed eucalypt stands were a possibility in about 44% of the area (about 6460 ha).

Annual stand biometric variables (e.g., hdom, dbh, volumes and biomasses), for all possible combinations of management unit and prescription, were simulated with the StandSIM-MD module (Barreiro et al. 2016), for eucalypt (GLOBULUS 3.0; Tomé et al. 2006), maritime pine (PINASTER; Nunes et al. 2011) and cork oak(SUBER; Tomé 2004; Paulo 2011; Paulo et al. 2011). An online simulation tool (https://manuelar.shinyapps.io/Quercusrobur_SimGaliza/) was used for pedunculate oak (Gómez-García et al. 2015, 2016) and yield tables were used for chestnut stands (CASTANEA; Patrício 2006; Filipe 2019). Stand-specific characteristics that affect yield levels were considered in all simulations,assigning a site index to each combination ofmanagement unit and prescription. The site indexes were determined for each stand according to the available forest inventory data from the study area (i.e.,height of dominant trees at a species-specific age).

Table 1 Silvicultural models characteristics and growth models applied in simulation

Climate scenarios

Considering a 90-year planning horizon, from 2017 to 2106, annual climate variables for the region were estimated with the CliPick online tool (Palma 2015, 2017;http://www.isa.ulisboa.pt/proj/clipick/), based on KNMIRACMO22E models (van Meijgaard et al. 2012).

Local climate normals (1981–2010), including annual precipitation and average air temperature, were repeated throughout the planning horizon to depict a “business as usual” scenario (BAU). The other scenario assumed local conditions will follow RCP8.5 global climate change projections, with increases of 2.36°C and 193 mm, in mean annual temperature and annual precipitation, respectively (IPCC 2014b). Annual average values of temperature and precipitation are represented in Fig. 1 for both scenarios.

To our knowledge,process-based models are not available for most of the studied tree species, so climate change influence on forest growth was considered by adjusting stand biometric variables to the expected yield changes. Specifically, according to Santos and Miranda(2006), a temperature increases of 2.36°C under the REF scenario will produce a cork yield increase of 3.52%, and a 4.2% increase in timber yield, along the 90-year planning horizon.

Soil erosion estimates

The Revised Universal Soil Loss Equation (RUSLE;Renard et al.1991)was used to estimate potential annual soil loss by surface runoff for each management unit under each possible prescription and climate scenario.This widely used long-term soil erosion modelling tool(Guo et al. 2019), calculates soil losses (A, Mg·ha-1·yr-1)as the multiplication of six factors (Eq. 1):

Fig.1 Annual average air temperature(dashed lines)and accumulated precipitation(solid lines)in the Vale do Sousa along the planning horizon,for current(black lines,BAU)and locally adapted RCP8.5 scenario(grey lines,REF)

R is the erosivity factor, estimated according to Constantino and Coutinho (2001) equations for the region(MJ·mm·h-1·ha-1·yr-1). The erodibility factor, K(Mg·ha-1MJ·mm-1), was calculated using GIS techniques (ESRI 2018), combining each management unit location with national soil map information(Agroconsultores and Geometral 1991a, 1991b) and literature-cited values for each soil type (Constantino and Coutinho 2001; Pimenta 1998a, 1998b). The LS factor, the slope length-steepness factor (dimensionless),was also determined using GIS techniques, following the USPED (Unit Stream Power Erosion and Deposition)method, by extracting the flow accumulation and slope from the landscape DEM layer (Digital Elevation Model)and applying Eq. 2 on the map algebra function of the software.

where: λAis the area of upland flow of the slope; θ is the slope angle, in degrees; m and n are variables related to the rill to interrill erosion ratio and soil susceptibility to erosion, being considered equal to 0.4 and 1.4, respectively, for the present study conditions (Mitasova et al. 1996; Pelton et al. 2012; Kim 2014).

The cover-management factor, C (dimensionless), was estimated as a function of tree cover (Rodrigues et al.2020), considering literature-cited minimum and maximum values of 0.001 and 0.3, respectively, for forest covers in the studied region (Carvalho-Santos et al.2014, 2016; Panagos et al. 2015b). Having no further information on local soil conservation practices, the support practice factor P (dimensionless) was considered equal to 1 for the present study purposes.

Linear programming model

Model variables computations assumed that all silvicultural practices occur at the end of each planning year,biometric, production and erosion values being reported afterwards. For each climate scenario, a total of 26,513 prescriptions - nearly 20 management alternatives per stand, on average - were used to build a resources capability model (RCM) as follows:

where: xijare the decision variables, representing the percentage of management unit i area under prescription j; M is the number of management units (1373); Niis the number of possible prescriptions for each stand i;aiis the area of management unit i; S is the number of forest species (6); SPECIESsis the set of prescriptions j that correspond to species s; erosionijtis the potential soil loss in Mg in period t, resulting from assigning to stand i prescription j; T is the number of planning periods (9); pineijtis the maritime pine volume standing resulting from prescribing j to stand i, at the end of period t; eucalyptijtis the eucalypt standing volume,resulting from prescribing j to stand i, at the end of period t; chestnutijtis the chestnut standing volume,resulting from prescribing j to stand i, at the end of period t; poakijtis the pedunculate oak standing volume,resulting from prescribing j to stand i, at the end of period t; coakijtis the cork oak standing volume, resulting from prescribing j to stand i, at the end of period t;wpineijtis the volume of maritime pine harvested in period t resulting from prescribing j to stand i; weucalyptijtis the eucalypt harvested volume, resulting from prescribing j to stand i, along period t; wchestnutijtis the chestnut harvested volume, resulting from prescribing j to stand i, along period t; poakijtis the pedunculate oak harvested volume, resulting from prescribing j to stand i,along period t; and coakijtis the cork oak harvested volume, resulting from prescribing j to stand i, along period t.

Equation 3 defines that the sum of the areas assigned to all chosen prescriptions for each management unit must be equal to the stand area (ai). Equations 4 to 21 are used to estimate the value of bookkeeping variables and reflect the problem resource capability model (Davis et al. 2001). The area assigned to each forest species(ASPECIESs) are defined by Equation set 4. Soil potential losses for each stand i, under each prescription j and period t, are expressed by accounting variables Erosiont(Eq. 5). Equation sets from 6 to 11 account for the landscape standing volume for each tree species and period t, while Eq. 12 calculates total landscape standing volume at the end of each planning period (VEIt). Similarly,equation sets from 13 to 17 define the volume of each species harvested within each time period, and eq. 18 represents total landscape timber yield for each t(Woodt). Accounting variables defined by eqs. 19 and 20 express the landscape total wood production (Wood)and potential soil losses (Erosion), respectively, from the landscape along the 90-year planning horizon. Inequalities 21 set the model’ non-negativity constraints.

Both climate scenario models were also tested for optimizing soil protection while demanding for a minimum timber supply, defined with a 9 million cubic meters threshold for harvested wood along the 90-year planning horizon (Eq. 22). In Addition, maximum 10% fluctuations between consecutive planning periods’ timber harvesting were included(Eq.23).

All problems were read and solved using CPLEX Interactive Optimizer(IBM 2013).

The model was tested for the two contrasting objectives of timber provision (only without constraints;Eq. 24) and soil protection (Eq. 25).

The RCM was built and solved independently for the two climate scenarios, corresponding to current (BAU)and changing (REF) conditions. Therefore, the RCM coefficients - erosionijt, pineijt, eucalyptijt, chestnutijt, poakijt,coakijt, wpineijt, weucalyptijt, wchestnutijt, wpoakijt, wcoakijt- were calculated separately.

Table 2 summarizes the problems addressed, the respective LP-models and constraints analyzed.

Results

Regarding tree species distribution on the landscape,model solutions for current (BAU) and reference climate change (REF) scenarios were very similar (Fig. 2). Yet, a deeper analysis of the results revealed that the prescriptions assigned to each management unit were not always identical (data not showed). Specifically, in the REF solution, rotation lengths were changed (either extended or shortened) for management units where oak species (i.e.,pedunculated or cork oak) were not an option due to land aptitude limitations.

Results show that, for timber maximization purposes,a bigger proportion of the landscape should be occupied by maritime pine (49%), followed by eucalypt (43%) and residual chestnut and cork oak areas (3.7% and 3.5%, respectively), none of the area being ascribed to pedunculate oak plantations (Fig. 2). When minimum soil loss was set as the objective function, more than a half of the landscape was suggested to be occupied by maritime pine (53%), followed by 29% for cork oak, 16% for pedunculate oak and less than 1% (41 ha) for chestnut, with no area attributed to eucalypt plantations.

To cope with wood provision requirements, the model responded by reallocating 30% to 32% of the landscape area to alternative forest species, mainly by reducing maritime pine, cork and pedunculate oak areas, ascompared to the unconstrained solution (Fig. 2). An approximately even distribution of this reassigned area between chestnut and eucalypt plantations, was designated when a 9 million cubic meter minimum harvested volume was requested. Constraining the model with maximum 10% harvesting fluctuations between time periods,has further increased the area assigned to eucalypt plantations (up to 25%) and reassigned the remaining to pedunculate oak (21% of total area) and, to a less extent,chestnut (ca. 4%). Considering the combined effect of minimum timber production and fluctuations constraints, followed the same area transference pattern,towards a more even distribution of species in the landscape, as follows: about 36%–37% for maritime pine, 29%–30% for eucalypt, 14%–15% for cork oak,11%–13% for pedunculate oak and 7%–8% for chestnut stands.

Table 2 Problem description: linear programming models and constraints tested

Fig.2 Vale do Sousa optimal area distribution of each forest species,for the tested problems under current climate conditions(BAU),at the end of the 90-year planning horizon

Under current climate conditions (BAU), soil erosion from the Vale do Sousa landscape can be managed to reach a minimum average of 26.4 Mg·ha-1·year-1, along the 90-year planning horizon, corresponding to a 7.6-Mg reduction of potential annual soil loss per hectare (ca.22%), compared to the timber production maximization solution (Table 3). Accordingly, changing the management goal - from wood production to soil protection -has increased the volume of ending inventory (VEI) in 0.52×106m3. Adding minimum or even flow timber production constraints led to increases of about 1.5 and 2.4 Mg·ha-1·year-1(6% and 9%), respectively, to the landscape potential soil loss, while resulting VEI was reduced in 0.18 (12%) and 0.48 (31%) million m3, respectively. Considering both timber constrains added 2.9 Mg(11%) to the potential annual soil erosion per hectare,and reduced final standing trees volume in 0.38×106m3(25%), compared to the unconstrained solution.

Increasing annual precipitation under the reference climate change scenario (REF) produced a 25% increase in potential soil loss from the landscape, compared to the current climate conditions solution (Table 3). Yet, a 9.8 Mg·ha-1·year-1(about 23%) soil loss can be prevented and an extra 0.53×106m3standing volume can be attained at the end of the planning horizon, if the landscape management plan is based on soil protection concerns, instead of timber production maximization.Climate change model solutions followed similar trends to those described for current conditions when timber constraints were added. Specifically, increases of 5%, 9%and 11% in potential landscape annual soil erosion, were associated with, respectively, requiring timber production to be above 9×106m3, demanding for an even harvested timber flow, and considering both. Although changes to the solution VEI values were only reduced in less than 3% (40×103m3decrease) when minimum harvested volume constraint was added, when volume fluctuations constraints were considered, alone or together with the former, this accounting variable reduction was proportional to that registered in BAU scenario, respectively, 31% and 25% (0.49×106and 0.36×106m3lower,respectively).

The harvested timber volume constraint (greater or equal than 9 million cubic meters) was binding under both climate scenarios, with dual prices equal to 1318 and 1499 kg, for BAU and REF scenarios, respectively(Table 3). This information highlights that a decrease of a cubic meter in volume harvested may contribute to a decrease of the potential soil loss from the landscape. When both minimum timber production and maximum harvest fluctuations constraints were considered, the amount of soil displacement in the landscape as a result of harvesting the last cubic meter of wood under BAU was 1172 kg, whereas for REF it was 1103 kg.

Table 3 Optimal LP solutions for potential soil losses (Erosion),timber provision(Wood) and average volume of ending inventory(VEI)for the tested problems

Fig.3 Optimal LP solutions’ average annual soil loss,total volume harvested and the volume of ending inventory for each period’final year under current(BAU,left)and changing(REF,right)climate conditions,for the addressed problems and each 10-year period(1 to 9 corresponds to 2017–2026,2027–2036,2037–2046,2047–2056,2057–2066,2067–2076,2077–2086,2087–2096 and 2097–2106)

Potential soil erosion estimates for each 10-year time period reflect the same trends as the global solutions:climate change conditions and harvesting are associated with higher potential for soil displacement (Fig. 3).

When timber production is the only goal, the optimal management solution produces the highest potential soil losses in every planning period, ranging from a minimum of 28.5 or 34 Mg·ha-1·year-1in period four, to a maximum of 40.1 or 50.7 Mg·ha-1·year-1in period two,under BAU and REF respectively (Fig.3).Timber harvest fluctuations were also higher than for other tested problems, BAU varying from 0.75×106to 2.68×106m3, and REF from 0.77×106to 2.78×106m3, for the fourth and ninth periods,respectively.

Conversely, if soil protection is the sole management objective, the lowest values of soil erosion can be achieved for most planning periods, from 17.6 to 33.8 for BAU, and 22.4 to 44 Mg·ha-1·year-1for REF,minimum values corresponding to the last planning period, and maximum values for the second decade in the planning horizon (Fig. 3). Wood provision is more evenly distributed among periods, compared to the maximum wood production solution, ranging 0.17×106– 0.99×106m3in BAU, and 0.17–1.14 million m3in REF, per period. Forest standing volume by the end of each period was also less variable under both climate scenarios, ranging from 0.43×106m3in the sixth period, to 1.55×106and 1.59×106m3by the end of the 90-year planning horizon, for BAU and REF, respectively (Fig. 3).

Combining soil protection with minimum timber production goals determined soil erosion intensification in most study periods (+0.19 for BAU, and+0.21 Mg·ha-1·year-1for REF, on average), except for the 2077–2086 decade (-1.6 and-1.7 Mg·ha-1·year-1, in BAU and REF, respectively), compared to the unconstrained model solutions for both climate scenarios (Fig. 3). Harvested and standing volumes inter-period variability was enhanced by the 9 million cubic meters requirement,when compared to the unconstrained models for both climate scenarios. Specifically, minimum harvested volumes in one period were 0.55×106and 0.52×106m3,for BAU and REF, respectively, corresponding to the 7th period (2077–2086), while maximum values reached 2.21×106and 2.20×106m3, in the last decade of the planning horizon, under BAU and REF, respectively.Similarly, volume standing by the end of each planning decade ranged from a minimum of 0.57 (in period 6) to a maximum of 1.91×106m3in BAU, and from 0.66×106to 1.88×106m3in REF scenario.

The harvested timber inter-period fluctuations constraints conveyed balance to timber provision between periods, ranging from 0.73 to 1.04 and 0.76×106–1.18×106m3for BAU and REF, respectively, and standing forest volume, with minimum 0.54 (2nd period) and maximum 1.34×106m3(4th period), for both BAU and REF, but created additional 2.4 and 3.2 Mg potential soil movement per hectare and year,on average per period, under BAU and REF conditions, respectively, when compared to the unconstrained model solutions (Fig. 3).

When both minimum wood production and fluctuations were considered, potential landscape soil displacement was enlarged, in comparison to the unconstrained models for both BAU and REF scenarios, by an average of 3 and 4.3 Mg·ha-1·year-1, respectively (Fig. 3). Wood provision in each period varied less than observed in the unconstrained model, from 0.76 to 1.18 million cubic meters in BAU, and from 0.77×106to 1.18×106m3in REF, respectively by the end of the first and sixth time period. Forest standing volume by the end of each planning period ranged from minimums 0.68 and 0.65 (second period) to maximums 1.53×106and 1.49×106m3(fourth period),under BAU and REF, respectively.

Less-than-or-equal constraints, defining timber production in one time-period must not exceed that from the previous period in more than 10% (Eq. 23), were binding for periods 2, 4, 5, 6 and 9 under both climate scenarios, with dual prices ranging from 759.1 to 8385.2 kg of soil loss for BAU, and 1091.5 to 13148.8 kg for REF (Table 4). On planning periods 3, 7 and 8, these constraints were not binding,all presenting positive slacks.

Volume fluctuation constraints expressing that harvested timber from one period must always exceed 90%of that in the previous period,were not binding,generating negative surpluses in both climate scenarios’ models,which ranged from 32.01 and 19.66 in period 7 to 215.14 and 214.07 thousand cubic meters in period 6,for BAU and REF conditions, respectively (Table 4).

Analysing the reduced costs of the alternative decision variables from the constrained model under current climate conditions, suggests that most of the management prescriptions options involving eucalypt plantations are associated with potentially higher costs in terms of soil loss, when compared with the other species range of alternatives (Fig. 4). Contrariwise, pedunculate oak alternative prescriptions, are generally related to lower soil erosion costs. It is also noteworthy that the potential costs of choosing some of these alternative prescriptions,regarding soil losses, diverge significantly from the average value for all prescriptions of the corresponding species, particularly for eucalypt and maritime pine. This suggests some of the considered management alternatives may not be compliant with soil protection goals,for what they may be identified, re-evaluated and, if necessary, omitted or replaced in the original model.

Table 4 Slacks or surpluses and dual prices results for harvested volume fluctuation constraints

Discussion

The present study attempted to incorporate forest stakeholders concerns with soil protection from water erosion, into a typical timber production forest plan, at landscape level. By enabling the comparison of a large number of alternatives, including different objectives,resources limitations and possible outcomes, the proposed LP-based optimization approach was effective in analysing the complexity of the Vale do Sousa forest management problem, thus constituting a valuable support to the stakeholders’ decisionmaking process.

Fig.4 Boxplots for the reduced costs of alternative prescriptions attributable to all management units,corresponding to maritime pine(Pb,n=4730),eucalypt(Ec,n=19,219),chestnut(Ct,n=1004),pedunculate oak(Qr,n=681) or cork oak (Sb,n=86)plantations,for current climateconditions(BAU)soil erosion minimization model,including constraints for minimum timber production and maximum inter-period fluctuations

Complying with the restrictions of land units’ suitability, the five alternative forest species were distributed across the landscape according to their performance towards the considered management objective. Specifically,for wood maximization purposes the solutions comprised eucalypt and maritime pine in most of the area,choosing the most productive alternative Fagaceae species (i.e., chestnut and cork oak) only for residual areas,where eucalypt and pine where not suitable. Conversely,concerning soil erosion minimization goals, cork and pedunculate oak plantations were ascribed to all suitable management units, while maritime pine was chosen for stands were none of the alternative broadleaved species were an option. These results are in agreement with a stand-level study for the same area (Rodrigues et al.2020), stating that longer rotations, less intensive silvicultural practices and wider tree crowns may offer better conditions for soil protection from water-related erosion.This same trend is further supported by the analysis of the reduced costs in the final solution, which emphasise that, for the Vale do Sousa landscape biophysical conditions, soil conservation goals are most likely to be met by cork and pedunculate oak or chestnut silvicultural models, than by eucalypt prescriptions, some of which represent considerably high potential soil losses. Accordingly, the timber provision constraints added to the soil conservation optimization problem, have expressed the expected changes to the area attributed to each forest species. When only an overall minimum timber production is considered, rearranging the landscape enabled higher productivity species (mainly eucalypt) to be represented, but when restrictions to production fluctuations are included, the necessary intensification of harvesting events has compelled some recovery on the best soil cover option (i.e., pedunculate oak) for some of the area, in order to balance the overall landscape erosion risk. These results denote not only the model flexibility, but also its information proficiency, since without these kind of optimization tools, it would hardly become clear for forest managers to what extent should one management option be extended, and to which the other should be reduced, to achieve the best feasible management solution.

As expected, a negative relationship between timber provision and soil protection was evidenced by the model solutions, namely, when comparing the two contrasting objective functions - maximize wood or minimize soil loss -, and when analysing changes to the optimal solution for soil protection goals when timber constraints - minimum volume and/or volume fluctuations - were included. Furthermore, the substantial dual price values in every binding constraints for the timber volume harvested, emphasizes the existence of an important trade-off between wood provision and soil protection services, which is in agreement with several other studies reporting this same trade-offs in forest ecosystems (e.g., Baskent and Kücüker 2010; Kele? and Ba?kent 2011; Selkim?ki et al. 2020).

The constraints aiming for a constant timber flow,have further pushed the model optimal solution to higher erosion values (10%–11% increases), if compared to the minimum volume constraint (only 5%–6% increases), reflecting the negative effect of increasing harvesting events frequency, over the landscape potential for soil losses. In fact, logging is among the major disturbance causes of soil erosion, along with wildfires(Borrelli et al. 2016). However, by concentrating harvesting events in fewer planning periods, the minimum timber volume constraint led the model solution to estimate higher rates of potential soil loss in specific periods, particularly in those following harvesting peaks (periods 2 and 6). This observation is crucial in the light of the inherent uncertainty of forest management planning.Therefore, although the solution from the model considering only the minimum 9 million cubic meters constraint may seem appealing, regarding soil conservation purposes for the 90-year planning horizon, the fact that this same solution enhances the risk of extreme soil erosion events in specific periods, may led decision-makers to reconsider. Indeed, given the unpredictability of extreme rainfall events, particularly when associated with climate change scenarios (Panagos et al. 2015a),the need to maintain soil as uniformly covered as possible, avoiding extensive bare soil areas, becomes evident (Borrelli et al. 2016).

The highest harvest-related soil loss risks were associated to periods 4 and 5 for both climate scenarios, as shown by the considerably higher dual prices on harvested volume fluctuation constraints Wood4≤1.1×Wood3and Wood5≤1.1 ×Wood4(up to 8.4 and 13.1 Mg soil per m3wood, for BAU and REF respectively). These results are related to the lack of a no-harvest alternative in the model, which is making it struggle with accommodating Fagaceae species and maritime pine rotations(35 to 60 years) to a more even flow of wood. In fact, although the model can choose some of the Fagaceae and maritime pine stands to be harvested on periods 3 (between years 31 and 40, or 2037–2046), all the remaining area must be harvested on periods 4 or 5 (years 41 to 60, 2047–2066). Remarkably, it was not difficult to target higher wood flows, which was achieved by managing more area for eucalypt, but it was a challenge to prevent it from increasing during periods of concentrated harvest. The LP-based optimization approach confirmed the identification of critical time points in which forest management decisions can be decisive for the achievement of targeted goals, including the limitations derived from the standard silviculture regimes.

The proposed RCM design was broad enough to accommodate the required wood production levels, both in quantity and even distribution, while maintaining potential soil loss rates below 35 Mg·ha-1·year-1for all planning periods.

Climate change is likely to further challenge forest management planning, in which soil conservation issues may become critical (Cosofret and Bouriaud 2019). Our case study results show that, under RCP8.5-compatible local climate change conditions, forest growth enhancement following temperature rise, may not fully compensate the predicted rainfall erosivity increase. Therefore,soil erosion potential may increase to values as high as 46 Mg·ha-1·year-1in some planning periods, reinforcing the susceptibility of the study region to the predicted climate change conditions (Rodrigues et al. 2020), and urging the discussion on adaptive and mitigation strategies, among which silvicultural practices are crucial. It is noteworthy that the presented LP-based optimization approach was able to identify and incorporate heavier precipitation years, in the event of climate change, namely by adjusting the solution with different rotation lengths to ensure tree cover when erosion risk might be critical, reinforcing the usefulness of optimization tools in unravelling, solving and interpreting complex forest landscape management problems.Notwithstanding, by considering annual precipitation values to represent rainfall erosivity (R factor), we may be imposing a false level of accuracy to the model’ solutions, given the forecast nature of the available local climatic data. This possible limitation can be improved by a sensitivity analysis to variations in precipitation patterns, including the use of precipitation periodic values.

The relatively high values of potential soil displacement in the landscape must be seen under the scope of the limitations associated to the RCM building, so the actual soil may be considerably lower. Specifically, the present research main limitations can be listed as follows: there were no process-based models available (to our knowledge) to simulate the growth of the studied forest species in the study region; available growth and yield models are not sensitive to soil disturbances, such as soil erosion; the value of ending inventory was not considered and this may lead to the extension of second or third rotation ages so that the corresponding clearcut occurs after the end of the planning horizon; understorey vegetation cover was not accounted in erosion computation; despite the acknowledged seasonal variability of rainfall erosivity (Zema et al. 2020), our approach was based on the available annual precipitation data for the region. In addition, since all stands are even aged, most management prescriptions (except for cork and pedunculate oaks) include at least a clearcut within the planning horizon. Thus, as a consequence of the current inventory and the limited number of rotation ages alternatives, a certain degree of inter-period variability, both in wood extraction and soil loss, becomes inevitable. Results suggest that increasing management flexibility by considering additional rotation ages may contribute to decrease that variability. Nevertheless,this will increase computational costs.

This approach intended to depict how assessing a large number of forest management alternatives, at a landscape level,and finding the best possible spatially explicit solution for a 90-year planning horizon, regarding soil conservation and timber provision goals may be accomplished by LP-based optimization. It may be further improved if other relevant ES are added to the RCM (e.g.,vulnerability to wildfires, biodiversity conservation, carbon sequestration goals, cultural services), so that their specific trade-offs or synergies regarding soil protection can be analysed. The RCM in the linear program may be also extended to take into account the diversity of wood products as well as the corresponding monetary returns(e.g., net present value). Moreover, the policy model in the linear program, e.g., the set of product flow constraints as well as the objective function, may be adjusted to the ecological and socioeconomic context of the management planning problem.

Conclusions

The present research demonstrated that a LP-based optimization approach can describe the influence of forest management options on soil protection ecosystem services. The described approach enabled the analysis of trade-offs between the conflicting goals of soil conservation and timber provision.Results for the case study area indicated that it is possible to manage the Vale do Sousa forest landscape for soil conservation, without disregarding wood demand. The usefulness of optimization tools for forest landscape management, concerning the provision of biophysical conditions-dependent ES, is therefore emphasized. By providing a spatialized management solution, information on possible trade-offs or synergies between goals, as well as an overview on the impact of changes and uncertainties, this methods are undoubtedly valuable for the development of broader management planning systems to support forest managers’decision process.

Acknowledgements

Authors thank Associa??o Florestal do Vale do Sousa for assistance with study area characterization and forest inventory tasks.

Authors’contributions

All authors contributed to the present manuscript preparation. ARR developed the methodology, analysed and interpreted the data and wrote the paper.SM was involved in conceptualization, methodology development, results analysis and interpretation, manuscript writing and reviewing. BB contributed to conceptualization, data collection,methodologies development,manuscript writing and reviewing. MM assisted in data collection and computation,methodology development and implementation and manuscript reviewing. JGB was involved as research coordinator, project manager and in conceptualization, development of methodologies and manuscript writing and final reviewing.All authors read and approved the final manuscript.

Funding

This research was funded by BIOECOSYS project “Forest ecosystem management decision-making methods: an integrated bioeconomic approach to sustainability” (LISBOA-01-0145-FEDER-030391, PTDC/ASP-SIL/30391/2017),MODFIRE project “A multiple criteria approach to integrate wildfire behaviour in forest management planning with the reference” (PCIF/MOS/0217/2017),and NOBEL project “Novel business models to sustainably supply forest ecosystem services”, under the umbrella of ERA-NET cofund ForestValue (Horizon 2020 research and innovation programme grant agreement N°773324).Centro de Estudos Florestais, research unit is funded by Funda??o para a Ciência e a Tecnologia I.P. (FCT), Portugal within UID/AGR/00239/2020.

Availability of data and materials

The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

Declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Received: 3 November 2020 Accepted: 10 June 2021