Assessing the Effect of Disturbances on the Functionality of Direct Protection Forests

Forests provide direct protection to human settlements from hydrogeomorphic hazards. This paper proposes a method for assessing the effect of natural disturbances on the functionality of direct protection forests (DPFs) in order to prioritize management interventions. We georeferenced disturbance data for wildfires, wind and snow damage, avalanches, and insects and overlaid them to a region-wide DPF map. Within each disturbance polygon, we used a Landsat-5 TM image to identify DPFs with insufficient vegetation cover, by using a maximum likelihood classifier of 6 spectral bands plus 5 vegetation indices. For each disturbance agent, we fitted a generalized linear model of the probability of finding a forested pixel, as a function of topography, time since disturbance, distance from disturbance edge, summer precipitation, and drought in the disturbance year. DPFs covered almost half of total forest area in the study region. Disturbance by insects occurred in more than one sixth of all forests. Avalanche and wildfire occurred each in about one tenth of total forest area, and wind and snow disturbance in only 1%. In the last 50 years, disturbances had a recurrence rate of 3% every 10 years. Almost one sixth of DPFs are currently lacking sufficient forest cover. Wildfires resulted in the highest rate of nonforested pixels (42% of all DPFs), followed by avalanches (21%). Forest recovery was explained by time elapsed, distance from edge (for conifers), and aspect. Summer precipitation and drought had a mixed influence. Our approach to assessing the effect of disturbances on the functionality of DPFs is reproducible in all mountain regions using institutional or open-access geographic data and provides a tool to prioritize DPF management by indicating where restoration of protection is most urgent.


