Study area
The Vosges mountain range is located in the north-east of France (area: 8,800 km²). It has contrasting soil and climate conditions. Altitude ranges between 200 m in the eastern plains of Alsace to 1,424 m (Fig. 1A). Rainfall varies between 500 and 2,000 mm, and mean annual temperatures between 6°C and 11°C for the 1981-2010 period (Aurelhy model, Météo France). Geology mainly consists of acidic sandstones in the north and west of the study area, and magmatic rocks in the centre of the mountain range (Fig. 1B). Forests cover more than 60% of the territory, and provide important economic and social benefits². The region experienced important tree decline events during the 1916-1925, 1943-1951, or 1973-1981 periods. During this last period, combined climatic stress after the 1976 drought and air pollution that increased nutritional difficulties led to magnesium deficiencies (Landmann et al. 1987). Silver fir and Norway spruce have gradually been declining since 2016, questioning managers and policy makers about the future of these species in the Vosges massif.
Mapping of dieback using remote sensing data
We used three Sentinel-2 images level N2A downloaded from Theia website[1]. We selected images taken on July 24, 2019, during a period with little cloudiness and maximum vegetation cover to avoid confusion between the soil and dieback, and differentiate dead trees from living trees. These images were corrected for atmospheric, topographic and geometric effects. Dieback was identified through a supervised classification by excluding bare soil, crops, ground vegetation, urban areas, water, and forested areas. The classification was trained using 288 units presenting dieback, and validated with 1,814 others units including 285 units presenting dieback, collected by photo-interpretation. The Sentinel-2 images were classified with ArcGIS Pro 2.5.0, using maximum likelihood classification (MLC, Sun et al. (2013)).
Sampling strategy
To avoid considering mortality patterns from other species than silver fir and Norway spruce, we used the IGN BD forêt v2[1] database to select vegetation units recorded on aerial photographs that showed these two species. A systematic 100-m grid was built, and 1,000 plots were randomly selected in dieback-struck areas and another 1,000 plots in healthy places.
We used aerial photos at 20 cm resolution to separate plots according to the dominant species (spruce or fir) and to identify the stand composition on 15-m-radius plots (seven classes, Online resource 1). We used photos dating back to 2018, taken prior to the dieback event, to identify the stand characteristics. To avoid the risk of confusion, only the classes where the studied species was dominant or not in mixture with the other studied species were kept for each dataset to avoid confusion (classes 1, 2, 3 and 5 in Online resource 1; 872 plots were selected for silver fir, and 1,043 plots for Norway spruce). For each plot, the stand’s status (healthy or dieback) and 112 explanatory variables were collected in a database: 6 variables described the stand characteristics, 50 described the site conditions, and 56 described climate evolution.
Variables describing the stand characteristics
We distinguished four levels of species mixture (COMP, modalities 1, 2, 3 and 5, Online resource 1) within the stand characteristics variables (Table 1). DENS differentiated dense stands from sparsely forested stands, where the ground was visible on more than 10% of the plot area. To describe the age and structure effects, STRUCT was recorded in three classes: old even-aged forest stand, young even-aged forest stand, and heterogenous structure. TYP_EDGE identified the presence of edges due to interfaces between different land uses (absence of an edge, edge between 2 forest stands of different heights, and with an open environment), and EDGE_S identified edges facing south. These variables were collected using aerial photographs, on 15-m-radius plots (similar to classical forest inventory plots in France), for COMP, DENS, and STRUCT, and on 50-m-radius plots for TYP_EDGE and EDGE_S to consider data around the plot that could influence the stand health. Finally, we added data about the 1999 windstorm damage that had a strong impact on the Vosges mountains forests by rejuvenating some forest stands. IGN maps differentiated areas with less than 10% damage from those with 10 to 50% damage or more than 50% damage.
Variables describing the site conditions
Fifty variables described the site conditions. They concerned geology (n=1), topography (n=9), soil nutrition (n=2), waterlogging (n=2), energy for plants (n=12) and water availability (n=24) (table 2). A simplified geological map (GEOL) was drawn by merging similar substrates on the local maps at 1:50,000 scale elaborated by the BRGM. Topographical indices were calculated with the IGN 50-m DEM for altitude (ALT), aspect – cosine of aspect (COS_ASP), sine of aspect (SIN_ASP) –, and slope (SLO). Topographical position was represented by four indices: containment (CONT), the difference of altitude between the pixel and the average altitude within a 1,500-m radius, curvature (CURV) that describes the shape of the topographical surface in the direction of the maximum slope within a 1,500-m radius (convex, flat or concave), TOPO and REL_ALT that are the relative distance and altitude of the cell along the toposequence (the ratio of the distance or altitude between the cell and the nearest ridge to the total length or the difference of altitude found in the toposequence), respectively. We characterised soil nutrition and waterlogging according to 1-km resolution maps built using bioindication techniques (Gegout et al. 1998) for soil surface pH, nitrogen nutrition (C/N ratio, CN), temporary waterlogging (TW) and permanent waterlogging (PW).
GIS layers of monthly climate variables were used to estimate energy and water available for plants over the 2009-2019 period. Solar radiation (RAD) was calculated at 50 m resolution using Helios model (Piedallu and Gegout 2007). Monthly minimum, average and maximum temperatures and rainfall (TMIN, TMEAN, TMAX, and RAIN, respectively), were modelled and mapped at 1 km resolution (Piedallu et al. 2016). Climatic water balance (CWB) is the difference between RAIN and potential evapotranspiration according to Turc’s formula (Turc 1961). One-km resolution maps of soil depth (DEPTH) and soil water holding capacity (SWHC) (Piedallu et al. 2011) were used to evaluate the soil water content (SWC, available water in the soil) and the soil water deficit (SWD, the difference between potential evapotranspiration and actual evapotranspiration) according to Thornthwaite formula (Thornthwaite and Mather 1955). In addition, three indices characterising lateral fluxes in the soil were extracted from GIS layers. Kh and Kv represented the horizontal and vertical hydraulic conductivity at soil saturation, respectively, and FLUX estimated the balance between the input and output surface and subsurface fluxes (Ondo et al. 2017).
All the variables including climate were averaged over the ten years before data collection (2009-1019), per season (winter -wi, December to February-; spring -sp, March to May-; summer -su, June to August; and autumn -au, September to November), and over the whole year (an).
Variables describing climate evolution
Two sets of calculations – climate anomalies and climate trends – were used for TMIN, TMEAN, TMAX, RAIN, CWB, SWC and SWD. Climate anomalies represented the evolution between a reference period and a contemporary period. A period of ten years before data collection (2009-2019) was selected to represent the contemporary period because tree death can be triggered by the cumulative effects of many disturbances (Das et al. 2013). We selected 1961-1985 as the reference period because significant warming occurred in France after 1986 (Taccoen et al. 2019).
Climate trends were calculated for the same variables over the 1986-2019 period. For each variable and each pixel of 1-km resolution, the trends were regressed against time for the 32 studied years, and the linear regression coefficients were mapped. A positive slope indicated an increase of the climatic constraint for TMIN, TMEAN, TMAX and SWD, and a decrease for RAIN, CWB and SWC.
Modelling approach and vulnerability map calculation
We modelled silver fir and Norway spruce mortality with logistic regression (McCullagh and Nelder 1997). For each species, we evaluated the 113 explanatory variables through a forward stepwise procedure. At each step, the variable explaining the maximum deviance was selected if it was statistically significant (P < 0.01), and little correlated with the previously selected variables (Pearson’s r < 0.7). Because monotonic responses were expected, quadratic responses were not considered for the climate variables. However, monotonic and quadratic response curves were both evaluated for the other variables. The process ended when the selection of a new variable led to a d² increase lower than 0.01 or a p-value > 0.01.
The receiver operating characteristic (ROC) curve was calculated to evaluate the performances (Fielding and Bell 1997), and different metrics were used: the area under the ROC curve (AUC), the kappa index, success (the rate of good prediction), sensitivity (the rate of good prediction of presence), and specificity (the rate of good prediction of absence) (Manel et al. 2001). The relative importance was calculated for each variable using the drop contribution (Lehmann et al. 2002), i.e., the difference in explained deviance between the full model and a partial model when the variable was dropped. The ratio of the drop contribution of the variable over the sum of the drop contributions for all the variables gave the relative importance of each predictor.
The ROC curve also provided the optimal threshold for predicting the presence / absence of dieback. It was used to calculate the vulnerability map: upper values delineated the “high level of risk” class where the stands were predicted dead by the model, while lower values were divided in three classes with equal intervals of increasing risk. The vulnerability map was calculated using the GIS layers describing the environmental predictors selected in the model. It included four levels of risk, calculated from a type stand used across the whole area.
Validation dataset
A map of the volumes of harvested dead wood from exceptional sanitary cuts was provided by the ONF for public forests, which cover 2/3 of the forested area of the Vosges mountains. We selected areas where one of the studied species was dominant or not in mixture with the other studied species. For 7,355 management units identified for silver fir (11,218 for Norway spruce), the average volume of dead wood was compared for the different vulnerability levels. Our vulnerability maps were drawn using 2019 satellite imagery, but dieback may have occurred before or after this date in the most vulnerable areas. As a consequence, we considered both the volumes of dead trees recorded in 2019 and the sum of the volumes over the 2018-2020 period.
2 https://theia.cnes.fr/atdistrib/rocket/#/search?collection=SENTINEL2
3 https://inventaire-forestier.ign.fr/spip.php?article646