Integrating Field Inventory and Sentinel-2 Imagery to Assess Carbon Stock and Biomass Dynamics of Gymnosperms in Ayubia National Park, Pakistan

INTRODUCTION

Carbon management through forest carbon stocking has become a critical strategy for addressing global climate change driven by escalating greenhouse gas (GHG) emissions. Human activity is the dominant contributor to increasing atmospheric CO₂ concentrations, the primary driver of global warming. Forest ecosystems serve as natural carbon sinks, absorbing and storing large amounts of carbon dioxide (Barber et al., 2021). Covering nearly one-third of the Earth’s land surface, forests sequester about 2.4 gigatons of carbon annually. Alarmingly, this crucial carbon sink is under increasing pressure, with significant forest cover loss reported over the past three decades (FAO, 2020). Within forest ecosystems, carbon is stored in various components: soils (44%), living biomass (42%), deadwood (8%), and litter (5%) (Sun and Liu, 2019). Forests account for about 80% of above-ground and 40% of below-ground terrestrial carbon (Shin et al., 2007). Between 2000 and 2007, global forests sequestered nearly 4.1 gigatons of carbon per year. However, deforestation and degradation, primarily caused by human activities, have sharply reduced this capacity (Ali et al., 2023). Deforestation not only releases large amounts of stored carbon into the atmosphere but also diminishes future carbon storage potential. In essence, forests can either mitigate or exacerbate climate change depending on how they are managed (Raihan et al., 2022). Sustainable forest management, reforestation, and conservation are essential strategies to reverse this trend. Beyond their role in carbon sequestration, forests offer critical ecosystem services and socio-economic benefits (Stinson et al., 2022). It is estimated that land-use change, including deforestation, contributes about 10% of global CO₂ emissions. Coniferous forests, in particular, play a pivotal role in carbon storage, with their woody biomass holding nearly 50.8% of stored carbon (Thomas and Martin, 2012). Yet, only 13.25% of global forests are formally protected. To safeguard carbon-rich ecosystems, countries with high deforestation rates must expand conservation zones especially in tropical and subtropical regions, where biodiversity and carbon stocks are most at risk (Ali et al., 2025). Forest biomass, which encompasses the total mass of living organisms within an ecosystem, is a widely accepted proxy for carbon storage potential (Siddiq et al., 2021). It includes above-ground biomass (AGB) such as stems, branches, and foliage and below-ground biomass (BGB), such as roots. AGB, which is more readily measurable than BGB, is commonly used as a surrogate for total forest biomass and health (Goetz and Dubayah, 2011). The international REDD+ framework (Reducing Emissions from Deforestation and Forest Degradation) incentivizes developing nations to reduce deforestation and enhance carbon stocks by offering financial compensation. A key requirement for REDD+ participation is the establishment of robust forest monitoring systems, including reference emission levels and MRV (Measurement, Reporting, and Verification) protocols (Goetz and Dubayah, 2011; Angelsen et al., 2018; Kumar et al., 2022). In Pakistan, the MRV framework is operationalized through the National Forest Monitoring System (NFMS), which includes satellite monitoring (SLMS), a National Forest Inventory (NFI), and GHG accounting. While satellite imagery effectively detects deforestation, measuring forest degradation remains a challenge. Community-based forest inventories have proven valuable, enhancing both data quality and local participation (Holmes et al., 2017). Pakistan adopted UN-REDD in 2010, initiating its preparatory phase in 2012. Integrating satellite-based remote sensing with ground truthing is essential for achieving accurate and cost-effective biomass estimation, especially given the limited availability of local datasets (Hilmi et al., 2022). Recent advances in freely accessible remote sensing data, such as Sentinel-2 imagery have revolutionized carbon stock estimation. Traditional biomass estimation methods, including destructive (tree cutting) and non-destructive (measurement-based) approaches, are often costly and labor-intensive. Allometric equations now allow biomass estimation using tree DBH, height, and wood density, offering a practical alternative (Picard et al., 2012). Remote sensing complements this by providing large-scale, repeatable data. Sentinel-2, with its high-resolution multispectral sensors and unique red-edge and shortwave infrared bands, offers exceptional vegetation analysis capabilities (Ahmad et al., 2023). The red-edge band, in particular, enhances the accuracy of biomass modeling (Pertille et al., 2019). Given these advancements, this study aims to assess the carbon sequestration potential of three major gymnosperm species Pinus wallichiana, Abies pindrow, and Picea smithiana in Ayubia National Park. Specifically, it seeks to estimate AGB, BGB, and total biomass; also, to develop regression models linking field biomass data with Sentinel-2 imagery; map the spatial distribution of carbon across the park; and to identify priority zones for REDD+ implementation.

In particular, the Himalayan forests in Pakistan, dominated by long-lived gymnosperm species, have made significant contributions to the national carbon budget, but have not been adequately quantified as part of the national carbon budget. The rugged terrain and logistical difficulties of these mountain ecosystems are major obstacles to large-scale field-based inventories, which are often limited in scope and frequency (Shafique et al., 2021). This highlights the importance of integrating remote sensing. Although field inventories are indispensable and high-accuracy data are essential for calibration and validation, Sentinel-2 imagery can be used to map biomass changes across whole landscapes due to its spatial and temporal scalability (Fassnacht et al., 2021). The multispectral sensor on Sentinel-2 is particularly effective in heterogeneous mountainous terrain because of its high resolution (10 m) and red-edge bands, which are highly sensitive to chlorophyll content and canopy structure, enhancing the accuracy of biomass estimation in dense conifer forests (Chrysafis et al., 2017). Although protected areas such as Ayubia National Park play an important ecological role, there is a severe lack of spatially explicit evaluations of carbon stocks in these areas, which prevents the successful local adoption of national REDD+ policies and climate strategies (Siddiqui et al., 2020). Therefore, this research marks the first time a scientifically robust, integrated approach has been used to measure both above‑ and below‑ground carbon stocks of the park’s dominant gymnosperms in Ayubia National Park. By combining ground‑truth field measurements with Sentinel‑2 satellite imagery, the study validates field‑derived above‑ground biomass (AGB) estimates and shows that a step‑wise, multi‑index regression model dramatically improves carbon‑stock predictions compared with traditional single‑index or multi‑band techniques. The authors also delineate species‑specific biomass ranges, identify the most effective spectral predictors, and test how these results could fit into REDD+ carbon accounting. Altogether, the paper presents a scalable, precise methodology that could be applied to other high‑altitude coniferous forests across the Western Himalayas.

