Comparative study of forest biomass and carbon stocks of Margalla Hills National Park, Pakistan

Abstract Forests can play an important role in climate change mitigation. However, limited information is available worldwide regarding forest carbon and biomass stocks. Financial mechanisms such as ‘reducing emissions from deforestation and forest degradation and the role of conservation of forest carbon, sustainable management of forests and enhancement of forest carbon stocks’ (REDD+) also emphasize the quantification of forest biomass and carbon. This study aimed to estimate the forest biomass in two forests of Margalla Hills National Park (MHNP): Sub-tropical Chir Pine Forest (SCPF) and Sub-tropical Broadleaved Evergreen Forest (SBEF). For this, circular sampling plots of a 20 m radius were used for the collection of the variables, “diameter at breast height (DBH) and height”. Statistical analysis was done for exploring regression relationships between the variables. We found a mean Aboveground Carbon (AGC) of 73.36 ± 32.55 Mg C ha−1 in SCPF and a mean AGC of 16.88 ± 25.81 Mg C ha−1 in SBEF. The mean Aboveground Biomass (AGB) for SCPF was recorded as 146.73 ± 65.11 Mg ha−1, while for SBEF it was 33.77 ± 51.63 Mg ha−1. It was therefore concluded that the SCPF had higher mean AGB and mean AGC than the SBEF. Similar differences were also noticed in the structural characteristics of the two forests. These could be valuable information while designing sustainable management plans and afforestation programmes for the future and also for accessing nature-based funding such as REDD+.


