Climate Change Impacts the Protective Effect of Forests: A Case Study in Switzerland

Introduction

Forests provide a series of ecosystem services to human beings, such as water filtration, biodiversity, wood production, carbon sequestration, and protection from natural hazards (Costanza et al., 1998; Pagiola et al., 2002). In steep, mountainous terrain, the protective effect of forests against natural hazards is of particular importance (Sakals et al., 2006; Grêt-Regamey et al., 2008). Thanks to this so-called “nature-based solution,” costs for technical structures for the protection of settlements and infrastructure can be reduced or even avoided (Getzner et al., 2017).

Various studies have shown the high effectiveness of forests in reducing the negative effects of natural hazards (see for review, Moos et al., 2018a; Sebald et al., 2019; Ruangpan et al., 2020). These, however, mostly focused on a specific (static) forest state and did not take the temporal evolution of forests into account. The evaluation of the long-term efficacy of protection forests remains complicated since uncertainties in the protective capacity are high (Woltjer et al., 2008). Natural forest dynamics and disturbances can result in temporal variation of the protective effect of forests (Silins et al., 2009; Vacchiano et al., 2015). Furthermore, climate change is expected to induce shifts in the forest structure and extent (e.g., treeline) and in the distribution of their constitutive species (Kulakowski et al., 2010; Lindner et al., 2014; Albrich et al., 2020; Scherrer et al., 2020). This may, on the one hand, increase the protective effect of temperature limited mountain forests and treeline ecotones, where higher temperatures are expected to increase the forest extent (e.g., along altitude) and growth (Bebi et al., 2017; Scherrer et al., 2020). On the other hand, more frequent drought periods at lower elevations lead to a shift of the current species composition to more drought adapted species (Peñuelas and Boada, 2003; Urli et al., 2013). This shift ideally happens gradually. However, abrupt diebacks of contemporary species, followed by a slow reforestation with better adapted species can also be expected in the future (Elkin et al., 2013). This may result in temporary or even irreversible losses in the protective effect of forests due to a reduced stem density and diameters (Bebi et al., 2017; Irauschek et al., 2017).

Dynamic forest landscape models (Schumacher et al., 2004; Lischke et al., 2006; Seidl et al., 2012) enable to both understand and predict (sensu Shmueli, 2010) forest development over larger areas under past or future climate scenarios. They are a promising way to investigate the evolution of the forest structure and composition at broad temporal and spatial scale (Petter et al., 2020). So far, few attempts have been made to assess the future development of the protective effect of forests. Furthermore, previous studies have mostly used simple indicators to estimate the protection efficacy of the forest without calculating the explicit risk reduction provided by the forest (Elkin et al., 2013; Albrich et al., 2020). In contrast, a risk-based approach allows for the translation of the physical effect of trees on the natural hazard process into monetary terms (Teich and Bebi, 2009; Moos et al., 2018b). This can then serve as input for economic evaluation techniques, such as cost-benefit analyses (Olschewski et al., 2012). An economic evaluation, however, can only be realistic if it accounts for the long-term evolution of the protective effect, which becomes increasingly challenging in the face of climate change (Dale et al., 2001; Allen et al., 2010; Elkin et al., 2013).

In this study, we used the dynamic forest model TreeMig (Lischke et al., 2006) to project the future forest development of a regularly managed rockfall protection forest complex in the Swiss Alps under two different climate change scenarios. We first compared the model projections to the forest cover under current climatic conditions, including the regeneration in areas where forest interventions have taken place in the past 30 years. We then assessed how the protective effect of the forest against rockfall may change under climate change in the next 150 years based on a quantitative risk analysis. The combined analysis of current data on forest regeneration and future forest simulation then provides a basis to evaluate on the role of management in protection forests under climate change.

Materials and Methods

Study Area

