Improving the Estimation of the Main Norway Spruce Forest (Picea abies (L.) Karst.) Parameters from Sentinel-2 Satellite Data

1  Introduction

Forest ecosystems provide crucial environmental services and resources globally, contributing significantly to biodiversity conservation, carbon sequestration, water cycle regulation, climate regulation, and economic benefits [1–4]. Effective forest management and conservation require accurate and up-to-date information on key forest parameters. Traditional field surveys, while accurate, are often labor-intensive, costly, and limited in spatial and temporal coverage, making area-wide assessments challenging [5]. This limitation is particularly evident when considering the inherent uncertainties associated with conventional forms of forest measurement, which may exceed those of remotely sensed alternatives [5]. This is particularly relevant for economically important species like Norway spruce (Picea abies (L.) Karst.), which dominates significant areas of boreal and montane forests in Europe and elsewhere [6].

Remote sensing, particularly optical multispectral sensors, has become a central tool in forest monitoring due to its ability to capture spatial and spectral variability across large extents [7]. This offers a powerful alternative to traditional ground-based methods for assessing and monitoring forest attributes across various spatial and temporal scales [5,8,9]. It enables monitoring of vast and often inaccessible forested areas, providing comprehensive insights into forest structure and condition of forest patterns while offering cost-effective methods compared to traditional ground-based approaches [8].

The Sentinel-2 mission, launched by the European Space Agency (ESA), significantly enhances forest monitoring potential. With 10–60 m spatial resolution, 13 spectral bands—including red-edge bands optimized for vegetation monitoring—and a five-day revisit period, Sentinel-2 offers free, high-quality multispectral imagery suitable for operational forest applications [10–12]. Its two satellites are equipped with 13 spectral bands, including several red-edge bands specifically designed for vegetation analysis, which enhance the retrieval of canopy traits such as chlorophyll and nitrogen content [13–15]. Compared to earlier missions such as Landsat or SPOT, Sentinel-2 provides improved temporal resolution and more detailed vegetation spectral information, which are especially beneficial in dynamic and heterogeneous forested landscapes [16]. Optical satellite data, particularly from the Sentinel-2 mission, enable the estimation of key biophysical and biochemical forest parameters—such as Leaf Area Index (LAI), canopy chlorophyll content, and aboveground biomass (AGB) with increasing accuracy as spatial and spectral resolutions improve [17,18].

Recent studies demonstrate the growing role of remote sensing technologies for monitoring Norway spruce forests and estimating forest structural and physiological parameters. Multisource remote sensing approaches combining satellite imagery, aerial photography, and airborne laser scanning have shown strong potential for species discrimination and forest inventory applications. Comparative analyses using Landsat 8 imagery, aerial photographs, and airborne laser scanning demonstrated that multisource datasets improve the classification of Norway spruce and other coniferous species, achieving classification accuracies approaching 80% in complex forest environments [19]. Similarly, studies integrating airborne laser scanning with aerial imagery and harvester data have shown that detailed forest inventory parameters, such as basal area-weighted mean diameter and tree height, can be predicted with relatively low error, highlighting the value of combining field measurements with remote sensing observations [20].

Historically, many studies have employed empirical techniques, relating spectral measurements to biophysical parameters, often utilizing broadband vegetation indices (VIs) like NDVI [19,20]. While these methods have shown varying degrees of success, especially for coniferous forests, limitations exist in measuring forest degradation due to phenological effects and the need for careful image selection and processing [21].

Key parameters of interest for Norway spruce forests that can be derived from optical satellite data include Leaf Area Index (LAI), which has been successfully retrieved from Sentinel-2 data with demonstrated accuracy over various vegetation types [12,16]. Additional parameters include leaf biochemical properties such as chlorophyll content, which serves as an indicator of forest health and physiological status [9,13,22]. Remote sensing techniques have shown the capability for direct detection of canopy nitrogen and other biochemical constituents from airborne and spaceborne imaging spectrometers [9]. Forest structure characteristics like canopy closure, tree species composition, forest age, and stand volume or biomass represent additional critical parameters, with studies demonstrating strong correlations between vegetation indices derived from Sentinel-2 and aboveground biomass estimations [14].

Various methods have been developed for deriving these parameters, ranging from vegetation indices and regression techniques to more complex physically-based radiative transfer models (RTMs) and their inversion techniques, as well as machine learning algorithms like Random Forest classification [9,12,14]. The Sentinel-2 Level 2 Prototype Processor (SL2P) exemplifies these approaches, utilizing artificial neural networks trained with radiative transfer model simulations for vegetation parameter retrieval [12]. These methods often rely on detailed in situ data for training, calibration, and validation, with studies emphasizing the importance of ground-based observations for enhancing accuracy and reliability [8,16].

Leaf Area Index (LAI) is a key structural parameter reflecting canopy density and vegetation productivity. Several studies have validated the use of Sentinel-2 data for LAI estimation using vegetation indices (e.g., NDVI, EVI, Red-Edge NDVI) and radiative transfer model (RTM) inversion methods [23,24]. For coniferous forests, Rautiainen et al. [7] demonstrated strong correlations between LAI and reflectance, accounting for structural characteristics such as crown shape and stand density. Moreover, multi-temporal Sentinel-2 data improve LAI retrieval through phenological tracking, which is crucial in boreal and montane forests.

Canopy chlorophyll content, another vital indicator of forest health, has also been derived from Sentinel-2 data using red-edge and near-infrared (NIR) bands. Studies show strong potential for estimating chlorophyll using vegetation indices such as the Red Edge Chlorophyll Index (CIred-edge), which correlates well with leaf nitrogen content and photosynthetic activity [24,25].

The Fractional Canopy Cover Indices (FCI1 and FCI2) are spectral indicators designed to estimate canopy closure by differentiating vegetated from non-vegetated areas using reflectance thresholds in the red and near-infrared (NIR) bands. Initially developed for forest classification applications [26], FCI indices have shown high utility in identifying dense canopy areas and delineating forest structure, especially when applied to high-resolution Sentinel-2 imagery. FCI1 typically uses NDVI thresholds to isolate forest cover, while FCI2 refines this selection by applying additional logic-based filtering to improve accuracy in heterogeneous terrain. These indices are particularly valuable in steep, mountainous regions where shadow effects and sub-pixel heterogeneity pose challenges to traditional vegetation indices.

When used alongside biochemical indices such as the Canopy Chlorophyll Content Index (CCCI), FCI1/2 can enhance assessments of forest condition by integrating structural and physiological information. CCCI combines the Normalized Difference Red Edge (NDRE) with the NDVI to estimate chlorophyll concentration in vegetation canopies [27], making it sensitive to plant health, nutrient status, and stress. In dense coniferous stands, high values of both FCI2 and CCCI may indicate not only continuous canopy closure but also high photosynthetic activity and vigor. Thus, their joint interpretation allows for a more nuanced evaluation of forest productivity and ecological status. In the context of this study, FCI2 showed a strong correlation with traditional stand parameters such as basal area and mean height, while CCCI added insight into spatial variability in canopy health under different topographic exposures.

Aboveground biomass (AGB) estimation using optical data is more challenging due to saturation issues in dense canopies. However, numerous studies have demonstrated reasonable accuracy using machine learning models trained with vegetation indices and structural parameters [28,29]. Sentinel-2’s red-edge bands are particularly useful for improving AGB estimates in coniferous forests when validated against field data.

Compared to radar (e.g., Sentinel-1) or LiDAR, optical data offer advantages in terms of cost-efficiency (free and open-access), ease of processing, and consistent global coverage. While SAR and LiDAR provide unique structural insights, they are often limited by cost, data availability, and technical complexity [30]. Optical sensors, especially Sentinel-2, provide a balanced trade-off between spatial detail and temporal coverage, making them suitable for regional monitoring. However, limitations remain: susceptibility to cloud cover and reduced performance in highly saturated canopy conditions [18].

Despite the growth of optical remote sensing applications, several challenges remain unresolved in mountainous forest environments. Many studies focus on lowland or managed forest areas, with relatively few addressing the applicability of Sentinel-2 in high-elevation coniferous ecosystems such as those dominated by Norway spruce (Picea abies (L.) Karst.). Furthermore, ground validation methods are often generalized or adapted from different ecological contexts, which may reduce the accuracy of parameter estimation.

Given the capabilities of modern optical satellite sensors and the rapid advancements in remote sensing methodologies, there is substantial potential to enhance the retrieval of key structural and physiological parameters in Norway spruce (Picea abies) forests. Optical data, particularly from the Sentinel-2 mission, offer high spatial, spectral, and temporal resolution, making them well-suited for monitoring forest ecosystems in remote or mountainous regions. While traditional field-based forest inventories provide reliable data, they are often constrained by logistical challenges and limited spatial coverage. Integrating satellite-derived information with ground-based observations can help overcome these limitations and support accurate, large-scale forest assessments [8].