Introduction
Climate change is one of the greatest threats of recent times (Po skus and Zukauskien_ e 2017; Soutter and Mõttus 2020). It is defined as the shift in climate patterns, which is mainly caused by greenhouse gas (GHG) emissions (Fawzy et al. 2020). CO 2 has mainly received widespread attention, among other GHG emissions, due to its increasing levels which have resulted in higher temperatures, shifting precipitation patterns, unpredictable extreme events, etc. (Reichstein et al. 2013). Globally, the CO 2 concentration has risen from approximately 277 parts per million (ppm) in 1750, which was the beginning of the Industrial Era, to 409.85 ± 0.1 ppm in 2019 (Joos and Spahni 2008;Dlugokencky and Tans 2020). Anthropogenic activities have been the main culprit causing climate change (Xi-Liu and Qing-Xian 2018). Deforestation and forest degradation, for instance, have led to an increase in GHG emissions from 490 ± 180 Gt CO 2 in 1970 to 680 ± 300 Gt CO 2 in 2010 (IPCC 2007). The severity of climate change can also be seen from the fact that the global surface temperature was 1.09 C higher in 2011-2020 than in 1850-1900(IPCC 2023. Aboveground Biomass (AGB), an important global measurement standard for climate change mitigation activities, generally, accounts for 70% to 90% of total forest biomass (FAO 2010). AGB is "all the biomass of living vegetation, both woody and herbaceous, above the soil including stems, stumps, branches, bark, seeds, and foliage" (FAO 2020a). It is usually measured in metric tonnes of dry matter per hectare (e.g. t ha À1 or Mg ha À1 ) or in metric tonnes of carbon per hectare (e.g. t C ha À1 or Mg C ha À1 ) (Omar et al. 2017;Rodr ıguez-Veiga et al. 2017). It can only be measured directly and accurately through destructive sampling, where the trees are cut down and dry weights are calculated using oven-drying methods, which are regressed against the volumes estimated using volumetric methods (Wilkes et al. 2018). However, this is labor intensive, expensive, often impractical, and rarely adopted (Dumitras , cu et al. 2020). Destructive sampling is also not recommended for places where there is threatened flora and fauna (Kumar and Mutanga 2017). Nondestructive methods include the use of already available allometric equations for the estimation of AGB without cutting down the trees (Rodr ıguez-Veiga et al. 2017; Debastiani et al. 2019). The allometric equations are statistical models based on the forest inventory data such as diameter at breast height (DBH), height, wood density, and forest type (Chave et al. 2004). The forest inventory data is one of the most reliable sources for carbon estimation globally (Zhao et al. 2019). The choice of the allometric equation is, however, critical, where an inconsiderate selection of the equation may lead to errors of up to 79% (Clark et al. 2001). For reduction in uncertainties in AGB estimation, site and species-specific allometric equations have been emphasized (Soares and Schaeffer-Novelli 2005).
The other nondestructive method, involving forest inventory data, includes the use of the Biomass Expansion Factors (BEF) and Wood Density (WD), which are used to convert the tree volume to the AGB (Asrat et al. 2020;Hossain et al. 2020). BEF is unitless and is defined as "the ratio of oven-dried tree biomass to either the commercial or stem oven-dried biomass including bark" (Schepaschenko et al. 2018). BEF approach is not flexible in terms of the use of various variables such as the allometric equations. In the BEF approach, only the tree volume, BEF, and wood density are used for the estimation of AGB, as mentioned before. Other forest carbon pools include; Belowground Biomass (BGB), deadwood, litter, and soil carbon (Buchholz et al. 2014). BGB is "all the biomass of live roots. Fine roots of less than 2 mm diameter are excluded because these often cannot be distinguished empirically from soil organic matter or litter" (FAO 2020a). Deadwood is "all the non-living woody biomass not contained in the litter, either standing, lying on the ground, or in the soil. It includes wood lying on the surface, dead roots, and stumps larger than or equal to 10 cm in diameter or any other diameter used by the country" (FAO 2020a). Litter includes "all non-living biomass with a diameter less than the minimum diameter for dead wood (e.g. 10 cm), lying dead in various states of decomposition above the mineral or organic soil" (FAO 2020a). Soil carbon is the "organic carbon in mineral and organic soils (including peat) to a specified depth chosen by the country and applied consistently through the time series" (FAO 2020a).
Generally, forest biomass is an indicator of its carbon status since carbon forms half of the dry weight of a tree (Vicharnakorn et al. 2014;Raihan et al. 2021). Agriculture, Forestry and Other Land Use (AFOLU) offers one of the largest potentials for GHG emissions reduction through reduced deforestation and forest degradation, amounting to 0.4-5.8 Gt CO 2 -eq year À1 (IPCC 2019). However, according to the fourth assessment report of the Intergovernmental Panel on Climate Change (IPCC) (IPCC 2007), there is a limited amount of information available regarding forest biomass and carbon status at both national and regional levels. For tackling climate change, large-scale implementation of mitigation responses would be required (IPCC 2023). The United Nations Framework Convention on Climate Change (UNFCCC) has also recognized the AGB as an Essential Climate Variable (ECV) (Omar et al. 2017). ECVs are the physical, biological, and chemical variables or groups of variables that critically affect the Earth's climate (Bojinski et al. 2014). The 'reducing emissions from deforestation and forest degradation and the role of conservation of forest carbon, sustainable management of forests and enhancement of forest carbon stocks' (REDDþ) is also a resultbased financial incentive for developing countries to counteract deforestation (Maraseni et al. 2014;) since 98% of the deforestation is occurring in developing countries (FAO 2020b). REDD þ is also one of the cheapest climate mitigation options for tackling climate change (Golub et al. 2018;Raihan et al. 2019;Raihan and Said 2022). However, such global initiatives require quantification of forest biomass and carbon (Campbell 2009;Koch 2010).
Pakistan's forest covers 4.6% of its total land area, according to the FAO (2020). It is a low forest cover country (LFCC) (Aftab and Hickey 2010). Forest per capita in Pakistan is only 0.03 ha compared to the global average of 1 ha (Ali et al. 2020a). Despite low forestry resources, Pakistan faced one of the highest deforestation rates in the world losing 0.94% annually between 2010 to 2020 (FAO 2020b). One of the main pressures on forestry resources in Pakistan has been the high reliance of local communities on these resources often for their subsistence (Shahbaz et al. 2007;Hussain, Zhou, et al. 2019). Such types of pressures are further aggravated by the increasingly adverse effects of climate change . Like many developing countries, Pakistan has contributed only a meager amount of 0.8% to the world's GHG emissions (Hussain, Butt, et al. 2019;Sajid 2020).
Forest biomass quantification becomes essential for the areas, such as the study area for this research, Margalla Hills National Park (MHNP), which are prone to increasing local pressures. The case of MHNP even becomes more unique, since it is adjoined with the capital city of Pakistan. The location is, thus, exposing the MHNP to encroachment due to the increasing population of the city. This, along with other internal pressures i.e. local communities' increasing dependencies on forest resources, forest fires due to humans, limestone extraction, etc. within the MHNP, are threatening the unique ecological resources of the MHNP. Forest biomass assessments can therefore help in the identification of the carbon potential of the MHNP which can be used for accessing financial resources such as REDDþ. The adequate financial resources, in turn, can help to develop and implement suitable sustainable management strategies, that can involve local communities of the MHNP in participatory management, since they depend on forestry resources. Given the aforementioned, the study aimed at the quantification of forest biomass at the MHNP. Specific objectives included i) Assessment of the AGB and Aboveground Carbon (AGC) of the trees of two forest types of the MHNP: Sub-tropical Chir Pine Forest (SCPF) and Sub-tropical Broadleaved Evergreen Forest (SBEF) ii) Determination of the structural characteristics (i.e. volume, basal area, DBH, and height) of the two forests iii) Identification of the contribution of individual tree species to each forest type iv) Comparison of the AGC between the two forests v) Regression analysis of the main tree species of the two forest types involving following variables: AGB, volume, basal area, DBH, and height. Since there is a lack of allometric equations in Pakistan for local tree species for AGB estimation, the BEF approach was used for the AGB assessment (Ismail et al. 2018).

Study area
The MHNP has an area of 17,386 ha ( Figure 1) , which is spread between 73 7'3.32" E, 33 41'59.61" N, and 73 1'34. It has 27 villages located within its boundaries. The climate is semi-arid. The area lies in the monsoon belt, and it experiences two rainy seasons: Winter rains from January to March, and summer rains from July to September. During the winter season, the temperature ranges from 1-15 C, and during the summer season, the temperature ranges from 20-40 C. It has five seasons: winter (December-February), spring (March-April), summer (May-June), monsoon (July-September), and autumn (November) (Anwar and Chapman 2000). The average annual rainfall is 1,000 mm (Himalayan Wildlife Foundation 2007). The soil is mainly composed of limestone (Butt et al. 2015). The flora of the MHNP includes 101 families, 548 genera, and 608 species (Akbar 2012).
The MHNP is composed of two types of forests, SCPF and SBEF ). The SCPF is predominantly located above 900 m (Himalayan Wildlife Foundation 2007) and is mainly characterized by the Pinus roxbughii Sarg. that forms the main canopy. Other areas in Pakistan where this forest is found include Murre Hills, Swat, Shinkiari, etc. (Nizami 2010). The SBEF is classified as an arid forest which is usually composed of branchy trees and varies from a complete closure, owing to favorable conditions, to scattered trees/shrubs on the drier sites (Sheikh 1993). These are low forests with species that are mostly thorny and often with evergreen leaves (Sheikh 1993). The larger trees can be seen in valleys where deep soil and adequate water are available. Generally, these forests are encountered by noticeable erosion and are characterized by gullies and deep ravines, which makes the rocks and boulders common features. In other areas in Pakistan, these forests are generally found in Kala Chitta Mountains (Kherimurat), salt ranges, Malakand Agency, Parachinar, Loralai, Sohawa, etc. (Nizami 2010).