The study area (Figure 1) is located in the east-west oriented main valley of the canton of Valais (Switzerland) in the south of the town Martigny (46°05′32″N/7°04′09″E). It ranges from 500 to 1,000 m a.s.l. and the slope has an inclination between 0° and 70° (0° and 55° in the forested area). The climate of Martigny is dry (one of the driest regions in Switzerland) with an average yearly precipitation of 880 mm and a maximum mean temperature of 14°C (June) and –2°C (January) and minimum mean temperature of 3°C (June) and –10°C (January). Due to already low precipitation, climate change is expected to distinctly influence the forests in this region, and, consequently, also their protective effect against natural hazards. Since the forest above Martigny plays an important role in protecting subjacent infrastructure (road/railway) and houses from rockfall, it is particularly suitable to analyze climate change effects on the protection service. The rockfall area consists of small, dispersed rock cliffs with a total area of ∼2 ha and consists mainly of gneiss and mica slate. The forest is mainly dominated by Fagus sylvatica, Picea abies, Pinus sylvestris, and Abies alba and covers a total area of ∼77 ha. In the past 30 years, regular interventions have been made across forest with cable-based harvesting, where relatively large openings for regeneration were created (information from local forest manager and orthophotos).

Figure 1. Study area with sample plots in the intervention areas (orange squares), sample plots in the rockfall zone (green dots), and the forest area modelled in TreeMig (dark green).

Forest Scenarios

We modelled current and future forest development with the dynamic, distribution based, spatially explicit forest landscape model TreeMig (Lischke et al., 2006; Lischke, 2020). TreeMig simulates reproduction, tree establishment, growth, competition, and mortality of trees per height class and cells of 25 m × 25 m; the state variables are the numbers of individuals per height class, grid cell and species. These processes are controlled by the yearly bioclimatic drivers day-degree sum (sum of mean daily temperatures above 5.5°C), minimum winter temperature (mean temperature of the coldest month), and a drought stress index (based on monthly precipitation, temperature, and on soil water holding capacity). Because of the yearly resolution, also extreme events of drought or heat are taken into account. Trees of different cells interact spatially through seed dispersal. The tree density and light distribution in a cell is represented by a Poisson distribution, assuming trees to be randomly distributed. Currently, 30 European species are parameterized in TreeMig (Supplementary Tables 1, 2). Because only rarely sufficient area-covering forest data for all height classes and species are available and initializing point data (e.g., NFI data) transformed to such format tends to lead to artefacts in the first simulated years, TreeMig is usually initialized by running a spin-up under different disturbance regimes (also representing simplified management) and then choosing the final output of the scenario which fits best to current forest data. TreeMig is embedded in an R-wrapper which helps to prepare the input files; the model and support thereon can be obtained upon request from lischke@wsl.ch, and in future on GitHub.

Simulations were run in different steps: (a) spin-up under different disturbance scenarios, (b) choice of simulated current forest at the end of the spin-up period by comparison to the real forest, and (c) future simulations from the simulated current forest under different climate change scenarios. (a) We simulated forest development in a forest mask based on the existing forest cover derived from land cover data (swisstopo, 2020b). The current forest state was modelled based on a simulation from bare ground for a time span of 400 years (“spin-up” simulation). Climate data (Figure 2) were taken from measurements from 1930 to 2016 with randomly sampled values for the years before. Soil water holding capacity (“bucket size”) for the study area was estimated between 11 and 33 cm based on water retention potential and soil wetness data from the Swiss soil suitability map (BFS, 2000; Lischke et al., 2006). The simulations were conducted considering the disturbing effect of rockfall, and with and without a random “background disturbance.” The rockfall disturbance was derived from three-dimensional rockfall simulations on the unforested slope following Moos and Lischke (in revision) and held constant over time (rockfall disturbance does not change over the considered time period). The “background disturbance” (Dist) is a spatially randomly distributed additional mortality that can represent natural disturbances or-very much simplified–forest management, for example. We implemented two scenarios of a low and a high background disturbance (Table 1). In the low disturbance scenario, a cell is on average disturbed once in 1,000 years with an intensity of 0.8 (80% of the cell are destroyed). Under the high disturbance scenario, a cell is on average disturbed once in 20 years with an intensity of 0.8. The latter represents very regular and rather intense forest management and/or natural disturbances. For both disturbance probabilities, we further tested a lower disturbance intensity of 0.3 to get an estimate on the sensitivity of forest development to disturbance intensity. (b) The simulated current forest state was validated based on data from field samples in recent and former forest intervention areas and samples in the rockfall zone (see section “Field Investigations”). Based on this validation, we used the spin-up of the current forest state with low background disturbance as initial state for the simulations with future climate.

Figure 2. Bioclimatic variables used for the TreeMig simulations: day degree sum (A), mean minimum winter temperature (B), and drought stress index (C) from 1970 to the beginning of the stabilizing period (2120). The spin up of the simulations run from 1616 to 2016, with climate data sampled from interpolated measured data from 1930 to 2016.

