Quantifying Land Use Changes and Carbon Stock Degradation in Ayubia National Park Using Application of Markov Chain Model

1. Introduction

The alteration of the earth's surface caused by human activities that disrupt the natural ecosystem and environmental processes of the planet is known as Land Use/Land Cover Change(LULCC) (Hassan.Z,2012). Land use change analysis necessitates a deeper comprehension of the factors that lead to LULCC, which typically results intricately integrated sequence(Winkler 2021). The accurate information of LULCC is of a great importance for the forests management, urban planning, agricultural land development, environmental management, tourism, industries and transportation etc. LULCC affect forests at both local and global level. The loss of forest cover as a result of LULCC has led to environmental problems, global warming, habitat destruction for wildlife, and disturbance. It has also had a negative impact on the carbon content of forest biomass(Tewabe 2020). Remote sensing images are valuable source of spatial and temporal information to track the dynamics of LULC changes and forecast their effects on any environmental problems on a regional scale. Changes in land use and cover have become key for many different applications, including those in agriculture, the environment, ecology, forestry, geology, and hydrology, according to empirical investigations by researchers from a variety of disciplines(Hao, Zhu et al. 2021). At the same time, a significant project to study land use change has evolved as a worldwide initiative in recent decades and has acquired enormous momentum in its efforts to understand the processes driving land use change(Lambin,1997).These initiatives have simulated researcher interest in using a variety of methodologies to identify and further model environmental dynamics at various levels.

Land use planning involves strategies and processes aimed at managing and regulating the transformation of the Earth's surface due to human activities was given significance in 1970 because of its assorted relevance after the approach of multispectral satellites images(Gondwe 2021). An efficient and precise method for exploring terrestrial ecosystems on a large scale has been developed in recent decades as a result of advancements in remote sensing (RS) techniques(Abebe, 2021). Remote sensing(RS) provides continuous temporal and spatial analyses as RS data obtained at various resolutions from different sensors can provide ecological analysis at different scales. Satellite images are periodic and cost effective allowing the researchers to conduct precise LULCC monitoring(Kumar and Arya 2021). Spatial temporal analysis of remote sensing data is also used to see the temporal patterns and identify the key drivers of deforestation(Reis 2008).

The RS data combined with the field data for validation to get precise results for Land use change(LUC) and land detection and the researchers are able to carry out precise LULCC monitoring because landsat satellite images are available on a regular basis and are inexpensive but not all satellite images are inexpensive. The cost of satellite imagery can vary significantly depending on resolution, coverage and frequency (Hassan, Shabbir et al. 2016). The key drivers of deforestation can also be identified through spatial temporal analysis of remote sensing data(Qureshi 2012).Urbanization, population growth, socioeconomic development, and forest fires are all contributing factors to the recent decline in forest biomass carbon (Sun and Liu 2019). According to(Nandal, Yadav et al. 2023), forest fires pose the greatest threat to stored terrestrial carbon, and each year, an estimated 2-3 Pg C, or 3-4 million km2 of burned forests, are released into the atmosphere. Forest fires cause changes in the world's natural and actual qualities which throughout the long term influence carbon stock(Song, Yang et al. 2021).

According to several studies, LULC Change caused an estimated 14% of the total forest area and 5531.5 Mg C yr−1 to be lost across Pakistan's subtropical forests between 1992 and 2022 (FAO, 2015; Ahmad et al., 2018; Mannan et al., 2018). In a similar vein, Ahmad et al. (2018) discovered that an estimated 629 ha (ha) of forests were lost in the Pakistani foothills of the Himalayan Mountains between 1994 and 2016. An estimated 2.1% annual carbon loss in forest biomass was noted in a related study (FAO, 2009) covering the years 1990–2010. According to (UN Habitat, 2010), LULCC, mostly caused by deforestation, have raised the average temperature in the region by almost 3.5°C. The number of rainy days and duration of the rainy season in the region have decreased as a result of these changes.The amount of forest cover on Earth was 2.5 million hectares in 1990; this decreased to 2.1 million hectares in 2000 and 1.7 million hectares in 2010, with annual decreasing rates of 1.6% and 2.0% in the first and second decades, respectively (FAO, 2010a). Over 40% of the potential biomass on Earth is found in forests, and between 1990 and 2007, forests sequestered between 1.9 and 2.6 Pg C of this biomass annually, according to Saeed et al. (2016). Similarly, it was discovered that forest carbon accounted for the majority of terrestrial carbon cycles, with soils containing over 33% and plants accounting for 77% (Saeed et al., 2016).Within the past few decades, LULCC have become a major challenge to the foothill forests of the Pakistan's Himalayan Mountains (Ahmad et al., 2018). In the meantime, between 1990 and 2017, the subtropical forests of Pakistan are estimated to have lost 5.70 Gg C yr−1 of forest carbon. During the same period, local temperatures were trending upward, there were fewer rainy days and more rainfall, and landslides were occurring more frequently (Rukya et al., 2014). According to a related study (Lu et al., 2018), between 1700 and 2005, the United States lost an estimated 3.5 million km2 of forests, which translated into a loss of roughly 10607 TgC.

In our study area the temporal effects of LULCC on forest Biomass Carbon and future prediction has not been studied before. This work attempts to investigate LULCC and quantify its effects on forest biomass carbon in the western Himalayan region of Pakistan using time series data together with Markov Chain Model. This would be extremely helpful in the formulation of policies that might assist global and local efforts for sustainable forest management.It advances our knowledge of the factors influencing the dynamics of carbon in forests and supports international initiatives to reduce global warming by improving land management techniques. Policymakers, land managers, and environmentalists can benefit from the findings by learning how to maintain and improve forest carbon stock in the change of shifting land use patterns. The study underscores the importance of protecting forest biomass carbon as a vital part of the global carbon cycle and climate management. Through measuring the effects of LULCC on forest carbon stocks, the research offers crucial data for programs aimed at mitigating climate change and planning sustainable land use. In order to preserve the ecosystem services that forests provide, such as carbon sequestration, biodiversity preservation, and watershed protection, it highlights the necessity of integrated approaches that strike a balance between land use changes and environmental conservation.