Field inventory measurements
There were 46 sampling plots in the SCPF, and 31 in the SBEF, using a stratified random sampling approach. The sampling plots were circular with a radius of 20 m and an area of 0.126 ha (Pearson et al. 2005). The DBHs of all the trees above 5 cm, in the sampling plots, were measured with the help of the DBH tape above the ground level at 1.3 m (Laar and Akça 2007). The Vertex IV (Hagl€ of Sweden, AB), along with its transponder, was used for the estimation of tree heights in the sampling plots. The Vertex IV has a mean error of À0.66 m. Overall, there were 640 trees measured from the SCPF and 443 trees measured from the SBEF. The coordinates were also recorded from the center of the plots using handheld Global Positioning System (GPS) device (Garmin GPS 60 Marine). The GPS device can have an error of 3-4 m.

Structural and AGB analysis
The quadratic mean diameter (cm) for each plot was computed as follows (Eq. 1) (Iles and Wilson 1977); Where n is the total number of trees. The quadratic mean diameter in this research was referred to as the mean DBH. The tree volume (m 3 ) was recorded using the following formula (Eq. 2) (Philip 1994); Where H (m) is the height of a tree and FF (unitless) is the form factor. The form factor of a tree is defined as "the stem volume which is expressed as the proportion of the cylinder volume of the same height which has a diameter equal to the stem diameter at the selected reference point" (Laar and Akça 2007). For this research, a form factor of 0.45 was used for Pinus roxburghii, and for other broadleaved species, a form factor of 0.59 was used as suggested by (Gray 1956).
The stand basal area is defined as the sum of the cross-sectional area of the stems at breast height and expressed per unit of ground area (West 2009). The basal area (m 2 ) was calculated by using the following formula (Eq. 3) (K€ ohl et al. 2006); The density is the number of trees per unit area (West 2009). Initially, the AGB (Kg) was calculated as follows (Eq. 4) (Brown et al. 1999;Pandey et al. 2019); (4) The 'Wood Density (Kg m 3 )' is defined as "the amount of actual wood substance which is present in a unit volume of wood" (Zobel and Jett 1995). Wood densities for sampled tree species were sourced from literature (see Supplementary Table 1). The BEF of 1.51 was used for Pinus roxburghii and 1.59 was used for all other broadleaved tree species which were encountered during the sampling process (Haripriya 2000). The 'AGB (Kg)' was converted to 'AGB Mg (megagram)' by dividing the former by 1000. The 'AGB' was converted to 'Aboveground Carbon (AGC)' by using the carbon conversion factor of 0.5 (K€ ohl, Ehrhart, et al. 2020). The above-mentioned variables: Volume (m 3 ), Basal Area (m 2 ), Density, AGB (Mg), and AGC (Mg C) were scaled up to ha by dividing the variables by 0.1 (Torres et al. 2020).

Statistical analysis
Means and standard deviations were calculated for the variables (mean ± SD). The data was checked for the normality distribution, for which Shapiro-Wilk test was used. Since the data were not normally distributed, the Wilcoxon rank sum test was used for probing the statistical difference between the variables. The significance levels of 0.01, 0.05, and 0.1 were used for investigating the differences between the variables of the two forest types. The Spearman correlation coefficient, which is a non-parametric measure, was used for the correlation analysis of the variables since the data was not normally distributed (Best and Roberts 1975). For the Spearman correlation coefficient, which ranges between À1 to 1, the interpretation was done following Ratner (2009), for assisting in understanding the relationships based on different levels; 1) 0 indicating no relationship 2) þ1 indicating perfect positive linear relationship 3) À1 indicating perfect negative relationship 4) values between (0 and 0.3) indicating weak positive linear relationship and values between (-0.3 and 0) weak negative relationship 5) values between (0.3 and 0.7) indicating moderate positive linear relationship and values between (-0.3 and À0.7) indicates a moderate negative linear relationship 6) values between (0.7 and 1) indicating a strong positive linear relationship and values between (-0.7 and À1) indicating strong negative linear relationship. The Spearman correlation coefficient, Shapiro-Wilk test and Wilcoxon rank sum test were determined using the ggpubr package (Kassambara 2023) in the RStudio (R Core Team 2023). Regression analysis was also performed in RStudio, using the stats package in R (R Core Team 2023). Regression analysis was done for analyzing the relationship between every two variables of interest for the three dominant species (highest number of measured species) for the two forest types i.e. i) Pinus roxburghii Sarg. for SCPF ii) Bauhinia variegata L. and Grewia optiva J. R. Drumm. ex Burret for SBEF. Linear regression analysis was performed for linear relationship and a power curve was used in case of a non-linear relationship between the variables i.e. DBH and height; volume and basal area; DBH and AGB. Significance levels of 0.01, 0.05, and 0.1 were used to explore the statistical significance, and the Coefficient of Determination (R 2 ) was then used for assessing the strength of the statistical relationship between the two variables. The evaluation of the strengths of the R 2 was based on the following; 1) 0.19 (weak) 2) 0.33 (moderate) and 3) 0.67 (substantial/significant) (Chin 1998). The negative sign implies the negative statistical relationship between the variables.