MATERIALS AND METHODS

Study area: This study was conducted in Ayubia National Park, located in the Abbottabad District of Khyber Pakhtunkhwa, Pakistan. The region lies entirely within the Lesser Himalayas. Its northernmost extent is 34° 01′ N, and its southernmost point is 34° 03.8′ N. Longitudinally, it stretches from 73° 22.8′ E to 73° 27.1′ E. Elevations range from about 1 450 m up to 3 033 m above sea level, giving the area a steep, mountainous character typical of the lower Himalayan foothills (Afza et al., 2016). The park is within the western Himalayan temperate zone and is characterized by mixed coniferous forests. Dominant gymnosperm species include, Blue Pine (P. wallichiana), Silver Fir (A. pindrow), and Spruce (P. smithiana). The park serves as a significant carbon sink due to its dense forest cover, making it a priority site for climate change mitigation and carbon sequestration studies.

Figure 1. Map of Galiyat showing the research area from Ayubia to Mushkpuri Top

Research design: A total of 63 circular plots (0.1 ha each; 17.84 m radius) were systematically established across Ayubia National Park to capture variability in stand structure and species composition. Within each plot, tree diameter at breast height (DBH) and height were recorded for all individuals of the three dominant conifer species (P. wallichiana, A. pindrow, and P. smithiana). A synthesis-based approach was adopted, integrating field inventory data with satellite-based remote sensing to estimate and map carbon stocks. The methodology included establishing sample plots for field data collection, applying species-specific allometric equations for biomass calculation, using Sentinel-2 imagery for vegetation analysis, and conducting correlation analysis between field-derived carbon stocks and spectral vegetation indices (particularly NDVI). This integrative approach provides a cost-effective, scalable solution for estimating forest biomass across large and heterogeneous landscapes.

Remote sensing data acquisition and pre-processing: Cloud-free Level-1C Sentinel-2A multispectral images were obtained from the Copernicus Open Access Hub. Imagery acquisition targeted the peak growing season to capture maximum vegetative activity. Sentinel-2 data, with high spatial resolution and unique red-edge bands, is well suited for forest biomass studies (Khan et al., 2020). Images were atmospherically corrected to Level-2A (bottom-of-atmosphere reflectance) using the Sen2Cor processor in the SNAP software suite, minimizing atmospheric effects such as scattering and absorption. All selected bands were resampled to 10-meter resolution to ensure uniformity across datasets. This preprocessing is standard for quantitative analysis of satellite imagery (Ali et al., 2023).

Vegetation index calculation: The Normalized Difference Vegetation Index (NDVI) was calculated using the red (Band 4) and near-infrared (NIR; Band 8) bands of Sentinel-2. NDVI is a widely used proxy for vegetation health and density and has shown strong correlations with AGB in similar ecosystems (Ali et al., 2023).

The generated NDVI map was used to correlate satellite-derived vegetation indicators with field-based biomass measurements (Khan et al., 2020).

Field data collection: A stratified random sampling design was employed, following REDD+ guidelines for national forest inventories (REDD+ Pakistan, 2019). Forest types, such as coniferous and mixed, were identified using Sentinel-2 imagery and stratified accordingly. Each circular sample plot had a radius of 17.84 meters, covering an area of 0.1 ha area. Within each plot, all trees with a DBH greater than 5 cm at 1.3 meters above ground level were measured. A diameter tape and clinometer were used to record DBH and total tree height, respectively. GPS devices captured the geographic coordinates of plot centers. Wood density values for P. wallichiana, A. pindrow, and P. smithiana were derived from core samples using an increment borer (Anwar, 2015). These species-specific values were necessary for accurate biomass calculations (IPCC, 2006).

Biomass estimation: Destructive sampling for biomass estimation is often impractical, so species-specific allometric equations were used to estimate AGB using DBH and height (Ismail et al., 2024). Species-Specific Allometric Equations are summarized in Table 1.

Table 01. Species-specific allometric biomass and height equations used for estimating above-ground biomass (AGB) and tree height of dominant conifer species (P. wallichiana, P. smithiana, and A. pindrow) based on diameter at breast height (D) and total tree height (H).

Where: Y = biomass (kg), D = DBH (cm), H = tree height (m), ln = natural logarithm

Below Ground Biomass (BGB) was estimated using the IPCC default root-to-shoot ratio of 0.26 for temperate coniferous forests (Ali et al., 2023):

BGB = AGB × 0.26

BGB = AGB * 0.26

Carbon stock calculation: Carbon stock was calculated by multiplying the total biomass (AGB + BGB) by a carbon fraction of 0.47, as per standard practice (Khan et al., 2020).

Carbon Stock (t C/ha) = (AGB + BGB) × 0.47

This provided an estimate of carbon stocks per plot, later extrapolated to a per-hectare basis.

Relationship between NDVI and Carbon stock: Linear and nonlinear regression models were developed to quantify the relationship between NDVI values (from Sentinel-2) and field-derived carbon stock. Model performance was evaluated using R² (coefficient of determination) and RMSE (root mean square error), consistent with recent remote sensing literature (Bhatti et al., 2023).