Table 1. Time period, climate data and disturbance regimes [disturbance probability (P_Dist) and intensity] for the spin-up forest simulation and the future forest scenarios.

(c) From this initial state, future forest development was simulated under warming conditions for 2016–2100 followed by a stabilizing period (2101–2200). Bioclimatic variables were derived from average daily temperature and precipitation data for a moderate (RCP4.5) and an extreme (RCP8.5) climate change scenario (Figure 2; CH2018 Project Team, 2018). The CH2018 scenarios are derived from the EURO-CORDEX ensemble of regional climate simulations. In CH2018, multiple Global Climate Models (GCMs) were dynamically downscaled by different RCMs for three future emissions scenarios (i.e., RCP2.6, RCP4.5, and RCP8.5). We chose RCP4.5 (moderate climate mitigation) and RCP8.5 (no climate mitigation) as the two most likely trajectories. CH2018 scenarios currently provide the best available dataset of future climate change estimates in Switzerland (Sørland et al., 2020). Gridded projections from CH2018 at a 1-km resolution were downscaled to a 25-m resolution based on a lapse rate (Broennimann, 2020). We linearly interpolated temperature and precipitation data from the available projections for 2035, 2060, and 2085 and added a random inter-annual variation based on the reference period of 1986–2016. The climate variables of the stabilizing period were randomly sampled from the last 20 years of the climate projection and complemented with the same random inter-annual variation.

Field Investigations

The simulated forest scenarios after the spin-up were compared to and validated with ground-truth data from the current forest. We measured adult trees and regeneration in 21 sample plots in areas with recent interventions (last intervention ≤ 10 years) and areas that have not been intervened in the past 10 years (Figure 1). The time since the last intervention was estimated based on a series of orthophotos (swisstopo, 2020a) and information from local forest managers. At each sampling site, we measured the diameter at breast height (dbh) and the species of each tree with a dbh ≥ 12 cm in a 20 m × 20 m plot (measured on the slope), and the species, age and height (for plants with a height ≤ 3 m) or dbh (for plants with a height > 3 m) of the regeneration in a 5 m × 5 m subplot. The subplots were selected based on the criteria that they were representative for the regeneration in the large plot. The plots cover the total slope and go beyond the area simulated with TreeMig (Figure 1).

Furthermore, data on species and dbh distribution (trees with dbh ≥ 8 cm) were available from 19 sample plots in the main rockfall zone (simulation area) from a previous study (Moos et al., 2020). Based on the data from these plots, we calculated tree density, species and diameter distribution for homogeneous forest areas (determined based on the orthophoto). We then generated trees with random positions in each area based on the derived tree densities and diameter distributions. This “treefile” served for validation of the rockfall risk assessment based on the simulated forest scenarios (see section “Protective Effect of Forest”).

Protective Effect of Forest

We quantified the protective effect of the different forest scenarios based on the reduction of rockfall risk for people and infrastructure at the bottom of the slope. The risk reduction is the difference between the risk without forest and the risk with forest (i.e., avoided costs). Rockfall risk is defined as the yearly expected damage of all possible damaging events (Agliardi et al., 2009). This is a useful indicator for a long-term perspective, although the damages actually caused can vary strongly from year to year. Rockfall risk for buildings and roads below and along the slope was calculated spatially explicitly following Moos et al. (2018b). The rockfall risk R_i,_j of a specific element at risk i and for a rockfall scenario j is defined as the product of the onset probability of rockfall scenario j (P_onset,j), the probability that the block j reaches the element at risk i (propagation probability P_prop,i,_j), the value of the element at risk E_i, its presence probability P_presence,i and its vulnerability V(I)_i,_j. The latter depends on the intensity I of the rockfall impact (quantified as kinetic energy). The risk is finally averaged over all possible impacts n_i,_j.

The total risk R is the sum of the risks of all elements at risk e and for all expected rockfall scenarios v:

By stopping or deviating falling blocks, the forest can on the one hand influence their propagation probability, and on the other hand reduce their energy (and thus their intensity). These two effects were quantified based on rockfall simulations.

Rockfall Simulations