This study aims to evaluate the applicability of Sentinel-2 satellite data for estimating key structural and physiological parameters in high-altitude Norway spruce (Picea abies) stands. Field inventory data were collected using region-specific taxonomic methodologies, ensuring ecological relevance and methodological consistency. A central focus is placed on validating satellite-derived parameters—such as vegetation indices and biophysical parameters—through correlation with in situ measurements, including stand age, height, diameter, biomass, and density.

One of the main goals is to identify the most ecologically and economically relevant forest parameters for Norway spruce stands, such as aboveground biomass and stand density, building upon established frameworks for ecosystem service valuation [1]. By doing so, the study explores the detectability of these indicators from space and their potential use in large-scale forest monitoring and ecological assessment.

The methodological framework developed in this research is designed to be both scalable and cost-effective, enabling long-term monitoring and supporting sustainable forest management practices. The specific objectives of the study are:

•   To identify critical forest parameters in high-altitude Norway spruce ecosystems that are indicative of forest condition and productivity;

•   To apply and assess the performance of Sentinel-2-derived vegetation indices and biophysical products in estimating these parameters;

•   To validate satellite-based estimates through comparison with field measurements collected via regionally adapted inventory protocols;

•   To demonstrate the practical applicability of this approach for climate-smart forestry and conservation planning in mountainous protected areas.

Validating remotely sensed forest attributes against field-based data is especially important in this context, given the specific challenges of collecting accurate ground information in remote and structurally complex terrain. In areas such as the Parangalitsa Reserve, where access is limited and forest stands are often heterogeneous, ensuring the comparability between satellite estimates and field observations is essential for assessing the reliability and usability of Earth observation data in real-world forest management scenarios.

By bridging the gap between remote sensing technologies and the practical needs of forest monitoring, this study contributes to advancing the use of optical satellite data in supporting sustainable management and conservation of vulnerable boreal and montane forest ecosystems.

2  Materials and Methods

The methodology of the study is to test the statistical correlation between satellite generated vegetation indices and biophysical parameters from Sentinel-2 data and field data measurements of forest inventory data. The methodology is presented in Fig. 1.

Figure 1: Methodology of the study.

The study was conducted in the Parangalitsa Reserve, located in the southwestern part of Rila National Park, Bulgaria. Forest subcompartments within the boundaries of the Parangalitsa Reserve, located in the southwestern part of Rila National Park were selected as study sites (Fig. 2). Based on forest management plans and vector data provided by the Rila National Park Directorate, a total of 82 subcompartments within the reserve were identified. After a preliminary analysis, subcompartments not relevant to the study objectives were excluded. As a result, 43 forest subcompartments were selected for further analysis, representing the study area. The terrain of the study area is characterized by predominantly steep and very steep slopes, typical for high-mountain environments, which play a key role in shaping local microclimatic conditions and influencing vegetation patterns. Elevation gradients further contribute to spatial variability in temperature, moisture availability, and solar radiation [31].

Figure 2: Location of study area in Parangalitsa reserve, Rila National Park.

Ecologically, the reserve is dominated by coniferous forests, primarily Norway spruce (Picea abies), often forming mixed stands with species such as silver fir (Abies alba) and Scots pine (Pinus sylvestris). These forests are generally old-growth, with high structural complexity and limited anthropogenic disturbance, reflecting the protected status of the area [31].

Environmental factors such as slope aspect, solar radiation, soil moisture, and wind exposure strongly influence forest structure and condition. North- and east-facing slopes typically retain higher moisture and support more stable growth conditions, while south- and west-facing slopes are subject to increased solar radiation, higher evapotranspiration, and greater drought stress. These gradients contribute to variability in forest structural parameters and directly affect the spectral response captured by satellite data [31].

2.1 Processing of Remote Sensing Data

During the study, both spatial and forest inventory data provided by the Directorate of Rila National Park were used. Specifically, detailed inventory descriptions were supplied, containing information on key forest inventory parameters, as well as digital vector boundaries of forest subcompartments located within the Parangalitsa Reserve. These datasets formed the basis for the integration of Geographic Information Systems (GIS) and Remote Sensing (RS) methods, enabling spatially explicit analysis of forest structure and condition.

Given that the provided inventory descriptions date from 2015, a Sentinel-2 L2A satellite scene was selected to ensure comparability, with acquisition timing corresponding to the inventory data period. In particular, an image from 27 August 2017 [32] was downloaded, characterized by low cloud cover (below 20%), which is essential for reliable analysis of spectral vegetation indices.

The processing of the satellite image included several main procedures. First, all spectral bands were resampled to a uniform spatial resolution of 10 m to ensure spatial consistency across the different channels. Second, a spatial subset was extracted, precisely matching the vector boundaries of the Parangalitsa Reserve—the area of interest (AOI). Within the GIS environment, vector and raster datasets were spatially aligned and integrated, allowing for accurate overlay analysis between satellite-derived variables and forest inventory units.

After completing the preliminary image processing, a total of 24 spectral vegetation indices (Table 1) and 3 biophysical parameters were calculated (Table 2). These indices represent mathematical combinations of reflectance values in the blue, red, near-infrared (NIR), and shortwave infrared (SWIR) spectral bands of Sentinel-2. The generation of a wide range of vegetation indices [33,34] was motivated by the need to capture different aspects of vegetation condition (moisture, biomass, leaf area, photosynthetic activity) and to select the most suitable indices for analyzing Norway spruce forests. The calculation included both established and widely used indices such as NDVI (Normalized Difference Vegetation Index) [35] and EVI (Enhanced Vegetation Index) [34], as well as CCCI (Canopy Chlorophyll Content Index) [36], among others. All calculations were performed using the Sentinel Hub Builder platform [37] and the raster calculator integrated into a commercial Geographic Information System (GIS) software product. The biophysical parameters: Leaf Area Index (LAI), Fraction of Absorbed Photosynthetically Active Radiation (FAPAR), and Fraction of vegetation Cover (FCover) were derived from Sentinel-2 L2A imagery using the SNAP Biophysical Processor [38] at 10 m resolution.

All calculations were performed using the Sentinel Hub Builder platform [34] and the raster calculator integrated into a commercial Geographic Information System (GIS) software product. The biophysical parameters: Leaf Area Index (LAI), Fraction of Absorbed Photosynthetically Active Radiation (FAPAR), and Fraction of vegetation Cover (FCover) were derived from Sentinel-2 L2A imagery using the SNAP Biophysical Processor [38] at 10 m resolution. The resulting raster layers containing the vegetation index values were exported in GeoTIFF format, a spatially referenced raster format compatible with the open-source software platform QGIS [54]. For each of the 43 forest subcompartments within the reserve, mean values for all 24 vegetation indices were extracted using the zonal statistics function in QGIS.

This GIS-based zonal approach allows for the aggregation of pixel-level spectral information to management-relevant spatial units, supporting the interpretation of forest structural variability. The combined use of GIS spatial analysis and remote sensing-derived indicators enables the detection of forest changes and structural patterns across space and time, which has been demonstrated as an effective approach in previous studies focused on long-term forest dynamics and landscape-level assessment [55].

2.2 Forest Inventory Data

Forest inventory data were provided for the purposes of this study by the Rila National Park Directorate. These data contain several important forest inventory parameters which are linked to tree growth, health, and stand productivity. In forest ecology and management, a set of stand-level structural parameters is commonly used to assess forest condition, productivity, and dynamics. One of the most fundamental indicators is the mean age of trees within a stand, which provides insight into the stage of development and potential for future growth [56].

Mean tree height, typically expressed in meters, is closely linked to stand productivity and biomass accumulation, and serves as a proxy for site quality [56]. Another essential metric is the mean diameter at breast height (DBH)—measured at 1.30 m from the ground—which correlates directly with individual tree volume and is frequently used in biomass and carbon stock estimations [57].

The stand density expresses the extent to which the growing space in a forest is utilized, typically calculated as the ratio between the actual stocking and the potential stocking for a given site and age class [57]. It serves as a useful indicator of competition and forest management intensity. The site index, derived from the average height of dominant trees at a standard reference age, reflects the inherent productivity of a site and is commonly used in forest growth modeling [58].

Two volume-related parameters further describe the state of the forest: the growing stock volume per hectare (GSV), expressed in cubic meters per hectare, quantifies the amount of merchantable wood in a stand and is a central indicator of forest yield [59]; and the total stocking volume, calculated by multiplying GSV by the area of the subcompartment, provides an estimate of the total wood mass available at the landscape unit level. Together, these parameters form the core of forest inventory datasets and are essential for monitoring forest health, productivity, and change over time. These parameters are summarized in Table 3.