Statistical analysis: Descriptive statistics including means, standard deviations, and ranges, were calculated for DBH, height, biomass, and carbon stock. A subset of the field plots not used in model training was used for validation. The validated models were then applied to NDVI maps to generate spatially explicit carbon maps for Ayubia National Park, providing essential input for REDD+ planning and forest conservation strategies. One-way ANOVA, followed by post-hoc Tuckey's test was used to analyze the means of different groups. Data normality was confirmed by the KS test.

RESULTS

Plot wise topographic variables: Table 2 summarizes the topographic variables for the study plots. Elevation ranged from 1,490 m to 2,534 m, with a mean of 2,362.1 m and a SD of 64.37. Aspect values ranged from 1.72° to 361.43°, with a mean of 183.58° (SD = 22.23). Slope ranged from 10.72° to 49.74°, with a mean of 31.12° (SD = 2.09). These variations reflect the diverse terrain of Ayubia National Park.

 Table 2. Summary statistics of topographic variables (n = 63)

Tree attributes of Gymnosperm species in Ayubia National Park Galyat: As summarized in Table 3, P. wallichiana had the highest mean DBH (44.95 cm) and height (24.7 m), followed by A. pindrow and P. smithiana. The maximum DBH and height recorded for each species were 74 cm and 35.8 m (P. wallichiana), 72.13 cm and 34.65 m (A. pindrow), and 70.45 cm and 32 m (P. smithiana), respectively. A total of 8,591 stems were recorded for each species.

Table 3. Summary statistic of DBH, Height and Stem Number of P. wallichiana, A. pindrow and P. smithiana (n = 8591)

DBH= Diameter at breast height

Biomass and Carbon stock of P. wallichiana: Table 4 provides detailed statistics for P. wallichiana. The mean AGB was 67.2 t/ha (SD = 46.53), while AGC was 31.58 t/ha. BGB and BGC had mean values of 17.47 t/ha and 8.21 t/ha, respectively. The total biomass reached 247.12 t/ha, with a mean of 84.67 t/ha. CO₂ equivalent emissions ranged from 10.62 t/ha to 425.1 t/ha. Thus CO₂e was obtained by converting carbon stock using the standard molecular weight ratio of CO₂ to carbon (44/12 = 3.67). Accordingly, CO₂ equivalent emissions were calculated as:

CO₂e (t/ha) = Carbon stock (t C/ha) × 3.67

Table 4. Summary of plot wise biomass (t/ha) and carbon stocks (t/ha) of P. wallichiana

Tree attributes vs. AGB (P. wallichiana): Figure 2 presents the relationship between tree attributes and above-ground biomass (AGB) of Pinus wallichiana. The analysis revealed a moderate positive correlation between diameter at breast height (DBH) and AGB (r = 0.68), with a coefficient of determination (R² = 0.4593), indicating that approximately 45.9% of the variation in AGB is explained by DBH. In contrast, tree height showed a stronger positive correlation with AGB (r = 0.76) and a higher coefficient of determination (R² = 0.58), suggesting that about 58% of the variability in AGB is explained by tree height. These results indicate that tree height is a better predictor of above-ground biomass than DBH for Pinus wallichiana in the studied forest stands.

Figure 2. Relationship between tree attributes and AGB (P. wallichiana). The figure illustrates a moderate positive correlation between diameter at breast height (DBH) and AGB (r=0.68, R² = 0.4593), while a stronger relationship is observed between tree height and AGB (r=0.76, R² = 0.58). These results indicate that tree height serves as a more reliable predictor of aboveground biomass than DBH for P. wallichiana in the studied area.

Biomass and Carbon stock of A. pindrow: As presented in Table 5, A. pindrow had an AGB of 53.68 t/ha (SD = 38.28) and AGC of 25.23 t/ha. BGB and BGC were 13.96  t/ha and 6.56 t/ha, respectively. The mean total biomass was 67.64 t/ha, corresponding to an average carbon stock of 31.79 t/ha and a CO₂ equivalent of 116.44 t/ha.

Table 5. Plot wise biomass (t/ha) and carbon stocks (t/ha) of A. pindrow

Tree attributes vs. AGB (A. pindrow): The regression analysis of tree attributes diameter at breast height (DBH) and height (H)-against above-ground biomass (AGB) for Abies pindrow revealed positive correlations for both parameters (Figure 3). Tree height demonstrated a stronger predictive relationship with AGB (r = 0.74, R² = 0.55) compared to DBH (r = 0.69, R² = 0.48). The coefficient of determination (R²) indicates that height explains approximately 55% of the variability in AGB, whereas DBH accounts for 48%. This suggests that vertical growth is a more reliable indicator of biomass accumulation than radial growth in this species. This observation is also supported by the accompanying bar graphs, which indicate that changes in AGB are strongly related to changes in tree height among sampled individuals. Both DBH and height are valid predictors of AGB, but height has greater explanatory power and should be prioritized in biomass models. However, since neither attribute alone explains a significant proportion of variance, future models for AGB estimation of A. pindrow should incorporate multiple regression using both DBH and height to enhance predictive efficiency and minimize estimation errors.

Figure 3. Relationship between key tree structural attributes and above-ground biomass (AGB) of Abies pindrow. The figure illustrates how variations in tree attributes such as diameter at breast height (DBH), tree height, and crown dimensions’ influence AGB accumulation. Each plotted relationship highlights the strength and direction of association between individual tree characteristics and biomass, providing insights into growth patterns and biomass allocation strategies of A. pindrow. These relationships are critical for improving biomass estimation accuracy and understanding forest structure and productivity in alpine and sub-alpine ecosystems

Biomass and Carbon stock of Picea smithiana: Table 6 highlights that Picea smithiana had an AGB of 52.85 t/ha (SD = 34.99), AGC of 24.84 t/ha, BGB of 13.74 t/ha, and BGC of 6.46 t/ha. The total mean biomass and carbon stock were 66.59 and 31.3 t/ha, respectively. CO₂ emissions reached 378.5 t/ha.