Introduction
Mountainous areas are associated with a number of hydrogeomorphic hazards, including rockfall, landslides, snow and rock avalanches, and debris and earth flows (Korup and Clague 2009). Forests may provide direct protection to human settlements and activities from some of these hazards (Motta and Haudemand 2000;Dorren et al 2005;Covi 2008;van Beek et al 2008;Bebi et al 2009), by preventing their occurrence (eg by snowpack stabilization) or mitigating their intensity (eg by slowing or blocking falling rocks) . This type of forest is referred to as a direct protection forest (DPF).
The extent of forests managed for their protection services is included in the 1995 Montreal Process (Brand 1997), the Helsinki Process (MCPFE 2002), and the Forest Resources Assessment of the Food and Agriculture Organization (FAO 2000) as an indicator of sustainable forest management in mountain ecosystems. This service (Dorren et al 2004) has long been recognized in areas with exposure to hazards, such as the European Alps, and this has resulted in bans on timber harvesting since at least the 14th century (Sablonier 1995).
Currently, DPFs must be identified and undergo special management under national or regional forest regulations in Austria, France, Germany, Italy, and Switzerland (Weiss 2000;Brang et al 2006;Wehrli et al 2007;Brosinger and Tretter 2009), as well as Greece (Zagas et al 2011), Japan (Shimamura and Togari 2006), and Turkey (Aydin and Dorren 2010). In other parts of the world, the generic role of forests in soil protection or hazard mitigation is recognized, but formal identification of DPFs is surprisingly lacking, even in areas with a high risk of hydrogeomorphic hazards, such as winter resorts in North American mountain areas (Ives et al 1976).
The effectiveness of DPFs depends not only on their location, but also on elements of forest structure that maximize physical hazard mitigation, such as density, tree diameter, and canopy cover (Bebi et al 2001;Frehner et al 2005;Berretti et al 2006;O'Hara 2006;Teich et al 2012;Dorren et al 2015). Moreover, like all ecosystems, DPFs are subject to natural disturbances (White and Pickett 1985)-including wildfires, windstorms, herbivores, and wood-boring insects-that can temporarily disrupt their structure and therefore their functionality. Disturbances are increasing in magnitude and frequency across the Alpine region (Sass and Oberlechner 2012;Seidl et al 2014b). The response of ecosystem functions to singlepulse or interacting disturbances can follow various trajectories: short-term change and recovery, gradual continuous modification, or abrupt transition to a different state (Scheffer et al 2012). The speed, direction, and mode of ecosystem response depends on ecosystem resistance and resilience (Holling 1973), and on the type, frequency, time of occurrence, intensity, and extent of the disturbance, considered either as a single event (eg Vacchiano et al 2015b) or integrated over a disturbance regime (Sousa 1984;Turner et al 1993).
An undisturbed DPF is not automatically functional, because age, site, or other limitations can result in forest structures being temporarily or permanently unsuitable to mitigate hydrogeomorphic hazards. However, rapid assessment and urgent action are often needed following sudden structural changes by natural disturbance agents. Building on previous knowledge on optimal stand structure in DPFs, this study aimed to assess the effect of abiotic and biotic disturbance agents on the current functionality of DPFs, measured by their vegetation cover. In particular, we integrated spatially explicit information on DPFs, time-labeled perimeters of natural disturbances by wildfire, wind and snow damage, avalanches, and insects for variable time periods in the last 50 years, and classified remotely sensed land cover into forested, herbaceous, and nonvegetated areas. The overlay of such information allowed us to assess the temporal recovery of forest structure following different disturbances, model the effect of disturbance size, time of occurrence, and other environmental variables on DPF effectiveness, and prioritize DPF restoration.

Study area
The study area, the Aosta Valley region in northwest Italy (Figure 1), covers 3262 km 2 and has a topography shaped by an east-west valley with several north-south protrusions, legacies of the last ice age. Mean elevation is 2100 m above sea level (masl). More than 90% of the land has an elevation of at least 1000 masl, and 60% is above 2000 masl. Climate is continental with cold winters and hot summers; mean annual rainfall in Aosta city (45u269N, 7u119E) is 494 mm (Biancotti et al 1998), with a dry season extending from June to September. Winter precipitation usually comes as snow. The study area exhibits both crystalline (granite) and metamorphic bedrocks, but most of the landscape is covered by quaternary deposits of glacial, gravitative, or colluvial origin. The pedological background is dominated by shallow (Lithic, Umbric, and Dystric Leptosols) or eroded soils (Eutric and Calcaric Regosols), acid soils with organic matter, iron oxides and aluminum accumulation (Dystric Cambisols, Haplic Podzols, and Humic Umbrisols), or alluvial Eutric Fluvisols in the valley bottoms (Costantini et al 2004).
Hydrogeological protection is one of the most important ecosystem services provided by forests in the region. One fourth of total forest area is susceptible to hydrogeomorphic processes, including landslides (2700 ha), avalanches (4000 ha), and rockfall (14,600 ha) (INFC 2007). Most forests (80%) sit on slopes .40%; this rather conservative threshold, together with watershed delineation algorithms, has been previously used to map potential DPFs as a subset of total regional forest cover . Potential DPFs cover 40,557 ha-12.4% of the region's territory and 42.7% of total forest cover ( Figure 2). Larch is the most common forest cover type with a protective function; however, chestnut has the greatest share of its forests in DPFs, 58% .
Wildfires occur mainly in late winter and early spring and are 95% anthropogenic in origin (Vacchiano et al 2015a). Surface fires usually start at the bottom of the valley and spread upward, often transitioning into crown fires in conifers due to the low moisture content of live foliage during the dormant season. Most wildfires are small (average size 5 8 ha); however, 70% of the burned area is due to sporadic but relatively large (.20 ha) wildfires. This fire regime is common to several inner Alpine valleys (Moser et al 2010;Zumbrunnen et al 2011) and reflects the influence of fire suppression policies adopted throughout the Alpine region (Pezzatti et al 2013;Valese et al 2014).

Natural disturbances
We selected for this analysis the following disturbance agents: wildfires, wind and snow, avalanches, and 2 biotic agents, the nun moth (Lymantria monacha L.) and the larch tortrix (Zeiraphera griseana Hü bner), which have undergone multiyear or recurring outbreaks in the study region (Baltensweiler and Rubli 1999;Pasquettaz and Negro 2007;Bottero et al 2013). Other biotic agents were not investigated due to the episodic nature of the outbreaks or insufficient temporal data. The following data sources were accessed to obtain a spatially and temporally explicit data set on natural disturbances for all forests in the study region. N Wildfires: We georeferenced all fire perimeters from a regional forest fire archive for the years 1961-2010. Fire perimeters had been mapped visually or by GPS (Global Positioning System) by Regional Forest Service staff. Fires smaller than 6 ha were recorded only from 1987; fires smaller than 0.5 ha were recorded since 1995 as points and were georeferenced as such.
Beginning in 1987, additional information was available for each fire, such as date and time of the fire, duration, cause of ignition, land use, main forest species, elevation of the ignition point, fire behavior type (surface, crown, or underground), and number of rainless days before ignition.
N Wind and snow damage: We queried all reports of windor snow-induced damage (excluding avalanches) in the Regional Forest Service's annual forest health reports for 1990-2010. Reported variables included the location and date, type of event, extent of damaged area, average or minimum/maximum elevation, management unit identification number, mean age and diameter of the damaged stand, forest species involved, and disturbance severity (%). Uncharted events were assigned a circular shape centered on the most likely location of the event, using location, elevation, forest cover, and extent of the damage as proxy information.
N Avalanches: We queried all avalanche polygons in the regional avalanche registry (source: Ufficio Neve e Valanghe-RAVDA), which was established in 1970 and digitized in 2010 (Lunardi et al 2009). Polygons were outlined so as to reproduce the maximum run-out area based on Forest Service yearly reports, historical photographs (more than 31,000 images from 1971 to 2010), and historical newspaper articles (250 from 1800 to 1900).
N Insect outbreaks: We georeferenced available maps of L. monacha (1984-1990) and Z. griseana (1994-2006) outbreaks. Forest health reports were consulted for additional data for 1990-2010, in the same process reported above for wind and snow damage. For Z. griseana only, disturbance severity (%) was available and integrated in the data set.
We computed the following summary statistics (overall and for each disturbance agent): total and yearly disturbance frequency, disturbed area, disturbance rotation period (Agee 1993), and relative frequency of each disturbance agent in total forest area and in the DPF subset by elevation, slope, aspect class, and plant functional type (ie with or without broadleaves) according to 1990 Corine Land Cover level-3 data (European Environment Agency 2014). Areas disturbed more than once (eg by reburn) were repeatedly accounted for in overall area statistics. For each disturbance agent, we tested whether the frequency distribution of disturbed areas could be fitted by a negative power law (Kolmogoronov-Smirnov 1-sample test), as is expected under the "self-organized criticality" hypothesis (Malamud et al 1998;Turcotte and Rundle 2002;Espírito-Santo et al 2014). Finally, we assessed recurrence of each disturbance agent and interactions between disturbances by cross-tabulating the proportion of total forest area disturbed by any temporally ordered binary combination of disturbance agents, relative to the total forest area disturbed by each agent.

Residual protection function in disturbed forests
For legal purposes in Italy, any area designated as forest that would-at some point-lose its tree cover because it has been disturbed is still classified as forest. For this reason, and under the consideration that these sites may sustain a forest cover potentially functional for hydrogeologic hazard protection, we analyzed the functionality of mapped DPF based on their current status as forested versus nonforested. We used remotely sensed imagery to identify disturbed DPFs with insufficient vegetation cover, which needed urgent management action to restore their protective function. We obtained a recent, cloud-free Landsat 5 TM image of the study region (WRS-2 path/row: 195/28) taken during the vegetative period (image taken on 20 September 2011; scene ID: LT51950282011263MOR00). Image resolution was 30 m. This is larger than the minimum area considered conducive to avalanche formation in forests, that is, a width of 15 to 20 m (Schneebeli and Meyer-Grass 1993). It is insufficient to detect a snow gliding hazard (Viglietti et al 2013), but it is the highest resolution currently available for unrestricted-access, large-scale multispectral Earth imagery.
For each Landsat band, we converted digital numbers to at-sensor radiance and then to top-of-atmosphere reflectance by applying image-specific rescaling coefficients and illumination parameters. Top-ofatmosphere values were atmospherically corrected based on dark-object subtraction, which assumes that any radiance received at the sensor for a dark-object pixel is due to atmospheric path radiance (Chavez 1996). The reflectance threshold for dark objects was selected in each band as the highest top-of-atmosphere value with an image-wide frequency of ,1000 pixels. This value was then subtracted from the reflectance across the whole scene to reduce scattering influences. Topographical correction was achieved by the C-correction method (Teillet et al 1982) modified with an additional inverse-squared elevation factor to account for atmospheric thickness.
We computed 5 reflectance indices recommended by the literature for the analysis of post-disturbance vegetation recovery (Table 1). The Forest Disturbance Index was computed following a tasseled cap transformation of the 6 TM reflectance bands (Crist and Cicone 1984); preliminary to this calculation, the 3-tasseled cap indices were normalized using their respective mean and standard deviation, computed over all pixels contained within polygons of a preexisting map of regional forest cover .
After filtering out water-and snow-covered pixels, we classified current land cover into 4 classes: forested, meadow (ie mowed and fertilized), pasture, and rocks/bare soil. To do so, we visually delimited 3 to 5 regions of interest for each land-cover class on the corrected Landsat image, and used them to train a maximum likelihood classifier (Foody et al 1992) based on 6 Landsat bands (excluding thermal infrared) plus the 5 disturbance indices. Minimum probability for successful classification was set to 80%; the land cover map was validated by computing classification error for each class in a random sample of 200 points (50 per class) extracted outside of the selected regions of interest and compared against visual classification of the most recent Google Earth image into the same land-cover classes. Finally, we produced summary statistics and maps of current land cover in DPFs, and compared the relative frequency of land cover classes inside and outside disturbance polygons. We assumed that DPF pixels currently covered with meadow, pasture, or rocks and bare soil were not functional.

Effect of disturbances on DPFs
In order to assess the ability of DPFs to recover from disturbances at the pixel (30 m) scale, we modeled the probability of each DPF pixel being functional (1 5 forested) or nonfunctional (0 5 nonforested) by means of generalized linear models, as functions of disturbance size, time since disturbance, topography, and climate. The response variable was modeled by a binomial distribution (logit link); a correction for zero inflation was applied when necessary (Hall 2000). A total of 9 models were built, 1 for each combination of disturbance type (crown and surface fire were fitted in separate models) and plant functional type (with or without broadleaves), as regeneration follows fundamentally different dynamics in each of these situations (eg Pausas 1999;Chen et al 2009;Ilisson and Chen 2009;Marzano et al 2012;Bergeron 2014;Vacchiano et al 2014). Combinations resulting in insufficient data were pooled (eg there was only 1 model for wind and snow disturbance). Mapping of plant functional types based on Corine Land Cover data (European Environment Agency 2014) required the consideration of a slightly different number of forest pixels than mapping based on the most recent forest cover map . This was especially true for areas at the tree line (eg when analyzing disturbance by Zeiraphera that occurs on high-elevation larch forests). Variables included the following: 1. Elevation, slope, and aspect (linearized by cosine transformation), from a regional digital elevation model at 10 m resolution. 2. Time since disturbance as of 2011 (not analyzed for avalanches). 3. Mean summer precipitation for 1950-2000 at 30 arcsec resolution from the Worldclim climate grid (Hijmans et al 2005). Worldclim data are strongly correlated to elevation in mountainous terrain; in order to avoid collinearity between independent variables, we computed rainfall variation independent of elevation using the residuals of a linear model between the two variables. 4. Drought in the growing season of the disturbance year (not analyzed for avalanches), expressed by computing the Standardized Precipitation Evapotranspiration Index (Vicente-Serrano et al 2010) for April to August using historical monthly data from the Aosta meteorological station (Figure 1). 5. Minimum distance to the disturbance polygon edge (only for wildfires, wind and snow, and avalanches, as polygons for Lymantria and Zeiraphera had already been aggregated by the data provider from smaller disturbance units).
All variables were spatially resampled at a 30 m resolution, filtered to avoid collinearity (Pearson's R , 0.7), and standardized in order to compare effect sizes. Generalized linear model goodness-of-fit was measured by percent deviance explained. For each model, we computed both absolute and standardized effect sizes, and scrutinized residual patterns in order to ascertain the absence of nonlinear effects. All analyses were carried out using the package glmmADMB (Bolker et al 2012) from the R statistical package 3.1 (R Core Team 2015).

Results
Data for different disturbances were available for different but overlapping periods (Figure 3). The disturbance database included 1765 wildfires, of which 1195 were smaller than 0.5 ha, 52 wind and snow damage events, 2538 avalanches, 7 yearly totals for L. monacha outbreaks, and 56 individual outbreaks of Z. griseana. Avalanche data were not available broken down by year. During the years for which data are available, avalanches disturbed the greatest land area, 52,975 ha. The total land area disturbed by insects was 16,813 ha, by wildfire was 7470 ha, and by wind and snow was 1093 ha (Table 2). When scaled to the periods of data availability, and relative to the total forest cover of the region, these figures result in a disturbance rotation time of 1806 years for wind and snow damage, 659 years for wildfires, and 120 years for insects.
Forest disturbance frequencies varied a great deal over time, with exceptional years for wildfire damage and wind and snow damage, and pulses of insect outbreaks showing up with multiyear frequency (Figure 4). The largest wildfire covered about 400 ha; avalanche areas can be much larger, up to 1650 ha. However, most of the events were small, especially the wildfires. Except for insect outbreaks, for which only area-based estimates were available instead of single-event records, none of the disturbance agents could be fitted by a negative power law distribution of disturbed areas (Kolmogorov-Smirnov test, P . 0.05). The proportion of small events was higher than that expected ( Figure 5), indicating that large, infrequent disturbances may be missing from the data set.
Zeiraphera was the most frequently recurring disturbance agent in the study region: more than 21% of disturbed forests were hit more than once (Table 3). Avalanche data, which were not temporally explicit, could not be evaluated for this comparison. When normalized by the number of monitoring years, wildfires, wind and snow, and Lymantria outbreaks had similar recurrence rates (0.29, 0.23, and 0.27% per year, respectively). The interaction between wildfire and wind and snow disturbance was negligible; in contrast, a non-negligible 13% of areas disturbed by wind and snow were subsequently disturbed by Zeiraphera or Lymantria, while 5% of Lymantria sites were subsequently disturbed by wind or snow.
Disturbances in DPFs occurred according to definite gradients of elevation, slope, and aspect. Fire and wind/ snow followed opposite topographic trends: DPFs were more likely to be disturbed by fires on low-elevation, steep, south-exposed slopes, and by wind and snow on mid-to high-elevation, north-exposed, gently sloping sites. Avalanche incidence was much higher at the tree line (2200 masl and higher), and exhibited no clear trend related to slope or aspect. The highest incidence of insect outbreaks centered at ,1800 m of elevation, with Lymantria occurring almost exclusively on north-exposed slopes and Zeiraphera preferring high-elevation, steep sites. All disturbances except wildfires had a higher incidence in conifer than in broadleaf DPFs (Supplemental Material, Figures S1-S4; http://dx.doi.org/10.1659/MRD-JOURNAL-D-15-00075.S1).
Temporal data not available for avalanches.

b)
Data on insect outbreaks were summarized by yearly damage and did not necessarily correspond to individual disturbance events. Classification of the Landsat scene by maximum likelihood produced a 4-class land cover map that was successfully validated against the most recent Google Earth imagery. Although meadows and pastures were often confused with each other (Table 4), omission and commission errors for forested cover, rocks and bare soil, and meadows and pastures (pooled) were always , 12%. Based on these figures, we calculated a weighted Cohen's kappa (Cohen 1968) of 0.89. The resulting land cover classification showed that DPFs in 2011 were covered 13.6% by herbaceous vegetation and 0.8% by rocks or bare soil. Disturbed DPFs did not exhibit a significantly lower proportion of forested cover than nondisturbed ones (x 2 test, P , 0.05). However, proportions of bare soil and rock, meadow and pasture, and forested pixels in both DPFs and total forest area differed significantly among disturbance agents (P , 0.05). Wildfire resulted in the highest rate of nonforested DPFs ( Figure 6: 450 out of 630 ha, ie 42% of all DPFs), followed by avalanches (21%), Zeiraphera (20%), Lymantria (6%), and wind and snow (4%).

MountainResearch
Generalized linear models of current functionality of disturbed DPFs had a low explanatory power (Table 5; deviance explained). However, even in the context of a noisy data set, some variables had a significant and consistent effect (Table 5 and Supplemental Material, Figures S5-S13; http://dx.doi.org/10.1659/MRD-JOURNAL-D-15-00075.S1). In conifer forests, distance from the disturbance edge was always negatively related (by 1-37% for every 100 m) to the probability of a DPF pixel being forested. Time elapsed since disturbance was always positively related to the probability of forested cover (by 0.1-1% per year). Among topographic variables, lower elevations and northerly aspects, when significant, were responsible for a more successful forest recovery after disturbances of all types. After wildfire, probability of forested cover dropped 2-4% for every 100 m increase in elevation. Aspect was the most important variable in 6 out of 9 models. Slope had a positive influence on forest recovery after wind and snow, L. monacha, and avalanche disturbances (conifers only), and a negative influence after Z. griseana and fire disturbances (broadleaves only). A strong and positive influence of mean summer precipitation was found after crown fire in conifers, surface fire in broadleaves, and L. monacha disturbances. Post-disturbance water deficit (Standardized Precipitation Evapotranspiration Index) had a negative effect on recovery after fire (conifers only), but positive after wind and snow and L. monacha disturbance.

Discussion
Avalanches dominated the regional disturbance regime across all forests and in DPFs. Avalanches of all magnitudes from 1971 to 2010 covered 10% of regional forest area. Because regional avalanche data were not temporally explicit, it was impossible to assess return intervals and recurrence rates. Even assuming that annual avalanches represent a small fraction of the cumulative area disturbed, the high spatial fidelity of this disturbance (Patten and Knight 1994) implies that many thousands of hectares of DPFs are disturbed each year. The location of disturbed areas should be fairly constant from year to year.
The second most common disturbance agent was insects, in terms of both overall and yearly disturbed area. The period of analysis    of 3 years each (Baltensweiler et al 1977). Insect disturbance had a very low recurrence, as the insects moved to adjacent areas from year to year. However, the area disturbed could be overestimated, because the Forest Service recorded only the general perimeter of the affected stands, and no estimate of disturbance severity exists. Defoliators were an important forest disturbance agent in the study region; their influence on forest structure may vary depending on outbreak duration and severity, but in the southern Alps these are seldom standreplacing (Battisti 2008). Outbreaks react to climate change via the direct influence of temperature on insect reproduction and the indirect effects of the quantity and vigor of plant hosts (Esper et al 2007;Vanhanen et al 2007;Netherer and Schopf 2010;Haynes et al 2014), so that future dynamics are quite complex to predict. Wildfires and wind are usually regarded as the most severe disturbance agents, especially in DPFs, due to their potential to be stand-replacing and to increase exposure to other natural hazards in their aftermath (Angst and Volz 2002;, Bebi et al 2015. However, in our study the overall area disturbed by wind/snow and wildfires combined was lower than that for either avalanches or insects; moreover, even without a direct estimate of disturbance severity, we must assume that not all events were stand-replacing. While DPFs at risk from avalanches mostly protect roads (as people tend not to settle below avalanche-prone slopes), DPFs at risk from wildfires and wind protect buildings and towns. The higher vulnerability or value of the hazard target increases the risk posed by these disturbances even if they affect a smaller proportion of regional forests.
The nonconformity of the size distribution of the disturbances to a negative power law (Turcotte and Rundle 2002) suggests that large, infrequent disturbances (Turner 2010) are missing from the data set. However, the maximum area disturbed by wildfire was 400 ha, which is comparable to other heavily settled inner-Alpine regions with a similar fire regime (eg largest wildfire in Valais: 300 ha; Gimmi et al 2004) and is suggested to correspond to a 100-year return interval event (Moser et al 2010). These large fires occurred mostly after the year 2000. In the Alpine region, landscape fragmentation due to anthropogenic land use has prevented disturbance agents from spreading across larger areas, and disturbance suppression used to be put into practice for both wildfires and insect outbreaks (Faccoli et al 2010;Zumbrunnen et al 2012). During the study period, however, the combined effect of fuel buildup due to fire suppression, increased fuel continuity due to land abandonment (Piussi and Farrell 2000), and higher fire danger due to climate warming may have induced a shift toward an increased frequency of large fires (Castellnou and Miralles 2009;Valese et al 2014).
Conclusive evidence for disturbance regime shifts cannot be provided using our data set, which, despite being one of the longest available for multiagent disturbances in any European region (see also Hanewinkel et al 2008), has an insufficient time span to capture exceptional events with a return interval of 100 years or more. We think that this is why we found scarce evidence for disturbance recurrence-less than 0.3% of the disturbed area was disturbed by the same agent a second time in each of the subsequent years. (It should also be kept in mind that the literature lacks reference figures for "pristine" temperate mountain  ecosystems.) Z. griseana was an exception, with a recurrence rate of 21%, probably due to the multiyear nature of the outbreak cycles (Baltensweiler and Rubli 1999) and the coarse delineation of disturbed areas by the Forest Service. Interactions between different disturbances were significant for wind and snow and insects. However, while there is an established literature on interactions between bark beetles and wind (Peltonen 1999;Bouget and Duelli 2004;Wermelinger 2004;Hanewinkel et al 2008;Temperli et al 2013), we could find no plausible mechanism involving defoliators, so we assume that the 2 disturbances simply occurred preferentially on similar sites, such as monospecific, dense spruce stands (Bejer 1988;Peltola et al 1999). Finally, wind/snow and wildfire did not significantly interact with each other within the period of analysis. Contrary to common assumptions (eg Myers and van Lear 1988), our data suggest that, at the regional scale, wind and snow damage does not have the potential to increase burn area. Increased dead fuel loads might not be sufficient by themselves to increase fire danger, especially if they do not create a sustained increase (in quantity and time) of fine fuels (Cannon et al 2014). Moreover, natural canopy gaps reduce leaf litter abundance, decreasing fuel availability and continuity for subsequent fires (O'Brien et al 2008). All disturbances followed their expected gradients of topography (wildfires on drier sites and wind and snow on moister sites-Schindler et al 2009) and plant functional type (higher incidence of fire in broadleaf forests, which are closer to ignition sources and less subject to prolonged snow cover in winter, and of all other disturbances in conifer forests).
The analyzed variables may be used to parameterize past disturbance regimes for the forests of Aosta Valley. However, we feel that this is not yet possible for several reasons: Validation metrics made us quite confident in our current land cover classification; the size of Landsat pixels allowed us to detect unforested gaps of 0.09 ha, which are enough to reduce or cancel protection effectiveness (Schneebeli and Meyer-Grass 1993). However, we cannot rule out that such gaps were unforested even before disturbance, as disturbance perimeters were not sampled at such a fine spatial resolution. The effect of disturbances on residual protection function was higher for wildfire (42% of DPF land disturbed by avalanches was classified as currently nonforested) and avalanches (21%), which reflects the stand-replacing nature of such disturbance agents, and may be related to other effects on the ecosystem that delay forest regeneration, such as removal or sterilization of soil, postdisturbance erosion, increased evaporation, and distance from the seed source (Neary et al 1999;Ice et al 2004;Certini 2005;Vacchiano et al 2014). Wind and snow and insect disturbances, on the contrary, were less severe and left a smaller proportion of DPFs without sufficient cover, likely due to the beneficial effect of disturbance legacies and moderate severity effects (Peterson and Pickett 1990;Ulanova 2000;Kuuluvainen and Kalmari 2003;Waldron et al 2013;Seidl et al 2014a).
For all disturbances, the probability of forest recovery was strongly dependent on time elapsed since  (Viegas 2004). Recovery models had a low explanatory power, probably due to the noise inherent in the data set (measurement error, unknown pre-disturbance ecosystem state) and to confounding factors, such as disturbance suppression and salvage logging, which are routinely carried out in the region. However, model results provide useful indications of the combinations of disturbance agent, vegetation, and environmental conditions that are most likely to slow the recovery of a forest's protective function and therefore should be prioritized for assisted recovery of this important ecosystem service.

Conclusions
We provided a comprehensive overview of the recent history of different disturbance regimes and their interactions in forests of an Alpine region. The combination of this analysis with an approximation of DPF functionality, as has been done is this study, also provides insight into how DPFs are affected by different disturbance regimes and how they have recovered after disturbances. This paper demonstrates a simple and universally applicable algorithm for assessing the effectiveness of DPFs at a regional scale based on remote-sensing data. A similar analysis can be carried out in all parts of the Information on this variable was not available for this disturbance agent. *P , 0.05; **P , 0.01. world using institutional or open-access geographical data and publicly available Landsat imagery, and can be replicated at different spatial extents (eg the forest management unit) and points in time (eg immediately after a natural disturbance). Our method can provide land planners with a tool to prioritize the management of DPFs, by indicating where restoration of their functionality (whether by mitigation measures such as temporary hazard barriers or by active ecological restoration) is most urgent. This is ever more relevant under the forecasted increase in natural disturbance frequency and severity due to climate change in mountain areas of Europe.