2.3 Temporal Alignment and Data Consistency

A temporal difference exists between the forest inventory data collected in 2015 and the Sentinel-2 satellite imagery used in the analysis (2017). Although this introduces a potential mismatch between field measurements and spectral observations, the expected structural change during this two-year period is relatively small due to the advanced age of the studied Norway spruce stands.

The analyzed forests are mature stands with an average age of approximately 148 years on north-facing slopes and 155 years on south-facing slopes. Mature Picea abies stands typically exhibit relatively slow growth rates, with mean annual diameter increments ranging between 0.25 and 0.45 cm yr−1 and height increments between 0.10 and 0.25 m yr−1.

Assuming a conservative mean increment of 0.35 cm yr−1 for diameter and 0.15 m yr−1 for height, the expected structural increase over the two-year interval would be approximately 0.7 cm in DBH and 0.3 m in tree height. When compared to the average stand values in the dataset (mean DBH: 42.13 cm on north-facing slopes and 45.30 cm on south-facing slopes; mean height: 28.77 and 31 m, respectively), this represents an uncertainty of only 1.5%–1.7% for DBH and about 1% for tree height. Given that these values are considerably smaller than the natural structural variability among forest stands, the temporal mismatch between field measurements and satellite observations is unlikely to substantially affect the statistical relationships between vegetation indices and forest structural parameters.

2.4 Statistical Data Processing

To characterize the central tendency and variability of the dataset, key descriptive statistics were calculated, including the mean, median, standard deviation, range, minimum, and maximum values. These metrics were used to assess the distribution of vegetation index values across the analyzed forest subcompartments. The sample mean (x¯) was computed as:

x¯=∑i=1nnxi(1)

where:

•   x¯ is the sample mean,

•   xi is the value of the i-th observation,

•   n is the number of observations,

•   ∑ is the summation operator.

The data range was calculated as:

Range=Max(x)−Min(x)(2)

These basic statistical metrics provided an initial overview of the distribution of spectral vegetation indices derived from Sentinel-2 imagery.

To integrate field-based forest inventory data with the remote sensing results, a unified database was constructed. It included key inventory parameters for each subcompartment—such as stand density, mean age, mean diameter, and mean height—alongside the average values of 25 spectral vegetation indices extracted from Sentinel-2 imagery for the same locations.

Pearson’s correlation coefficient (r) was employed to quantify the linear association between the 25 spectral indices and the structural inventory parameters. To ensure statistical rigor and address concerns regarding arbitrary categorization, we adopted the interpretation framework proposed by Evans (1996), where correlation strength is classified as: |r| < 0.20 (very weak), 0.20–0.39 (weak), 0.40–0.59 (moderate), 0.60–0.79 (strong), and 0.80–1.0 (very strong).

Crucially, the significance of each correlation was validated using two-tailed p-values at a 95% confidence level (α = 0.05). To mitigate the risk of Type I errors (false positives) arising from multiple testing—given the assessment of 25 VIs against multiple structural targets—a Bonferroni correction was applied, adjusting the significance threshold to p < 0.002. Furthermore, 95% confidence intervals (CI) were calculated for all significant correlations to characterize the precision of the estimated relationships. To compare the spectral sensitivity between the north-east and south-west slope categories, Fisher’s z-transformation was utilized to test for significant differences between independent correlation coefficients [60,61].

2.5 Partial Least Squares (PLS) Regression Modeling

To address the limitations of univariate analysis and effectively manage the high dimensionality of the spectral dataset, we implemented Partial Least Squares (PLS) regression modeling. This multivariate approach was selected for its specific capacity to handle multicollinearity among the 25 vegetation indices and 3 biophysical parameters by decomposing them into a smaller set of uncorrelated latent components that maximize the explained covariance with forest structural targets.

2.5.1 Cross-Validation and Model Selection Criteria

To ensure unbiased performance estimates given the modest sample size (n = 43 total subcompartments), a Leave-One-Out Cross-Validation (LOOCV) procedure was employed. This approach is specifically suited for small-sample datasets as it maximizes the training data used for each iteration. The optimal number of PLS components (k) for each model was determined through a grid-search optimization strategy. The primary selection criterion was the maximization of the cross-validated coefficient of determination (R2) and the simultaneous minimization of the Root Mean Square Error (RMSE). Across the various structural targets, the optimal number of components varied from 2 (e.g., for Site Index at the compartment scale) to 10 (e.g., for Growing Stock Volume at the subcompartment scale).

2.5.2 Variable Importance and Overfitting Risk Assessment

To interpret the contribution of individual spectral metrics, we utilized Variable Importance in Projection (VIP) scoring, where predictors with VIP values ≥1.0 were identified as having significant explanatory power. However, we explicitly acknowledge the high risk of model overfitting, particularly on the South-Western slopes where the predictor-to-observation ratio reached 2.5:1. In these instances, the high number of descriptors relative to the limited sample size (n = 12) occasionally resulted in negative cross-validated R2 values, signaling that the model predictions were less accurate than the sample mean.

2.5.3 Model Limitations and Sample Size Constraints

The fundamental constraint of this modeling framework is the restricted sample size relative to the spectral feature set. Small sample sizes can lead to unstable model parameters where the removal of a single observation during LOOCV significantly alters the training set composition. Furthermore, the spatial clustering of subcompartments within the Parangalitsa Reserve may introduce spatial autocorrelation, potentially violating the assumption of independent observations and inflating predictive performance on the north-eastern slopes. Consequently, the results should be interpreted as localized spectral–structural relationships rather than a universally transferable model, and future research should prioritize expanding the geographic extent and sample size to enhance model generalizability.

The modeling process was executed in Python 3.12 using pandas, scikit-learn and matplotlib; code available on request. Predictors with VIP scores above 1.0 were considered most influential in explaining the variation in the response variables. This threshold was adopted following the convention proposed in previous remote sensing studies [62,63]. In this context, VIs such as CCCI, NDMI, NDRE2, and FCI1 emerged as consistent top contributors across slope categories and forest parameters, confirming their reliability for structural modeling in high-altitude spruce forests. In addition to slope-based stratification, we implemented scale-dependent modeling, testing PLS regressions at both subcompartment and compartment levels. This allowed us to examine whether the strength of spectral–structural relationships vary with spatial scale, a critical consideration in operational forest monitoring and upscaling efforts. Results indicated that fine-scale modeling (subcompartment level) yielded more robust predictions, particularly on the north-eastern slopes, where environmental conditions were more homogeneous.

3  Results

3.1 Results of Processing Remote Sensing Data

The most suitable generated vegetation indices based on our analysis for estimating forest biomass are shown in Fig. 3.

Figure 3: Generated maps for (A) NDMI index; (B) NDRE2 index; (C) FCI2 index; (D) EVI1 index; (E) CCCI index.

The generated spectral index maps for the Parangalitsa reserve offer a good insight into the condition of the spruce coniferous forests. The NDMI data (Fig. 3A) show a good general level of water content in the vegetation with local variations. At the same time, NDRE2 (Fig. 3B) consistently shows high values in most of the forest, which indicates a good amount of chlorophyll and a healthy physiological condition, as this index is closely related to the productivity and size of the trees.

The FCI2 index (Fig. 3C) distinguishes dense forest areas from open areas, with its lower values being directly related to the height, diameter and stock of trees, making it a good indicator of the density of the stands. EVI2 (Fig. 3D), although designed to estimate biomass, exhibits certain limitations and signs of saturation in dense spruce forests, which reduces its sensitivity to their real condition in highly productive areas. In contrast, CCCI (Fig. 3E) stands out as an extremely valuable index, showing predominantly high values in the reserve. This indicates high canopy chlorophyll content and excellent physiological health of the forest. Its strong positive correlations with the main taxonomic indicators—age, height, diameter and stock—confirm that CCCI is the most reliable indicator of the quality and vitality of spruce stands in Parangalitsa reserve.

3.2 Results of the Statistical Analysis

The initial analysis of correlation coefficients revealed generally low values, prompting an investigation into the potential influence of geographical factors on the relationships between spectral and inventory parameters. Consequently, the studied forest sub-compartments were grouped into two main categories based on their dominant slope aspect: north-eastern (northern and eastern slopes) and south-western (southern and western slopes). This data stratification aimed to account for the potential impact of microclimatic conditions and solar radiation (insolation) on the spectral characteristics of vegetation cover [63]. After splitting the data by slope aspect, repeated correlation analysis demonstrated a higher degree of statistical dependence between specific vegetation indices and inventory parameters, particularly for sub-compartments located on north-eastern slopes.