Table 6. Plot wise biomass (t/ha) and carbon stocks (t/ha) of Picea smithiana

Tree attributes vs. AGB (P. smithiana): The relationship between major tree structural attributes and above-ground biomass (AGB) of Picea smithiana is presented in Figure 4. The analysis revealed a strong positive association between AGB and key dendrometric parameters, including diameter at breast height (DBH), total tree height, and crown characteristics. Trees with larger DBH values showed a substantial increase in biomass accumulation, indicating that stem diameter is a primary determinant of above-ground biomass. As DBH increased, the estimated AGB also increased consistently, suggesting that mature trees contribute significantly to the overall biomass stock of the forest stand. Similarly, total tree height exhibited a positive relationship with AGB. Taller individuals accumulated greater biomass compared to shorter trees, reflecting the progressive growth and structural development of P. smithiana with age. This trend indicates that height growth plays an important role in biomass production and reflects the species’ capacity to occupy vertical canopy space in montane forest ecosystems. Crown characteristics also noticeably influenced biomass accumulation. Trees with broader, well-developed crowns were associated with higher AGB values. Larger crown dimensions may enhance photosynthetic capacity, leading to greater carbon assimilation and subsequent biomass formation. The relationship between crown attributes and AGB highlights the importance of canopy structure in determining the productivity and growth performance of P. smithiana. Overall, the observed patterns demonstrate that variations in individual tree structural attributes significantly influence biomass accumulation. Among the studied parameters, DBH showed the strongest relationship with AGB, followed by tree height and crown traits. These findings provide valuable insights for improving biomass estimation models and contribute to a better understanding of forest structure, productivity, and carbon storage potential in montane and sub-alpine ecosystems where P. smithiana is dominant.

Figure 4. Relationship between major tree structural attributes and above-ground biomass (AGB) of Picea smithiana. The figure depicts the associations between AGB and key dendrometric variables, including diameter at breast height (DBH), total tree height, and crown characteristics. The observed trends demonstrate how variations in individual tree attributes contribute to biomass accumulation, reflecting species-specific growth behavior and biomass allocation patterns of P. smithiana. These relationships provide an empirical basis for improving biomass estimation models and for understanding forest structure, productivity, and carbon storage potential in montane and sub-alpine forest ecosystems.

Total biomass and Carbon stock estimation: Tables 7 and 8 summarize total biomass and carbon stock across Ayubia National Park. The mean total biomass ranged from 201.72 to 208.31 t/ha, with carbon stock averaging between 94.81 and 97.83 t/ha. High variability was observed across plots, with some exceeding 393.61 t/ha in biomass and 185 t/ha in carbon. The total carbon stock across the park was estimated at 5,878.08 tons.

Table 7. Summary statistics of above-ground biomass (AGB), below-ground biomass (BGB), total biomass, and associated carbon stocks across all sampled gymnosperm plots in Ayubia National Park

The table presents mean values, variability measures (standard error, standard deviation), and distribution metrics (minimum, maximum, range), along with total cumulative biomass and carbon (t/ha) derived from 63 field plots.

 Table 8. Summary of plot wise total biomass (t/ha) and total carbon stocks (t/ha) of Ayubia national park

Temporal changes in biomass and Carbon (1988-2023): Table 9 reveals a substantial decrease in carbon stocks of coniferous species from 1988 to 2023. For example, Abies pindrow carbon stocks declined by nearly 40%. The overall forest area decreased from 3,320 ha to 2,945 ha (department data), with a corresponding loss observed in Landsat imagery (Table 10). Carbon density also declined from 827 t/ha to 684 t/ha.

 Table 9.  Species-wise carbon stock estimation for Galies Forests in 1988 and 2023 based on wood density (WD), biomass expansion factor (BEF), and stand volume

Parameters include wood density (WD), biomass expansion factor (BEF), stand volume (m³), total biomass (t), and carbon stocks (t) for fir, kail, spruce, and broad-leaved species. The table illustrates historical changes in forest structure and carbon sequestration capacity.

The analysis of carbon stocks associated with forest cover change in Ayubia National Park from 1988 to 2023 indicates a substantial decline in carbon storage due to deforestation (Table 10). In 1988, the Landsat-derived forest area was 2,052 ha, corresponding to an estimated 2,748,441 tons of carbon stock with an average carbon density of 827 t/ha. By 2023, the forest area had decreased to 1,681 ha, reducing the total carbon stock to 2,023,296 tons and the average carbon density of 684 t/ha. Overall, the decline in forest cover led to a loss of approximately 725,145 tons of carbon stock over the 35-year period. The reduction in carbon density per hectare further indicates degradation of forest biomass and structural composition. The Landsat-based analysis shows a forest area loss of 371 ha, which directly contributed to the decrease in carbon storage capacity. These findings highlight the significant role of deforestation on reducing the carbon sequestration potential of Ayubia National Park. The loss of forest biomass affects ecosystem stability and contributes to increased atmospheric carbon emissions. Restoring degraded forest areas and implementing REDD+ strategies could play a crucial role in enhancing carbon storage and mitigating climate change impacts within the park.

 Table 10. Carbon emissions from Deforestation (1988-2023)

Sentinel-2 spectral indices and biomass estimation: Regression results (Figure 5) showed that ARVI had the highest individual correlation with AGB (R² = 0.6994), followed by NDVI (R² = 0.622) and MSAVI (R² = 0.6089). Individual spectral bands (Figure 6) performed poorly, with Band 4 being the most useful (R² = 0.21). The correlation matrix (Figure 7) indicated strong multicollinearity between indices like NDVI, SAVI, and MSAVI. Table 11 indicates that the stepwise multi-index regression model achieved the highest accuracy (RMSE = 20.22), while models using individual bands or backward regression performed poorly.

