Addressing soil protection concerns in forest ecosystem management under climate change

Climate change may strongly influence soil erosion risk, namely through variations in the precipitation pattern. Forests may contribute to mitigate the impacts of climate change on soil erosion and forest managers are thus challenged by the need to define strategies that may protect the soil while addressing the demand for other ecosystem services. Our emphasis is on the development of an approach to assess the impact of silvicultural practices and forest management models on soil erosion risks under climate change. Specifically, we consider the annual variation of the cover-management factor (C) in the Revised Universal Soil Loss Equation over a range of alternative forest management models to estimate the corresponding annual soil losses, under both current and changing climate conditions. We report and discuss results of an application of this approach to a forest area in Northwestern Portugal where erosion control is the most relevant water-related ecosystem service. Local climate change scenarios will contribute to water erosion processes, mostly by rainfall erosivity increase. Different forest management models provide varying levels of soil protection by trees, resulting in distinct soil loss potential. Results confirm the suitability of the proposed approach to address soil erosion concerns in forest management planning. This approach may help foresters assess management models and the corresponding silvicultural practices according to the water-related services they provide.


Background
Soil erosion, that is, the process that transforms soil into sediments, is one of the major and most widely spread forms of land degradation (Lal 2014;Weil and Brady 2017). It encompasses the destruction of the physical structure that supports the development of plant roots. Moreover, surface soil removal may result in substantial nutrient and water losses, as well as in the decrease of productivity and the increase of pollution of surface waterways. Soil erosion impacts thus the sustainability of ecosystems and the provision of ecosystem services. Soil conservation efforts address concerns with these impacts and meet the increasing needs for food and raw materials (Hurni et al. 2008).
Soil erosion by water is linked to desertification processes. Its severity is prone to increase as a consequence of changes in the amount of precipitation as well as in its temporal and spatial distribution under prospective climate scenarios (IPCC 2014a). This will exert further pressure on ecosystems water balance and calls thus for adequate soil protection and conservation practices in the framework of ecosystems management (Coutinho and Antunes 2006;Jones et al. 2011;Panagos et al. 2015bPanagos et al. , 2015cAnaya-Romero et al. 2016;Seidl et al. 2016).
Trees are widely known to impact the ecosystem hydrological cycle and resultant water availability and quality (Brown and Binkley 1994;Marc and Robinson 2007;Keenan and Van Dijk 2010;Carvalho-Santos et al. 2014). As vegetation cover plays a crucial role in erosion and runoff rates, afforestation is considered among the best options for soil conservation (Durán Zuazo and Rodríguez Pleguezuelo 2008;Lu et al. 2004;Gyssels et al. 2005;Panagos et al. 2015b;Ganasri and Ramesh 2016). Water-related forest ecosystem services include the provision, filtration and regulation of water, along with stream ecosystem support and water-related hazards control, e.g., soil protection from erosion and runoff (Bredemeier 2011). In this context, forest management practices that involve vegetation cover modifications may have a substantial impact on the provision of waterrelated ecosystem services (Ellison et al. 2012;Panagos et al. 2015b). Moreover, forest ecosystems interactions with water and energy cycles have been highlighted as the foundations for carbon storage, water resources distribution and terrestrial temperature balancing. Forest management may thus play a key role to meet climate change mitigation goals (Ellison et al. 2017).
The cover-management factor (C-factor) within the Revised Universal Soil Loss Equation (RUSLE) is used as an indicator of soil protection by different landuses and management options (Renard et al. 1991). Yet, few studies have addressed its potential as a dynamic tool for erosion control (Panagos et al. 2015b). Experimentally determined values for the C-factor for most land uses and management systems are easily found in the literature (e.g., Pimenta 1998a). Moreover, both remote sensing and geographical information systems (GIS) techniques can be efficiently used to estimate the C-factor at landscape level (Wang et al. 2003;Lu et al. 2004;Durigon et al. 2014). Nevertheless, the literature does not report the use of the C-factor to address impacts of vegetation density changes over time under the same land use or management type. This provided the motivation for this research.
We aim at assessing the impacts of forest ecosystem management practices (e.g., selection of tree species, harvesting) on soil protection, as its planning schedule impacts soil erosion over the long-term (Lu et al. 2004;Panagos et al. 2014Panagos et al. , 2015b. Our research examines how management practices contribute to change the vegetation cover over time. It further encapsulates these changes within the RUSLE, by determining the corresponding C-factor. Seven stand-level forest management models (sFMM), i.e., sequences of management practices, with species-specific rotations, over a 90-year time span, are used for testing purposes. Specifically, we assess and compare sFMM according to their potential for the provision of water-related ecosystem services under two climate scenarios.

Study area and climate scenarios
For testing purpose, we considered a 14,837-ha study area located in the northwestern region of mainland Portugal (Lat: 41.1343, Lon: − 8.2951; Fig. 1). This area was selected according to its representativeness of the country's forest management context (e.g., species composition, ownership structure; Novais and Canadas 2010;ICNF 2013).
We considered further two local climate scenarios over a temporal period extending from 2017 to 2106. The first scenario assumes that current conditions will stay the same (mean annual temperature equal to 13.8°C and annual precipitation of 1194 mm; local climate normals from 1981 to 2010). As an example of predicted global climate changes, a second scenario was considered, consisting of a local adaptation of the RCP8.5 global climate change scenario (IPCC 2014b). It encompasses increases of 2.36°C and 193 mm, in mean annual temperature and annual precipitation, respectively. The values of the annual climate variables in the region over the 90-year planning horizon were estimated using the CliPick online tool (Palma 2015(Palma , 2017; http://www.isa.ulisboa.pt/ proj/clipick/) and its KNMI-RACMO22E models (van Meijgaard et al. 2012).

Stand-level forest management models
Four currently used stand-level forest management models (sFMM) were identified as prevalent in the study area. A mixture between maritime pine (Pinus pinaster Ait.) and eucalyptus (Eucalyptus spp.) characterizes the first two: sFMM1, where maritime pine is dominant (73% coverage) with fewer eucalyptus (27%); and sFMM2, where eucalyptus is the predominant species, with less maritime pine coverage (67% and 33%, respectively). sFMM3 consists of pure chestnut stands (Castanea sativa Mill.) and sFMM4 are pure eucalyptus plantations.
Recent concerns with wildfire risk and the provision of several ecosystem services triggered the development of three alternative stand-level forest management models in a participatory approach involving local stakeholders (e.g., forest owners, pulp and paper industry agents, municipalities and forest authorities) (see Marques et al. 2020). A fifth model (sFMM5) was proposed that involves the plantation of pure maritime pine stands with lower densities than the currently used (ca. 2200 trees·ha − 1 ), while a sixth (sFMM6) and a seventh (sFMM7) model proposed the plantation of pure pedunculated oak stands (Quercus robur L.) and pure cork oak stands (Quercus suber L.), respectively. The full list of sFMMs (Table 1) was developed in consultation with the stakeholders. This participatory process highlighted their concern with the provision of water-related ecosystem and triggered the development of an approach to compare the potential for soil erosion regulation by each sFMM.
The values of each species biometric variables, for each sFMM, were obtained using the StandSIM-MD simulation module (Barreiro et al. 2016). The simulation considered a 90-year temporal horizon. It considered further three site indexes (SI), representing higher, lower and intermediate (hereafter referred as site index 5, 1 and 3, respectively) site quality conditions in the study area for each tree species. Moreover, the impact of climate change on the trees growth was assessed according to results reported by Santos and Miranda (2006) for the region, by adjusting linearly the simulated standing volume (and respective biometric variables) to the expected temperature increase along the 90-year planning horizon. Specifically, the cork oak yield is expected to increase by 3.52%, while timber yields in the case of other forest species are expected to increase by 4.2%, over the 90-years temporal horizon.

Potential soil erosion estimates by RUSLE
To provide an approximate and preliminary ranking of both current and alternative sFMM water services provision, the Revised Universal Soil Loss Equation (RUSLE) was applied to estimate the potential annual soil loss by surface runoff under each sFMM. The RUSLE method is the most widely used long-term soil erosion prediction tool (Renard et al. 1991;Guo et al. 2019) and can be expressed as (Eq. 1): where A is the estimated soil loss (Mg·ha − 1 ·yr − 1 ); R is the rainfall-runoff erosivity factor (MJ·mm·h − 1 ·ha − 1 · yr − 1 ); K is the soil erodibility factor (Mg·h·MJ − 1 ·mm − 1 yr − 1 ); L is the slope length factor; S is the slope steepness factor; C is a cover-management factor; and P is a conservation practice factor.

K, LS and P factors
The soil erodibility factor (K) value was estimated as 0.1265 Mg·h·MJ − 1 ·mm − 1 ·year − 1 using GIS techniques (ArcGIS, version 10.6) to merge national soil map information Geometral 1991a, 1991b) and reference K values for the study area available in the literature (Pimenta 1998a(Pimenta , 1998bConstantino and Coutinho 2001). Soils in the study area are mostly Umbric Leptosols and Leptic Regosols, developed over The topographic factor (LS) was calculated considering a reference area of 100 m × 100 m, with a 2% slope. Using L (Eq. 2) and S (Eq. 3) equations by Renard et al. (1991), a LS value of 0.45 was estimated.
where λ is the length horizontal projection of the slope; m is a variable exponent related to the rill to interrill erosion ratio, considered 0.4 for our study conditions (Pelton et al. 2012;Kim 2014); θ is the slope angle; and s is the slope in percentage. Since no specific information was available regarding the implementation of erosion control support practices in the study area, the conservation practice factor (P) was considered as equal to 1 when estimating the impact on erosion by all sFMM (Panagos et al. 2015c).

Erosivity factor (R)
The mean annual precipitation (i.e., precipitation normal, given by the previous 30-year average) is lower than 1425 mm in the case of the two local climate scenarios. Thus, the erosivity factor in each year was estimated by the linear relation between annual precipitation (P) and erosivity (R) as proposed by Constantino and Coutinho (2001) for the Portuguese northern region (Eq. 4).
Cover-management factor (C) The cover-management factor associated to specific land use conditions can be determined experimentally (Renard et al. 1991). However, experimental methods are generally high demanding in both time and resources. The literature reports approximate values that can be adapted to estimate soil erosion. For example, Panagos et al. (2015b) suggested the use of an equation (Eq. 5) to encompass the great variability of literaturecited values for similar non-arable land uses.
where C land use stands for the values of land use covermanagement factors reported in the literature; and F cover is the proportion of soil covered by vegetation, varying between 0 (no vegetation cover) and 1 (soil fully covered by vegetation). Based on an extensive literature research on similar forest species and conditions (Pimenta 1998a;Märker et al. 2008;Jones et al. 2011;Vieira et al. 2014;Carvalho-Santos et al. 2014Panagos et al. 2015b;Fernández and Vega 2016), we have selected 0.001 and 0.3 as minimum and maximum values, respectively, of the land use cover-management factor. We have estimated the proportion of covered soil (F cover ) using the tree crown width (CW). Specifically, this proportion in each year of the temporal horizon was estimated taking into account the stand density and considering that each tree would cover a circular area defined by its width (Shaw 2003;Pretzsch 2009). CW was determined for each tree species as suggested by Condés and Sterba (2005) where dbh and h are the diameter at breast height and the height, respectively, in each year of the temporal horizon; and a 0 , a 1 and a 2 are the species-specific coefficients by Condés and Sterba (2005) (see Table 2).

Rain erosivity under climate change
The value of the R factor ranged from 351 to 1900 MJ·mm·ha − 1 ·h − 1 and from 585 to 2550 MJ·mm·ha − 1 ·h − 1 in the case of, respectively, current climate and the RCP8.5-compatible conditions (climate change scenario).
In the case of the former the average value is 1189 MJ·mm·ha − 1 ·h − 1 , while in the latter it is 1498 MJ·mm·ha − 1 ·h − 1 . The expected precipitation increase under RCP8.5 will lead to an overall rain erosivity factor increase of about 309 MJ·mm·ha − 1 ·h − 1 , on average (Fig. 2).

A dynamic C-factor
The value of the cover-management factor reflected the impacts of both tree growth and silvicultural practices on the proportion of soil covered by vegetation under each sFMM over the 90-year temporal horizon (Fig. 3).
Harvests (e.g. thinning, clearcut) produced high C-factor values in the case of all sFMMs. Values were generally lower for higher site indexes within the same sFMM. Climate change (RCP8.5 scenario) had little impact over the cover-management factor values. All Fagaceae species, namely chestnut, pedunculated oak and cork oak, have reached total soil surface coverage, where the minimum C-factor value (0.001) was computed. However, the time needed to reach such soil cover conditions, as well as the ability to keep it for long time periods, was notoriously dependent on both, the site index and the silviculture model (Fig. 3).
The C-factor reaches its maximum values in all sFMM in years when the stands are clearcuted. Thinning or partial harvest operations (in mixed stands, sFMM1 and sFMM2) had comparatively lower impacts over the Cfactor. Consequently, mixed eucalypt and maritime pine stands (sFMM1 and sFMM2) showed lower average Cfactors (0.18 to 0.26, respectively), when compared with the same species pure stands (sFMM4 and sFMM5; 0.23 to 0.26, respectively). Moreover, eucalypt-based sFMM (1 and 4) showed lower average C-factor value when compared with maritime pine-based ones (sFMM2 and 5; Fig. 4).
The highest average C-factor, within the 90-year temporal horizon was determined for pure maritime pine stands (sFMM5, 0.16 minimum), while the lowest was estimated for pedunculated oak system (sFMM6; Fig. 4).

Soil erosion estimates
Climate change leads to a 24% to 46% increase in annual soil erosion potential over the 90-year temporal horizon (Fig. 5). Average annual soil losses were consistently higher in the case of lower site indexes. The impact of the site index is greater in the case of sFMM3.
In the case of current sFMM, pure eucalypt stands (sFMM4) and pure chestnut stands (sFMM3) are associated to the highest and the lowest average potential soil losses, respectively, over the 90-years temporal horizon (Fig. 5). Nevertheless, the average potential soil erosion of the pure maritime pine alternative (sFMM5) is the highest among all sFMM, with predicted annual soil losses being 3% to 420% higher than in the case of any other sFMM.
Erosion control goals seem to be best met by Quercus robur L. stands (sFMM6). The average potential soil loss over the 90-year period is reduced by 35% to 86%, when compared to other sFMMs.
The introduction of cork oak (sFMM7) may also provide considerably better soil protection than eucalypt and maritime pine models, both pure and mixed. When compared with the current Fagaceae species (chestnut, sFMM3), sFMM7 annual soil erosion potential was slightly lower in the case of site index 1, but considerably higher in the case of SI 3 and 5.

Discussion
Forest managers are challenged by the need to safeguard ecological values while extracting tangible economic products (e.g. timber) from forest ecosystems. The effectiveness of forest management planning depends on the provision of information about the trade-offs between timber and other ecosystem services (Keleş and Başkent 2011;McKay 2011;Seidl et al. 2016). This research aimed at clarifying how forest stands management impacts soil physical protection. This is influential to evaluate its potential to provide water-related ecosystem services such as soil erosion risk reduction.

Rain erosivity increase under climate change
The expected increments in mean annual temperature and annual rainfall, under climate change scenario RCP8.5, in the north-western Portugal case study area, will certainly impact in the provision of water-related ecosystems services by forests and aggravate soil erosion risks. Our results show that such effect can be mainly attributable to the overall increase of rain erosivity (R), while the increase of soil coverage due to higher productivities under RCP8.5 do not offset the impact of R. The values of R are in line with data from the study area, reported by other authors (e.g., Brandão et al. 2006;Ferreira 2013;Meneses 2014;Panagos et al. 2015a). Over the 90-years temporal horizon, average erosivity values under current conditions were above 1000 MJ·mm·h − 1 ha − 1 , confirming these areas' high potential for soil losses by water erosion and run-off. Moreover, under the local climate change scenario, erosivity may increase, on average, by more than 300 MJ·mm·h − 1 ·ha − 1 year − 1 , while in exceptionally high precipitation years R values under RCP8.5 conditions may exceed those under current conditions by more than 1600 units. These contrasts highlight the importance of considering interannual climate variability when estimating long-term soil erosion. They further highlight the complexity of erosion processes, which must not be overlooked when modelling, interpreting and up-scaling results.

Managing C-factor for soil erosion control
Our results have demonstrated that management can strongly influence soil coverage by trees and may thus have a substantial impact on potential erosion, within the same forest land-use system. Soil erosion potential increased with harvesting and thinning operations, and decreased with tree growth, as a consequence of the cover-management factor (C) variability. Furthermore, clearcuts had higher impacts than thinning or partial harvests (in mixed stands). This suggests that the proposed dynamic C-factor method is effective, reproducible and may be easily applied in the framework of practicable long-term soil erosion evaluation procedures such as RUSLE.
In this research, we assumed that the minimum and the maximum C-factor values did not vary across forest species. Thus, the estimated soil protection potential associated with each sFMM depends only on the corresponding species crown widening pattern. This research does not consider the influence of the species-specific leaf area (LAI) or their deciduous characteristics. Results should thus be interpreted with caution as we did not address the canopies capacity to intercept, store and break the mechanical energy of rainfall (Coutinho and Antunes 2006). However, our results do reflect the general differences between canopy structure and geometry (Keenan and Van Dijk 2010). For example, sFMM2, where a conifer species is dominant (maritime pine), had the highest soil erosion estimates, compared to other current systems, all of which included broadleaved species (pure and mixed eucalypt, or pure chestnut). The same trend was observed for the alternative sFMM5, consisting of pure maritime pine with reduced tree Species-specific canopy characteristics and growth patterns have also clearly distinguished Fagaceae-based sFMMs (3, 6 and 7) from those comprising only eucalypt and/or maritime pine. The simulation of the development of chestnut, pedunculate and cork oak canopies highlighted these species efficiency in achieving total soil coverage (i.e., minimum C-factor values) at some point within the temporal horizon, thus resulting in the lowest risks of soil erosion. However, the impacts of these sFMM on soil erosion over the 90-year time Fig. 3 Cover-management factor (C) in the period from 2017 to 2106 for each sFMM and Site Index (SI) under current and RCP8.5-compatible climate conditions. sFMM1 mixed eucalypt and maritime pine; sFMM2 mixed maritime pine and eucalypt; sFMM3 chestnut; sFMM4 eucalypt; sFMM5 maritime pine; sFMM6 pedunculated oak; sFMM7 cork oak span, depend also on the corresponding rotation lengths and thinning regimes, as these have determined the duration of the undisturbed periods where the C-factor was low. Moreover, as pointed out, species-specific leaf area (LAI) and deciduous characteristics may also impact the contrasts between sFMM.
Fire hazard concerns by the stakeholders led to the proposed reduction of tree density in the case of the pure maritime pine alternative sFMM5. Yet, results have shown that the associated lower soil coverage can increase soil erosion risk substantially. However, as wildfire hazards pose pressure on forest management, the strict correlation between fire risk and biomass load highlights the importance of taking into account the trade-offs between soil cover and forest resistance to wildfire (Garcia-Gonzalo et al. 2012;Botequim et al. 2019). Nevertheless, these trade-offs should be analysed at the landscape-level (e.g. Borges et al. 2014Borges et al. , 2017Marques et al. 2017), in order to provide information to support the spatial distribution of sFMM over the landscape mosaic.
Mixing eucalypt and maritime pine seemed a better alternative for soil protection than investing in pure forest Fig. 4 Boxplots for the C-factor values grouped by sFMM and Site Index (SI) along the temporal horizon (2017-2106), under current climate conditions. sFMM1 mixed eucalypt and maritime pine; sFMM2 mixed maritime pine and eucalypt; sFMM3 chestnut; sFMM4 eucalypt; sFMM5 maritime pine; sFMM6 pedunculated oak; sFMM7 cork oak horizon, under the current (solid colors) and RCP8.5-compatible (white) climate conditions. sFMM1 mixed eucalypt and maritime pine; sFMM2 mixed maritime pine and eucalypt; sFMM3 chestnut; sFMM4 eucalypt; sFMM5 maritime pine; sFMM6 pedunculated oak; ssFMM7 cork oak stands, as it enabled the reduction of potential soil losses by 0.3 to 1.3 Mg·ha − 1 ·year − 1 . Yet, mixed species sFMM require more frequent harvests, and that must be considered in the light of rain events variability and uncertainty. The combination of heavy rainfall with bare soil intensifies potential soil losses, so minimizing such erosion-prone conditions is usually advocated (Kort et al. 1998;McKay 2011).

Applicability, limitations and future studies
This research considered a constant slope and a onehectare squared area, in order to assess and compare the impacts of sFMM on soil erosion potential. In the case study area, slopes range from 1.3% to 82%, the average being close to 20%, so the impacts of the sFMM on soil erosion potential may vary substantially between management units. In this research, we did not consider the use of process-based models to project tree growth under climate change as they are not available for most forest species in the region. Nor did we consider the impact of extreme climate events (e.g. rain events) on forest growth and productivity and soil erosion. Nevertheless, our approach may easily be implemented to check the impacts of sFMM across the full range of slopes. Moreover, it may also be implemented taking advantage of tree growth models sensitive to climate change when they become available. This will be influential to develop effective forest-level approaches to allocate sFMM over the landscape.
Our approach may help forest managers address erosion concerns when developing their plans. It provides information needed to assess trade-offs between soil erosion and other ecosystem services when checking alternative spatial distributions of sFMM over the landscape mosaic. For example, it may be used to support the development of fuel management strategies to address concerns with wildfire and erosion risk. Future research will address the analysis of trade-offs between wildfire risk and soil erosion management criteria at the landscape-level spatial scale. Also, further research may monitor the implementation of the sFMM in these plans in order to enhance soil erosion modelling, its accuracy and applicability (Ferreira 2013).
Finally, it must not be overlooked that this research did not address the influence of the understorey vegetation layer or of the LAI and deciduous features. Spontaneous herbaceous and shrub vegetation can develop relatively fast after harvests, thus protecting the soil in the early stages of forest stand development and effectively reducing potential soil losses (Pimenta 1998a;Durán Zuazo and Rodríguez Pleguezuelo 2008). Future studies should focus on depicting the combined impact of trees and understorey vegetation on the provision, of forest water-related ecosystems services.

Conclusions
Our research demonstrated that a dynamic C-factor approach can be useful to improve estimates of long-term potential soil losses by water erosion, under alternative stand-level forest management models. Broadleaf species and longer rotation systems contribute to increase soil coverage and decrease erosion risks. This information may be used to support the analysis of trade-offs between water related ecosystem services and other forest products and services. Results under local climate change conditions have highlighted the case study area susceptibility to rain erosivity, stressing the importance of considering soil protection services when scheduling stand-level management options. In summary, this research developed an approach that may be used by forest managers in order to address soil protection concerns when developing forest management plans.