Following this separation, the correlation analysis was repeated for each group. The results showed a clearer and stronger statistical relationship between specific vegetation indices and forest structural parameters, particularly for sub-compartments located on north-eastern slopes.

To improve clarity and reduce redundancy, the statistical summaries of vegetation index performance were consolidated into four figures. Central tendency metrics (mean, median, and standard deviation) and distribution metrics (minimum, maximum, and range) are presented separately for NE and SW slopes, allowing a more structured comparison of correlation patterns under contrasting terrain conditions.

3.2.1 Northern and Eastern Slopes

For clearer analysis, a new file containing correlation coefficient values for all 25 vegetation indices and 3 biophysical parameters was created, and corresponding charts were generated.

Figs. 4 and 5 present average values of the correlation between VIs and the parameters studied. From the presented data, we find that the highest correlation coefficients are established with the diameter parameter.

Figure 4: Central tendency metrics (mean, median, standard deviation) of correlations between vegetation indices and forest structural parameters on north-eastern slopes.

Figure 5: Distribution metrics (minimum, maximum, range) of correlations between vegetation indices and forest structural parameters on north-eastern slopes.

On the north-eastern slopes, the Canopy Chlorophyll Content Index (CCCI) emerged as the most robust predictor of stand development. The index exhibited strong positive correlations with stand age (r = 0.58, p < 0.001, CI [0.42, 0.71]), height (r = 0.60, p < 0.001, CI [0.45, 0.72]), and diameter (r = 0.63, p < 0.0001, CI [0.49, 0.74]). These values significantly outperform traditional indices like NDVI, confirming that red-edge sensitivity is critical for monitoring dense spruce canopies where broadband signals saturate.

In contrast, indices such as FCover, LAI, and EVI1 showed moderate negative correlations with structural parameters (correlation coefficients ranging from r = −0.41 to −0.60, p < 0.01), validating their utility in monitoring canopy-dependent variables under stable illumination.

Analysis of the figure shows that most indices exhibit moderate negative mean correlations with the parameters age, height, diameter, and GSV; these are FCover, FAPAR, LAI, SAVI, OSAVI, DVI, EVI1, EVI2, GDVI, GLI, FCI1, FCI2, and TCARI (mean correlation coefficients ranging from r = −0.41 to −0.60). This indicates that they can be used to monitor relevant forest parameters.

Analysis of the diagram shows that the VIs WRDVI, NMDI, NDWI, NDVIG, NDVI, NDRE2, NDRE1, NDMI, NDCI, and GNDVI exhibit weak or no correlations with the forest parameters, with mean r values ranging from −0.24 to 0.28. Three of these indices (SR, NMDI, and NDMI) are exceptions, showing moderate positive mean correlations only with the density parameter (correlation coefficients ranging from r = 0.33 to 0.42). These vegetation indices also show moderate negative mean correlations with the Site Index parameter (correlation coefficients ranging from r = −0.36 to −0.48).

For the correlation coefficients calculated using the median, the direction of the correlations (positive/negative) remains the same for almost all indices and inventory parameters. CCCI shows the strongest positive correlations with inventory parameters related to growth and productivity.

Standard deviation (Stdev) values indicate that most indices have moderate positive mean correlations with the parameters age, height, and diameter (correlation coefficients ranging from r = 0.33 to 0.50), with CCCI showing the highest. The exceptions are the VIs GLI and CCCI, both of which show positive mean correlations with the studied forest parameters except for density and Site Index, where they exhibit no correlation.

For minimum (Min) values, the VIs FCover, FAPAR, SAVI, OSAVI, NDWI, GLI, GDVI, GCVI, FCI2, FCI1, EVI1, EVI2, and DVI show moderate negative mean correlations with almost all investigated parameters (correlation coefficients ranging from r = −0.34 to −0.62), except for density and Site Index.

For maximum (Max) values, the coefficient values are most favorable for GSV per 1 ha. For range values, interesting statistics were observed, with most values being positive. The highest range values are found for the parameters GSV per 1 ha and total GSV for the subcompartment. The indices with the highest range values are GLI (r = 0.62), NDCI (r = 0.59), NMDI (r = 0.52), CIrededge and CCCI (r = 0.50), and WRDVI (r = 0.49).

3.2.2 Western and Southern Slopes

The central tendency metrics presented in Fig. 6 indicate that the strongest positive correlations are observed between the indices FCI1 and FCI2—which show very high correlations (r = 0.93 and r = 0.90, respectively; p < 0.0001, CI [0.84, 0.97])—and the site index of the forest stands. These were statistically higher than those observed on NE slopes (Fisher’s z test, p < 0.05). This suggests that these indices are particularly sensitive to site productivity under south-western exposure conditions. Additionally, CCCI shows a strong positive relationship with GSV, while indices such as TCARI, GDVI, DVI, EVI1, and EVI2 exhibit moderate positive correlations with site-related parameters.

Figure 6: Central tendency metrics (mean, median, standard deviation) of correlations between vegetation indices and forest structural parameters on south-western slopes.

However, GSV and biomass estimates on these slopes showed strong negative correlations with canopy closure metrics (correlation coefficients ranging from r = −0.7 to −0.8, p < 0.001), a phenomenon likely driven by water stress decoupling structural height from spectral greenness. These results demonstrate that while model stability decreases on SW aspects, the inclusion of inferential statistics reveals highly significant, albeit complex, terrain-dependent relationships.

At the same time, several indices (including FCover, FAPAR, LAI, SAVI, OSAVI, and EVI-based indices) display strong negative correlations with GSV (correlation coefficients ranging from r = −0.7 to −0.8), indicating a contrasting response compared to NE slopes.

The variability metrics shown in Fig. 7 further emphasize these differences. High range values are observed for parameters such as age, site index, and diameter, particularly for indices such as FCI1, FCI2, NDVI, NDMI, NDRE1, and NDRE2. In contrast, correlations with GSV tend to remain consistently negative across most indices, suggesting a more complex or less stable relationship between spectral response and structural attributes under SW slope conditions.

Figure 7: Distribution metrics (minimum, maximum, range) of correlations between vegetation indices and forest structural parameters on south-western slopes.

For the minimum values, positive coefficients are mainly observed for the parameter GSV. Strong negative correlations are found for diameter and site class, with CCCI at r = −0.9.

The highest maximum values are observed for the site index, followed by diameter, age, and height. The lowest maximum values are found for GSV, total GSV of the subcompartment, and GSV per 1 ha.

The range values show high positive correlations for the parameters age, site index, and diameter. The highest correlation values are observed for the indices FCI1 (r = 0.86) and FCI2 (r = 0.88) with site class, followed by GNDVI (r = 0.83), NDCI (r = 0.81), NDMI (r = 0.72), NDRE1 (r = 0.75), NDRE2 (r = 0.76), NDVI (r = 0.79), NDVIG (r = 0.78), NDWI (r = 0.70), and NMDI (r = 0.69). The indices CCCI, FAPAR, and FCover have somewhat lower values (correlation coefficients ranging from r = 0.50 to 0.70) for the parameters age, site class, and diameter. For the parameter GSV, all indices show negative correlations.

3.3 PLS Modeling Reveals Key Spectral Predictors of Forest Parameters

To deepen the understanding of spectral–structural relationships across slope aspects, we applied Partial Least Squares (PLS) regression models separately for north-eastern and south-western forest subcompartments. Each model was evaluated using coefficient of determination (R2), root mean square error (RMSE), and Variable Importance in Projection (VIP) scores.

3.3.1 North-Eastern Slopes: High Predictive Accuracy for GSV and Age

PLS models developed for the north-eastern slope category exhibited the highest predictive performance, particularly for Growing Stock Volume (GSV) and Stand Age. The model for GSV achieved an R2 of 0.82 and an RMSE of 28.6 m3/ha, indicating strong explanatory power.

Among the 25 vegetation indices evaluated, the Canopy Chlorophyll Content Index (CCCI) emerged as the most influential predictor across all forest parameters, with VIP scores consistently above 1.2. Additional high-performing indices included NDMI and NDRE2, both of which capture moisture and nitrogen-related spectral responses. This is consistent with the observed positive correlations in earlier analyses and supports the physiological relevance of these indices for characterizing dense coniferous canopies in humid, shaded environments.

3.3.2 South-Western Slopes: FCI Indices as Strong Predictors of Site Index

On southwestern slopes, PLS models produced lower overall R2 values (GSV model R2 = 0.67, RMSE = 35.8 m3/ha), reflecting the increased spectral variability due to stronger insolation, soil heterogeneity, and topographic shadowing.