Rockfall simulations were conducted with the three-dimensional rockfall module RockyFor3D (Dorren, 2016). This “probabilistic process-based” rockfall model simulates flying, bouncing and rolling blocks depending on terrain characteristics and standing trees. After each impact of a block on a tree, the energy reduction of the block is calculated depending on the tree species, the tree diameter and the impact height (Dorren, 2016). The rockfall release area was determined based on slope inclination and orthophotos. Rock density was set to 2,700 kg m^–3. Soil types and roughness were determined for homogenous areas based on land cover data (swisstopo, 2020b), the geological map (swisstopo, 2019) and the slope inclination and verified in the field. We simulated 100 blocks per block volume scenario and start cell (in total 684,000 simulations per block volume). Blocks were simulated with a cubic form with uniform edge length. The simulations were conducted with a digital elevation model with a resolution of 2 m. Buildings were considered as “rockfall nets” in the simulations with an energy absorption capacity of 600 kJ and their actual height. The rockfall simulations were conducted for the current forest cover (based on field samples and based on the forest simulations; see sections “Forest Scenarios” and “Field Investigations”) and for the future forest scenarios (after 50, 100, and 150 years of forest simulations; Figure 3). For this purpose, the cell-based output of the TreeMig simulations was translated into single tree positions by assuming a random tree distribution per 25 m × 25 m cell.

Figure 3. Schematic overview of the methodological approach of this study.

Rockfall Risk Calculation

We calculated rockfall risk per forest scenario for four different rockfall release scenarios with return periods of 10, 30, 100, and 300 years (Table 2). The respective block volumes were determined by expert estimation and validated based on rockfall deposits measured in the forest sample plots in the rockfall zone (Moos et al., 2020). For each scenario, we estimated a varying number of blocks potentially releasing per event (“event type” ET, Table 2). An ET > 1 means that multiple blocks of this block volume class release in the time period of the rockfall scenario (e.g., 20 blocks in 10 years). The resulting block volume distribution was generally in good agreement with the distribution of the blocks assessed in the field (power law coefficient of estimated volume-frequency distribution: –1.80; power law coefficient of block deposits: –1.38). We further estimated the overall onset frequency of blocks ≥0.05 m³ based on tree damages below the release area, following Trappmann and Stoffel (2013). We counted all visible fresh tree damages (younger than approximately 1 year) in a strip of 10 m width 20 m below the cliff (Moos et al., 2020). We then used the “conditional impact probability” (CIP) concept (Moya et al., 2010) resulting in a rockfall frequency of blocks ≥0.05 m³ of ∼30 year^–1 below the cliff. This is well in line with the derived rockfall release scenarios used in the risk analysis (P_onset (v ≥ 0.05 m³) = 34 year^–1 and P_onset (v ≥ 0.1 m³) = 2 year^–1).

Table 2. Rockfall scenarios considered in the risk analysis.

The varying number of blocks per rockfall scenario was considered probabilistically in the rockfall risk analysis based on a Monte-Carlo simulation. The risk calculation was repeated 50 times and for each repetition, the number of released blocks per scenario was sampled randomly. We finally calculated the mean risk for each rockfall scenario. Varying the number of releasing blocks allows at least for a partial representation of the uncertainty associated with the release frequency.

The propagation probability (P_prop,i,_j) of a block j at the element at risk i was determined based on the rockfall simulations (section “Rockfall Simulations”). It is defined as proportion of the potential rockfall trajectories of volume j through i (n_i,_j), and the total number of simulated trajectories of volume j (n_j,tot):

The energy of each block reaching an element at risk in the simulations was registered and classified according to the Swiss risk concept as low (0–30 kJ), medium (30–300 kJ) and high (>300 kJ) intensities (Bründl et al., 2009). We considered roads (with a width > 4 m), buildings and the railway below the slope as elements at risk. Object categories and occupancy were determined based on national dwelling data (FOS, 2020) and values and vulnerabilities of objects and persons were derived for the different intensity classes based on standard values of the Swiss risk concept (FOEN, 2020; Supplementary Tables 3, 4). We calculated risk only for direct impacts and infrastructure damage, and did not consider road/railway closure. Mean daily traffic (MDT) was determined based on census data from the cantonal authorities and estimations based on GoogleMaps.

The risk reduction (in CHFyear^–1) provided by a forest was finally calculated as the difference of the risk without forest and the risk with forest. It was evaluated for the current forest cover based on field data and the forest simulations, and for the future forest scenarios after 50, 100, and 150 years of simulated forest development. We only considered changes in forest cover under future climate, and no other potential (climate related) changes in risk, such as a shift in the release frequency and magnitude of rockfall.