2. Materials and Methods

2.1. Study Area

The study area Ayubia National Park is Pakistan's sole temperate coniferous forest, and it supports a wide variety of species of fragile plants and animals. In the park area, there are roughly 10 different types of Gymnosperm trees, 200 species of herbs and shrubs, and other plants. It occupies an area of 34.09 sq.km hectares and is located in Abbottabad's Galiyat Forest division between 34° 1' 45.1164" N latitude and 73° 24' 9.1836"E longitude as shown in Figure 1-1. On April 17, 1984, the region was designated as a national park and situated at the northwest corner of Murree on a series of hills that stretch north to south 1220 to 2865 meters above sea level. The highest peaks are Mukshpuri (2865m) and Mirangani 2228m (Waseem,2005). In addition to receiving precipitation in the form of substantial winter snowfall, the average annual rainfall is above 1,500 mm, and the average annual temperature is 21 C with 66% a relative humidity. Kundala, Toheedabad, Mallach, Lahurkas, Kalabun, Derwaza, Mominabad, Ram kot, Raila, and Pasala are significant settlements around the park. The substantial depletion of the vegetation in Ayubia National Park is a result of the intense population pressure from the neighboring communities. The minimum elevation of the study area is 1242 and maximum lies at 2518 m.

2.1. Data Acquisition

United States Geological Survey (USGS) were used to download Landsat 5 Thematic Mapper (TM) images for the time periods of 1992, 2002, and 2012 as well as Landsat 8 Operational Land Imager (OLI) for 2022. Landsat 8 (OLI) includes nine bands, of which bands 1 through 7 and band 9 have a spatial resolution of 30 m, while band 8 is a panchromatic band with a spatial resolution of 15 m. The Landsat 5 TM has seven spectral bands with a spatial resolution of 30 m. Table 1-1 lists the specifics of the data sources and their specifications.

The images with the least amount of clouds (0–10%) were downloaded. Table 1 provides information about the sensor specifications.

2.1. Land-Use Classification

With the use of a GPS Magellan Explorist 110, which represented a homogeneous class, a number of training signature polygons were first produced for each land use to be classified. During the field survey, a GPS device was utilized to precisely capture the geographic coordinates (latitude and longitude) of certain places that reflect various classifications of land cover. Using the polygon feature in ArcGIS Pro, training samples were gathered for this investigation. To set each class apart from the others, they were each given a distinct identity and given a certain color. By defining polygons around typical locations, training samples were chosen for each of the preset land cover/use types. Pixels contained by these polygons were used to record the spectral signatures for the various land cover types acquired from the satellite imagery. The spectral signature that guarantees "minimum confusion" among the land covers that need to be mapped were considered satisfactory. The previous pattern of LULC of the studied region was determined using supervised classification algorithms with Landsat satellite data for the years 1992, 2002, 2012 and 2022. The area was categorized into five major land use/land cover classes as shown in Figure 2-1. For each LULCC class, the land cover area (ha) was computed. Results were produced for each year using the kappa coefficients and error matrix. The kappa coefficients, producer accuracy, total accuracy, and user accuracy for the chosen years were calculated.

Table 2.-1. Land use land change categories.

Serial no	LULC Classes	Descrption
1	Bare Land	A land with no vegetation or grasses
2	Conifer forest	Trees that grow needles instead of leaves and cones instead of flowers conifers tend to be evergreen they bear needles all year long.
3	Built-up/Settlements	Includes residential areas like town, roads,villages,strip transportation and commercial areas
4	Grassland	A grassland is an area where the vegetation is dominated by grasses with flat large area
5	Mixed Forest	Mixed forest is a vegetation type dominated by a mixture of broadleaf trees and shadow conifers.

Table 2.-1. Land use land change categories.

Serial no	LULC Classes	Descrption
1	Bare Land	A land with no vegetation or grasses
2	Conifer forest	Trees that grow needles instead of leaves and cones instead of flowers conifers tend to be evergreen they bear needles all year long.
3	Built-up/Settlements	Includes residential areas like town, roads,villages,strip transportation and commercial areas
4	Grassland	A grassland is an area where the vegetation is dominated by grasses with flat large area
5	Mixed Forest	Mixed forest is a vegetation type dominated by a mixture of broadleaf trees and shadow conifers.

2.1. Image Pre-Processing

Using ENVI 5.4, the satellite images for 1992, 2002, 2012, and 2022 years were radio-metrically corrected. After that we applied the FLAASH atmospheric correction. The top of the atmosphere (TOA) technique was used to convert the digital numbers (DN) into radiance and reflectance. For image correction, the equations spectral radiance(1-1) and TOA(1-2), were applied. Based on MODTRAN, FLAASH is an advanced method that can more precisely account for atmospheric phenomena. The average elevation of the study region, the scene center's coordinates, type of sensor, the date and time of the flight, and details on aerosol dispersion, visibility, and water vapor conditions are all input into the FLAASH module. The FLAASH atmospheric correction is the basis of our classification.

$${\mathbf{L}}_{\mathsf{\lambda }}=\mathbf{G}\mathbf{a}\mathbf{i}\mathbf{n}\times \mathbf{p}\mathbf{i}\mathbf{x}\mathbf{e}\mathbf{l} \mathbf{v}\mathbf{a}\mathbf{l}\mathbf{u}\mathbf{e}+\mathbf{o}\mathbf{f}\mathbf{f}\mathbf{s}\mathbf{e}\mathbf{t}$$

(1- 1)

$${\mathbf{p}}_{\mathsf{\lambda }}={\displaystyle \frac{\mathsf{\pi }\times {\mathbf{L}}_{\mathsf{\lambda }}\times {\mathbf{d}}^2}{{\mathbf{E}\mathbf{S}\mathbf{U}\mathbf{N}}_{\mathsf{\lambda }}\times \mathbf{s}\mathbf{i}\mathbf{n}\mathsf{\theta }}}$$

(1- 2)