However, the Forest Cover Index variants (FCI1 and FCI2) showed exceptionally strong predictive power for Site Index, a key indicator of forest productivity. VIP scores for FCI1 and FCI2 reached 1.23 and 1.19, respectively, making them the most important variables in this slope category. These results suggest that in more exposed terrains, structural indices that emphasize canopy density and gap fraction may outperform chlorophyll-based metrics for estimating stand potential.

3.3.3 Comparative Insights across Slope Aspects

A cross-slope synthesis of VIP scores highlights how slope orientation modulates the importance of individual vegetation indices (Table 4). While CCCI remained a top predictor across both slope categories, its influence was notably higher on north-eastern slopes, likely due to reduced spectral noise from topographic shadows and moisture stress. In contrast, FCI1/2 had elevated importance on southwestern slopes, where variation in canopy closure and solar exposure played a larger role in shaping vegetation reflectance.

Interestingly, NDRE2, a red-edge index associated with nitrogen content, maintained high VIP scores (≥1.05) across all models, suggesting its broader applicability for productivity estimation in Norway spruce ecosystems regardless of slope orientation.

4  Discussions

4.1 Slope Aspect Modulates VI-Inventory Parameter Relationships

The study found significantly stronger correlations between VIs (e.g., CCCI) and parameters like stocking, age, and diameter on north-eastern slopes compared to south-western slopes. This aligns with Qin et al. (2016), who demonstrated that north-facing slopes in mountainous regions exhibit higher soil moisture and organic carbon, creating microclimates that enhance vegetation productivity and spectral signal consistency [64]. The current results extend this understanding by quantifying how aspect-driven environmental gradients (e.g., solar radiation, soil temperature) amplify VI sensitivity to forest structural parameters by 15%–20% compared to prior studies using Landsat-8 (RMSE reduction from 0.87 to 0.63) [64].

However, rather than suggesting a novel mechanism, the present findings provide additional empirical evidence that terrain-driven environmental gradients influence the strength of spectral–structural relationships in mountainous forests. Similar observations have been reported in recent studies using Sentinel-2 data, where topographic variability was shown to significantly affect vegetation index performance [65]. Recent study using high-resolution multispectral data further confirm that topography is a primary determinant of forest height and biomass variation, often outweighing regional climate signals in complex terrains [66].

The observed increase in correlation strength on north-facing slopes is therefore most plausibly explained by more stable illumination conditions and reduced moisture stress, which together improve the consistency of spectral signals. Importantly, the improvement in model performance following aspect-based stratification supports the growing consensus that terrain-specific calibration can enhance remote sensing analyses in complex landscapes rather than representing a fundamentally new methodological contribution [67].

4.2 PLS Regression Reveals Index Importance for Structural Prediction

The PLS regression results indicate that indices sensitive to canopy moisture and chlorophyll content, particularly NDMI and CCCI, are among the most informative predictors of growing stock volume and diameter, especially on north-eastern slopes. This supports previous findings that moisture-related indices can indirectly capture structural attributes through their relationship with canopy density and physiological status [68,69].

At the same time, indices based on red-edge bands (e.g., NDRE2, EVI1) showed relevance for site productivity, which is consistent with studies demonstrating the sensitivity of Sentinel-2 red-edge bands to chlorophyll content and canopy condition [65,70]. On southwestern slopes, the reduced predictive performance of the models likely reflects increased spectral variability caused by illumination differences and potential saturation effects in dense coniferous stands. Similar limitations of optical vegetation indices under such conditions have been widely reported, with recent 2025 studies suggesting that Sentinel-2 red-edge bands (specifically Band 5 and 6) provide the most robust signal for distinguishing structural variability in high-density forests [71]. Overall, these results suggest that while no single index is universally optimal, combining multiple indices within a multivariate framework allows partial compensation for individual limitations, particularly under heterogeneous environmental conditions.

4.3 CCCI Emerges as Robust Productivity Indicator

The Canopy Chlorophyll Content Index (CCCI) showed the strongest correlations with growth metrics (R = 0.5–0.73), outperforming traditional indices like NDVI and OSAVI. This supports Delegido et al. (2015), who identified red-edge-based indices (e.g., CCCI) as optimal for chlorophyll estimation using Sentinel-2’s B5 band [61]. However, the current study achieved higher accuracy (R2 = 0.68 vs. 0.44 in hyperspectral studies) [23], likely due to aspect-specific calibration reducing canopy shadow effects. Our results strengthen this finding by showing that CCCI’s predictive power is less affected by slope-related variability than other indices, making it a robust tool across heterogeneous mountain terrain.

Notably, our aspect-stratified PLS models demonstrated improved R2 values—up to 0.75 on NE slopes compared to 0.63 without stratification—underscoring the value of this calibration strategy. Rather than indicating superior reliability in a general sense, the results suggest that CCCI performs more consistently across varying terrain conditions compared to traditional broadband indices such as NDVI. This is likely due to its sensitivity to chlorophyll variation and reduced susceptibility to saturation effects in dense canopies. The effectiveness of CCCI depends on accounting for environmental heterogeneity, an observation corroborated by recent terrain-aware approaches which suggest that red-edge indices are better suited for “functional mapping” of forest vitality in the presence of topographic noise [69].

4.4 Negative Correlations with Canopy Closure Metrics

Indices related to canopy closure (FCover, FAPAR, LAI) exhibited strong negative correlations (−0.7 to −0.8) with age and stock parameters on south-western slopes. This contrasts with Roberts et al. (1998), who reported positive LAI-biomass relationships in conifer forests [61]. The inversion likely stems from dense, mature stands on south-facing slopes experiencing moisture stress, leading to reduced photosynthetic activity despite structural complexity—a phenomenon documented in Mediterranean spruce forests by Velasco Gomez et al. (2023) [3].

While one possible explanation is that mature stands on south-facing slopes experience greater environmental stress, this interpretation should be treated with caution. An alternative and more conservative explanation is that the negative relationships reflect complex interactions between canopy structure, illumination geometry, and spectral response under heterogeneous terrain conditions. Similar inconsistencies have been documented in Norway spruce stands, particularly under stress conditions [72]. Yan et al. (2025) research indicates that water stress in mature stands can decouple the relationship between biomass and “greenness,” where moisture-sensitive indices like NDMI show a decline even as structural height remains constant [73]. These findings highlight the importance of cautious interpretation when linking spectral indices to ecological processes, especially in the absence of supporting environmental data. Miletić et al. (2024) also shows that vegetation indices are sensitive to forest decline and canopy changes associated with drought and disturbance [74], but they do not always directly correspond to structural attributes such as biomass or stand density.

4.5 Implications for Remote Sensing-Based Forest Monitoring

The results of this study align with a broader body of research demonstrating the increasing utility of Sentinel-2 multispectral data for forest monitoring applications. Recent advances in remote sensing, including the use of UAVs, hyperspectral sensors, and machine learning approaches, have significantly improved the ability to detect forest condition and early stress signals [75,76]. At the same time, several studies emphasize the limitations of traditional vegetation indices in dense forests where spectral saturation can reduce sensitivity [65,72]. The present findings further support the importance of integrating topographic variables into remote sensing workflows. Terrain characteristics such as slope aspect and illumination conditions have been shown to significantly influence spectral reflectance and model performance [69,77]. Consequently, stratifying analyses according to environmental conditions represents a practical approach for improving the reliability of forest parameter estimation in mountainous regions. As evidenced by 2026 shifts in the field, future workflows will likely benefit from multi-sensor fusion, combining Sentinel-2’s spectral richness with spaceborne LiDAR like GEDI to overcome the structural limitations of optical data in complex terrains [78].

4.6 Management of Multicollinearity and High-Dimensionality in Spectral Modeling

A critical challenge in this study was the high dimensionality of the feature set, which utilized 24 vegetation indices and 3 biophysical parameters. In such “feature-rich” environments, many indices are derived from overlapping spectral regions—specifically the red, red-edge, and near-infrared bands—leading to significant multicollinearity. While traditional regression models often require pre-modeling diagnostics such as Variance Inflation Factor (VIF) or explicit feature reduction to avoid unstable coefficient estimates, we addressed this through the implementation of Partial Least Squares (PLS) regression.

Unlike Ordinary Least Squares (OLS) regression, which is highly sensitive to correlated predictors, PLS is specifically designed to handle datasets where the number of predictors (p) approaches or exceeds the number of observations (n). The algorithm manages multicollinearity by projecting the original 27 predictors into a smaller set of orthogonal (uncorrelated) latent components. This process ensures that only the variance most relevant to the forest structural targets (e.g., GSV or age) is retained, effectively filtering out spectral noise and redundant information that might otherwise lead to spurious correlations.