Figure 5. Stepwise regression analysis of above-ground biomass (AGB) using Sentinel-2–derived spectral indices. The figure presents the results of a stepwise multiple regression model illustrating the contribution of selected spectral indices to AGB estimation. Indices derived from key Sentinel-2 bands (including visible, red-edge, and near-infrared regions) were progressively included based on their statistical significance, highlighting the most influential predictors of biomass variability. The model output demonstrates the relative importance of individual spectral indices and their combined explanatory power, underscoring the effectiveness of multispectral remote sensing data for accurate, large-scale estimation of forest AGB and carbon stocks.

The figure below, as shown in Figure 6, indicates that simple linear regressions between AGB and individual Sentinel-2 spectral bands were mostly weak predictors across the spectrum. Among the twelve bands, Band 4 (Red) was most closely related to AGB, with a correlation of about 21 percent (R² = 0.21). Other bands with slight but relatively low correlations included Band 1 (Coastal Aerosol, R² = 0.1686), Band 9 (Water Vapour, R² = 0.1502) and Band 8 (Near-Infrared, R² = 0.1499). Contrastingly, the remaining bands performed very poorly with the shortwave-infrared (SWIR) bands, B11 and B12, being the least effective predictors, showing negligible correlations with AGB (R² = 0.0157 and R² = 0.0104, respectively). These findings show that a single spectral band cannot adequately estimate AGB; instead several spectral bands need to be considered to achieve statistically significant results.

Figure 6. Regression analysis of above-ground biomass (AGB) with Sentinel-2 spectral bands. The scatter plots show the relationship between AGB (kg/ha) and the Sentinel-2 spectral bands (B1 to B12). Each plot includes the regression equation and R² values, indicating the strength of the correlation between AGB and the spectral bands.

As seen in Figure 7 the correlation matrix indicates the interconnections between Total Above-Ground Biomass (TAGB) and various vegetation indicators. Several key indices are strongly and positively correlated with TAGB, with the strongest correlations observed for the Normalized Difference Vegetation Index (NDVI), Soil-Adjusted Vegetation Index (SAVI), Modified Soil-Adjusted Vegetation Index (MSAVI) and Atmospherically Resistant Vegetation Index (ARVI). The big dark blue circles indicate these relationships and imply high, statistically significant correlation coefficients. Moderate positive correlations are also indicated in other indices, including the Difference Vegetation Index (DVI) and Ratio Vegetation Index (RVI). On the other hand, indices such as the Red Edge Inflection Point (REIP) were observed to be insignificant indicators of TAGB, as indicated by cross marks. There is also a very high level of multicollinearity between the indices themselves; e.g. NDVI, SAVI and MSAVI are all extremely strongly correlated. This suggests that although several of these indices are effective standalone predictors of biomass, careful variable selection is necessary in regression analysis to avoid redundancy.

Figure 7. Correlation of AGB with various vegetation indices. The heatmap shows the correlation between AGB (TAGB t ha^-1) and different vegetation indices (NDVI, SAVI, MSAVI, etc.), with the color intensity representing the strength of the correlation. The symbols (X) indicate the indices that show no significant correlation with AGB, while the size of the circles reflects the degree of correlation. The color bar on the right shows the correlation scale from -1 (negative correlation) to 1 (positive correlation).

Table 11 showed that the analysis of different regression models for estimating The AGB showed a huge correlation in terms of predictability, as identified by the Root Mean Square Error (RMSE). The stepwise linear regression, which combined vegetation indices such as ARVI, MSAVI, NDVI, and DVI, was the most effective, followed by the standard linear regression model. The RMSE for the stepwise linear regression model is 20.22, indicating the highest accuracy among all approaches tested. In contrast the other models, did not perform well. The linear regression model had a much larger error (RMSE = 39.712). The models using spectral bands alone, as well as the complex backward stepwise regression using indices, had the worst errors, with RMSE values of 79.933 and 87.14 respectively. These findings clearly indicate that a curated multiple linear regression model using vegetation indices is the most reliable method for AGB estimation in this study.

 Table 11. Regression models for estimating above-ground biomass

An Analysis of Variance (Table 12) was done to establish whether there were significant differences in group means for various measured parameters. The findings showed that a number of parameters were greatly influenced. In particular, statistically significant differences (p < 0.05) were observed in Biomass 1 (F(3, 6) = 1.33, p = 0.035), SAVI (F(3, 6) = 23.86, p = 0.001) and in DVI (F(3, 6) = 0.19, p = 0.039), PVI (F(3, 6) = 2.44, p = 0.01 The post-hoc test results are in the form of the superscript letters (e.g., d in case of SAVI, or bc in the case of RENDVI, REIP, and RVI), meaning that the means of the groups that do not have a common letter differ significantly between each other. On the other hand, NDVI, MSAVI, ARVI, RENDVI1, RENDVI2, and S2REP did not show statistically significant differences at p > 0.053), 0.060), 0.080), 0.066), and 0.050 respectively. In the case of a parameter such as NDVI and S2REP, the post-hoc letter a means that the means of all groups were statistically equal.

Table 12. Analysis of variance (ANOVA) of Sentinel-2 derived spectral indices in relation to biomass estimation