Where ${\mathbf{E}\mathbf{S}\mathbf{U}\mathbf{N}}_{\mathsf{\lambda }}$ is the solar irradiance, d is the Earth-Sun distance, ${\mathbf{L}}_{\mathsf{\lambda }}$ is the radiance, ${\mathbf{p}}_{\mathsf{\lambda }}$ is the TOA reflectance, and is the Sun elevation in degrees. Using Fast Line-of-Sight Atmospheric analysis of hypercubes (FLAASH), which lessens the impacts of humidity, haze, water vapors, aerosols, and other variables, the top of the atmosphere reflectance was then transformed to surface reflectance. The spectral signatures vary depending on the atmospheric conditions on different days. Before detecting land-use change, the conversion of TOA reflectance to surface reflectance was completed using the correct atmospheric calibration approach, as explained by (Lambin,et al. 2003, Luo,et al. 2022).

2.1. Image Classification and Accuracy

In ArcGIS 10.5, the shape file for the research region was created and the area was clipped using the topographic sheet (1:25000) from the Survey of Pakistan (SOP) as a reference. Using the Random Forest classifier approach in ArcGIS Pro 10.4, supervised classification was used to categorize various land uses because a Random Forest has a higher training time than other classifers and we increase the number of trees in a random forest, the time to each train sample also increases. It also reduces the overfitting risk and providing high level of precision. We collected 80 training samples for each class which was enough for good classification. The different land cover classes of the area were represented using a stratified random approach for the accuracy evaluation of land cover maps that were separated from satellite images. Based on ground truth data and visual interpretation, the accuracy was surveyed using 250 points. An error confusion matrix was used to compare the classification results with the reference data. As it represents diagonal elements and all the components in the confusion matrix, a non-parametric Kappa coefficient was also used to quantify the degree of classification accuracy (Musetsho, Chitakira et al. 2021). Kappa is a statistic that assesses the degree to which user ratings and specified producer ratings concide. Equation 1-4 is the calculation for this:

K=P(a)-P(e)/1-p(e)

(1-4)

Where P (a) is the number of time the k raters agree, and P (e) is the number of time the k raters are expected to agree only by chance (Rijal, Rimal et al. 2021).

The application of a bootstrap method to measure overall decision tree categorization accuracy and confidence was reviewed, along with its description and application to land use sampling strategies. Consistent with previous findings, there was a substantial correlation between classification accuracy and confidence for overall accuracy and the number of fields assessed within a region. The association between sample size, accuracy, and confidence was not as significant for the accuracy of particular land use groups. While some classes displayed bigger confidence intervals and lower accuracies even when higher numbers of samples were used, others demonstrated high estimated accuracy and confidence even when relatively limited samples were utilized.

2.6. Markov Chain Model

A Markov Chain is a stochastic process that evolves over time in discrete steps, where the probability of transitioning from one state to another depends only on the current state(Weng, 2002). According to the Markovian system, a land use system's future state is predicted based on the current situation is the transformation of a system from one state to another, and the chance of this state shift is known as the Transition probability(Araya, 2009).The Markov chain is defined by the state space and the corresponding transition probabilities (Damjan, 2009). Set of states have been described in the Markov Chain model as shown in equation 2-1.

$${\mathbf{S}}_{\mathbf{t}}=\{\mathbf{S}0,\mathbf{S}1,\mathbf{S}2,\dots \dots \mathbf{S}\mathbf{n}\}$$

(2- 1)

Where S_t is the current state, while in the next step it changed into S_j with transitional probabilities p_ij, the state S in the Markov Chain model could also be used for the determination of state S_t+1 in the system by using the formula given in equation 2-2, 2-3 and 2-4.

$${\sum }_{\mathit{j}=1}^{\mathit{m}}{\mathit{P}}_{\mathit{i}\mathit{j}}=1,\mathit{i},\mathit{j}=\mathrm{1,2},\dots ,\mathit{m}$$

(2- 2)

$$\mathit{P}=\left(\mathit{p}\mathit{i}\mathit{j}\right)=\left[\begin{array}{ccc}\mathit{P}11& \mathit{P}12\cdots & \mathit{P}1\mathit{m}\\ \mathit{P}21& \mathit{P}12& \mathit{P}2\mathit{m}\\ \mathit{P}\mathit{m}1& \mathit{P}\mathit{n}2& \mathit{P}\mathit{m}\mathit{m}\end{array}\right]$$

(2- 3)

Where P_ij is the transitional probability from one land use type to another, however the value of P_ij is within the range from$1\le Pij\le 0$, where m is the land use type studied.

$${\mathbf{L}}_{\mathbf{t}+1}=\mathbf{P}\mathbf{i}\mathbf{j}\times {\mathbf{L}}_{\mathbf{t}}$$

(2- 4)

Here L_t is the current land use status and t; t+1 is the time. In order to obtain both the transition and probability matrix of land use type, this study applied the Markov Chain analysis during two periods including 2002-2012 and 2022-2040. The CA-Markov model resulted in three steps of prediction for the possible changes in LULCC: (1) calculating transition matrices from LULCC map of 2022; (2) calculating the LULC transition potential maps and (3) using the CA model to anticipate the geographical distribution of LULC by applying it to the variations matrices as well as the potential transition maps.

2.7. Land Use Change Analysis

To compute Land Use/Land Cover Change (LULCC) from satellite imagery, a systematic approach was employed involving data acquisition, pre-processing, classification, and change detection. Firstly, multiple-period satellite images were collected and processed to normalize radiometric values and correct geometric corrections. Then applied a supervised classification approaches, which classify pixels into different land cover classes including forest land, barren land,grassland and urban areas using training samples taken from field surveys and high-resolution image interpretation. Changes in land cover are revealed by post-classification comparisons of classified images collected at different times using either object-based or pixel-based change detection techniques. Validation data was used to evaluate accuracy in order to improve classifications and guarantee consistency. In order to quantify change areas and comprehend the causes and effects of land use and land cover change over time, the findings were evaluated using statistical analysis and mapping. Land use change was determined using equation 2-5.