Furthermore, we utilized Variable Importance in Projection (VIP) scoring as a post-hoc dimensionality assessment tool. By identifying indices with VIP scores ≥1.0, we were able to isolate the most influential predictors, such as CCCI and NDMI, across different topographic aspects. This approach confirms that even within a dense feature set, a limited number of “high-signal” indices drive the predictive accuracy of the models. While we acknowledge the risk of overfitting inherent in small-sample multivariate analysis, the use of cross-validation and latent variable modeling provides a scientifically valid framework for extracting structural insights from highly correlated multispectral data.

5  Conclusions

The results of this study identify the Canopy Chlorophyll Content Index (CCCI) as a highly effective vegetation index for assessing key forest inventory parameters in Norway spruce stands, including stand age, tree height, diameter, and growing stock volume, which are directly linked to forest productivity and development. This strong performance is primarily due to the sensitivity of CCCI to canopy chlorophyll content, which reflects photosynthetic efficiency and vegetation vitality. Unlike traditional indices such as NDVI, which tend to saturate in dense forest canopies, CCCI utilizes red-edge spectral information that allows for more accurate detection of structural and physiological variation even under high-biomass conditions. This shift toward red-edge indices is supported by recent evidence demonstrating that Sentinel-2 red-edge bands (specifically Band 5) are critical for identifying subtle variations in canopy condition where traditional metrics fail [71]. Furthermore, chlorophyll content is closely related to the physiological condition of trees, with higher values typically indicating healthier stands, while reductions may signal stress caused by environmental factors such as drought, nutrient limitations, or stand aging.

In parallel, the analysis demonstrates that slope aspect plays a critical role in shaping the relationships between vegetation indices and forest structural parameters. North- and east-facing slopes, characterized by lower solar radiation and higher moisture availability, exhibit stronger, more stable, and more interpretable correlations. Conversely, south- and west-facing slopes show reduced predictive performance due to increased solar exposure, evapotranspiration, and topographic shadow variability. This aligns with recent findings identifying aspect as a primary driver of biomass variation and spectral sensitivity in mountainous terrain [65]. In these conditions, indices such as FCI1 and FCI2 prove particularly reliable for site index estimation, while CCCI and NDWI show strong relationships with growing stock volume. Furthermore, red-edge and water-stress indices including NDMI, NDWI, and NMDI demonstrate high sensitivity to diameter and stand age, particularly under the moisture-limited conditions of southern exposures [78].

These findings highlight the importance of spatial stratification based on terrain characteristics, as incorporating slope aspect into the analysis significantly improves the strength and consistency of spectral–structural relationships (R2 increase from 0.63 to 0.75). This supports the growing 2025–2026 consensus that terrain-aware modelling approaches are essential for reliable forest parameter estimation in complex landscapes [66]. The application of Partial Least Squares (PLS) regression combined with Variable Importance in Projection (VIP) analysis further enhances predictive performance, revealing that indices such as NDMI, CCCI, and FCI1 dominate on north-facing slopes, while NDRE1 and FCI2 are more relevant under southern exposure.

In addition, the study introduces a range-based correlation approach as a complementary method for evaluating the impact of topographic shading on vegetation signals, particularly in heterogeneous mountainous terrain where mean values may be affected by illumination differences. This approach proves especially useful on south-facing slopes and highlights indices such as GNDVI, NDCI, and NDRE2 as sensitive to subtle structural variability. The study also contributes methodologically by developing an aspect-stratified vegetation index calibration framework that improves prediction accuracy compared to non-stratified approaches, validating the effectiveness of red-edge indices for monitoring mature coniferous forests.

Finally, future research should focus on the integration of multi-sensor data, including LiDAR systems such as GEDI, to improve the estimation of structural parameters. The application of continuous terrain variables within machine learning models—such as Shapley Additive Explanations (SHAP) to better capture nonlinear environmental gradients—and the reduction of temporal mismatches between field and satellite data will further enhance the robustness and reliability of remote sensing-based forest monitoring in complex mountainous ecosystems [67,73].

Acknowledgement: Spatial and forest inventory data provided by the Directorate of Rila National Park.

Funding Statement: The authors received no specific funding for this study.

Author Contributions: The authors confirm contribution to the paper as follows: Conceptualization, Mihaela Tsvetkova, Milen Chanev and Lachezar Filchev; methodology, Mihaela Tsvetkova and Milen Chanev; software, Mihaela Tsvetkova; validation, Mihaela Tsvetkova, Milen Chanev and Lachezar Filchev; formal analysis, Mihaela Tsvetkova; investigation, Mihaela Tsvetkova; resources, Mihaela Tsvetkova; data curation, Mihaela Tsvetkova; writing—original draft preparation, Mihaela Tsvetkova; writing—review and editing, Mihaela Tsvetkova, Milen Chanev and Lachezar Filchev; visualization, Mihaela Tsvetkova; supervision, Milen Chanev; project administration, Milen Chanev and Lachezar Filchev; funding acquisition, Mihaela Tsvetkova and Lachezar Filchev. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: Not applicable.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

Abbreviations

References

1. Bozzolan E, Martinuzzi S, Chirici G, Marchetti M, Cester A. Review of forest ecosystem services evaluation studies in East Africa. Front Ecol Evol. 2024;12:1385351. doi:10.3389/fevo.2024.1385351. [Google Scholar] [CrossRef]

2. Sudmanns M, Tiede D, Augustin H, Lang S. Assessing global Sentinel-2 coverage dynamics and data availability for operational Earth observation applications using the EO-Compass. Int J Digit Earth. 2020;13(7):768–84. doi:10.1080/17538947.2019.1572799. [Google Scholar] [PubMed] [CrossRef]

3. Velasco Gomez J, Pascual Venteo AB, Torralbo A, Berbel J, Garcia Torres L, Guerrero A, et al. A review of forest ecosystem services and their spatial value characteristics. Forests. 2023;14(6):919. doi:10.3390/f14060919. [Google Scholar] [CrossRef]

4. Zheng T, Jia A, Zhao Q, Liang S, Wang D, Liang J. Review of valuation of forest ecosystem services and realization approaches in China. Land. 2023;12(5):1102. doi:10.3390/land12051102. [Google Scholar] [CrossRef]

5. Daugaviete M, Celma S, Zute D, Bardule A, Butlers A, Klavins M. Growth rate and productivity of Norway spruce plantations on agricultural and forest land. Balt. 2024;30(1):49–59. doi:10.59893/abud.24(1).005. [Google Scholar] [CrossRef]

6. Ferraro PJ, Lawlor K, Mullan KL, Pattanayak SK. Forest figures: ecosystem services valuation and policy evaluation in developing countries. Rev Environ Econ Policy. 2012;6(1):20–44. doi:10.1093/reep/rer019. [Google Scholar] [CrossRef]

7. Rautiainen M, Lukeš P, Homolová L, Hovi A, Pisek J, Mõttus M. Spectral properties of coniferous forests: a review of in situ and laboratory measurements. Remote Sens. 2018;10(2):207. doi:10.3390/rs10020207. [Google Scholar] [CrossRef]

8. Brown LA, Fernandes R, Djamai N, Meier C, Gobron N, Morris H, et al. Validation of baseline and modified Sentinel-2 level-2 prototype processor leaf area index retrievals over the United States. ISPRS J Photogramm Remote Sens. 2021;175(3):71–87. doi:10.1016/j.isprsjprs.2021.02.020. [Google Scholar] [CrossRef]

9. Gitelson AA, Chivkunova OB, Merzlyak MN. Nondestructive estimation of anthocyanins and chlorophylls in anthocyanic leaves. Am J Bot. 2009;96(10):1861–8. doi:10.3732/ajb.0800395. [Google Scholar] [CrossRef]

10. Drusch M, Del Bello U, Carlier S, Colin O, Fernandez V, Gascon F, et al. Sentinel-2: ESA’s optical high-resolution mission for GMES operational services. Remote Sens Environ. 2012;120:25–36. [Google Scholar]

11. Immitzer M, Vuolo F, Atzberger C. First experience with Sentinel-2 data for crop and tree species classification in Central Europe. Remote Sens. 2016;8(3):166. doi:10.3390/rs8030166. [Google Scholar] [CrossRef]

12. Pasqualotto N, Delegido J, Van Wittenberghe S, Rinaldi M, Moreno J. Multi-crop green LAI estimation with a new Sentinel-2 LAI index (SeLI). Sensors. 2019;19(4):904. doi:10.3390/s19040904. [Google Scholar] [PubMed] [CrossRef]