Structural characteristics of the two forests of the MHNP
The mean AGC for SCPF was higher than the SBEF. The mean AGC for SCPF was 73.36 ± 32.55 Mg C ha À1 and the mean AGC for SBEF was 16.88 ± 25.81 Mg C ha À1 (Table 1). The mean AGB of the SCPF was also higher than the mean AGB of the SBEF. The mean AGB for SCPF was 146.73 ± 65.11 Mg ha À1 and the mean AGB for SBEF was 33.77 ± 51.63 Mg ha À1 . High significant differences were recorded between the AGB and AGC of the two forests (W ¼ 43, p < 0.01). SCPF was also greater than SBEF for all the other structural variables, except for the mean density ha À1 (Table 1).

Correlation analysis of structural variables of two forests with the AGB
Spearman correlation coefficients were calculated between the AGB and other structural characteristics of the SCPF (Table 2). Strong positive and significant correlations were recorded between AGB and volume, basal area, mean DBH, and mean tree height.
Spearman correlation coefficients were also calculated between AGB and other structural characteristics of the SBEF (Table 3). Strong positive and significant correlations were recorded between AGB and volume, basal area, mean DBH, and mean tree height.

DBH distribution of the two forests of MHNP
The DBH interval with the highest number of stems was recorded as 35-40 cm for the SCPF (Figure 2). It was followed by the DBH interval of 25-30 cm and the DBH interval of 20-25 cm. The total number of 114 stems was recorded in the DBH interval class 35-40 cm, which formed 17.81% of the total stems of the SCPF. Overall, there were 640 stems recorded for the SCPF. The DBH interval 25-30 cm recorded 107 stems, which formed 16.72% of the total stems, for the DBH interval 20-25 cm, 85 stems were recorded, forming 13.28% of the total stems of the SCPF. The maximum DBH was recorded as 67 cm, for Pinus roxburghii, and the minimum DBH as 10 cm, also for Pinus roxburghii, for the SCPF.
For SBEF, the DBH interval of 10-15 cm contained the highest number of stems ( Figure 2). It included 167 stems, which formed 37.70% of the total stems of the SBEF. There were overall 443 stems recorded for the SBEF. The DBH interval 15-20 cm recorded the second highest number of trees with 102 trees, constituting 23.02% of the total number of trees. The DBH interval of 5-10 cm formed the interval with the third highest number of trees of 80, constituting 18.06% of the total stems. The highest DBH was recorded as 87 cm for Ficus palmata subsp. virgata and the lowest DBH was recorded as 5.5 cm for Grewia optiva for the SBEF.
AGB (Mg ha 21 ) against mean DBH (cm) and mean tree height (m) The AGB (Mg ha À1 ) of the SCPF increases with the increasing mean DBH (cm) and mean tree height as shown in the 3D graph ( Figure 3).   The AGB (Mg ha À1 ) of the SBEF has shown an erratic pattern against the mean DBH (cm) and the mean tree height (m) as can be seen in the 3D graph (Figure 4).

AGC of the different trees species in the two forests of the MHNP
Pinus roxburghii contributed the highest percentage of 98.85%, which constitutes 3336.13 Mg C ha À1 , of the overall total AGC of 3374.80 Mg C ha À1 of the SCPF (Table 4). The Bombax cieba formed 0.81%, with 27.01 Mg C ha À1 , and the Bauhinia variegata contributed 0.34%, with 11.64 Mg C ha À1 , to the total AGC of the SCPF of MHNP.
For the SBEF, Eucalyptus camaldulensis formed the highest contribution of 33.65%, with 176.14 Mg C ha À1 , in the total AGC (Table 5). Ficus palmata subsp. virgata and Bauhinia variegata formed the second and third highest contributions with 94.24 Mg C ha À1 , 18%, and 73.35 Mg C ha À1 , 14.01%, respectively, in the total AGC of the SBEF of the MHNP. Similarly, Acacia modesta, Bombax ceiba, and Grewia optiva followed with 42.29 Mg C ha À1 (8.08%), 36.94 Mg C ha À1 (7.06%), and 30.84 Mg C ha À1 (5.89%). The lowest contributions were recorded for Ficus benghalensis with 0.7%, constituting 3.67 Mg C ha À1 , and Flacourtia indica, with 0.1%, forming 0.51 Mg C ha À1 , of the total AGC of the SBEF of the MHNP.