$$\mathbf{S}\mathbf{i}={\displaystyle \frac{{\mathbf{L}\mathbf{U}}_{(\mathbf{i},\mathbf{t}\mathbf{2})}-{\mathbf{L}\mathbf{U}\mathbf{A}}_{\mathbf{i}}}{{\mathbf{L}\mathbf{U}}_{(\mathbf{i},\mathbf{t}\mathbf{1})}}}\times {\displaystyle \frac{\mathbf{1}}{{\mathbf{t}}_{\mathbf{2}}-{\mathbf{t}}_{\mathbf{1}}}}\times \mathbf{100}\mathbf{\%}$$

(2- 5)

Si stands for change in land use, LU(i,t1) for land use change at a earlier date, LU(i,t2) for land use change at a later date, and LUAi represents areas with no change at all.

2.7. Biomass Carbon Change

Biomass carbon was assessed by field inventory in temperate coniferous forest (TCF) and mixed forest(MF). A total 60 sample plots of 30m × 30m were taken randomly in the study area (30 plots in each forest type). In each sample plot, the diameter of all trees greater than 4 cm was measured at breast height (DBH). We measured the height of 25 trees per species, which was enough to create a DBH/height function for a particular plot. The height of the trees were measured using Abney's level. According to equation 2-6 to determine the volume of all trees present in the sampling sites calculations were carried out (Williams, Phalan et al. 2018). Total biomass kgha-1 was estimated by dividing biomass by the wood density kg m-3 of a specific species in accordance with the recommendations suggested by the Intergovernmental Panel on Climate Change (IPCC) in 2006.

$$\mathbf{B}\mathbf{A}\left({\mathbf{m}}^{\mathbf{2}}{\mathbf{h}}^{-\mathbf{1}}\right)={\displaystyle \frac{\mathbf{D}\mathbf{e}\mathbf{n}\mathbf{s}\mathbf{i}\mathbf{t}\mathbf{y}}{\mathbf{h}\mathbf{a}}}\times \mathbf{D}\mathbf{C}\times \mathsf{\pi }{\mathbf{r}}^{\mathbf{2}}$$

(2- 6)

Tree volume (TV) was measured by multiplying BA and height (h) and form factor (f.f) , as shown in eq 2-11, formula(Woltz, Peneva-Reed et al. 2022).

$$\mathbf{T}\mathbf{V}\left({\mathbf{m}}^{\mathbf{3}}{\mathbf{h}}^{-\mathbf{1}}\right)=\mathbf{B}\mathbf{A}\left({\mathbf{m}}^{\mathbf{2}}{\mathbf{h}}^{-\mathbf{1}}\right)\times \mathbf{h}\left(\mathbf{m}\right)\times \mathbf{f}.\mathbf{f}$$

(2- 7)

Here f.f is the volume of cylinder to actual tree form(Zeng, et al. 2022) and f.f, 0.37 was used for conifers and 0.6 for broad leaved species in sub-continent as explained by(Zhou, Hartemink et al. 2019). Wood density (W.D) ton per m³ of particular species was used to calculate stem biomass (BM_S) by using equation 2-8. The W.D of each species was determined according to (Olorunfemi, Komolafe et al. 2019).

$${\mathbf{B}\mathbf{M}}_{\mathbf{S}}\left(\mathbf{M}\mathbf{g}{\mathbf{h}\mathbf{a}}^{-\mathbf{1}}\right)=\mathbf{S}\mathbf{V}\left({\mathbf{m}}^{\mathbf{3}}{\mathbf{h}}^{-\mathbf{1}}\right)\times \mathbf{W}\mathbf{D}\left(\mathbf{M}\mathbf{g}{\mathbf{m}}^{-\mathbf{3}}\right)$$

(2- 8)

Above Ground Biomass (BM_ABG) $Mg{ha}^{-1}$was measured by using equation 2-9. The Biomass expansion factor (BEF) as described by (Kauffman, Hughes et al. 2009) is 1.51 for temperate forests of Pakistan.

AGB = Volume × Wood Density × BEF

(2- 8)

In each sample plot, a carbon values of below-ground vegetation (BGV) sub-plot measuring 3 m by 3 m was created in order to calculate the biomass. Each plot's vegetation was first harvested, their fresh weight was measured, and the collected items were all placed into the labelled bags for subsequent examination. After 48 hours of drying in an oven at 72 °C, the samples' dried weights were calculated for biomass. Similar to litter, dead wood, and cones in the sub-plot were gathered, they were likewise subject to the same USV for biomass carbon assessment process described above. Using a soil auger and a core sampler, three replications of the soil samples were taken from the chosen plots at depths of 0–15 cm and 15–30 cm. In the field, each soil sample was weighed and labelled. The bulk density of each soil sample was determined using a core sampler with a volume of 0.0001256 m3 (diameter = 4 cm, height = 10 cm). The soil samples were taken to the lab for further examination. The technique of soil oxidized organic carbon described by (Tiwari, Singh et al. 2019) was used to measure the amount of organic material in the soil. Utilizing a total plant biomass convertible factor, which represents the average biomass content of the plant, carbon was calculated. The carbon content of all plant biomass, as determined by this conversion factor, was estimated to be 50% (Rashid, Bhat et al. 2017). According to (Morreale, Thompson et al. 2021) the bulk density, carbon concentration, and layer thickness were multiplied to get the soil's (MgC ha-1) carbon values.

2.9. Carbon Loss Assessment

The deforestation resulted in the conversion of forestland into grassland and built-up. Change in carbon is determined by Equation 2-10 (Gandhi and Sundarapandian 2017).

$$∆\mathbf{C}={\mathbf{T}\mathbf{o}\mathbf{t}\mathbf{a}\mathbf{l} \mathbf{F}\mathbf{o}\mathbf{r}\mathbf{e}\mathbf{s}\mathbf{t} \mathbf{C}\mathbf{a}\mathbf{r}\mathbf{b}\mathbf{o}\mathbf{n}}_{\mathbf{t}2}-{\mathbf{T}\mathbf{o}\mathbf{t}\mathbf{a}\mathbf{l} \mathbf{F}\mathbf{o}\mathbf{r}\mathbf{e}\mathbf{s}\mathbf{t} \mathbf{C}\mathbf{a}\mathbf{r}\mathbf{b}\mathbf{o}\mathbf{n}}_{\mathbf{t}1}$$