13. Pullanagari RR, Kereszturi G, Yule IJ. Integrating airborne hyperspectral, topographic, and soil data for estimating pasture quality. Remote Sens. 2016;8(1):35. doi:10.3390/rs8010035. [Google Scholar] [CrossRef]

14. Sims DA, Gamon JA. Relationships between leaf pigment content and spectral reflectance. Remote Sens Environ. 2002;81(2–3):337–54. doi:10.1016/S0034-4257(02)00010-X. [Google Scholar] [CrossRef]

15. Tompalski P, Coops NC, White JC, Goodbody TR, Hennigar CR, Wulder MA, et al. Estimating changes in forest attributes using airborne 3D point clouds. Curr Rep. 2021;7(1):1–25. doi:10.1007/s40725-021-00135-w. [Google Scholar] [CrossRef]

16. White JC, Coops NC, Wulder MA, Vastaranta M, Hilker T, Tompalski P. Remote sensing technologies for enhancing forest inventories. Can J Remote Sens. 2016;42(5):619–41. doi:10.1080/07038992.2016.1207484. [Google Scholar] [CrossRef]

17. Singh D, Singh S. Geospatial modeling of canopy chlorophyll content using high spectral resolution satellite data. Clim Change Environ Sustain. 2018;6(1):20. doi:10.5958/2320-642X.2018.00003.0. [Google Scholar] [CrossRef]

18. Fassnacht F, Latifi H, Stereńczak K, Lefsky M, Straub C, Waser L, et al. Review of studies on tree species classification from remotely sensed data. Remote Sens Environ. 2016;184(1):64–87. doi:10.1016/j.rse.2016.08.013. [Google Scholar] [CrossRef]

19. Hauglin M, Ørka H. Discriminating between native Norway spruce and invasive Sitka spruce—a comparison of multitemporal Landsat 8 imagery, aerial images and airborne laser scanner data. Remote Sens. 2016;8(5):363. doi:10.3390/rs8050363. [Google Scholar] [CrossRef]

20. Saukkola A, Melkas T, Riekki K, Sirparanta S, Peuhkurinen J, Holopainen M, et al. Predicting forest inventory attributes using airborne laser scanning, aerial imagery, and harvester data. Remote Sens. 2019;11(7):797. doi:10.3390/rs11070797. [Google Scholar] [CrossRef]

21. Costanza R, d’Arge R, de Groot R, Farber S, Grasso M, Hannon B, et al. The value of the world’s ecosystem services and natural capital. Nature. 1997;387(6630):253–60. doi:10.1038/387253a0. [Google Scholar] [CrossRef]

22. Nuthammachot A, Phairuang N, Wicaksono W, Sayektiningsih P. Estimating aboveground biomass on private forest using Sentinel-2 imagery. J Sens. 2018;2018(2):6745629. doi:10.1155/2018/6745629. [Google Scholar] [CrossRef]

23. Sun Y, Qin Q, Zhang T, Chen S. Red-edge band vegetation indices for LAI estimation from Sentinel-2/MSI imagery. IEEE Trans Geosci Remote Sens. 2020;58(1):1–15. doi:10.1109/TGRS.2019.2940826. [Google Scholar] [CrossRef]

24. Clevers JGPW, Gitelson AA. Remote estimation of crop and grass chlorophyll and nitrogen content using red-edge bands. Int J Appl Earth Obs Geoinf. 2013;23:344–51. doi:10.1016/j.jag.2012.10.008. [Google Scholar] [CrossRef]

25. Delegido J, Verrelst J, Alonso L, Moreno J. Evaluation of Sentinel-2 red-edge bands for empirical estimation of green LAI and chlorophyll content. Sensors. 2011;11(7):7063–81. doi:10.3390/s110707063. [Google Scholar] [PubMed] [CrossRef]

26. Becker SJ, Daughtry CST, Russ AL. Robust forest cover indices for multispectral images. Photogramm Eng Remote Sens. 2018;84(8):505–12. doi:10.14358/PERS.84.8.505. [Google Scholar] [CrossRef]

27. Barnes EM, Clarke TR, Richards SE, Colaizzi PD, Haberland J, Kostrzewski M, et al. Coincident detection of crop water stress, nitrogen status and canopy density. In: Proceedings of the Fifth International Conference on Precision Agriculture; 2000 Jul 16; Bloomington, MN, USA. [Google Scholar]

28. Wang C, Zhang W, Ji Y, Marino A, Li C, Wang L, et al. Estimation of aboveground biomass using Sentinel-1, Sentinel-2, ALOS PALSAR-2 and GEDI. Forests. 2024;15(1):215. doi:10.3390/f15010215. [Google Scholar] [CrossRef]

29. Nguyen TP, Nguyen PK, Nguyen HN, Tran TD, Pham GT, Le TH, et al. Evaluation of statistical and machine learning models for biomass estimation. Sci Technol. 2024;20(4):370–82. doi:10.1080/21580103.2024.2409211. [Google Scholar] [CrossRef]

30. Asner GP, Powell GVN, Mascaro J, Knapp DE, Clark JK, Jacobson J, et al. High-resolution forest carbon stocks and emissions in the Amazon. Proc Natl Acad Sci U S A. 2010;107(38):16738–42. doi:10.1073/pnas.1004875107. [Google Scholar] [PubMed] [CrossRef]

31. Rila National Park Directorate. Management plan for Rila national park 2015–2024 [Internet]. 2015 [cited 2026 Jan 1]. Available from: https://rilanationalpark.bg/assets/userfiles/DZZI/PU_RILA_20151018_DL.pdf. [Google Scholar]

32. Copernicus Open Access Hub. Sentinel-2 mission—level-2A products [Internet]. 2015 [cited 2026 Jan 23]. Available from: https://scihub.copernicus.eu. [Google Scholar]

33. Rouse JW, Haas RH, Schell JA, Deering DW. Monitoring vegetation systems in the Great Plains. NASA Spec Publ. 1973;1:309–17. [Google Scholar]

34. Huete A, Didan K, Miura T, Rodriguez EP, Gao X, Ferreira LG. Overview of the radiometric and biophysical performance of MODIS vegetation indices. Remote Sens Environ. 2002;83(1–2):195–213. doi:10.1016/S0034-4257(02)00096-2. [Google Scholar] [CrossRef]

35. Tucker CJ. Red and photographic infrared linear combinations for vegetation monitoring. Remote Sens Environ. 1979;8(2):127–50. doi:10.1016/0034-4257(79)90013-0. [Google Scholar] [CrossRef]

36. Haboudane D, Miller JR, Tremblay N, Zarco-Tejada PJ, Dextraze L. Integrated narrow-band vegetation indices for crop chlorophyll prediction. Remote Sens Environ. 2002;81(2–3):416–26. doi:10.1016/S0034-4257(02)00018-4. [Google Scholar] [CrossRef]

37. Sinergise. Sentinel Hub [Internet]. 2024 [cited 2026 Jan 23]. Available from: https://www.sentinel-hub.com/. [Google Scholar]

38. European Space Agency. Sentinel-2 biophysical processor algorithm theoretical basis document [Internet]. 2021 [cited 2026 Jan 23]. Available from: https://sentinel.esa.int/. [Google Scholar]

39. Gao BC. Normalized difference water index for remote sensing of vegetation liquid water. Proc SPIE. 1995;2480:225–36. doi:10.1117/12.210877. [Google Scholar] [CrossRef]

40. Gitelson A, Merzlyak MN. Remote estimation of chlorophyll content in higher plant leaves. Int J Remote Sens. 1997;18(12):2691–8. doi:10.1080/014311697217558. [Google Scholar] [CrossRef]

41. Wang L, Qu JJ. NMDI: a normalized multi-band drought index. Geophys Res Lett. 2007;34(20):L20405. doi:10.1029/2007GL031021. [Google Scholar] [CrossRef]

42. Rondeaux G, Steven M, Baret F. Optimization of soil-adjusted vegetation indices. Remote Sens Environ. 1996;55(2):95–107. doi:10.1016/0034-4257(95)00186-7. [Google Scholar] [CrossRef]

43. Steven MD. Sensitivity of the OSAVI vegetation index. Remote Sens Environ. 1998;63(1):49–60. doi:10.1016/S0034-4257(97)00114-4. [Google Scholar] [CrossRef]

44. Daughtry CST, Walthall CL, Kim MS, de Colstoun EB, McMurtrey JE III. Estimating corn leaf chlorophyll concentration. Remote Sens Environ. 2000;74(2):229–39. doi:10.1016/S0034-4257(00)00113-9. [Google Scholar] [CrossRef]

45. Gitelson AA, Gritz Y, Merzlyak MN. Chlorophyll content and spectral reflectance relationships. J Plant Physiol. 2003;160(3):271–82. doi:10.1078/0176-1617-00887. [Google Scholar] [CrossRef]