Spatial biomass mapping: The spatial distribution of AGB across the study area is illustrated in Figure 8. The results derived from multiple remote sensing datasets, including Sentinel-2, Landsat-8, WorldView-2 imagery, LiDAR canopy height data, and Hansen Global Forest Cover datasets, revealed clear spatial variability in biomass distribution. Higher AGB values were predominantly observed in the northern valleys and on north-facing slopes, where dense vegetation cover and well-developed forest stands contributed to increased biomass accumulation. In contrast, lower AGB values were detected in areas with sparse vegetation and reduced canopy cover. The spatial patterns derived from Sentinel-2 and Landsat-8 imagery demonstrated similar distribution trends, indicating the reliability of medium-resolution satellite data for biomass estimation. The WorldView-2 imagery provided more detailed spatial information and highlighted localized variations in biomass distribution within forest patches. The LiDAR-derived forest canopy height model supported these findings, showing a strong correspondence between taller canopy structures and higher AGB values. Regions with greater canopy height were associated with higher biomass density, confirming the importance of vertical forest structure in biomass estimation. Additionally, the Hansen Global Forest Cover dataset showed a spatial pattern consistent with the biomass distribution derived from other datasets. Areas with dense forest cover corresponded with zones of elevated AGB, validating the accuracy of the remote sensing-based biomass mapping approach. Overall, integrating multiple remote sensing datasets enabled a comprehensive assessment of biomass distribution across the landscape. The observed spatial patterns indicate that topographic factors, particularly slope orientation and valley locations, significantly influence biomass accumulation in the study area. These results demonstrate the effectiveness of multi-source remote sensing data for accurate large-scale biomass estimation and for assessing forest carbon storage potential.

Figure 8. Spatial mapping of AGB in the study area. The maps show AGB distribution from different data sources: Sentinel-2, LiDAR Forest Canopy Height, Landsat-8, Hansen Global Forests Cover 2019, and WorldView-2, with a Google Earth image for reference. AGB values are represented by a color scale, with red indicating high AGB and green indicating low AGB.

Deforestation trends and REDD+ implications: The Landsat-derived NDVI analysis revealed a clear decline in the reserved forest area of Ayubia National Park between 1988 and 2023 (Figure 9). The spatial distribution of forest cover shows that the reserved forest area decreased from 2,052 ha in 1988 to 1,681 ha in 2023, indicating a net loss of 371 ha over 35 years. This corresponds to an approximate 18.1% reduction in forest cover within the study area. The comparative maps illustrate that forested regions (represented in green) were more continuous and dense in 1988, whereas the 2023 map shows increased fragmentation and the emergence of non-forest patches across several parts. The reduction in forest cover is particularly noticeable in the central and southern portions of the park, where previously continuous forest areas have become fragmented. These newly identified blank or degraded areas represent potential sites for forest restoration and conservation initiatives. From a management perspective, these zones can be considered priority areas for REDD+ (Reducing Emissions from Deforestation and Forest Degradation) interventions, as restoring forest cover here could support carbon sequestration and biodiversity conservation. The spatial identification of these areas aligns with the forest department’s designated working circles, which aim to promote sustainable forest management within the park.

Figure 9. Changes in reserved forest area of Ayubia National Park (1988–2023). The maps show the reserved forest areas in Ayubia National Park for 1988 (left) and 2023 (right). The green areas represent reserved forests, while the blank areas are potential sites for REDD+ implementation. The forest area decreased by 371 hectares, from 2052 ha in 1988 to 1681 ha in 2023.

DISCUSSION

The mean tree density of 560 stems/ha in this study aligns with previous findings in Himalayan gymnosperm forests (Khan et al., 2024; Amir et al., 2018). Differences in tree density likely reflect site-specific forest composition, past management, and DBH cutoff thresholds (Nizami, 2012). Maximum DBH and height values for P. wallichiana, A. pindrow, and P. smithiana closely match those reported by Shaheen et al. (2016). The highest AGB values observed (312 t/ha) compare favorably with regional records (Singh et al., 2012; Sharma et al., 2016). NDVI was the best single predictor of biomass (R² = 0.622), consistent with results by Pandit et al., (2018) and others. However, the stepwise multi-index model substantially improved accuracy (R² = 0.915), confirming findings by Lu and Kang (2008), Lu (2005), Salako et al., (2025), and Imran et al., (2020). The limitations of NDVI in dense forests were addressed by combining it with MSAVI, ARVI, and DVI. This study validates the integration of Sentinel-2 imagery with regression models for biomass estimation in complex terrain. Though advanced tools like LiDAR or machine learning (e.g. Random Forests) were not employed here, they present a future opportunity (Imran et al., 2021; Khan et al., 2023).

As a result, this structural maturity is indicated by high above-ground biomass (AGB) values, with a maximum of about 312 t/ha. These high carbon stocks are characteristic of forests with large, long-lived gymnosperms, which serve as important long-term carbon sinks (Singh et al., 2012; Sharma et al., 2016). The findings highlight the immense importance of the park as a means of conserving biodiversity and as a nationally important carbon reservoir, whose protection is crucial for achieving Pakistan’s climate mitigation goals.

From a methodological perspective, this research offers critical information on biomass estimation using Sentinel-2 data in complex mountainous surfaces. Although NDVI was the best single predictor (R² = 0.622), NDVI cannot be utilized in densely canopied areas because the signal becomes saturated and does not respond to changes in high-biomass regions (Mutanga and Skidmore, 2004). The most valuable finding is the substantial increase in predictive power achieved by the stepwise multi-index model (R² = 0.915). This improvement results from the synergistic incorporation of other indices that address NDVI’s limitations. As an example, the Modified Soil-Adjusted Vegetation Index (MSAVI) reduces background soil reflectance in smaller dense areas and the Atmospherically Resistant Vegetation Index (ARVI) corrects for aerosol scattering prevalent in mountainous ridges, enhancing the vegetation signal (Haboudane et al., 2002). The high accuracy of the model confirms that using multiple indices tailored to certain environmental and canopy-related challenges offers a strong and economical alternative to single index models for carbon accounting, as supported by analogous multi-index evaluations in other forested environments (Salako et al., 2025; Chrysafis et al., 2017; Imran et al., 2020).