(2- 10)

The annual rate of carbon loss is determined by Equation 4-6:

$$\mathsf{\Delta }\mathbf{C}={\displaystyle \frac{{(\mathbf{T}\mathbf{o}\mathbf{t}\mathbf{a}\mathbf{l} \mathbf{F}\mathbf{o}\mathbf{r}\mathbf{e}\mathbf{s}\mathbf{t} \mathbf{C}\mathbf{a}\mathbf{r}\mathbf{b}\mathbf{o}\mathbf{n}}_{\mathbf{t}\mathbf{2}}-{\mathbf{T}\mathbf{o}\mathbf{t}\mathbf{a}\mathbf{l} \mathbf{F}\mathbf{o}\mathbf{r}\mathbf{e}\mathbf{s}\mathbf{t} \mathbf{C}\mathbf{a}\mathbf{r}\mathbf{b}\mathbf{o}\mathbf{n}}_{\mathbf{t}\mathbf{1}})}{{(\mathbf{t}}_{\mathbf{2}}-{\mathbf{t}}_{\mathbf{1}})}}$$

(2- 11)

ΔC = Carbon loss or gain

TFC_t2 = Total forest carbon at time t₂

TFC_t1 = Total forest carbon at time t₁

3. Results and Analysis

3.1. Land Use Land Change Analysis

The results of LULCC showed that the mixed forest decreased from 1253.79 (ha) in 1992 to 943.21 (ha) in 2002 with -0.75% loss. In 2022 conifer forest was reduced from 1510.2 (ha) to 1161.09 ha or approximately -7.68%. Similarly, the bare land was reduced from 601.65 ha in 2002 to 499.79 ha in 2012. In 2022 it was increased to 689.13 (ha) In 2022 conifer forests were reduced from 1510.2 (ha) to 1161.09 ha and for it is projected to 1295.45 (ha) in 2040. Similarly, the mixed forest slighlty increased from 943.21 ha in 2002 to 1152.74 ha in 2012. On the other hand, from 1129.71 (ha) in 2022, mixed forest is projected to reduce to 1019.56 (ha) in 2040 see Figure 2-1, Figure 2-1a. The coniferous forest were reduced from 1510.20 ha to 1161.09 ha and grassland was slighly increased to 65.79 ha from 31.6 ha in 1992 then again showed a downward trend and reduced to 46.71 ha. The main cause of the gradual loss of coniferous forests, bare land and grasslands over time was urbanization and deforestation in the study area which is the result of numerous human activities that modify the natural environment. First of all, deforestation damages the fragile ecosystems that these woods support in addition to removing trees. Coniferous forests and grasslands are frequently transformed into urban areas to meet the demands of expanding inhabitants and infrastructure. The ongoing construction surrounding the national park causes biodiversity loss, habitat fragmentation, and changes to regional climatic patterns. Table 3-1 shows that the builtup area was increased continuously throughout the last decades and on the other hand it clearly states that other land use cover classes were reduced between the respective years 1992 and 2022. The built up area was recorded a increment up to 368.23 ha in 2022 and further increased to 416.29 ha in 2040. The actual landcover map of 2022 shows that the builtup area was increased significantly throughtout the study period and bare land also showing a upward trend with area increased from 493,56 ha to 689.13 due to the deforestation in the study area. While conifer forest and mixed forest decreased from 1992-2022.The projected and predicted land cover map was generated using Markov chain model which shows that the projected area in 2022 mixed forest,conifer forest,bare land,builtup area and grassland was 1102.9 ha,1182.45,432.2,398.76 and 34.67 respectively.The builtup area will be expanded as predicted in 2040 up to 416.29 ha due to which the other land use types may reduced more and cause deforestation in the area. The time series graph in Figure 2-1c states clearly more about the downward and upward trend of different land use land change classess.

The results showed that during the period 1992-2002, CF were decreased by 24.7% (0.247% yr^-1), during 2002-2022, mixed forest was increased by 10.495 (0.104% yr^-1,) and 8.13% (0.183% yr^-1) respectively as shown in Table 4-1. Apart from a decreased and slightly increased in area under mixed forest, major changes were also observed for barren land, conifer forest, grassland and settlements (Table 4-1 and Figure 2-1, Figure 2-2). Our results showed that a total of 371.94 ha forest land (CF, MF) had been converted into other land uses from 1992 to 2022 and Indicating the increasing scale of anthropogenic pressure on forest resources of the Park. During the predicted period 2022-2040 the built up area will increase more due to population growth and economic factors.

Table 4.-1. (%) Statistics of LULCC per decade and annual.

Land use categories	1992-2002 change (%)	Annual change (%)	2002-2012 change (%)	Annual change (%)	2012-2022 change (%)	Annual change (%)	2022-2040 change (%)	Annual change (%)
Mixed forest	-16.7	-0.167	+10.49	+0.104	+18.3	+0.183	+16.47	+0.164
Conifer forest	+3.98	+0.039	-0.65	-0.0065	-25.57	-0.255	+11.57	+0.11
Bare land	+25.7	+0.257	-14.6	-0.146	+30.07	+0.30	-59.4	- 0.59
Built-up/settlements	+0.92	+0.09	-10.12	-0.101	+13.80	+0.138	+34.18	+0.341
Grassland	+65.3	+0.653	+25.62	+0.256	-29.00	-0.290	-41.04	-0.410

Table 4.-1. (%) Statistics of LULCC per decade and annual.

Land use categories	1992-2002 change (%)	Annual change (%)	2002-2012 change (%)	Annual change (%)	2012-2022 change (%)	Annual change (%)	2022-2040 change (%)	Annual change (%)
Mixed forest	-16.7	-0.167	+10.49	+0.104	+18.3	+0.183	+16.47	+0.164
Conifer forest	+3.98	+0.039	-0.65	-0.0065	-25.57	-0.255	+11.57	+0.11
Bare land	+25.7	+0.257	-14.6	-0.146	+30.07	+0.30	-59.4	- 0.59
Built-up/settlements	+0.92	+0.09	-10.12	-0.101	+13.80	+0.138	+34.18	+0.341
Grassland	+65.3	+0.653	+25.62	+0.256	-29.00	-0.290	-41.04	-0.410