Results

Forest Today: Modelling Results Compared to Field Data

The initial forest as generated by the TreeMig model after spin up had a mean basal area of 50.1 m² ha^–1 [standard deviation (sd.) = 15.4 m² ha^–1] with low background disturbance (Dist) and a mean basal area of 43.7 m² ha^–1 (sd. = 15.9 m² ha^–1) with high Dist. For both disturbance scenarios, highest basal area was found in the lower (<700 m a.s.l.) and upper (900–1,100 m a.s.l.) parts of the slope, whereas it was slightly lower at medium elevation. This is well in line with the actual measures of the current forest, which has a mean basal area up to 60 m² ha^–1 in the upper part of the slope and a slightly smaller basal area (45–50 m² ha^–1) at medium altitude, where the stand is generally denser (800–900 trees ha^–1) but with lower diameters (dbh) (mean dbh = 25–28 cm). In the upper part of the forest, the current forest has a tree density of 550–600 trees ha^–1 and a mean dbh of ∼35 cm. In contrast, TreeMig predicted a slightly denser forest (on average 870 trees ha^–1 with low Dist; 1,270 trees ha^–1 with high Dist), but with a lower dbh (mean dbh = 23.0 cm with low Dist and 20.5 cm with high Dist). Both the modelled forest and the actual forest are mainly dominated by F. sylvatica, mixed with P. abies, Acer pseudoplatanus, Fraxinus excelsior, and Tilia cordata (Figure 4). However, the model predicted a lower abundance of Larix decidua and A. alba than observed in the sample plots. In contrast, there were more Quercus pubescens in the modelled forest.

Figure 4. Mean basal area and standard deviation of the species measured throughout the study area and/or modelled by TreeMig for the spin-up simulation with background disturbance (“modelled_Dist”; red), without background disturbance (“modelled_noDist”, blue), and based on field measurements in sample plots in recent and older intervention areas (“field”; orange) and sample plots in the rockfall zone (“field2”; green), where deciduous tree species with a low abundance were summarized as “deciduous”. The standard deviation represents the variation in the basal area between the sample plots (field data) and the cells (modelled data), respectively.

The regeneration in the modelled forest as well as found in the field samples was also clearly dominated by F. sylvatica, followed by F. excelsior (Figure 5). The regeneration density of P. abies, Acer sp., and of P. sylvestris and Q. pubescens (in simulations with Dist) were higher in the simulations than in the field data. The regeneration density in the present forest clearly decreased with increasing age of the intervention and density of mature trees (see Supplementary Figure 2).

Figure 5. Mean density and standard deviation of the regeneration of the species found in the field samples in recent and older intervention areas (“field”; orange), and modelled in TreeMig with background disturbance (“modelled_Dist”; red) and without background disturbance (“modelled_noDist”; blue).

Forest in Future: Modelling Results

The forest simulations predicted a clear decrease in basal area of most of the currently present species under climate scenario RCP4.5 and RCP8.5 after 150 years (Figure 6). The decline started after ∼20–30 years of simulation (∼2040–2050; Figure 7). The forest finally turned into a Q. pubescens dominated forest, for both scenarios, and mixed with P. sylvestris in RCP4.5 (Figures 6, 8). The total basal area of the forest declined to 20–50 m² ha^–1 (mean BA = 27.5 m² ha^–1) in RCP4.5 and to 10–50 m² ha^–1 (mean BA = 19.8 m² ha^–1) in RCP8.5, whereas the strongest decrease can be observed in the steep part of the slope (500–900 m a.s.l.; Figure 8), where a relatively dense forest of mainly Q. pubescens with low diameters was predicted (Figure 8). F. sylvatica completely disappeared under the extreme climate scenario RCP8.5, whereas the model predicted an average basal area of ∼5 m² ha^–1 for the climate scenario RCP4.5. A. alba and P. abies disappeared in both scenarios. A small abundance of L. decidua remained in RCP8.5, while a slight increase is predicted in RCP4.5. In both scenarios, the basal area of F. sylvatica (RCP4.5) and Q. pubescens is slightly smaller with enhanced disturbance (Figures 6, 7). No difference or even a slight increase was found for P. sylvestris and L. decidua (Figure 6, 7). In contrast, an increase in the regeneration density with increasing disturbance could be observed for Q. pubescens, P. sylvestris, L. decidua, and B. pendula after 50 years of simulation (Figure 9). The comparison to simulations with the same disturbance probabilities, but a low disturbance intensity (0.3) did not reveal a significant difference in the predicted basal area (not shown).