Regression analysis of the main tree species of the two forests
All the equations were highly significant ( Table 6). The R 2 values for all the equations ranged between 0.36 to 0.98. The Bauhinia variegata recorded a moderate statistical association between the DBH and height (R 2 ¼ 0.36). The highest R 2 value (0.98) was recorded for the Pinus roxburghii, for its power curve regression, between the DBH and AGB.  Amir et al. (2018), 264.5 Mg C ha À1 , for the SCPF at Murree Hills, Pakistan, is also higher than this study. The study by Shaheen et al. (2016) for Muzaffarabad District, AJK, also recorded higher AGC for SCPF compared to this study. In three of these studies by Amir et al. (2018), Nizami (2012), and Mannan et al. (2019), the belowground carbon (BGC) was also included in the mean AGC. The BGC generally constitutes 20% of the AGC (Cairns et al. 1997). The result of this study is also in agreement with the estimates provided by Singh et al. (1985) who reported AGC of 35 to 75.2 Mg C ha À1 for the degraded forests of the Central Himalayan, India. The lower AGC for the SCPF can be attributed to the lower AGB yield of the study area (Segura and Kanninen 2005). This low yield may also be due to anthropogenic activities and the natural disturbances occurring in the area (Bhadwal and Singh 2002;Rossi et al. 2007). The depletion in the AGC of the forest ecosystems is mainly attributed to forest degradation in developing countries (Pant and Tewari 2014).

AGB and AGC of the SCPF
The AGB estimation enables the quantification of the AGC, which the forest stores (Naveenkumar et al. 2017). AGB estimation, therefore, is essential for the quantification of carbon budgets and for the national development planning for overall forest management, which makes it also more important given the challenges posed by climate change (Devagiri et al. 2013). The mean AGB reported for this study, 146.73 ± 65.11 Mg ha À1 , is lower  Himalayan, India, which is again higher than the estimate of this study. These disagreements may mainly be due to the differences in sites and the sampling approaches.