Table 5.-1. LUC statistics from 1992 to 2022.

Land use type	Land use changes to other land uses	Area change (ha) (1992-2022)	Total area change (ha)
Conifer forest	Bare land	104.94	472.95
	Mixed forest	241.92
	Settlement	126.09
Mix forest	Bare land	7.740	61.7
	Settlements	5.490
	Conifer forest	48.51
Built-up/Settlements	Mixed forest	0.81	29.34
	Bare land	22.86
	Grassland	5.67
Grassland	Conifer forest	2.340	25.11
	Bare land	14.40
	Settlements	8.37

Table 5.-1. LUC statistics from 1992 to 2022.

Land use type	Land use changes to other land uses	Area change (ha) (1992-2022)	Total area change (ha)
Conifer forest	Bare land	104.94	472.95
	Mixed forest	241.92
	Settlement	126.09
Mix forest	Bare land	7.740	61.7
	Settlements	5.490
	Conifer forest	48.51
Built-up/Settlements	Mixed forest	0.81	29.34
	Bare land	22.86
	Grassland	5.67
Grassland	Conifer forest	2.340	25.11
	Bare land	14.40
	Settlements	8.37

3.1. Accuracy Assessment of LULCC

Accuracy assessment is an important step in the for processing land use change analysis. It contains the information value of the resulting data to a user. Ground reference data were obtained through the use of stratified random sampling. Every plot had a design measuring 30 by 30 meters. Based on the information gathered in the field, a LULCC class identification number was assigned to each plot. The classifiers were trained using roughly 70% of the reference data, with the remaining 30% being utilized to evaluate the accuracy assessment. Several vector and grid data types, including roadways, population centers, and DEM 30m, which was utilized to create slopes, aspects, and elevation to create the vector map layers. LULCC classification was applied in this investigation using five classes. Using confusion matrices, incorporating commission and omission errors, overall accuracy, and kappa value, classification accuracy processes were established. The results showed kappa coefficient 0.845, 0.862, 0.873 and 0.928 for the year 1992, 2002, 2012 and 2022 respectively shown in Table 5-2. The simulated results showed kappa statistics values to range from 0.7044 to 0.9086. while transitional matrix statistics for 2022 are shown in Table 5-3.The use of a CA-Markov model resulted in three steps of prediction for the upcoming change in LULC: (1) calculating transition matrices from LULCC maps for 2022; (2) calculating the LULC transition potential maps and (3) using the CA model to anticipate the geographical distribution of LULC by applying it to the variations matrices as well as the potential transition maps.

The results from model validation also showed the total agreement was 0.93, while the simulation error was 0.82 as shown in Table 5-2a, while the probability transitional matrix statistics for 2022 are shown in Table 5-3.

Table 5.-2. Accuracy Assessment for 1992-2022.

Classified map	Overall Accuracy	Kappa Coefficient
LULC map 1992	86.81 %	0.845
LULC map 2002	89.43%	0.862
LULC map 2012	90.43%	0.873
LULC map 2022	94.08%	0.928

Table 5.-2. Accuracy Assessment for 1992-2022.

Classified map	Overall Accuracy	Kappa Coefficient
LULC map 1992	86.81 %	0.845
LULC map 2002	89.43%	0.862
LULC map 2012	90.43%	0.873
LULC map 2022	94.08%	0.928

Table 5.-2a. Validation results and Agreement/Disagreement.

Agreement/Disagreement (%)
Allocation disagreement	7.45
Quantity disagreement	5.16
Allocation agreement	14.95
Quantity agreement	62.10
Chance agreement	10.34

Table 5.-2a. Validation results and Agreement/Disagreement.

Agreement/Disagreement (%)
Allocation disagreement	7.45
Quantity disagreement	5.16
Allocation agreement	14.95
Quantity agreement	62.10
Chance agreement	10.34

Table 5.-3. Transition matrix, simulation for LULCC in year 2022.

	Probability to change from-to
	Mixed forest	Conifer forest	Bare land	Built-up/Settlements	Grassland
Mixed forest	0.8836	0.1014	0.0108	0.0030	0.0012
Conifer forest	0.2726	0.6831	0.0367	0.0006	0.0071
Bare land	0.1178	0.2279	0.5990	0.0082	0.0471
Builtup/Settlements	0.0012	0.4710	0.1765	0.2796	0.0718
Grassland	0.0000	0.0000	0.7272	0.0000	0.2728

Table 5.-3. Transition matrix, simulation for LULCC in year 2022.

	Probability to change from-to
	Mixed forest	Conifer forest	Bare land	Built-up/Settlements	Grassland
Mixed forest	0.8836	0.1014	0.0108	0.0030	0.0012
Conifer forest	0.2726	0.6831	0.0367	0.0006	0.0071
Bare land	0.1178	0.2279	0.5990	0.0082	0.0471
Builtup/Settlements	0.0012	0.4710	0.1765	0.2796	0.0718
Grassland	0.0000	0.0000	0.7272	0.0000	0.2728

3.3. Carbon Dynamics

The estimated carbon value recorded in TCF was 135.19 ±9.74 MgC ha^-1, whereas it was 86.43 ±8.25 MgC ha^-1 in MF (Table 6-1). In each forest, the higher value of carbon was recorded for Above ground vegetation (ABGV) followed by soil, below ground vegetation and litter dead wood. The results of the total carbon stock (Table 6-2) and fig 2-2 showed that in 1992 the total forest carbon stock was 527 GgC. Similarly, the total forest carbon stock was reduced in 2002, 2012 and 2022 i.e. 522 GgC, 489 GgC, and 435 GgC respectively. These findings revealed the decreasing trend in carbon value associated with LULC.

4. Discussion and Conclusion

4.1. Land Use Land Cover and Carbon Change