Figure 6. Basal area per main species after 150 years of simulation with climate scenario RCP4.5 (above) and RCP8.5 (below) and for low (blue) and high (red) background disturbance (Dist) as well as the basal area of the current forest (orange).

Figure 7. Simulated evolution of the basal area of the main species in Martigny under climate scenario RCP4.5 (yellow) and RCP8.5 (orange) and for low (points) and high (triangle) background disturbance (Dist).

Figure 8. Current (“today”; based on field data) and future (“RCP4.5” and RCP8.5” in 50 years) mean basal area per elevation level and main species. Species abbreviations: AA, Abies alba; AP, Acer pseudoplatanus; FE, Fraxinus excelsior; FS, Fagus sylvatica; LD, Larix decidua; PA, Picea abies; PS, Pinus sylvestris; QP, Quercus pubescens; TC, Tilia cordata.

Figure 9. Simulated mean density and standard deviation of the regeneration of the main species under the future climate scenarios RCP4.5 (above) and RCP8.5 (below) in 50 years and for high background disturbance (“future_Dist”; red) and low background disturbance (“future_noDist”; blue).

Protective Effect of Forest

The current forest substantially reduces rockfall risk for subjacent buildings, roads and the railway. Today, a total risk of 19,100 CHFyear^–1 is expected for people and infrastructure (current forest based on field samples; Figure 10 left and Supplementary Figure 4). With the current forest simulated with TreeMig, risk was slightly higher. Without forest, the risk would be five times as high. The current forest consequently provides a risk reduction of ∼1,350 CHF per ha forest and year. The strong increase under non-forested conditions was mainly caused by small block volumes (≤1 m³), which are currently stopped by trees and rarely reach the roads or houses (see also Supplementary Figure 3).

Figure 10. Mean risk for all considered elements at risk for the current forest based on field data (top; “field_data”), based on TreeMig simulations (tm) with (above; “tm_today_Dist”) and without (left, “tm_today”) disturbance, and for future forest scenarios in 50 (tm_50yr), 100 (tm_100yr), and 150 years (tm_150yr) under climate scenario RCP4.5 (middle) and RCP8.5 (right) (only shown with background disturbance).

With climate warming, a clear increase in risk is expected for both climate change scenarios (Figure 10 middle and right). In ∼150 years, an increase in risk to 59,250 CHFyear^–1 (two times the risk of the spin-up forest) and 70,100 CHFyear^–1 was predicted (2.5 times the risk of the spin-up forest). While highest risk under scenario RCP4.5 was expected after 100 years (69,700 CHFyear^–1), followed by a decline, a steady increase in risk could be observed for RCP8.5. The risk increased mainly due to an increase in the occurrence frequency of blocks. Currently, a block with a volume between 0.1 and 0.5 m³ is expected to reach the railway about once in 10 years. In future, its occurrence frequency could increase to once in 3 years (Figure 11).

Figure 11. Predicted development of the occurrence frequency of the considered block volume scenarios at the level of the railway until 2200 for the future forest scenarios under climate scenario RCP4.5 (above) and RCP8.5 (below).

Discussion

The results of our case study suggest that the protective effect of forests against rockfall may substantially decrease under future climate change at sites where significant shifts in tree species distribution due to increased drought stress can be expected. The reduced protection effect observed in this study is mainly due to a decrease in the basal area in the middle part of the mountain slopes (at mid-elevations), where drought stress will be particularly pronounced due to the steep slopes with shallow soils. Under both climate scenarios, our simulations predicted a shift towards a Q. pubescens dominated forest, mixed with P. sylvestris, leading to significant decrease in tree diameters, partially also a reduced stem number under the extreme climate scenario. This would result in a substantial increase in risk, and thus in the expected damages due to rockfall for subjacent houses and infrastructure in the next decades. In the medium term (i.e., 50–100 years), the results of this study suggest a stabilization of the risk, or even a slight decline in the long-term (i.e., >100 years) since the Q. pubescens dominated forest would increase in biomass.