The direct implications of this validated methodology affect forest management and policy. This method will enable the strategic implementation of REDD+ and other payment-for-ecosystem-services programs, as it is based on accurately quantifying biomass and determining regions that have been previously degraded. Although we have a highly useful regression model, further studies can improve its accuracy. By combining our spectral data with structural data from active sensors, such as LiDAR or GEDI mission data, we could better estimate biomass by directly measuring the vertical structure of forests (Imran et al., 2021). Similarly, using machine learning models such as Random Forests, which capture complex non-linear relationships between spectral data and field measurements, is another promising approach for controlling carbon in the area (Khan et al., 2023).

Conclusion: This study establishes that Ayubia National Park is a significant carbon source, with above-ground biomass reaching up to 312 t/ha in its gymnosperm-dominated forests. Our results show that a stepwise multi-index regression framework based on Sentinel-2 data (R² = 0.915) is an excellent predictor of biomass and serves as a highly valid and useful tool for carbon accounting in ecosystems. Despite this high capacity for carbon sequestration, there was a huge loss of forest cover around 18 percent between 1988 and 2023. Notably, the areas where forest cover was lost are those where REDD + interventions would be most feasible. The study confirms that combining field inventory with multispectral satellite data offers a replicable and affordable model to enhance Pakistan’s national forest monitoring capacity. In the long term, it is essential to scale-up this tested methodology across various forest ecosystems in Pakistan to monitor carbon dynamics accurately and promptly. Future research should incorporate data from more sophisticated sensors, such as LiDAR, to further refine biomass models and strengthen the country’s contribution to global climate change mitigation.

Acknowledgements: NA

Author Contributions: Conceptualization, MG, AU; Data curation, MG, AU; Formal analysis, MG, MNK; Investigation, MG; Methodology, MG; Project administration, AU; Resources, MG, MNK; Software, MG, MNK; Validation, MG, AU; Writing – original draft, MG; Writing-review & editing, AU. All authors contributed significantly, have read and agreed to the published version of the manuscript.

Funding: NA

Data Availability: All data is available either in the article or in the literature cited.

Conflicts of Interest: All authors declare that they have no conflicts of interest.

REFERENCES

Ahmed, M., M. Wahab, N. Khan, J. Palmer, K. Nazim, M. U. Khan and M. F. Siddiqui (2010). Some preliminary results of climatic studies based on two Pine tree species of Himalayan area of Pakistan. Pak. J. Bot. 42(2): 731–738. http://www.pakbs.org/pjbot/PDFs/42(2)/PJB42(2)0731.pdf

Ali, S.R. and N. Mujahid (2025). Agricultural productivity under climate change vulnerability: does carbon reduction paths matter for sustainable agriculture. Environ. Dev. Sustain, 1–25. https://doi.org/10.1007/s10668-025-06076-9

Afza, R., H. Ahmad, Z. Saqib, N. Ali and J. Khan (2016). Phytodiversity of Ayubia national park Pakistan: Conservation and management issues. J. Biodiversity Environ. Sci. 8(2): 327-336.

Barber, G., A. Edwards, K. Zander (2023). Fire, rain and CO₂: potential drivers of tropical savanna vegetation change, with implications for carbon crediting. Fire 6(12): 465. https://doi.org/10.3390/fire6120465

Salako, J., N. Millar, A. Kendall, B. Basso (2025). Assessing tree root distributions using ground penetrating radar and machine learning algorithms. Agrosystems, Geosciences & Environment. 8(4): 1-22. https://doi.org/10.1002/agg2.70217.

Bhatti, U.A., Z. Yu, A. Hasnain, S. A. Nawaz, L. Yuan, L. Wen, M. A. Bhatti (2022). Evaluating the impact of roads on the diversity pattern and density of trees to improve the conservation of species. Environ. Sci. Pollut. Res. 29(10): 14780–14790. https://doi.org/10.1007/s11356-021-16627-y

Chrysafis, I., G. Mallinis, I. Gitas, M. Tsakiri-Strati (2017). Assessing the relationships between growing stock volume and Sentinel-2 imagery in a Mediterranean forest ecosystem. Remote. Sens. Appl. Soc. Environ. 8: 95–105. https://doi.org/10.1080/2150704X.2017.1295479

Dar, J. A., and S. Sundarapandian (2015). Variation of biomass and carbon pools with forest type in temperate forests of Kashmir Himalaya, India. Environ. Monit. Assess. 187(2): 55. https://doi.org/10.1007/s10661-015-4299-7

Din, I.U., S. Muhammad, I.U. Rehman and C. Tokatli (2023). Spatial distribution of potentially toxic elements contaminations and risk indices of water and sediments in the Darband and Samana streams, Pakistan. Environ. Monit. Assess. 195(11): 1343. https://doi.org/10.1007/s10661-023-11914-2

Goetz, S. and R. Dubayah (2011). Advances in remote sensing technology and implications for measuring and monitoring forest carbon stocks and change. Carbon Manag. 2(3): 231–244. https://doi.org/10.4155/cmt.11.18

Hilmi, E., R. Dewi, E. Sudiana, A. Mahdiana, L.K. Sari and T.N. Cahyo (2022). The clustering and distribution of heavy metal accumulation and translocation as an ability of mangrove vegetation to reduce impact of pollution. 30 March 2022, PREPRINT (Version 1) available at Research Square. https://doi.org/10.21203/rs.3.rs-1445862/v1

Holmes, I., K.R. Kirby and C. Potvin (2017). Agroforestry within REDD+: experiences of an indigenous Emberá community in Panama. Agrofor. Syst. 91(6): 1181–1197. https://doi.org/10.1007/s10457-016-0003-3

Imran, A.B., S. Ahmed, W. Ahmed, M. Zia-ur-Rehman, A. Iqbal, N. Ahmad and I. Ullah (2021). Integration of Sentinel-2 derived spectral indices and in-situ forest inventory to predict forest biomass. Pak. J. Sci. Ind. Res. Ser A. Phys. Sci. 64(2): 119–130.