In Pakistan's temperate and subtropical areas, where LULCC and deforestation are ongoing processes, natural forests are primarily found. The current study approaches the land use land cover changes effect on forest biomass and their drivers who is responsible for the carbon value due to LULCC. The rise in Settlements(SM) and deforestation between 1992 and 2022 reveals the extent of LULCC in the park. In comparison to our base year 1992, we discovered that SM had grown by 25.9% during the past 30 years. According to Figure 2-3, the SM for Ayubia was increased from 315.06 to 368.23 between 1992 and 2022. Additionally, the number of SM in and around the park during 2022-2040 will increase more due to settlement growth (Waseem, Mohammad et al. 2005). The predicted LULCC map for the year 2040 is generated using Markov Chain Model with earlier and later classified images of 1992 and 2022 years. The procedure of model validation combines the contents of two data sources and takes into consideration the features of actual and simulated LULCC data. The LULC for 2040 is predicted using the LULCC data between 1992 and 2022, the spatial variables selected for simulation, and the transition probability matrix, and a kappa value of 0.87 was attained.

For instance, the Conifer forests(CF) were decreased from 1510.2 hectares in the base year 1992 to 1161.09 ha in 2022. This downward trend in CF was comparable to that noted by Malik. A (2003; 48.4% in 1994, 54.8% in 2000). The uncontrolled and illegal logging by local villages for fuelwood, livestock grazing, and forest fires may be to blame for the transformation of CF into BL in and around the park. Fires can alter vegetation structure and composition, leading to changes in land cover and biodiversity and overgrazing by livestock can lead to soil erosion, compaction, and degradation of vegetation cover. This reduces the ability of ecosystems to sequester carbon and maintain soil fertility. Carbon loss from deforestation, land degradation, and changes in land use practices contributes to greenhouse gas emissions. This exacerbates climate change impacts and affects ecosystem services such as carbon sequestration, water regulation, and biodiversity conservation. Increased carbon emissions can amplify climate change effects, potentially leading to more frequent and severe forest fires, altered precipitation patterns, and shifts in vegetation dynamics that further impact carbon storage and ecosystem resilience. Integrated land management approaches are crucial for mitigating these impacts. Strategies such as sustainable forestry practices, controlled grazing management, urban planning for resilience, and carbon-friendly policies can help reduce carbon loss, enhance ecosystem resilience, and support sustainable development goals. This can influence subsequent land use decisions and practices. As is the case elsewhere in the world, forest fires frequently occur in ANP, particularly during the dry summer months(Ahmad and Javed 2008). These fires usually damage the CF more severely as compared to MF. The continuous reduction in the area under CF and increase in areas under bare land and built up showed CF to have been mainly converted into bare land, grassland and built up areas . Land use and Land Cover Change (LULCC) due to urbanization, agriculture, and forestry practices have profound impacts on landscapes, ecosystems, and socio-economic dynamics. Urban expansion was a major driver of land cover change in the study area. It leads to the loss of vegetative cover, increase in impervious surfaces, and modification of natural ecosystems. The spatial extent and intensity of urbanization are often quantified through satellite imagery and urban growth models. Expansion of agricultural land often involves deforestation or conversion of grasslands and wetlands. Changes in land use patterns, crop types, irrigation practices, and use of agrochemicals significantly impact terrestrial ecosystems and their services. Logging and deforestation contribute directly to changes in forest cover and carbon dynamics. Afforestation and reforestation efforts can mitigate some impacts by restoring forest cover and enhancing carbon sequestration capacities.

Carbon was computed with values ranging from 86.43 ± 8.25 in MF to 135.19 ± 9.74 in TCF. Similar findings were reported for the carbon store in the temperate woods of Nathiagali (Amir, Muhammad et al. 2022). Similarly, the stated figure by (Gandhi and Sundarapandian 2017) from the Indian Himalayan area conflicts with our current findings. The loss of carbon due to the conversion of forestland to other land uses decreased from 527 GgC in 1992 to 435 GgC in 2022. The same findings were previously published and demonstrated a loss of 32.7 PgC from Brazil's Amazonian protected areas in 2014 (MohanRajan, Loganathan et al. 2020). Similar to this, more than 3 million km² of US protected areas' forest, grassland, and shrub land were converted to agricultural land between 1700 and 2005, resulting in a 10,607 TgC loss (Núñez-Regueiro, Siddiqui et al. 2021) . This decrease in carbon corroborates the idea that forest land's capacity to store carbon is diminished as a result of conversion to other land uses. Forest conversion brought about by deforestation, agricultural development, and settlement growth linked to population increase led to a drop in carbon stock. Over these different periods, the highest rate of carbon losses was recorded for the period 2012-2022 as shown in Table 6-3. This period also coincided with the maximum of population growth and human migration within the study area. Although this LULC affects carbon loss, the intensity is much higher where CF is converted into SM, or BL (Ohler, Schreiner et al. 2023)

Between 1992 and 2022, Ayubia's population grew from 18,000 to nearly 55,000. In addition to other environmental and socioeconomic issues, these demographic characteristics are the main causes of LUCC in and around ANP (PKS 2019). Additionally, it was predicted that 16,000 people lived in nearby villages in 1992, but by 2022, this number had grown to almost 55,000, spread across 12 communities (Iqbal and Khan 2014). The viability of these natural environments is also seriously threatened by this demographic strain (Li, Yang et al. 2018, Ohler, Schreiner et al. 2023). The growth of built-up regions and the extraction of wood (fuel wood, artisanal, and commercial logging) are also influenced by this demographic pressure and contribute to the alteration of the forest cover and the deterioration of natural ecosystems(Jiyuan, Mingliang et al. 2002, Justice, Gutman et al. 2015). Anthropogenic activities pose serious hazards to the world's forest cover, according to several research (Rokityanskiy, Benítez et al. 2007, Sasmito, Taillardat et al. 2019). Historically, this increase in population and other underlying factors particularly in developing countries have been observed to be the most important driving force of LULC (Rosero-Bixby and Palloni 1998). According to (Tan, Ge et al. 2023), these underlying drivers of LULCC are poverty, population growth, and socio-economic situation. Increased population often results in the migration of people in search of farmlands and forests are easier to convert to agricultural land (Luo,et al. 2022). Rural areas around the ANP have experienced unchecked urban growth, which has forced these areas' rural areas into the region's wooded areas. Because National Parks feature high-quality limestone and sandstone needed for cement manufacture, stone mining for this purpose also contributes to carbon loss (Avtar, Tsusaka et al. 2020). Around the ANP, there are over 300 stone quarries, which have a negative impact on the flora, fauna, and 58 kinds of animals and over 900 plant species (The News, july 13, 2014). In a related research,(Khan, Khan et al. 2021) found that the primary causes of deforestation in the ANP region were uncontrolled harvesting, forest disease, and the removal of wood for commercial and domestic use.