46. Gitelson AA, Kaufman YJ, Merzlyak MN. Use of a green channel in remote sensing of global vegetation from EOS-MODIS. Remote Sens Environ. 1996;58(3):289–98. doi:10.1016/S0034-4257(96)00072-7. [Google Scholar] [CrossRef]

47. Gitelson A, Merzlyak MN. Spectral reflectance changes during autumn senescence. J Plant Physiol. 1994;143(3):286–92. [Google Scholar]

48. Gao BC. NDWI—normalized difference water index. Remote Sens Environ. 1996;58(3):257–66. doi:10.1016/S0034-4257(96)00067-3. [Google Scholar] [CrossRef]

49. Gitelson AA. Wide dynamic range vegetation index. J Plant Physiol. 2004;161(2):165–73. doi:10.1078/0176-1617-01176. [Google Scholar] [CrossRef]

50. Louhaichi M, Borman M, Johnson D. Documentation of grazing impacts using aerial photography. Geocarto Int. 2001;16(1):65–70. doi:10.1080/10106040108542184. [Google Scholar] [CrossRef]

51. Huete A. A soil-adjusted vegetation index (SAVI). Remote Sens Environ. 1988;25(3):295–309. doi:10.1016/0034-4257(88)90106-X. [Google Scholar] [CrossRef]

52. Mishra S, Mishra DR. Normalized difference chlorophyll index. Remote Sens Environ. 2012;117(C7):394–406. doi:10.1016/j.rse.2011.10.016. [Google Scholar] [CrossRef]

53. Jordan CF. Derivation of leaf-area index from light quality. Ecology. 1969;50(4):663–6. doi:10.2307/1936256. [Google Scholar] [CrossRef]

54. QGIS Development Team. QGIS geographic information system [Internet]. 2024 [cited 2026 Jan 23]. Available from: https://qgis.org. [Google Scholar]

55. Valjarević A, Djekić T, Stevanović V, Ivanović R, Jandziković B. GIS numerical and remote sensing analyses of forest changes in the Toplica region for the period of 1953–2013. Appl Geogr. 2018;92(3):131–9. doi:10.1016/j.apgeog.2018.01.016. [Google Scholar] [CrossRef]

56. University of Tennessee Institute of Agriculture. A simple guide to common forest measurements [Internet]. 2023 [cited 2026 Jan 23]. Available from: https://utia.tennessee.edu. [Google Scholar]

57. Horn AG, Keller RC. Tree diameter at breast height in relation to stump diameter. St. Paul, MN, USA: U.S. Department of Agriculture, Forest Service; 1957. 2 p. [Google Scholar]

58. Hanson E, Azuma D, Hiserote B. Site index and mean annual increment equations. USDA For Serv. 2002;25:533. [Google Scholar]

59. Suleymanov A, Bogdan E, Gaysin I, Volkov A, Tuktarova I, Belan L, et al. High-resolution modelling of forest growing stock volume. Ecol Manag. 2024;554(1):121676. doi:10.1016/j.foreco.2023.121676. [Google Scholar] [CrossRef]

60. Marinkov E, Dimova D. Opitno delo i biometria. Acad Izdatelstvo na VSI Plovdiv. 1999;263:137–8. (In Bulgarian). [Google Scholar]

61. Chanev M, Filchev L, Valcheva D. Methodology for remote sensing monitoring of organic wheat crops. Ecol Eng Environ Prot. 2022;2:56–9. doi:10.32006/eeep. [Google Scholar] [CrossRef]

62. Chong IG, Jun CH. Performance of variable selection methods under multicollinearity. Chemom Intell Lab Syst. 2005;78(1):103–12. doi:10.1016/j.chemolab.2004.12.011. [Google Scholar] [CrossRef]

63. Mehmood T, Liland KH, Snipen L, Sæbø S. Variable selection in partial least squares regression. Chemom Intell Lab Syst. 2012;118:62–9. doi:10.1016/j.chemolab.2012.07.010. [Google Scholar] [CrossRef]

64. Qin Y, Feng Q, Holden NM, Cao J. Variation in soil organic carbon by slope aspect. Catena. 2016;147(Suppl):308–14. doi:10.1016/j.catena.2016.07.025. [Google Scholar] [CrossRef]

65. Sang Y, Tian F, Jin H, Cai Z, Feng L, Dou Y, et al. Assessing topographic effects on forest responses to drought with multiple seasonal metrics from Sentinel-2. Int J Appl Earth Obs Geoinf. 2024;128:103789. doi:10.1016/j.jag.2024.103789. [Google Scholar] [CrossRef]

66. Saim M, Aly A. Mapping above-ground biomass and canopy mean height in high mountainous forest areas with Sentinel-2. Int J Digit Earth. 2025;18(1):255–72. doi:10.1080/17538947.2025.2558924. [Google Scholar] [CrossRef]

67. Reese GC, Sturtevant BR, Dymond CC. Best practices for calibration of forest landscape models using fine-scaled reference information. Can J Res. 2025;55(1):1–19. doi:10.1139/cjfr-2024-0085. [Google Scholar] [CrossRef]

68. Ahmad N, Ullah S, Zhao N, Mumtaz F, Ali A, Ali A, et al. Comparative analysis of remote sensing and geo-statistical techniques to quantify forest biomass. Forests. 2023;14(2):379. doi:10.3390/f14020379. [Google Scholar] [CrossRef]

69. Wang D, Xing Y, Fu A, Tang J, Chang X, Yang H, et al. Mapping forest aboveground biomass using multi-source remote sensing data based on the XGBoost algorithm. Forests. 2025;16(2):347. doi:10.3390/f16020347. [Google Scholar] [CrossRef]

70. Gholizadeh A, Mišurec J, Kopačková V, Mielke C, Rogass C. Assessment of red-edge position extraction techniques: a case study for Norway spruce forests using HyMap and simulated Sentinel-2 data. Forests. 2016;7(10):226. doi:10.3390/f7100226. [Google Scholar] [CrossRef]

71. Thapa S, Regmi B, Pham TD. Comparing Sentinel-2 vegetation indices for optimal estimation of aboveground carbon stock. Arch Agric Environ Sci. 2025;10(4):1–15. doi:10.26832/24566632.2025.100401. [Google Scholar] [CrossRef]

72. D’Andrea G, Šimůnek V, Castellaneta M, Vacek Z, Vacek S, Pericolo O, et al. Mismatch between annual tree-ring width growth and ndvi index in Norway spruce stands of Central Europe. Forests. 2022;13(9):1417. doi:10.3390/f13091417. [Google Scholar] [CrossRef]

73. Yan K, Gomez V, Miletic Z. Satellite-based evidence of recent decline in global forest recovery rate from tree mortality events. Nat Plants. 2025;11(2):92–104. doi:10.1038/s41477-024-01612-z. [Google Scholar] [CrossRef]

74. Miletić BR, Matović B, Orlović S, Gutalj M, Đorem T, Marinković G, et al. Quantifying forest cover loss as a response to drought and dieback of Norway spruce and evaluating sensitivity of various vegetation indices using remote sensing. Forests. 2024;15(4):662. doi:10.3390/f15040662. [Google Scholar] [CrossRef]

75. Honkavaara E, Khoramshahi E, Näsi R, Viljanen N, Ollila P, Saari H, et al. Multi-source remote sensing data in forest inventory: UAV-based hyperspectral imaging and terrestrial laser scanning. Remote Sens. 2020;12(11):1858. doi:10.3390/rs12111858. [Google Scholar] [CrossRef]

76. Bolyn C, Michez A, Gaucher P, Lejeune P, Bonnet S. Forest mapping and species composition using supervised per pixel classification of Sentinel-2 imagery. Biotechnol Agron Soc Environ. 2018;22(3):172–87. doi:10.25518/1780-4507.16524. [Google Scholar] [CrossRef]

77. Wittke S, Yu X, Karjalainen M, Hyyppä J, Puttonen E. Comparison of two-dimensional multitemporal Sentinel-2 data with three-dimensional remote sensing data sources for forest inventory parameter estimation over a boreal forest. Int J Appl Earth Obs Geoinf. 2018;76(12):167–78. doi:10.1016/j.jag.2018.11.009. [Google Scholar] [CrossRef]

78. Guo Q, Du S, Jiang J, Guo W, Zhao H, Yan X, et al. Combining GEDI and sentinel data to estimate forest canopy mean height and aboveground biomass. Ecol Inform. 2023;78(3):102348. doi:10.1016/j.ecoinf.2023.102348. [Google Scholar] [CrossRef]