Imran, A.B., K. Khan, N. Ali, N. Ahmad, A. Ali and K. Shah (2020). Narrow band based and broadband derived vegetation indices using Sentinel-2 imagery to estimate vegetation biomass. Glob. J. Environ. Sci. Manage. 6(1): 97–108. https://doi.org/10.22034/gjesm.2020.01.08

Ismail, S.N.H., N. Rusli and N. Jamirsah (2024). Optimizing carbon capture in urban forests through strategic tree selection. (November 22, 2024). Available at SSRN: http://dx.doi.org/10.2139/ssrn.4872009.

Khan, M.N.I., M.R. Islam, A. Rahman, M.S. Azad, A.S. Mollick, M. Kamruzzaman and A. Knohl (2020). Allometric relationships of stand level carbon stocks to basal area, tree height and wood density of nine tree species in Bangladesh. Glob. Ecol. Conserv. 22: e01025. https://doi.org/10.1016/j.gecco.2020.e01025

Khan, P., K. Ullah, B. M. Hung, A. Ahmad, S. Y. Lee, L. D. Hai and R. N. Muhammad (2024). Assessing the biomass and carbon productivity of pure Pinus wallichiana forest in Kumrat Valley, Pakistan. Biol. Bull. 51: 1863–1876. https://doi.org/10.1134/S1062359024609133

Lu, C. and J. Kang (2008). Generation of transgenic plants of a potential oilseed crop Camelina sativa by Agrobacterium-mediated transformation. Plant. Cell. Rep. 27(2): 273–278. https://doi.org/10.1007/s00299-007-0454-0

Lu, H. (2005). A novel high-order tree for secure multicast key management. IEEE. Trans. Comput. 54(2): 214–224. https://doi.org/10.1109/TC.2005.15

Pandit, S., S. Tsuyuki and T. Dube (2018). Estimating above-ground biomass in sub-tropical buffer zone community forests, Nepal, using Sentinel-2 data. Remote. Sens. 10(4): 601. https://doi.org/10.3390/rs10040601

Pertille, R.H., M.R. Sachet, M.T. Guerrezi and I. Citadin (2019). An R package to quantify different chilling and heat models for temperate fruit trees. Comput. Electron. Agric. 167: 105067. https://doi.org/10.1016/j.compag.2019.105067

Picard, N., L. Saint-André and M. Henry (2012). Manual for building tree volume and biomass allometric equations: from field measurement to prediction. FAO. pp 1-215. http://www.fao.org/docrep/018/i3058e/i3058e.pdf

Raihan, A. and A. Tuspekova (2022). Dynamic impacts of economic growth, energy use, urbanization, tourism, agricultural value-added, and forested area on carbon dioxide emissions in Brazil. J. Environ. Stud. Sci. 12: 794–814. https://doi.org/10.1007/s13412-022-00782-w

Sun, W. and X. Liu. (2019). Review on carbon storage estimation of forest ecosystem and applications in China. Forest Ecosys. 7: 4. https://doi.org/10.1186/s40663-019-0210-2

Fassnacht, F.E., J. Poblete-Olivares, L. Rivero, J. Lopatin, A. Ceballos-Comisso, and M. Galleguillos (2021). Using Sentinel-2 and canopy height models to derive a landscape-level biomass map covering multiple vegetation types. Int. J. Appl. Earth Observ. Geoinform 94: 102236. https://doi.org/10.1016/j.jag.2020.102236

Shafique, T., M.H. Zuberi and Z.I. Shams (2021). Geospatial assessment of carbon stock inventory by vegetation indices in Pai Forest, Sindh, Pakistan. EQA – Int. J. Environ. Qual. 43: 47–64. https://doi.org/10.6092/issn.2281-4485/12203

Shaheen, H., R. W. A. Khan, K. Hussain, T. S. Ullah, M. Nasir and A. Mehmood (2016). Carbon stocks assessment in subtropical forest types of Kashmir Himalayas. Pak. J. Bot. 48(6): 2351–2357. http://www.pakbs.org/pjbot/PDFs/48(6)/20.pdf

Sharma, R. P., Z. Vacek and S. Vacek (2016). Individual tree crown width models for Norway spruce and European beech in Czech Republic. For. Ecol. Manage. 366: 208–220. https://doi.org/10.1016/j.foreco.2016.01.040

Shin, M. Y., M. D. Miah and K. H. Lee (2007). Potential contribution of the forestry sector in Bangladesh to carbon sequestration. J. Environ. Manage. 82(2): 260–276. https://doi.org/10.1016/j.jenvman.2005.12.025

Siddiq, Z., M. U. Hayyat, A. U. Khan, R. Mahmood, L. Shahzad, R. Ghaffar and K. F. Cao (2021). Models to estimate the above and below ground carbon stocks from a subtropical scrub forest of Pakistan. Glob. Ecol. Conserv. 27: e01539. https://doi.org/10.1016/j.gecco.2021.e01539

Siddiqui, M. F., M. Ahmed, M. Iqbal and I. Khan (2020). Forest cover change and carbon stock assessment in the protected areas of Khyber Pakhtunkhwa, Pakistan. Pak. J. Agric. Sci. 57(3): 825–835. http://dx.doi.org/10.19045/bspab.2017.60089

Singh, I. S., O. P. Awasthi, R. S. Singh, T. A. More and S. R. Meena (2012). Changes in soil properties under tree species. Ind. J. Agric. Sci. 82(2): 146.

Stinson, B., A. Laguerre and E. T. Gall (2022). Per-person and whole-building VOC emission factors in an occupied school with gas-phase air cleaning. Environ. Sci. Technol. 56(6): 3354–3364. https://doi.org/10.1021/acs.est.1c06767

Thomas, S.C., and A.R. Martin (2012). Carbon content of tree tissues: a synthesis. Forests 3(2): 332–352. https://doi.org/10.3390/f3020332