Monitoring montane-subalpine forest ecotone in the Pyrenees through sequential forest inventories and Landsat imagery
Annals of Forest Science volume 80, Article number: 32 (2023)
Abies alba Mill.–Pinus uncinata Ramond. ecotone dynamics are examined along both altitudinal and protection level gradients by combining field inventories and Landsat data. An upward expansion of A. alba to the subalpine belt is observed in the last decades as a result of stand maturation after logging cessation.
High-mountain forests constitute sensitive locations to monitor the impacts of global change on tree-species composition and ecotone dynamics. In this study, we focus on the Spanish Pyrenees where silver fir (Abies alba Mill.) coexists with mountain pine (Pinus uncinata Ramond.) forming montane-subalpine ecotones.
The main goal of this study is to assess the spatiotemporal dynamics of the silver fir–mountain pine ecotone and its underlying driving factors.
We reconstructed the spatial distribution and dynamics of the species by combining remote sensing imagery and field plot data from 1989 to 2015, employing support vector machine techniques for image classification. Using variance analysis and mixed effects models, we then analyzed the evolution of basal area and replacement index, a measure of relative change in species composition, over time and altitude range. Additionally, we explored their relationship with site factors and protection level (National Park vs. protection buffer zone).
Silver fir has expanded its distribution in both the National Park and the protection buffer zone, whereas mountain pine has remained stable. Both species exhibit increased basal area associated with stand maturation and a higher level of protection. The replacement index indicates a rise in silver fir in the understory on North-facing slopes, attributed to stand densification. These findings are particularly noticeable in the area with the highest level of protection.
The cessation of traditional land uses has led to ongoing stand densification, promoting succession and favoring the increased abundance of silver fir at its uppermost locations, where this species outcompetes mountain pine.
The legacy of past forest management together with climate change is leading to shifts in forest composition, modifying forest structure, dynamics, and related ecosystem services all over the world (Lu et al. 2020). Climate warming might induce species range shifts in mountain ecosystems by alleviating the constraints of cold on growth and regeneration at the upper limit (Lenoir and Svenning 2014) or by inducing shifts in forest composition (Gehrig-Fasel et al. 2007; Fadrique et al. 2018). The cessation of forestry exploitation activities such as logging and traditional land uses such as firewood extraction or grazing have triggered an increase in forest cover and density, as well as stand maturation in European mountains (Gehrig-Fasel et al. 2007; Camarero et al. 2011; Bodin et al. 2013; Kouba et al. 2012), favoring a transition to more shade-tolerant species (Knott et al. 2019).
Ecotones like the alpine treeline constitute sensitive locations to monitor the effects of global change on forest ecosystems. The spatiotemporal changes in ecotones can reflect forest dynamics and the response to drivers of global change, alleviating or aggravating the consequences at local to regional scales. There is a considerable body of research dedicated to disentangling the impact of global change on treeline dynamics (e.g. Batllori et al. 2009; Lu et al. 2020; Bader et al. 2021), but few studies have focused on subalpine-montane forest ecotones located at lower altitudes (but see Foster and D’Amato 2015).
Mediterranean basin mountain forests such as those found in the Spanish Pyrenees have been historically subjected to important human demographic fluctuations and related changes in logging pressure, as well as varying climate conditions. As a result, Iberian mountains represent sensitive areas for monitoring the effects of global change on the distribution and structure of mountain forests (Ameztegui et al. 2016). The reduction in human pressure due to the massive migration of the rural population to the cities in the 1960s, the increase of protected areas, and the steady rise in temperatures since the 1970s (Macias et al. 2006) have modified forest dynamics in the Pyrenees over recent decades (Serra-Maluquer et al. 2021). This has led to forest densification and expansion (Gracia et al. 2011; Ameztegui et al. 2021), a decline in certain tree species in low-altitude populations under warmer-drier conditions (Camarero et al. 2011; Hernández et al. 2019), altitudinal upward shifts (Sangüesa-Barreda et al. 2020; Hernández et al. 2014), or tree species replacement (Ameztegui and Coll 2011). Silver fir (Abies alba Mill.) reaches its South-Western margins of distribution in the Spanish Pyrenees. Although a general decline in low-altitude mixed forests of silver fir has been reported in areas affected by drought-induced dieback in the Central and Western Pyrenees (Camarero et al. 2011; Sangüesa-Barreda et al. 2015), this species seems to have expanded into the subalpine Pyrenean vegetation belts and towards the East, where dieback is rare (Ameztegui and Coll 2011; Ameztegui et al. 2015; Hernández et al. 2019).
Remote sensing data provide valuable sources for tracking forest dynamics over time (Gómez et al. 2012; Hermosilla et al. 2015). By applying multitemporal or time-series approaches to remote sensing imagery, trends and discrete changes in land cover (Gómez et al. 2016; Hermosilla et al. 2018; White et al. 2018; Wulder et al. 2020) and species dynamics (e.g. Ameztegui et al. 2016; Sánchez de Dios et al. 2020) can be described. The Landsat program in particular has provided a global record for monitoring forest dynamics since 1972 (Wulder et al. 2012). This temporal range makes Landsat archives an essential source of data for the characterization of forests over time, which combined with geo-referenced data from periodic field inventories provide information on changes in forest structure and dynamics (White et al. 2016). However, the identification of species continues to be an issue, especially in coniferous forests (Rautiainen et al. 2018).
Although within the Pyrenean area several studies have been carried out at a larger scale, these were mainly focused on dieback processes in montane Pinus sylvestris L.–Fagus sylvatica L.–Abies alba Mill. (Scots pine-European beech-silver fir) mixed forests. There are no studies focused on tree species dynamics at Pyrenean montane-subalpine ecotones integrating field data and satellite imagery. The “Aigüestortes i Estany de Sant Maurici” National Park and surrounding buffer area with different level of protection provides a unique scenery to study silver fir–mountain pine (Pinus uncinata Ramond.) dynamics at the montane-subalpine ecotone combining historic and intensive forest monitoring information with Landsat data. For this, we apply a two-scale approach during the period 1989–2015. Integrating field data and remotely sensed information with machine learning techniques, we analyze three main indicators: (i) changes in silver fir and mountain pine basal area, (ii) changes in a specific replacement index which assesses the composition differences between the regeneration strata and the dominant strata, and (iii) changes in the altitudinal distribution of the populations of both species, paying special attention to the underlying causes that may be fostering these dynamics.
2 Material and methods
2.1 Study area
Our study focuses on the “Aigüestortes i Estany de Sant Maurici National Park” and its protection buffer zone, situated in the Eastern Spanish Pyrenees (Fig. 1). The area includes different levels of protection: logging has been excluded in the National Park since 1993 whereas forest management limitations were established in the buffer zone (Memoria. Parc Nacional d’Aigüestortes i Estany de Sant Maurici 2018). The study area covers an extent of 40,852 ha (National Park: 14,119 ha, buffer zone: 26,733 ha). The altitude of the study area ranges from 1000 m to 3033 m asl with valleys of different orientations. N-W valleys have an Oceanic influence and S-E valleys are more influenced by Mediterranean continental climatic conditions, which determine the distribution of the different forest ecosystems. The average annual precipitation is 1500 mm per year, and the average temperature is 4.7 °C at 1332 m with monthly temperatures dropping below 0 °C from November to February. The soils are predominantly acidic, mostly developed on granite and slate substrates (Memoria. Parc Nacional d’Aigüestortes i Estany de Sant Maurici 2018).
The montane belt (located at about 1000–1700 m), comprises European beech and birch (Betula sp.) stands, other broad-leaved mixed forests, mixed forests of coniferous and broad-leaved species, and forests of Scots pine (mainly in low altitude and South-facing slopes) or silver fir. Upwards, silver fir and mountain pine coexist to form the montane-subalpine (1700–2100 m) ecotone, while the subalpine belt is almost entirely dominated by mountain pine. Silver fir and mountain pine are the only species forming high mountain forests in this area. Silver fir covers a smaller area than mountain pine, occupying the moister, cooler, shaded slopes or lower valleys. Mountain pine, in turn, covers a large area ranging from cooler, shaded slopes to sunny, windy mountain ridges. Mountain pine forms the alpine treeline at an altitude of 2300–2400 m.
Mountain pine is a pioneer, shade-intolerant species which can reach 800–1000 years in the study area (Camarero et al. 2016), whereas silver fir is more shade-tolerant (Dobrowolska et al. 2017) and can reach ages of 500–600 years (Camarero et al. 2016).
2.2 Data collection
2.2.1 Field data
Since 1961, a systematic National Forest Inventory (NFI) has been implemented in Spain over the entire national territory (Alberdi et al. 2017). Each inventory comprises a cycle of approximately 10 years and measurements in every plot have been replicated in each cycle since the second NFI (NFI2), after 1989. The fourth NFI (NFI4) was conducted in the study area in 2015 (MITECO 2021). Additionally, the National Parks Autonomous Agency (OAPN) has carried out three independent field inventories over the last three decades (NPFI), measuring plots in the study area (Fig. 1 and Table 1).
Field plots of both inventories, NFI and NPFI, consist of four concentric circular areas with radii of 5, 10, 15, and 25 m, where trees are measured depending on the diameter at breast height (dbh). In the 5-m, 10-m, 15-m, and 25-m radius plots, trees with dbh ≥ 7.5 cm, dbh ≥ 12.5 cm, dbh ≥ 22.5 cm, and dbh ≥ 42.5 cm are measured, respectively. Trees are identified to the species level, and their dbh and height are measured; thereafter, the basal area (BA) is estimated using a scale factor that varies according to the plot radius (Alberdi et al. 2016, 2017). Therefore, both NFI and NPFI plots have the same design, making them comparable and complementary.
For this study, we used three cycles of NFI data (NFI2, NFI3 and NFI4) and three cycles of NPFI data (NPFI1, NPFI2, NPFI3), dating from 1989 to 2015. We merged NFI and NPFI data within a time-lapse of 5 years over each of the three reference years: 1991 for NFI2 and NPFI1, 2002 for NFI3 and NPFI2, and 2015 for NFI4 and NPFI3 (Table 1). We created two databases: one considering only permanent plots (plots remeasured in the three reference years) and the other considering all measured plots (permanent and temporary plots).
To analyze silver fir–mountain pine ecotone dynamics we classified tree species capitalizing on information on species presence/absence from every NFI and NPFI 25 m radius plot. To compare basal area and species replacement index (RE) between 1991, 2002, and 2015 and to analyze the underlying drivers of change, we included only the permanent plots. Although permanent plots are placed in the same location across successive inventories, the different plot radii regarding the tree dbh may result in changes in the sampling area for the same tree. This prevents the analysis of the rate of change by period analysis at a plot level; therefore, BA and RE were analyzed for the whole sampling area (Fig. 1).
2.2.2 Remotely sensed data
The Landsat scene identified by path/row 199/030 covers the entire study area with a pixel resolution of 30 m. Landsat data have been proven to be a powerful tool for forest species classification thanks to its temporal and spatial resolution and have been applied in diverse studies all over the globe (e.g. Meng et al. 2020; Thompson et al. 2015). Landsat holdings in the United States Geological Survey (USGS) archive were explored, and the availability of Landsat collection 1 level 2 surface reflectance products was analyzed. Landsat collection 1 level 2 data are radiometrically, geometrically, and atmospherically corrected, and the terrain distortions are rectified. Epochal or multi-year image composites with target date 15 July were built (White et al. 2014) averaging pixel values from the reference or target year (1991, 2002, and 2015) and pixel values from the nearest previous and posterior years with available data (Table 2). Only images with cloud cover ≤ 70% were considered. Landsat 7-ETM + SLC-off images were also disregarded. By employing multi-year image composites, inter-annual variation effects on vegetation reflectance and subsequent classification errors are minimized. A forest mask was built from the Spanish Forest Map 1:25.000 (Vallejo-Bombín 2005) and applied to the epochal composites. All image processing was implemented using ProLand, a software package developed in Matlab® by the ICIFOR—INIA, (CSIC) team.
Both field inventory data and remote sensing-derived data are freely available at CSIC digital repository (Aulló-Maestro et al. 2023).
2.3 Metrics derived from field data
2.3.1 Species basal area analysis
To better understand forest ecotone dynamics, we compared the BA of mountain pine and silver fir by altitude intervals over the three periods employing data from the permanent NFI and NPFI plots. In addition, we studied the differences in BA changes between the National Park and the buffer zone area. As our data do not fulfill the assumptions for normal distribution, we applied a Friedman test and the post hoc Wilcoxon test for non-parametric repeated measures analysis (R Core Team 2020; Kassambara 2021) to assess differences in BA between years using paired plots (Table 1) grouped by altitudinal intervals of 200 m. Therefore, only those plots measured in the three referenced years (permanent plots) were considered (Table 3).
2.3.2 Evaluation of the species replacement index (RE)
In order to characterize the dynamics between the two tree species at the stand level, we calculated the species replacement index (RE) defined in Ledo et al. (2009) by using the dbh of all trees across the diameter range from 7.5 cm to the dominant trees in each permanent plot (Table 1).
The RE is based on the r12(Δdbh) function (Ledo et al. 2009) which summarizes, by dbh-difference intervals, the number of tree pairs where the smaller tree belongs to species 2 (mountain pine) minus the number of tree pairs where the smaller tree belongs to species 1 (silver fir):
Δdbh is the between-trees dbh difference lags. For each pair of individuals i,j, Iij(Δdbh -δ, Δdbh + δ) is 1 if the absolute difference between the dbh of trees i and j is within the interval [Δdbh -δ, Δdbh + δ)] and 0 otherwise. f(i,j) is 1 if the smaller tree between i and j belongs to species 2 (mountain pine) and is -1 if the smaller tree between i and j belongs to species 1 (silver fir).
The RE is the weighted sum of the accumulated r12(Δdbh) values for all the dbh difference intervals, being the weight for each interval the dbh difference for that interval. In this way, those pairs of trees with greater dbh differences (juvenile vs. dominant trees) have more weight in the RE index:
where max(Δdbh) is the dbh difference between the largest and the smallest trees in the plot and min(Δdbh) is the dbh difference between the most similarly sized trees.
Positive values of the RE index indicate a greater proportion of mountain pine in juvenile classes, 0 indicates coexistence between both species, and negative values indicate a greater proportion of silver fir in juvenile classes. The replacement index was calculated using a software package developed in Matlab® by the ICIFOR-INIA, (CSIC) team.
As in BA evaluation, we applied a Friedman test and the post hoc Wilcoxon test for non-parametric repeated measures to assess changes in RE over time (1991, 2002, and 2015) within 200 m altitudinal classes and within National Park and buffer zone. Only permanent plots were considered (Table 1).
2.4 Statistical analysis
2.4.1 Mixed-effects modeling
To evaluate the site and temporal effects over total basal area (BAT) and RE index, we applied a linear mixed-effect model with a normal distribution of errors. For modeling, we studied the inclusion of altitude, aspect, level of protection (National Park vs. protection buffer zone), reference year (categorical) and all possible interactions between them as fixed effects, and plot as random effect. For RE modeling, we considered as well BAT as fixed effect. To identify the best-supported model we analyzed the Akaike information criterion (AIC) within all the possible models using the maximum likelihood (ML) fitting. The final models were fitted with restricted maximum likelihood (REML) (Zuur et al. 2009), and the normality of residuals was checked with DHARMa’s (Hartig 2022) Kolmogorov–Smirnov test. For the generation of the linear-mixed effects, we used MuMIn R package for AIC analysis (Barton 2022) and lme4 for the model fitting (Bates et al. 2015).
2.4.2 Tree species classification
Support vector machine (SVM) binary learners were applied to retrieve the presence/absence areas of both mountain pine and silver fir species at reference years (1991, 2002, and 2015), deriving six models in total. Predictor variables included values of reflectance from the six Landsat optical bands (0.45–0.52 μm, 0.52–0.60 μm, 0.63–0.69 μm, 0.76–0.90 μm, 1.55–1.75 μm, 2.09–2.35 μm). The SVM classifier was chosen as it is specifically designed for one-class binary classification, being a suitable classifier for single-species distribution. Moreover, SVM achieves higher overall accuracy, and it has been proven suitable for classifying forest attributes when using remote sensing imagery (Khatami et al. 2016).
We used a presence/absence indicator derived from field inventory data as the dependent variable in each model. This indicator was determined by the species' relative basal area as a percentage of BAT. For the mountain pine indicator variable, a value of 1 (presence) was assigned to plots where mountain pine’s relative basal area exceeded a cutoff threshold, and 0 (absence) otherwise. Similarly, the silver fir indicator variable was assigned values of 1 or 0 based on its relative basal area. To identify the optimal relative basal area threshold for the SVM classifier, we constructed a receiver operating characteristic curve (ROC). This curve shows the relationship between the sensitivity, measured by the true positive ratio (TPR), and 1 – specificity, represented by the false positive ratio (FPR). We trained the SVM classifier for 11 different relative basal area thresholds for each species, ranging from 0 to 100% with increments of 10%. For each threshold, we calculated the mean sensitivity, specificity, and total accuracy of the two species classifications. The threshold selected for the presence/absence indicator in the final SVM model was determined by the optimal cut-off point of the mean ROC curve.
NFI and NPFI plots corresponding to each reference year served as training data (Table 1). Reflectance values of the Landsat epochal composites obtained for the reference years (1991, 2002, 2015) were extracted from 3 × 3 pixel windows centered on each plot (that is, covering the 25-m radius plot). This resulted in a training sample of 1557, 1701, and 1395 pixels for the respective reference years 1991, 2002, and 2015.
SVM-trained classification models were employed to predict the presence of the species above the chosen threshold across the entire study area. Reflectance values served as predictors for each pixel, specifically focusing on forested areas according to the Spanish Forest Map (Vallejo-Bombín 2005). To optimize the SVM’s Kernel scale parameter, we utilized the sequential minimal optimization algorithm (Fan et al. 2005) through the fitcsvm function in Matlab® with 50 iterations. The cost parameter utilized default values from the fitcsvm function (1 for all binary learner classifications). Subsequently, using the resulting presence/absence maps, we examined the influence of altitude and aspect on species expansion, retraction, and stable areas.
The accuracy of the SVM results was first evaluated through confusion matrices obtained from 10-fold cross-validation. Omission-commission errors and total accuracy were derived from each confusion matrix. Additionally, we evaluated performance using the area under the curve for the receiver-operator characteristic curve (AUCROC). The receiver-operator characteristic curve was constructed by applying thresholds to the class scores, following a common approach in machine learning methods (Valavi et al. 2022). Furthermore, we utilized the Mathews correlation coefficient (MCC, Matthews 1975) to evaluate the correlation between observed and predicted values. MCC was normalized within the range of 0 to 1 and provides accurate performance assessment for binary classifications (Chicco & Jurman 2023).
3.1 Changes in species basal area
The post hoc analysis revealed a progressive increase in the basal area (BA) of both mountain pine and silver fir within the study area’s forests from 1989 to 2015.
Despite no significant difference evidenced in silver fir BA in plots located at lower altitudes (1400–1600 m), its BA significantly increased at higher altitudes (1600–2000 m) from 1991 to 2015, whereas stability was detected in the altitude class 2000–2200 m.
Similarly, mountain pine BA showed a non-significant difference at its lower altitude (< 1800 m), but a significant increasing trend was observed in the 1800–2000 m altitude class from 1991 to 2015 and 2000–2200 m class from 1991 to 2002 (Fig. 2 and Table 4).
The Friedman and post hoc analysis considering levels of protection (National Park vs. buffer zone) revealed that the buffer zone has only shown BA increments in 2002–2015 for silver fir. More intense changes started even before within the National Park, where we observed that both silver fir and mountain pine displayed a significant increase in BA since 1989 (Table 4).
3.2 Changes in the species replacement index
At the lowermost altitudes (1400–1600 m), the null RE index value indicates stability between both species since 2002. Negative values for altitudes between 1600 and 1800 m reflect a predominance of silver fir in the understory at these altitudes. No significant changes in the RE index were observed at altitudes from 1400 to 1800 m. According to Friedman’s test results, between 1800 and 2000 m RE values were positive, meaning that mountain pine was the predominant species in the juvenile strata. However, these positive RE values significantly decreased during the last reference year (even becoming negative), denoting a slight increase in silver fir in the understory at the expense of mountain pine. Above 2000 m, positive RE values denoted a clear predominance of mountain pine in the understory, but a significant drop was observed within the entire studied period (1991 – 2015), especially during the second period (2002–2015), indicating an increase of silver fir (Fig. 3 and Table 4).
The analysis of changes between reference years in the National Park, where logging cessations had taken place earlier, and the protection buffer zone, reveals that the buffer zone has not been subject to significant changes throughout the study period. In contrast, we observed differences in the National Park area over the entire studied period (1991–2015), especially during the second period (2002–2015) (Table 4).
3.3 Factors influencing species abundance and replacement
The best-supported BAT model in terms of AIC included altitude, aspect, level of protection, and year of reference as fixed effects and year of reference × altitude interaction, and plot as random effect. On the other hand, for RE, the best-supported model included altitude, aspect, level of protection, year of reference, and BAT as fixed effects and BAT × altitude interaction and plot as random effect. The difference between marginal and conditional R2 suggests that plot effect explains part of the variance linked to factors not included in the model such as climatic or soil conditions. DHARMa test results show normality of residuals for both models (Kolmogorov–Smirnov p-value = 0.18 for BAT and 0.15 for RE) (Table 5).
We found that BAT is mainly affected by the level of protection. BAT increases in those plots inside the National Park area and slightly in South-facing slopes. Regarding temporal effects, BAT increased in both periods, but this effect was greater in the second period (2002–2015). In RE, all variables included in the selected model were significant except the first period (1991–2002) and the level of protection. Altitude and BAT × altitude interaction got the greatest significance. As negative values of RE indicate a predominance of silver fir in juvenile classes, a negative correlation in models reflects a positive effect in silver fir abundance in juvenile classes and a positive correlation indicates a positive effect in mountain pine abundance in juvenile classes. Silver fir recruitment increased in dense areas and on North-facing slopes while mountain pine recruitment is associated with lower BAT and higher altitudes. The BAT × altitude interaction reflects changes in the effect of stand density on RE within altitude due to the different altitudinal ranges of both species, RE being less sensible to BAT fluctuations as altitude increases.
3.4 Changes in tree species distribution from SVM classification
The presence/absence indicator was obtained for each species by applying the 40% threshold of relative basal area. This threshold was identified as the optimal cut-off point derived from the ROC curve analysis (the point of the ROC curve which maximizes TPR – FPR difference) with a mean global accuracy of 77.52%. The silver fir and mountain pine indicators were calculated for the NFI and NPFI as training data for SVM classifiers, using Landsat epochal mean composite reflectance as predictors. The Kernel scale parameters of the fitted SVM models were 0.401, 0.556, and 0.990 for the reference years 1991, 2002, and 2015, respectively, for silver fir, and 1.595, 0.889, and 2.324 for the reference years 1991, 2002, and 2015, for mountain pine.
Silver fir classification showed accuracies of 83.0% (1991), 82.1% (2002), and 79.5% (2015). Mountain pine SVM classification models had accuracies of 72.1% (1991), 75.4% (2002), and 76.7% (2015) (Fig. 4). Silver fir SVMs showed greater specificity, whereas the mountain pine SVMs showed greater sensitivity. MCC ranged from 0.69 for silver fir in 2015 (AUCROC = 0.80) to 0.74 for mountain pine in 2002 (AUCROC = 0.82), indicating the reliability of the SVM classifications (Table 6).
Analyzing changes in the predicted area we observed differences between silver fir and mountain pine. Silver fir registered greater changes within the National Park area, while in mountain pine populations, changes are balanced between both areas. From 1991 to 2015, the area where silver fir BA is > 40% has shown an increment of 20.14% in extent compared to the first reference year (1991), with most changes occurring from 2002 to 2015. In contrast, the area where mountain pine BA is > 40% has experienced a minimal reduction of 0.68% from 1991 to 2015 (Table 7, Fig. 5). Regarding the level of protection, silver fir registered greater changes within the National Park area. From 1991 to 2002, silver fir predicted area experienced a slight reduction in the National Park area and an increase in the buffer zone; however, this tendency reverts in the second period (2002–2015) where there is a great increase in the area occupied within the protected area, mainly at 1800–2000 on North-facing slopes. The changes in mountain pine populations are evenly distributed between the National Park and the buffer zone. Losses in mountain pine distribution occurred at lower altitudes (1200–1800 m) and on North-facing slopes, while an increase in its predicted distribution area was observed at higher altitudes (1800–2400 m) on NE–NW facing slopes (Table 7, Figs. 6 and 7).
This work analyzed the change in basal area, tree replacement dynamics, and species spatial distribution ranges in the montane-subalpine ecotone of the Eastern Spanish Pyrenees. The use of different sources of forest information enabled the analysis of species dynamics and a discussion about the underlying drivers at multiple spatial scales.
The complementarity of NFI and NPFI forest inventories’ spatial distribution and the results derived regarding the succession in the composition of forest along the altitudinal gradient supports the observations based on Landsat image classification. These findings highlight the complementarity of remote sensing and field plot data for modeling species distribution (Pasquarella et al. 2018; Sánchez de Dios et al. 2020) and monitoring forest dynamics (Rautiainen et al. 2018). Remote sensing facilitates the assessment of canopy changes, and field data provides information about changes in forest structure and understory composition. The average epochal composites enable increasing the stability of species classification smoothing the effects of inter-annual climatic variations on vegetation reflectance (Hermosilla et al. 2015). The temporal continuity, free accessibility, and spatial resolution of Landsat imagery make it a valuable source of data for the characterization of forest dynamics (White et al. 2016; Gómez et al. 2019). Thus, changes in the stand dynamics may be monitored by remote sensing analysis and can be used as early-warning signals of impending changes in forest composition.
The confusion matrices show differences in sensitivity and specificity of the SVM classification of both species (Fig. 4), underestimating silver fir presence and overestimating mountain pine. Since a similar pattern affects the three reference years, the species distribution would be consistent in the studied period and would not affect the primary objective of our analysis, which is the assessment of change. The underestimation of silver fir may hide a greater expansion since some areas where this species increases over 40% of BA may not be represented. However, the bias in the classifications should be considered in the case of species distribution modeling. SVM performance could improve with a weighted-by-species abundance approach (Valavi et al. 2022). Although we tested this option (results not shown), weighting did not improve our classification performance.
During the first period analyzed (1991–2002), our findings show a slight increase in silver fir basal area and predicted distribution area in all the altitudinal ranges of the species. Regarding mountain pine, there was an upward expansion of the predicted distribution area and a basal area increase above 1800 m, remaining stable at lower altitudes. RE index reflected moderate changes in stand dynamics in this period. During the second period (2002 to 2015), there was a notable expansion of the silver fir predicted distribution area linked to a BA increase. This increase was more pronounced on 1800–2000 m − North-facing slopes within the area of highest protection level. At the same time, mountain pine experimented a range contraction driven by an upward shift of its trailing edge. RE index decreased at 1600–2400 m indicating an increase of silver fir in the understory in this period. These changes give rise to an increase in areas where both species coexist between 1800 and 2000 m.
From 1993, logging within the National Park was completely excluded, while in the buffer zone, only certain limitations for forestry were established (Memoria. Parc Nacional d’Aigüestortes i Estany de Sant Maurici 2018). Our results show the ecological succession of montane-subalpine ecotones in the study area over the analyzed period, characterized by a gradual and slow replacement of mountain pine by silver fir. This replacement is more pronounced within the National Park, especially at altitudes above 1800 m. Management legacies of historical logging may have limited silver fir expansion in the past (Camarero et al. 2011; Camarero 2017). The reduction of the logging pressure throughout the area since the implementation of the 1993 and 2003 National Park’s management plan has favored current forest densification, especially inside the National Park area where these plans were more restrictive. This circumstance benefited the regeneration of silver fir, a shade-tolerant species well adapted to lower altitudes and mesic conditions, thus increasing the competitiveness of fir (Ameztegui and Coll 2011; Hernández et al. 2019) and increasing its dominance in the silver fir–mountain pine ecotone, especially in North-facing slopes. However, these results contrast with the retraction reported in mixed silver fir, beech, and Scot’s pine forests in the Western and Central Spanish Pyrenees at lower altitudinal ranges (Sangüesa-Barreda et al. 2015). Denser mixed forests of silver fir and beech and Scots pine at lower altitudes endure drought and the spread of pathogens such as Heterobasidion spp (Oliva and Colinas 2010). This, together with changes in the historical management (logging), has led to the decline and dieback of some montane, low-altitude silver fir forests in the more xeric areas (Oliva and Colinas 2007; Camarero et al. 2011). Therefore, conditions between Western-Central Pyrenees, where silver fir populations reach their lowermost altitudes and encounter more xeric conditions, differ with Eastern silver fir forests at higher altitudes. Protecting mountain pine forests by precluding silvicultural practices in Eastern Pyrenees may have increased stand density, facilitating the gradual establishment of silver fir in the understory due to the greater tolerance to competition of this species (Dobrowolska et al. 2017).
Given the relationship between altitude and temperature, the shift of both species’ distribution peaks towards higher altitudes may also indicate the effect of a rise in temperatures as a consequence of climate change on mountain species dynamics. Although these changes in species distribution are based on SVM predictions and, therefore, should be interpreted with caution, our findings provide an example of the contemporaneous redistribution of biota, and in particular of tree species as a plausible consequence of accelerating global change drivers such as climate warming and land-use modifications over recent decades (Batllori et al. 2020). Furthermore, the results support the predicted upward range shift of silver fir in the Pyrenees towards subalpine forests and the contraction of the low-altitude mountain pine range under forecasted warmer conditions and low-intensity management scenarios (Benito Garzón et al. 2008). The shift of silver fir towards North-facing slopes of the subalpine belt in search of suitable climatic conditions may be favored by the densification of montane pinewoods. This is compatible with the expected increase in climate sensitivity of lower-altitude montane silver fir populations (Hernández et al. 2019), although this upward shift may be hindered at higher altitudes by the shallow soils and summit conditions of high mountains (Camarero and Gutiérrez 2017).
Our study highlights the influence of historical land-use policies on forest dynamics. The abandonment of traditional land uses has led to denser forests, resulting in increased BA for both silver fir and mountain pine. This shift has facilitated the presence and regeneration of silver fir, particularly at mountain pine lower limits. In the altitude range of 1600–1800 m, where both species coexist, we observed an increase in the abundance of silver fir in the juvenile strata, while mountain pine populations were observed to be retracting. These dynamics are more intense within the National Park, where logging has been completely excluded since 1993, compared to the protection buffer zone, only subjected to forest management regulation.
The presence/absence classifications derived from Landsat epochal mean composites revealed an expansion and a dominance increase of both silver fir and mountain pine towards higher altitudes, indicating a range shift in the montane-subalpine ecotone and active successional dynamics. These dynamics became more evident during the period 2002–2015 compared to 1991–2001, reflecting the extended time required for the response of long-lived tree species to land use or climate changes to become apparent.
The integration of remote sensing techniques and forest field inventory data has proven to be a complementary and effective approach for studying forest dynamics at multiple spatial scales within the context of global change.
Availability of data and materials
The datasets generated and analyzed during the current study are available in the CSIC repository, https://doi.org/10.20350/digitalCSIC/15332.
Alberdi I, Hernández L, Condés S, Vallejo R, Cañellas I (2016) Spain. In: Vidal C, Alberdi I, Hernández Mateo L, Redmond J (eds) National Forest Inventories. Springer, Cham
Alberdi I, Vallejo R, Álvarez-González JG, Condés S, González-Ferreiro E, Guerrero S, Hernández L, Martínez-Jauregui M, Montes F, Oliveira N, Pasalodos-Tato M, Robla M, Ruiz-González AD, Sánchez-González M, Sandoval V, San Miguel A, Sixto H, Cañellas I (2017) The multi-objective Spanish National Forest Inventory. For Syst 26:e04S
Ameztegui A, Coll L, Brotons L, Ninot J (2016) Land-use legacies rather than climate change are driving the recent upward shift of the mountain tree line in the Pyrenees. Glob Ecol Biogeogr 25:263–273
Ameztegui A, Coll L (2011) Tree dynamics and co-existence in the montane-sub-alpine ecotone: the role of different light-induced strategies. J Veg Sci 22:1049–1061
Ameztegui A, Coll L, Messier C (2015) Modelling the effect of climate-induced changes in recruitment and juvenile growth on mixed-forest dynamics: the case of montane-subalpine Pyrenean ecotones. Ecol Modell 313:84–93
Ameztegui A, Morán-Ordóñez A, Márquez A, Blázquez-Casado Á, Pla M, Villero D, García MB, Errea MP, Coll L (2021) Forest expansion in mountain protected areas: trends and consequences for the landscape. Landsc Urban Plan 216:104240
Aulló-Maestro I, Gómez C, Hernández L, Camarero JJ, Sánchez-González M, Cañellas I, Vázquez De La Cueva A, Montes F (2023) Monitoring montane-subalpine forest ecotone in the Pyrenees through sequential forest inventories and landsat imagery - DATASET [Dataset]; DIGITAL.CSIC; https://doi.org/10.20350/digitalCSIC/15332
Bader MY, Llambí LD, Case BS et al (2021) A global framework for linking alpine-treeline ecotone patterns to underlying processes. Ecography (cop) 44:265–292
Barton K (2022) MuMIn: multi-model inference
Bates D, Maechler M, Bolker B, Walker S (2015) Fitting linear mixed-effects models using lme4. J Stat Softw 67(1):1–48
Batllori E, Blanco-Moreno J, Ninot J, Gutiérrez E, Carrillo E (2009) Vegetation patterns at the alpine treeline ecotone: the influence of tree cover on abrupt change in species composition of alpine communities. J Veg Sci 20:814–825
Batllori E, Lloret F, Aakala T, Lloret F, Aakala T, Anderegg WRL, Aynekulu E, Bendixsen DP, Bentouati A, Bigler C, Burk CJ, Camarero JJ, Colangelo M, Coop JD, Fensham R, Floyd ML, Galiano L, Ganey JL, Gonzalez P, Jacobsen AL, Kane JM, Kitzberger T, Linares JC, Marchetti SB, Matusick G, Michaelian M, Navarro-Cerrillo RM, Pratt RB, Redmond MD, Rigling A, Ripullone F, Sangüesa-Barreda G, Sasal Y, Saura-Mas S, Suarez ML, Veblen TT, Vilà-Cabrera A, Vincke C, Zeeman B (2020) Forest and woodland replacement patterns following drought-related mortality. Proc Natl Acad Sci U S A 117:29720–29729
Benito Garzón M, Sánchez De Dios R, Sainz Ollero H (2008) Effects of climate change on the distribution of Iberian tree species. Appl Veg Sci 11:169–178
Bodin J, Badeau V, Bruno E, Cluzeau C, Moisselin JM, Walther GR, Dupouey JL (2013) Shifts of forest species along an elevational gradient in Southeast France: climate change or stand maturation? J Veg Sci 24:269–283
Camarero JJ (2017) The multiple factors explaining decline in mountain forests: historical logging and warming-related drought stress is causing silver-fir dieback in the Aragón Pyrenees. In: Catalan J, Ninot JM, Aniz MM (eds) High Mountain Conservation in a Changing World. Springer International Publishing, Cham, pp 131–154
Camarero JJ, Gutiérrez E (2017) Wood density of silver fir reflects drought and cold stress across climatic and biogeographic gradients. Dendrochronologia 45:101–112
Camarero J, Bigler C, Linares J, Gil-Pelegrín E (2011) Synergistic effects of past historical logging and drought on the decline of Pyrenean silver fir forests. For Ecol Manage 262:759–769
Camarero J, Gutiérrez E, Sangüesa-Barreda G, Galván D (2016) El Parque Nacional de Aigüestortes i Estany de Sant Maurici como santuariode bosques y árboles viejos de pino negro. In: En: La investigació al Parc Nacional d’Aigüestortes i Estany de Sant Maurici, X Jornades sobre Recerca al Parc Nacional d’Aigüestortes i Estany de Sant Maurici, Espot (Pallars Sobirà, Catalunya). p. 89–95.
Caudullo G, Welk E, San-Miguel-Ayanz J (2019) Chorological data for the main European woody species, DATA IN BRIEF, 2017, ISSN 2352–3409, 12, p. 662–666, JRC104492
Chicco D, Jurman G (2023) The Matthews correlation coefficient (MCC) should replace the ROC AUC as the standard metric for assessing binary classification. BioData Mining 16:4
Dobrowolska D, Bončina A, Klumpp R (2017) Ecology and silviculture of silver fir (Abiesalba Mill.): a review. J for Res 22:326–335
Fadrique B, Báez S, Duque Á, Malizia A, Blundo C, Carilla J, Osinaga-Acosta O, Malizia L, Silman M, Farfán-Ríos W, Malhi Y, Young KR, Francisco CC, Homeier J, Peralvo M, Pinto E, Jadan O, Aguirre N, Aguirre Z, Feeley KJ (2018) Widespread but heterogeneous responses of Andean forests to climate change. Nature 564:207–212
Fan RE, Chen PH, Lin CJ (2005) Working set selection using second order information for training support vector machines. J Mach Learn Res 6:1889–1918
Foster JR, D’Amato AW (2015) Montane forest ecotones moved downslope in northeastern USA in spite of warming between 1984 and 2011. Glob Chang Biol 21:4497–4507
Gehrig-Fasel J, Guisan A, Zimmermann N (2007) Tree line shifts in the Swiss Alps: climate change or land abandonment? J Veg Sci 18:571–582
Gómez C, Wulder MA, Montes F, Delgado JA (2012) Modeling forest structural parameters in the mediterranean pines of central Spain using QuickBird-2 imagery and classification and regression tree analysis (CART). Remote Sens 4:135–159
Gómez C, White JC, Wulder MA (2016) Optical remotely sensed time series data for land cover classification: a review. ISPRS J Photogramm Remote Sens 116:55–72
Gómez C, Alejandro P, Hermosilla T, Montes F, Pascual C, Ruiz LA, Álvarez-Taboada F, Tanase MA, Valbuena R (2019) Remote sensing for the Spanish forests in the 21stcentury: a review of advances, needs, and opportunities. For Syst 28:1–33
Gracia M, Meghelli N, Comas L, Retana J (2011) Land-cover changes in and around a National Park in a mountain landscape in the Pyrenees. Reg Environ Chang 11:349–358
Hartig F (2022) DHARMa: residual diagnostics for hierarchical (multi-level/mixed) regression models. R package version 0.4.6. https://CRAN.R-project.org/package=DHARMa
Hermosilla T, Wulder MA, White JC, Coops NC, Hobart GW (2015) An integrated Landsat time series protocol for change detection and generation of annual gap-free surface reflectance composites. Remote Sens Environ 158:220–234
Hermosilla T, Wulder MA, White JC, Coops NC, Hobart GW (2018) Disturbance-informed annual land cover classification maps of Canada’s forested ecosystems for a 29-year Landsat time series. Can J Remote Sens 44:67–87
Hernández L, Cañellas I, Alberdi I, Torres I, Montes F (2014) Assessing changes in species distribution from sequential large-scale forest inventories. Ann for Sci 71:161–171
Hernández L, Camarero JJ, Gil-Peregrín E, Saz-Sánchez MA, Cañellas I, Montes F (2019) Biotic factors and increasing aridity shape the altitudinal shifts of marginal Pyrenean silver fir populations in Europe. For Ecol Manage 432:558–567
Kassambara A (2021) rstatix: Pipe-friendly framework for basic statistical tests
Khatami R, Mountrakis G, Stehman SV (2016) A meta-analysis of remote sensing research on supervised pixel-based land-cover image classification processes: general guidelines for practitioners and future research. Remote Sens Environ 177:89–100
Knott JA, Desprez JM, Oswalt CM, Fei S (2019) Shifts in forest composition in the eastern United States. For Ecol Manage 433:176–183
Kouba Y, Camarero J, Alados C (2012) Roles of land-use and climate change on the establishment and regeneration dynamics of Mediterranean semi-deciduous oak forests. For Ecol Manage 274:143–150
Ledo A, Montes F, Condes S (2009) Species dynamics in a montane cloud forest: identifying factors involved in changes in tree diversity and functional characteristics. For Ecol Manage 258:S75
Lenoir J, Svenning JC (2014) Climate-related range shifts - a global multidimensional synthesis and new research directions. Ecography (cop) 38:15–28
Lu X, Liang E, Wang Y, Babst F, Camarero JJ (2020) Mountain treelines climb slowly despite rapid climate warming. Glob Ecol Biogeogr 30:305–315
Macias M, Andreu L, Bosch O et al (2006) Increasing aridity is enhancing silver fir (Abies alba Mill.) water stress in its south-western distribution limit. Clim Change 79:289–313
Matthews B (1975) Comparison of the predicted and observed secondary structure of T4 phage lysozyme. Biochim Biophys Acta - Protein Struct 405(2):442–451
Memoria. Parc Nacional d’Aigüestortes i Estany de Sant Maurici (2018) Generalitat de Catalunya, Departament de Territoti i Sostenibilitat
Meng Y, Banghua C, Mao P, Dong C, Xidong C, Qi L, Wang M, Wu Y (2020) Tree species distribution change study in Mount Tai based on Landsat remote sensing image data. Forests 11(2):130
MITECO (2021) Web: Ministerio para la Transición ecológica y el Reto demográfico. The National Forest Inventory. https://www.miteco.gob.es/es/biodiversidad/servicios/banco-datos-naturaleza/informacion-disponible/cartografia_informacion_disp.aspx
MITECO (2022) Ministerio para la Transición Ecológica y el Reto Demográfico. The Spanish Forest Map. https://www.miteco.gob.es/es/cartografia-y-sig/ide/descargas/biodiversidad/mfe.aspx. Accessed 3 Jan 2022
Oliva J, Colinas C (2007) Decline of silver fir (Abies alba Mill.) stands in the Spanish Pyrenees: role of management, historic dynamics and pathogens. For Ecol Manage 252:84–97
Oliva J, Colinas C (2010) Epidemiology of Heterobasidion abietinum and Viscum album on silver fir (Abies alba) stands of the Pyrenees. For Path 40:19–32
Pasquarella VJ, Holden CE, Woodcock CE (2018) Improved mapping of forest type using spectral-temporal Landsat features. Remote Sens Environ 210:193–207
R Core Team (2020) R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna
Rautiainen M, Lukeš P, Homolová L, Hovi A, Pisek J, Mõttus M (2018) Spectral properties of coniferous forests: a review of in situ and laboratory measurements. Remote Sens 10(2):207
Sánchez de Dios R, Gómez C, Aulló I, Cañellas I, Gea-Izquierdo G, Montes F, Sainz-Ollero H, Velázquez JC, Hernández L (2020) Fagussylvatica L. peripheral populations in the Mediterranean Iberian peninsula: climatic or anthropic relicts? Ecosystems 24:211–226
Sangüesa-Barreda G, Camarero JJ, Oliva J, Montes F, Gazol F (2015) Past logging, drought and pathogens interact and contribute to forest dieback. Agric for Meteorol 208:85–94
Sangüesa-Barreda G, Esper J, Büntgen U, Camarero JJ, Di Filippo A, Baliva M, Piovesan G (2020) Climate–human interactions contributed to historical forest recruitment dynamics in Mediterranean subalpine ecosystems. Glob Chang Biol 26:4988–4997
Serra-Maluquer X, Gazol A, Igual JM, Camarero JJ (2021) Silver fir growth responses to drought depend on interactions between tree characteristics, soil and neighbourhood features. For Ecol Manage 480:118625
Thompson SD, Nelson TA, White JC, Wulder MA (2015) Mapping dominant tree species over large forested areas using Landsat best-available-pixel image composites. Can J Remote Sens 41(3):203–218
Valavi R, Guillera-Arroita G, Lahoz-Monfort JJ, Elith J (2022) Predictive performance of presence-only species distribution models: a benchmark study with reproducible code. Ecol Monogr 92(1):e01486
Vallejo-Bombín R (2005) El mapa forestal de España escala 1:50.000 (MFE50) como base del tercer inventario forestal nacional. Cuad La Soc Española Ciencias for 19:205–210
White JC, Wulder MA, Hobart GW, Luther JE, Hermosilla T, Griffiths P, Coops NC, Hall RJ, Hostert P, Dyk A, Guindon L (2014) Pixel-based image compositing for large-area dense time series applications and science. Can J Remote Sens 40(3):192–212
White JC, Coops NC, Wulder MA, Vastaranta M, Hilker T, Tompalski P (2016) Remote sensing technologies for enhancing forest inventories: a review. Can J Remote Sens 42:619–641
White JC, Saarinen N, Kankare V, Wulder MA, Hermosilla T, Cops NC, Pickell PD, Holopainen M, Hyyppä J, Vastaranta M (2018) Confirmation of post-harvest spectral recovery from Landsat time series using measures of forest cover and height derived from airborne laser scanning data. Remote Sens Environ
Wulder MA, Masek JG, Cohen WB, Loveland TR, Woodcock CE (2012) Opening the archive: how free data has enabled the science and monitoring promise of Landsat. Remote Sens Environ 122:2–10
Wulder MA, Hermosilla T, Stinson G, Gougeon FA, White JC, Hill DA, Smiley BP (2020) Satellite-based time series land cover and change information to map forest area consistent with national and international reporting requirements. Forestry 93:331–343
Zuur AF, Ieno EN, Walker NJ, Saveliev AA, Smith GM (2009) Mixed effects models and extensions in ecology with R. In: Gail M, Samet JM (eds) Statistics for Biology and Health. Springer, New York, NY
We thank the Ministry for Ecological Transition and Demographic Challenge (MITECO) for providing data from the National Forest Inventory (NFI). We thank María Merced Aniz Montes, manager of the “Aigüestortes i Estany de Sant Maurici” National Park, and the National Park staff for their support in the fieldwork, and Adam Collins for the language revision.
The custom code and/or software application generated during and/or analyzed during the current study are available from the corresponding author upon reasonable request.
This work was supported by the Spanish Ministry of Science and Innovation (formerly Ministry of Economy, Industry, and Competitiveness) through the FPI program (BES-2017–081606) and the AGL2016-76769-C2-1-R and PID2020-119204RB-C21 project and by the National Parks Autonomous Agency (Spanish Ministry for the Ecological Transition and the Demographic Challenge) through the project 2481S/2017 OLDFORES. We acknowledge support of the publication fee by the CSIC Open Access Publication Support Initiative through its Unit of Information Resources for Research (URICI).
Ethics approval and consent to participate
Consent for publication
All authors gave their informed consent to this publication and its content.
The authors declare that they have no competing interests.
Handling Editor: Matteo Vizzarri
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Aulló-Maestro, I., Gómez, C., Hernández, L. et al. Monitoring montane-subalpine forest ecotone in the Pyrenees through sequential forest inventories and Landsat imagery. Annals of Forest Science 80, 32 (2023). https://doi.org/10.1186/s13595-023-01198-4