With average carbon stored, the computed value of carbon was varied from 135.19 ±9.74 MgC ha-1 in TCF to 86.43 ±8.25 MgC ha-1 in MF. Comparable outcomes for the carbon stock in the subtropical/chir pine forests of Ghoragali, Pakistan, were documented (Nizami 2012). In a similar results, our current estimations agree with the figure published by Rana et al. (1989) from the Indian Himalayan area. According to Tao et al. (2015) and Zhu et al. (2021), the specific change process involved improving the regional carbon storage by increasing the amount of LULCC type with high-density carbon and decreasing the amount of LULC type low-density carbon. Furthermore, our analysis revealed that over the previous 30 years, the area covered by built-up areas rose from 315 to 368 hectares, resulting in a loss of 87 GgC in carbon storage. However, as built-up area expanded and conifer forests declined, the amount of low-density carbon land cover and high-density carbon land cover decreased, reducing the region's capacity to store carbon (Tao et al., 2015, Zhu et al., 2021). While the growth of built-up land in the National Park area has been directly associated with social and economic development, it has also had the large effect of reducing carbon storage. This is consistent with the findings of other scholars. Li et al (2020) takes Wuhan City as an example, the study showed that urban expansion led to a decreased of 2.62 Tg of carbon storage in Wuhan from 2000 to 2015(Li et al., 2020). Xiang et al (2022) similarly found that the encroachment of forest land by built-up land from 2000 to 2020 was the main reason for the decrease in carbon storage in Chongqing (–5.78 Tg)(Xiang et al., 2022). The increase of built up area as one of the main sources of deforestation in ayubia national park. This conclusion was consistent with that of (Sajjad et al.,2018) who observed that the primary cause of deforestation in Tehsil Barawal, Dir Upper, Pakistan, was the increase of settlements. According to a study by Ali et al., which was endorsed by the Forest Department, the primary causes of carbon loss was deforestation in Basho Valley, Skardu, are ineffective administration and illegal timber harvesting for commercial purpose.

Due to falling dry needles from conifer forests throughout the summer and dry season around May, forest fires frequently happen in the ANP. According to (de Groot, et al. 2009), the amount of carbon lost each year as a result of forest fires ranges from 50.2% to 57.6%, depending on the intensity of the fire. The majority of forest fires that started in Pakistan's temperate woods were unintentional as a result of carelessness and leisure activities. In Pakistan's temperate woods, surface fires have been seen to be trending upward(Amir, Muhammad et al. 2022). According to (Kauffman, Hughes et al. 2009), forest fires are known to decrease soil carbon, soil erosion, understory vegetation, and soil bareness. Livestock grazing and trampling have a negative impact on plant regeneration and contribute to LULCC(Ahmad and Nizami 2015).The research also covered the main causes of land use change and carbon loss in ANP, such as the exploitation of fuel wood for commercial purposes, population pressure and settlement growth, which have a negative impact on the park's environment. Additionally, decreased forest area is a result of the construction of new infrastructure. These elements, when combined with lack governmental regulations and poor land use management strategies, encourage a decline in the amount of carbon stored in the ANP's forest biomass. This deforestation process may be slowed down, stopped, and even reversed by better forest management, sensible forest conservation laws, the creation of land use management plans, and urbanization controls.

5. Conclusion

Interpretation of multi temporal satellite imagery was found to be a practical method for determining the dynamics of land use change. For the years 1992, 2002, 2012, and 2022, Geographic Information System (GIS) tools were used to analyze the land use change patterns of Ayubia National Park (ANP) present in the study area. The results showed that the area under mixed forest was increased slightly while the area under conifer forests were decreased continuously. According to the carbon shift dynamics of ANP, the area covered by conifer forests were decreased by 25.57 percent (or 0.25 percent per year) from 1992 to 2022.It is predicted that the builtup area will increased to 416.79 ha in 2040 which likely decreased the forestland from 1510.20 in intial period to 1295.45 ha in 2040 and mixed forest from 1253.79 to 1019.56 in 2040.Similarly, settlements were increased by 13.80% (0.13% yr^-1) respectively, during the same period. To estimate confidence in classification model accuracy bootstrap approach was used. This bootstrap approach is relatively easy to implement and can be applied to any automated classification method used regardless of the specific classes within a region. Moreover, it provides a robust, non-parametric approach to defining class confidence that can enable users to gain a better understanding of classification errors. The conversion of forests into other land uses were resulted in the reduction of carbon value highlighting that land use changes significantly altered carbon dynamics of the Park. The findings demonstrates that the anticipated temporal and spatial dynamics of LULCC changes and their impact on forest carbon throughout the park between 1992 and 2022. Furthermore the land changes were predicted for the year 2040 which can be linked to a number of causes such as resource exploitation, urbanization, and intensification of agriculture. Overall, this study highlighted the possible difficulties brought on by the LULCC on the forest carbon in the national park, and the findings highlight the necessity of making the best use of available resources, implementing wise land management techniques, and keeping an eye on the rapid urbanization that is occurring. Highlights methodological limitations encountered, such as challenges in integrating different data sources, modeling uncertainties, or assumptions made during analysis. For future implementations research encourages that integrates socio-economic factors (e.g., demographic trends, policy interventions) with biophysical models apply to get the better understanding of the drivers and consequences of LULCC. Overall, the findings contribute to the broader understanding of LULCC changes in rapidly developing regions and emphasize the importance of incorporating sustainable land management strategies for the long-term well-being of the Ayubia National Park and other similar regions across the world.