Structural characteristics of the SCPF
The forest structure defines the arrangement and relationship between different tree components: trunks, branches, and leaves (Ackermann 2015). It consists of various attributes including type, size, shape, and spatial distribution of forest (Pommerening 2002;McElhinny et al. 2005). The mean basal area reported by Nizami (2012) for the SCPF in Ghoragli and Lethrar, which were 30.38 m 2 ha À1 and 26.11 m 2 ha À1 respectively, are higher than the estimate of this study, 14.48 ± 4.80 m 2 ha À1 . The basal area is among the common variables such as DBH, total height, volume, and crown diameter, which are used for the estimation of the AGB . The result of this study is also not in agreement with Amir et al. (2018) who provided a mean basal area of 32.3 m 2 ha À1 for the SCPF of Murree Hills, Pakistan. Sharma et al. (2020) provided an estimate of 32.1 m 2 ha À1 for a plantation site of Pinus roxburghii in Kathmandu, Nepal, which is also not in agreement with this study. The estimate for this study however is higher than the minimum mean basal area provided by Kumar et al. (2019) which was 11.12 m 2 ha À1 for the SCPF of Garhwal Himalaya, India. The mean tree volume of 198.75 ± 87.76 m 3 ha À1 reported by this study is on the lower side compared to Nizami (2012) who reported 243 m 3 ha À1 for the SCPF of Ghoragli, but slightly higher than the 197 m 3 ha À1 for Lethrar area. The result for the mean tree volume for this study is more than the lower end of the mean tree volume of 83 m 3 ha À1 , which was for the young stands of SCPF, reported by the Amir et al. (2018) for Murree Hills, Pakistan. The higher-end estimate of the study, 549.4 m 3 ha À1 , by Amir et al. (2018) for SCPF for Murree Hills, Pakistan, is higher than the estimate provided by this study. The mean tree volume estimate provided by Kumar et al. (2019) for SCPF of Garhwal Himalaya, India, which was 214.88 m 3 ha À1 , is also higher than the result of this study. The mean tree density reported by Nizami (2012) for SCPF was 878 trees ha À1 for Ghoragali and 776 trees ha À1 for Lethrar areas. These estimates are on the higher side compared to the mean tree density of 139.13 ± 22.88 trees ha À1 for this study. The mean tree density estimate for this study is also lower than the estimate provided by Kumar et al. (2019) for SCPF, which was 244.34 trees ha À1 for Garhwal Himalaya, India. The tree density is also indicative of deforestation pressures on the forest ecosystems since people are dependent on forests for their livelihood and this can further be aggregated by additional demands for construction and industrial demands (Wani et al. 2010;Shaheen et al. 2016). In addition, lower tree densities may also portray lower forest carbon stock values (Kindermann et al. 2006). The mean tree density of this study is higher than the estimate provided by  who reported 103-104 trees ha À1 for the SCPF of Murree Forests, Pakistan, and thereby have attributed their lower densities to the anthropogenic pressures such as lopping and high tourism. The DBH distribution for the SCPF has shown a hump-shaped (unimodal type) (Subedi et al. 2018), which shows that there is a lesser number of stems in the lower DBH classes, and the numbers of stem increase with the increasing DBH classes until the classes 25-30 cm and 35-40 cm and after which the number of stems starts to decrease again (Figure 2). These two classes contained 16.72% and 17.81% of the total stems and were also the two highest classes in the DBH distribution. Pariyar et al. (2019) also reported similar results for the mixed Pinus roxburghii forest, of the Parbat District of Western Nepal, where most of the stems had a DBH of 30 cm compared to the pure stands of Pinus roxburghii which had higher representations from higher DBH classes. This type of distribution has also been observed elsewhere in different areas of the Himalayan region (Kunwar and Sharma 1970;Vetaas 2000;Sharma et al. 2008). This unimodal distribution is not a good sign for the forest ecosystem, since it is attributed to anthropogenic and natural disturbances (Wangda and Ohsawa 2006;Bin et al. 2012). Furthermore, looking at the classification of forest structure for Pinus roxburghii based on DBH; < 10 cm (sapling), 10-30 cm (young), 30-50 cm (middle-aged), 50-70 cm (mature), > 70 cm (overmature) (FRI 2018), it was also found that highest percentage of the stems belonged to the young and middle-aged DBH classes with 39.84% and 52.18%, respectively. This distribution represents imbalance and poses further challenges for the sustainable management of the forest (FRI 2018).
Regression analysis between DBH and the height of the Pinus roxburghii showed a significant association (Table 6). Tree heights are essential parameters for forest inventory variables and carbon budget models (Vanclay 1994). Tree heights are difficult to collect in the forests and are also costlier compared to DBH data collection (Peng et al. 2001).
Mostly, height-DBH equations are developed for predicting the missing heights, which were skipped during the field expeditions (Peng et al. 2001). In this study, all the heights for the tree stem for both forests were measured during the field expedition, 640 for the SCPF and 443 for SBEF. The highest height recorded for Pinus roxburghii was 41.7 m and the lowest was 13.9 m. Pinus roxburghii is a fast-growing pioneer species, which can reach a height of 51 m and DBH of 100 cm (Jackson 1994 Nizami et al. (2009), andSubedi et al. (2018). Tree DBH and height can also directly influence biomass production, where lower DBH and heights can result in lower biomass and carbon stock values for the forests (Feldpausch et al. 2012).
The regression analysis between the basal area and volume of the Pinus roxburghii has shown a significant linear association (Table 6). The result of the linear regression analysis is not in agreement with Amir et al. (2018) who have shown a polynomial relationship between the two variables. The result, however, is in agreement with Nizami et al. (2009). The regression analysis between DBH and AGB showed a significant non-linear relationship (Table 6). This is in agreement with Ali et al. (2020a) and Solomon et al. (2017), who also reported a significant non-linear relationship between the two variables. Many studies have shown that DBH is a reliable predictor of AGB (Giday et al. 2013;Hasen-Yusuf et al. 2013;Solomon et al. 2017).

AGB and AGC of the SBEF
The carbon stock assessment in the sub-tropical forests is also of utmost importance because they have been poorly studied in comparison to the other forest types (Rosenfield and Souza 2013). Table 8 shows AGC estimates from different studies, including this study. The result of this study is on the very low side compared to the study by Mannan et al. (2019) who reported 73.81 Mg C ha À1 for the AGB of the SBEF which also included the BGC. The result for this study is also on the lower side compared to the mean AGC reported by Ali et al. (2019), for SBEF, which was 49.54 Mg C ha À1 for the sites in Swat, Dir Upper, Dir Lower, and Buner Districts in the KPK, Pakistan. The result of this study is also on the lower side compared to Ali et al. (2018) who reported a mean AGC of 39.02 Mg C ha À1 for SBEF at a site in Khanpur Range, KPK, Pakistan. The estimate for SBEF from this study is also lower than the estimates by Nizami (2012) who had reported mean AGC (including BGC) of 25.54 Mg C ha À1 and 20.23 Mg C ha À1 for SBEF at the sites of Kherimurat and Sohawa respectively. The result of this study is higher than the estimate reported by Ali et al. (2020b) which was 4.52 Mg C ha À1 for the SBEF for the entire KPK, Pakistan.
The mean AGB reported by Nizami (2012) which was 50.93 Mg ha À1 and 40.43 Mg ha À1 (both estimates also included the BGB) for the SBEF for the sites of Kherimurat and Sohawa respectively are on the higher side comparing to this study, which is 33.77 ± 51.63 Mg ha À1 . The mean AGB reported by Ali et al. (2018) which was 83 Mg ha À1 for the SBEF for Khanpur Division, which is also on the higher side compared to the result of this study; they attributed their result to the high tree density of SBEF located at their study sites.

Structural characteristics of the SBEF
The mean tree volume recorded for the SBEF in this study, 37.82 ± 50.38 m 3 ha À1 , is on the higher side compared to the estimate provided by Nizami (2012) which was 12.86 m 3 ha À1 and 11.40 m 3 ha À1 for SBEF for the sites of Kherimurat and Sohawa, respectively. Mannan et al. (2019) reported the mean tree volume for SBEF as 115.71 m 3 ha À1 , which is on the higher side compared to the mean tree volume estimate of this study. Ali et al. (2018) recorded a mean tree volume of 57.66 m 3 ha À1 for SBEF for Khanpur Range, KPK, Pakistan, which is also on the higher side compared to this study. The mean tree basal area for this study was recorded as 4.28 ± 3.43 m 2 ha À1 , which is slightly higher than the estimates of 3.06 m 2 ha À1 and 2.65 m 2 ha À1 reported by Nizami (2012) for SBEF for the sites of Kherimut and Sohawa, respectively. The estimate for this study is less than the estimate of 13.1 m 2 ha À1 reported by Mannan et al. (2019) for SBEF for the sites at Margalla and Murree Forest Division. The estimate reported from this study is also close to the estimate of 6.69 m 2 ha À1 reported by Ali et al. (2019). The mean tree density of 142.90 trees ha À1 was recorded for this study, which is lower than the mean densities of 197 trees ha À1 and 179 trees ha À1 reported by Nizami (2012) for SBEF for Kherimurat and Sohawa, respectively. The mean tree density of this study is also on the lower side compared to Mannan et al. (2019) who reported a mean tree density of 761.18 trees ha À1 for SBEF. Ideally, the DBH distribution of a forest should follow an inverse J shape (Shimano 2000). The DBH distribution of the SBEF is positively skewed which means tree density increases till the DBH class is 10-50 cm and after which it starts to decrease. The DBH class 10-15 cm also contained the highest percentage of the stems, 37.70%. The highest number of stems in the lower DBH classes shows a good regeneration potential (Ballabha et al. 2013). It also highlights that a forest is in evolving stage (Campbell et al. 1992). But at the same time, the absence of trees in the higher DBH classes suggests historical disturbances which must also be not neglected while undertaking the overall planning for management purposes (Nath et al. 2017). The height distribution of the SBEF also showed a positively skewed distribution, with most of the stems being present in the two classes, 6-9 m, and 9-12 m.
The regression analysis between DBH and the height of the Bauhinia variegata had shown a moderate linear association (Table 6). The DBH range recorded for Bauhinia variegata was 5.6-35 cm. The height range was 5.4-18.1 m. The regression analysis between basal area and tree volume of Bauhinia variegata had shown a significant linear association ( Table  6). The mean basal area for Bauhinia variegata was recorded as 1.16 ± 0.77 m 2 ha À1 and the mean tree volume was recorded as 6.88 ± 4.77 m 3 ha À1 . A significant non-linear regression association was recorded between DBH and AGB of Bauhinia variegata ( Table 6). The result of this study was a non-linear significant regression association between DBH and AGB of Bauhinia variegata, which is in agreement with Kanwar (2019) who also recorded a strong non-linear relationship between DBH and AGB. The mean AGB of Bauhinia variegata, 7.33 ± 5.02 Mg ha À1 , from this study, is not in agreement with Kumar and Sharma (2015) who reported a mean AGB of 25.82 ± 48.57 Mg ha À1 for the Balganga Range of the Garhwal Himalaya region, India.
The regression analysis between DBH and the height of the Grewia optiva showed a moderate linear association ( Table 6). The DBH for Grewia optiva ranged from 5.5 cm to 22 cm and the height ranged from 5 m to 15 m. The regression analysis between the basal area and volume of the Grewia optiva had shown a significant linear association ( Table 6). The mean basal area and the mean tree volume were recorded as 0.83 ± 0.54 m 2 ha À1 and 4.75 ± 3.30 m 3 ha À1 , respectively. The regression analysis between DBH and AGB showed a moderate non-linear association ( Table 6). The result for the moderate non-linear association of this study between the DBH and AGB of Grewia optiva is not in agreement with Kanwar (2019) who had reported a strong non-linear association between the two variables at the Balganga Range of the Garhwal Himalaya region, India. The mean AGB for Grewia optiva, in this study, was recorded as 5.14 ± 3.57 Mg ha À1 . This mean AGB for Grewia optiva is on the lower side compared to Rana et al. (2020) who reported a mean AGB of 35.9 Mg ha À1 for Grewia optiva in Garhwal Himalaya, India. The result, however, is slightly higher than the mean AGB estimate of 4.91 Mg ha À1 for Grewia optiva reported by Kumar et al. (2012) for Dungripanth Village, Srinagar. Grewia optiva is a very important agroforestry species, which can also provide fodder for livestock during the winter season (Verma et al. 2014). It is also a very good fuelwood species since it has a calorific value of 4920 k cal kg À1 (Joshi and Dhiman 1992).

Comparison of AGC between SCPF and SBEF
The SCPF recorded higher AGC compared to the SBEF, which has also resulted in a highly significant statistical difference between the two. This can be attributed to the higher mean DBH and higher mean height recorded for the SCPF which were 36.37 ± 6.51 cm and 27.34 ± 3.67 respectively ( Table 1). The carbon storage in the forests is the function of size (girth), wood density, and height (Chave et al. 2005;Mac ıas et al. 2017;Moussa et al. 2019). The small-sized trees store less carbon compared to large-sized trees (Nero et al. 2018;Xu et al. 2018;Dimobe et al. 2019;Portela et al. 2020). Generally, therefore, the mature trees, which tend to have bigger and larger sizes such as DBHs, in a forest ecosystem store more carbon (Gandhi and Sundarapandian 2017;K€ ohl et al. 2017). In addition, larger trees have been reported to store more carbon despite lesser density and richness (Lutz et al. 2018;Mensah et al. 2020). The SCPF recorded 7.81% of the trees in the mature category; DBH 50-70 cm (Figure 2). At the same time, SCPF recorded 39.84% of the stems in the small-sized young-aged DBH category < 30 cm. This category, therefore, can also be not neglected and hence can play a critical and important role in future carbon sequestration (Gandhi and Sundarapandian 2017). Comparing all the plots, the SCPF had higher AGC in all the plots compared to SBEF. The difference between the AGC of the two forests can therefore be attributed to the higher DBHs and heights recorded for the trees in the SCPF, as already mentioned above. The difference could also be attributed to the species composition since it is also considered an important factor in modulating forest carbon stocks (Reyes-Palomeque et al. 2021). Negi et al. (2003), for instance, reported that the maximum AGC was stored in the following order: pines > deciduous > evergreen > bamboo. A similar, pattern was also reported by Kaushal and Baishya (2021). The result of this study is in agreement with Nizami (2012), Mannan et al. (2019), and Sharma et al. (2010) who also reported higher AGC for the SCPF over the SBEF. Gogoi et al. (2020), however, reported higher AGC in the SBEF in comparison to the Pine Forests which were 140.4 Mg C ha À1 and 133.6 Mg C ha À1 for SBEF against the 74.7 Mg C ha À1 and 63 Mg C ha À1 for Pine Forests of the Meghalaya District, India. They also attributed the higher carbon stocks to the presence of overmature trees, with DBH > 66 cm, and the relatively higher density of younger trees in the SBEF.

Species-wise AGB contribution
The Pinus roxburghii contained the highest AGC, with 98.85% of the total, for the SCPF (Table 4). The result is in agreement with Pariyar et al. (2019) and Shaheen et al. (2016) showing higher contributions of AGC by the Pinus roxburghii. Pinus roxburghii is usually a dominant species and mostly grows in pure stands throughout Pakistan, India, Nepal, and Bhutan (Gupta and Dass 2007;Khan et al. 2014). The species is commonly known as 'Chir Pine' (Anwar et al. 2019). It is spread over Bhutan, China, India, Nepal, and Pakistan from 400 to 2300 m a.s.l. (Polunin and Stainton 1984). It is native to the foothills of the Himalayan Hindu Kush Region (Aryal et al. 2018). The species provides various economical and ecological benefits, due to which nearby communities leaving in or around its zone rely on it for various purposes (Gupta and Dass 2007;Khan et al. 2014).
The Eucalyptus sp. is native to Australia (Afzal et al. 2018) and was introduced in Pakistan in the late 19s (Shah et al. 1991). Amongst the 40 sp. introduced in Pakistan, Eucalyptus camaldulensis has been the most widespread and adaptable (Bilal et al. 2014). It is a fast-growing species (Afzal et al. 2018) and is locally known as 'Sufeda' (Khan et al. 2017). The DBH and height for Eucalyptus camaldulensis in this study ranged between 38-86 cm and 22-35.5 m, respectively. The Eucalyptus camaldulensis recorded a total AGC of 176.14 Mg C ha À1 in this study, contributing 33.65% of the total AGC for the SBEF (Table 5). The result of this study is in agreement with Arif et al. (2017) who reported the highest AGC contribution from the Eucalyptus camaldulensis toward the Chichawatni Plantation in Pakistan. Mandal et al. (2017) also recorded the highest AGC contribution by the Eucalyptus camaldulensis in a mixed plantation in the Mahottari District of Nepal.
The Ficus palmata subsp. virgata was the species with the second highest AGC contribution to the SBEF. It contributed AGC of 94.24 Mg C ha À1 which was 18% of the total AGC. The DBH ranged from 8.2 cm to 87 cm and the height of the Ficus palmata subsp. virgata for this study ranged from 6.2 m to 18.6 m. The geographical distribution of the Ficus palmata subsp. virgata includes; Afghanistan, Arabian Peninsula, Ethiopia, Egypt, Iran, Nepal, Pakistan, Somalia, Sudan, and India (Chaudhary et al. 2012). The species is well known for its medicinal properties (Saqib et al. 2014;Salehi et al. 2021). The result of this study is not in agreement with Shaheen et al. (2016) who reported other tree species with higher AGC concentrations than the Ficus palmata subsp. virgata in AJK. The study is also not in agreement with Ali et al. (2019) who reported lower AGC contribution from Ficus palmata subsp. virgata., compared to other species, in four districts of Buner, Dir Upper, Dir Lower, and Swat of Khyber Pakhtunkhwa, Pakistan.

Conclusion
In this study, it was found that SCPF had a higher mean AGB and mean AGC compared to the SBEF. There were highly significant statistical differences recorded between the mean AGB and mean AGC of the two forests. The higher mean AGB and mean AGC of the SCPF in comparison to the SBEF can be attributed to the larger sizes of the trees i.e. DBHs and heights, and species compositions. Pinus roxburghii recorded higher AGC compared to other tree species in the SCPF. Eucalyptus camaldulensis recorded the highest AGC, amongst all the tree species, for the. The equations presented (Table 6), based on their statistical significance and relationship strength, will allow for future estimation of the investigated variables. The unimodal DBH distribution of the two forests suggests they are experiencing deforestation and degradation. Further comparison with other similar forest types also hinted toward the lower AGB and AGC in the study site. The lower AGB and AGC in a forest ecosystem are mainly associated with human or natural disturbances (Lugo and Brown 1992). Adequate forest management, therefore, is necessary for conserving and enhancing not only the AGB and AGC but also for safeguarding the associated biodiversity and the local communities which depend upon the forestry resources for their sustenance (Feldpausch et al. 2005). For this purpose, sustainable forest management plans should be developed for MHNP which should also include local communities for creating a sense of ownership. Adequate mechanisms should be developed for the implementation of environmental laws for the conservation of forestry resources. Afforestation programmes using Assisted Natural Regeneration (ANR) technique should be initiated to address the imbalance in the structures of the two forests. These afforestation programmes could also contribute to the Paris Climate Agreement, the Trillion Trees initiative, and the Bonn Challenge -which aims to restore 350 million ha of degraded and deforested areas by 2030. The site and species-specific allometric equations should be developed for future AGB and AGC studies.