Monitoring montane-subalpine forest ecotone in the Pyrenees through sequential forest inventories and Landsat imagery

Abstract

Key message

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.

Context

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.

Aims

The main goal of this study is to assess the spatiotemporal dynamics of the silver fir–mountain pine ecotone and its underlying driving factors.

Methods

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).

Results

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.

Conclusion

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.

1 Introduction

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).

Fig. 1

Study area. a Distribution of silver fir and mountain pine in Europe (source: Caudullo et al. 2019) and location of the study area (star). b Location of the Inventory plots in the study area: Spanish National Forest Inventory (NFI, triangles) and Aigüestortes I Estani de Sant Maurici National Park Forest Inventory (NPFI, diamonds); distribution of silver fir and mountain pine in the study area (sources: The Spanish Forest Map (MITECO 2022); Vallejo-Bombín, 2005); and National Park and buffer zone limits

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).

Table 1 Forest inventory data employed in this study. NFI refers to the National Forest Inventories. NPFI refers to the Inventories carried out by the National Parks Autonomous Agency

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.

Table 2 Landsat images employed in this work. Target years in italics

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).

Table 3 Silver fir and mountain pine basal area (m²ha⁻¹) and standard deviations (in parenthesis) by altitudinal classes for the three referenced years in permanent plots

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 r₁₂(Δ_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):

$${r}_{12}\left({\Delta }_{dbh}\right)=\frac{\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}f(i, j){\mathrm{I}}_{ij}({\Delta }_{dbh}-\delta ,{\Delta }_{dbh}+\delta )}{\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}{\mathrm{I}}_{ij}({\Delta }_{dbh}-\delta ,{\Delta }_{dbh}+\delta )}$$

(1)

Δ_dbh is the between-trees dbh difference lags. For each pair of individuals i,j, I_ij(Δ_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 r₁₂(Δ_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:

$$RE=\sum\limits_{\min(\Delta_{dbh})}^{\max(\Delta_{dbh})}\Delta_{dbh}\cdot r_{12}(\Delta_{dbh})$$

(2)

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 (BA_T) 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 BA_T 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 BA_T. 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 (AUC_ROC). 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 Results

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).

Fig. 2

Mean basal area (BA) of silver fir and mountain pine by 200 m altitude intervals in the three reference years: 1991, 2002, and 2015. Error bars represent standard deviations. Significant changes are marked with asterisks (*)

Table 4 Friedman test (p-value) results among years for silver fir basal area, mountain pine basal area, and replacement index within paired plots, grouped by altitudinal intervals of 200 m and level of protection (National Park and buffer zone). Post hoc Wilcoxon results are only given when significant (a, b, c)

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).

Fig. 3

Species replacement index (RE) evolution by 200 m altitude intervals in the three reference years 1991, 2002, and 2015: positive and negative RE values indicate the dominance in the lower diameter classes of mountain pine or silver fir respectively, whereas RE values close to 0 indicate species equilibrium across diameter classes. Error bars represent standard deviations. Significant changes are marked with asterisks (*)

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 BA_T 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 BA_T as fixed effects and BA_T × altitude interaction and plot as random effect. The difference between marginal and conditional R² 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 BA_T and 0.15 for RE) (Table 5).

Table 5 Goodness of fit and coefficients (standard errors in parentheses) for the best linear mixed-effects models of total basal area and replacement index. Significance levels *P < 0.05, **P < 0.01, and ***P < 0.001, respectively. Predictors have been standardized previously to fit the models

We found that BA_T is mainly affected by the level of protection. BA_T increases in those plots inside the National Park area and slightly in South-facing slopes. Regarding temporal effects, BA_T 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 BA_T × 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 BA_T and higher altitudes. The BA_T × 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 BA_T 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 (AUC_ROC = 0.80) to 0.74 for mountain pine in 2002 (AUC_ROC = 0.82), indicating the reliability of the SVM classifications (Table 6).

Fig. 4

Confusion matrices obtained through 10 k-fold cross-validation of the 6 SVM classifiers for silver fir and mountain pine classifications in a 1991, b 2002, and c 2015. 0 indicates absence and 1 presence. For each species for each year, the upper left data indicates the true positive classifications, the upper right—false positive, lower left—false negative, and lower right—true negative. For the total values: the first column indicates the true positive rate (TPR) and the second column the false positive rate (FPR). Bold data indicates the total accuracy

Table 6 Area under the curve for the receiver-operator curve (AUC_ROC) and Mathews correlation coefficient (MCC) estimated through 10 k-fold cross-validation of silver fir and mountain pine SVM classifiers for 1991, 2002, and 2015

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).

Table 7 Predicted area (ha) occupied at the three reference years (1991, 2002, 2015) and the changes in area during the period 1991–2002 (% compared to 1991), 2002–2015 (% compared to 2002), and the percentage of total changes (1991–2015) for both silver fir and mountain pine within the National Park and buffer zone

Fig. 5

Evolution from 1991 to 2015 of expanded, retracted, and stable predicted populations of silver fir and mountain pine in the National Park and buffer zone

Fig. 6

Total predicted area (ha) occupied by silver fir and mountain pine where the basal area of each species was above the 40% threshold by 200 m altitude intervals in the three different reference years (1991, 2002, and 2015)

Fig. 7

Expansion (gray), retraction (orange), and stable area (green) (ha) regarding altitudes (x-axes) and aspect (y-axes) during the studied period for silver fir (a, c, e) and mountain pine (b, d, f)

4 Discussion

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).

5 Conclusion

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.