The comparison of detailed field data on the species and diameter distribution, the tree density and the regeneration of the current forest with the TreeMig simulations allow concluding that the overall accuracy of the model was satisfactory. The slightly higher abundance of F. sylvatica and the lower abundance of L. decidua and P. sylvestris in the spin-up simulations indicate that the current forest corresponds to an earlier stage in the forest simulations (see Supplementary Figure 1). This is possibly an effect of local forest management, since regular forest interventions in the past decades (or even centuries; Farquet and Métral, 2004) accelerated a continuous regeneration, without reaching the equilibrium state of the forest.

Limitations of the Study

The predictions of future forest cover are, however, subject to several uncertainties. Firstly, while species parameters in the TreeMig model are primarily based on general ecological knowledge, partly they are based on observational data sources, which might be of low precision. For example, in the Swiss NFI, only few plots are located in dry conditions, so that an accurate estimation of the dry distribution limit of the species is difficult. Furthermore, the observed forests might be influenced by other factors (e.g., management history). Secondly, they are based on species behavior under past and current climate, while some species may adapt to changing climatic conditions (Kurparinen et al., 2010; Lefèvre et al., 2013). Thirdly, the soil conditions, strongly affecting the drought experienced by the trees in both model and reality, are spatially variable and hard to cover by measurement based empirical models such as the underlying soil suitability map. Fourthly, the TreeMig model is particularly sensitive to extreme drought events compared to other forest landscape models, due to a different formulation of mortality, whereas it is not clear which formulation is more realistic (Petter et al., 2020). Finally, climate change scenarios used in this study are highly uncertain (CH2018 and 2018) and the downscaling of the data added a potential bias. In general, our simulation results are in line with previous modelling studies, predicting a decline in biomass at lower elevations in the Valais or comparable regions as well as a shift from F. sylvatica to P. sylvestris (Schumacher and Bugmann, 2006; Speich et al., 2020) and/or an expansion of Q. pubescens (Vacchiano and Motta, 2015). The analysis of long-term data further revealed an increasing abundance of Q. pubescens in low-elevation forests in the Valais in the past few years (Weber et al., 2007; Rigling et al., 2012). Finally, observations of the local forest managers confirmed increasing problems particularly of F. sylvatica in drought periods in the past few years (R. Métral, personal communication).

The rockfall risk with the forest simulated in the spin-up (initial state for future simulations) is higher than the risk with the forest based on field samples, since the basal area of the latter is slightly higher than the basal area simulated with TreeMig. This suggests that risk under future scenarios might also be somewhat overestimated due to a potentially general underestimation of the basal area by the model. Additionally, TreeMig is a cell-based forest model and thus remains limited in representing local differences in forest structure potentially relevant for the risk reducing effect of the forest. However, the current forest does not represent true positions of trees, but is constructed based on forest samples. Consequently, its protective effect might also be biased because of local spatial variations in the forest structure (e.g., insufficient representation of gaps). Combining field data with tree positions derived from remote sensing data may increase the accuracy of the representation of the current forest cover (Mathys et al., 2004; Engler et al., 2013). Overall, the calculated rockfall risk of today seems plausible. Based on our risk assessment, a block with a volume between 0.1 and 0.5 m³ reaches the bottom of the slope about once in 10 years. Although hardly any events have been registered in the past few years, this is quite realistic based on the observations of the local forest engineer, as well as our own observation of rockfall deposits at the bottom of the slope, demonstrating that blocks regularly reach the lower part of the slope. Nevertheless, the risk assessment faces several sources of uncertainties, including, among others, the rockfall onset frequency, which was based on a scarce data basis, and the input parameters of the rockfall simulations. Additionally, the parameters for the quantification of the elements at risk are highly disputable due to missing observational data (Agliardi et al., 2009; Daniell et al., 2015). These parameters are also likely to change in future due to changing economic and social dynamics (Fuchs et al., 2013). Moreover, rockfall release frequency and magnitude may be affected by climate warming and thus change in future. Recent observations indicate, however, that climate change will mainly affect high mountain areas, above the permafrost limit, but not significantly alter rockfall frequency at lower elevations (Sass and Oberlechner, 2012). Since these potential additional changes are independent from the evolution of the forest, the relative differences between the forest scenarios do not change and, thus, the simplified assumptions regarding the risk assessment underlying this study remain reasonable.

Generalization of the Results

The results of our study are not only relevant for our case study site, but can provide important information for other sites with similar forest and rockfall conditions in the Alps (i.e., sites with similar species distribution and climate predictions). By combining dynamic forest modelling with rockfall risk assessment, it provides a methodological framework that can be applied at other case study sites, where climate change has not yet been considered in the forest scenarios (Bigot et al., 2009; Moos et al., 2019a) or even be transferred to regional scale (Scheidl et al., 2020b). Our results also indicate that regular disturbances can promote regeneration, and therefore favor the transition of current species to more drought resistant ones, which may positively influence the protective function of the forest in long-term (Scheidl et al., 2020a). Consequently, forest management may play a crucial role in accelerating the regeneration and immigration of drought adapted species and thus increase the resistance and resilience of the protection forest analyzed in this study under a changing climate (Brang, 2001). At the same time, forest interventions as well as natural disturbances can cause temporary declines in the protective effect of a forest due to the removal of trees (Bigot et al., 2009; Irauschek et al., 2017). Temporary measures, such as lying logs in openings, can at least partially bridge such periods of reduced protection. Furthermore, a sudden decline in the protection capacity of a forest can mostly be avoided through well planned and site-adapted management interventions (Woltjer et al., 2008). According to our analysis of current regeneration as a function of intervention age, an intervention interval of ∼30 years is indicated to maintain a constant level of regeneration (see also Supplementary Material 4). This may, however, strongly be influenced by the prevailing browsing pressure (Cailleret et al., 2014; Kupferschmid et al., 2019). Additionally, large disturbances, such as fire, insect outbreaks or wind throw, are likely to become more frequent and/or severe with climate change and might have even more drastic effects on species composition and structure than the direct effects of climate change as analyzed in this study (Bebi et al., 2017; Seidl et al., 2017; Wohlgemuth et al., 2017; Taeroe et al., 2019).

Future Avenues

Future research on the long-term efficacy and efficiency of forest management is needed to finally conclude, under which conditions, where and when, targeted forest interventions are appropriate, and in which form (e.g., planting of drought adapted species, assisted migration), to improve the protection function of forests. Furthermore, future studies should quantify the influence of changing natural disturbance regimes, which was only included as random background disturbance with varying probabilities and intensities in this study. In general, higher frequencies and intensities of large-scale disturbances can be expected (Seidl et al., 2017). In the specific context of our study area, an increasing risk of forest fires has to be reckoned (Schumacher and Bugmann, 2006; Wastl et al., 2012). Furthermore, it is probable that new and invasive species, such as Ailanthus altissima or Paulownia tomentosa (of which single exemplars can already be found on the slope or nearby), will immigrate and maybe even suppress some current tree species. Recent studies, however, indicate that these species should not prevail on larger scales (Knüsel et al., 2017; Wunder et al., 2018), and that the protection capacity of the single tree is comparable to that of other deciduous trees (Moos et al., 2019b). Large-scale modelling of forest development including species’ natural or assisted migration, e.g., with TreeMig, would shed light on the role of immigration of new species.

Conclusion

This case study evidences that climate change may lead to a decline in the protective effect of forests against rockfall. By combining dynamic forest modelling with quantitative rockfall risk assessment, we could show that increased drought at lower and medium elevations in the Swiss Alps potentially could cause critical periods of a decreased protection service of forests due to shifts in species distribution. This means that the risk for subjacent houses and infrastructure may substantially increase in future. Hence the need for engineered protection measures, such as rockfall nets, may arise, and cause additional costs. Our results further indicate that forest management can potentially promote the transition to more drought-adapted species and thus dampen the strong decrease in the protective effect of the forest. However, future research should address the long-term effect of different forest management strategies on the protection service of forests under climate change. The study was conducted on local scale, but the underlying methods as well as the results can be transferred to other sites with similar conditions and serve decision-makers as an important basis in natural hazard and ecosystem management. Upscaling the methods used in this study to a larger scale will further facilitate detecting areas where a particularly strong decrease in the protection service of the forest can be expected and thus facilitate prioritization in forest management.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

All authors contributed to the conception and method of the study. CM organized the field investigations and conducted the data analysis, the simulations and the risk analysis, and wrote the first draft of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

The funding of this research was provided by the “Interdisciplinary Centre for Mountain Research” (CIRM) of the University of Lausanne as part of a postdoctoral project on the temporal evolution of the protective effect of forests against rockfall.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank R. Métral for the information on local conditions and the valuable discussions about this study, as well as N. Khelidj and Z. Bontognali for helping in the field.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/ffgc.2021.682923/full#supplementary-material

References