Monitoring Carbon Stock and Land-Use Change in 5000-Year-Old Juniper Forest Stand of Ziarat, Balochistan, through a Synergistic Approach

: The Juniper forest reserve of Ziarat is one of the biggest Juniperus forests in the world. This study assessed the land-use changes and carbon stock of Ziarat. Different types of carbon pools were quantiﬁed in terms of storage in the study area in tons/ha i.e., above ground, soil, shrubs and litter. The Juniper species of this forest is putatively called Juniperus excelsa Beiberstein. To estimate aboveground biomass, different allometric equations were applied. Average above ground carbon stock of the forest was estimated as 8.34 ton/ha, 7.79 ton/ha and 8.4 ton/ha using each equation. Average carbon stock in soil, shrubs and litter was calculated as 24.35 ton/ha, 0.05 ton/ha and 1.52 ton/ha, respectively. Based on our results, soil carbon stock in the Juniper forest of Ziarat came out to be higher than the living biomass. Furthermore, the spatio-temporal classiﬁed maps for Ziarat showed that forest area has signiﬁcantly decreased, while agricultural and barren lands increased from 1988 to 2018. This was supported by the fact that estimated carbon stock also showed a decreasing pattern between the evaluation periods of 1988 to 2018. Furthermore, the trend for land use and carbon stock was estimated post 2018 using a linear prediction model. The results corroborate the assumption that under a business as usual scenario, it is highly likely that the Juniperus forest will severely decline.


Introduction
The present forest area of the world is 4.06 billion ha and it has lost 178 million ha of forest since 1990, presenting a decrease of 4.2% [1]. It is estimated that the world's forest area has decreased from 31.6% to 31% of the total Earth's land area. With approximately 5% of the land area under forest, Pakistan finds its place among the low-forested countries in the world [2][3][4]. Major forest types in Pakistan include coastal mangroves, riverine forests, sub-tropical scrub forests, moist temperate conifer forests, dry temperate conifer forests, and irrigated plantations including linear plantations [5]. Based on the 'Forestry Sector Master Plan of 1992 , 4,200,000 ha (4.8%) of Pakistan has natural forest [6]. Compared to that, in 2015, FAO (Food and Agriculture Organization) reported 1.9% or about 1,472,000 ha of forest area in Pakistan. This translates to a loss 2,728,000 ha of the forest cover, which has also been documented by the Pakistan Forestry Outlook study in 2009 and the Forest Resource Assessment in 2015. Evidently, there is a declining trend in the forest cover while agricultural and urban areas are expanding [7].

Monitoring Land-Use and Land-Cover (LULC) Changes
The archived Landsat imagery provide a unique opportunity to identify land use changes over the past 30 years [31]. Because of that, for the detection of the land cover for District Ziarat, images of Landsat 4-5 TM (Thematic Mapper) and Landsat 8 OLI (Operational Land Imager) with no or minimum cloud cover (0-10%) were selected and downloaded for the years of 1988, 1998, 2008, and 2018 from the USGS (United States Geological Survey) Earth Explorer website (https://earthexplorer.usgs.gov/). Two tiles (path row number 152/39 and 153/39) were downloaded for each period to cover the study area. The downloaded images belong to the month of September (Table 1) since it was found very feasible with almost no clouds present.

Sensor
Year

Monitoring Land-Use and Land-Cover (LULC) Changes
The archived Landsat imagery provide a unique opportunity to identify land use changes over the past 30 years [31]. Because of that, for the detection of the land cover for District Ziarat, images of Landsat 4-5 TM (Thematic Mapper) and Landsat 8 OLI (Operational Land Imager) with no or minimum cloud cover (0-10%) were selected and downloaded for the years of 1988, 1998, 2008, and 2018 from the USGS (United States Geological Survey) Earth Explorer website (https://earthexplorer.usgs.gov/). Two tiles (path row number 152/39 and 153/39) were downloaded for each period to cover the study area. The downloaded images belong to the month of September (Table 1) since it was found very feasible with almost no clouds present. The images were automatically atmospherically corrected using Semi-Automated Classification plugin in QGIS 3.6.2 (QGIS Development Team 2002). Images along with the metadata file were uploaded in the pre-processing Table DOS1 (Dark Object Subtraction) algorithm option was checked okay before running the processing chain. Atmospherically corrected images were opened in Arc-Map 10.3 (Esri 2014, Redlands, CA, USA), where the bands were stacked using the composite band tab in the Image Analysis window. The stacked images were mosaiced followed by extraction of study area using extract by mask algorithm. The study area was thoroughly examined on Google Earth maps. Both true and false color of Landsat images were observed to ensure that the visual interpretation was done correctly. Once the area was studied, pixels were assigned to the classes and classification was carried out via maximum likelihood classification. The area for each class was calculated in ArcMap and the data were further analyzed by using MS Excel. Line graphs and pie charts were formulated for easy analysis of the data.

Carbon Stock Assessment
This study estimated carbon from all the pools namely above ground tree biomass, shrubs, litter and soil C. Below ground biomass, however, was not considered in this study due to limited resources.
The field work was carried out in the month of February 2018. Six random clusters were selected throughout the study area, each with one primary and four secondary sampling units. The number of total plots were 30, however, due to absence of trees in some areas, the number of plots were reduced to 21 plots as depicted in Figure 2 below. The images were automatically atmospherically corrected using Semi-Automated Classification plugin in QGIS 3.6.2 (QGIS Development Team 2002). Images along with the metadata file were uploaded in the pre-processing Table DOS1 (Dark Object Subtraction) algorithm option was checked okay before running the processing chain. Atmospherically corrected images were opened in Arc-Map 10.3 (Esri 2014, Redlands, CA, USA), where the bands were stacked using the composite band tab in the Image Analysis window. The stacked images were mosaiced followed by extraction of study area using extract by mask algorithm. The study area was thoroughly examined on Google Earth maps. Both true and false color of Landsat images were observed to ensure that the visual interpretation was done correctly. Once the area was studied, pixels were assigned to the classes and classification was carried out via maximum likelihood classification. The area for each class was calculated in ArcMap and the data were further analyzed by using MS Excel. Line graphs and pie charts were formulated for easy analysis of the data.

Carbon Stock Assessment
This study estimated carbon from all the pools namely above ground tree biomass, shrubs, litter and soil C. Below ground biomass, however, was not considered in this study due to limited resources.
The field work was carried out in the month of February 2018. Six random clusters were selected throughout the study area, each with one primary and four secondary sampling units. The number of total plots were 30, however, due to absence of trees in some areas, the number of plots were reduced to 21 plots as depicted in Figure 2 below.  The above-ground biomass was calculated for all the plots, however, only one soil sample was collected from the primary plot of each cluster thus giving us six soil samples. Plots for above-ground biomass, shrubs, litter and soil had a radius of 17.8 m, 5.64 m, 0.56 m respectively, keeping in mind that soil and litter plots have the same radius. Shrubs and litter samples were collected as per their availability, and their location data points are displayed in Figure 2.
Tree height was measured using a clinometer and diameter was measured at breast height of 1.3 m using a caliper. All these data including elevation, site coordinates, time and date of sampling were also recorded. Living biomass data was entered into excel sheets to run the allometric equations. The allometric equations (see Table 2) used for biomass estimation were taken from Jenkin et al. [15], Ali [32] and Chave et al. [33]. The Jenkin et al. [15] equation has been used by the forest department of Gilgit Baltistan [14]. The equation used by Ali [32] and Chave et al. [33] has also been used by WWF (World Wide Fund) for the National Forest Inventory of Pakistan. Fraction of carbon in the above ground living biomass can be assessed using the following equation adopted from FAO [34].
Total Carbon in above ground biomass (kg) = 0.475 × Above ground Biomass (kg) (1) where the above ground biomass is assessed using the allometric equations listed in Table  2 The above ground living biomass volume formula was derived from the stem biomass formula of FAO [35] as mentioned below: where V s is stem volume and WD is wood density. Similarly, tree volume can be calculated by using following equation: where V t is tree volume and WD is wood density. The carbon stock of each carbon pool in tons was converted into tons per hectare using the following formula: C in ton/ha pool = C in tons/Area of the plot pool (5) Total above ground carbon (million tons) for the total forest area, in each respective year of 1988, 1998, 2008 and 2018 was estimated using the following equation: Total forest Above Ground Carbon (million tons) = C in ton/ha × Area of forest in ha for each year (6) Total forest carbon for the years before 2018 were estimated using average carbon stock of the three allometric equations used. Additionally, total above carbon for the year 2028, 2038 and 2048 was also estimated using the above formula. Area in hectares of these future years were inferred using the linear forecast model in Microsoft excel worksheet 2013 (Supplementary Material). Total carbon of the forest (million tons) for each year of past and future is presented in results.
The method for estimating total carbon stock of the forest by taking the product of average carbon ha −1 and the total forest cover (ha) at a specific temporal period was also used by Mannan et al. [18].

Linear Forecast Model
Linear forecast function uses a linear regression method to predict future values based on historical figures. It is a method of defining the relationship between two or more variables in a way that changes in dependent variable can be accounted by changes in the independent variable. In our study carbon stock and forest area are the dependent variables and time (against year) is the independent variable.

Soil Carbon Stocks
Soil samples were collected at three depths; 0 to10 cm, 10 to 20 cm and 20 to 30 cm. An auger was used for the collection of soil samples. A total of Six soil samples were collected, one from each cluster. Bulk density of each soil sample was calculated in g/cm 3 . For the determination of carbon in soil, the Walkley-Black titration procedure was applied. To find soil C in grams, the following equation was utilized.
where SC Soil Carbon, BD is the Bulk Density, SOC is the Soil Organic Carbon, and HT is the Horizon Thickness. Soil carbon in grams was converted into tons and further into ton/hectare by dividing it by the plot area.

LULC Changes 1988-2018
Land cover maps ( Figure 3) showed drastic changes in the land use pattern of Ziarat District. Forest area, which is the major focus of the study decreased from 21.5% in 1988 to 15.5% in 2018 showing a decrease of 6% in the total area. This represents significant deforestation over the past decades. Similarly, agriculture area increased from 1.5% in 1988 to 3.5% in 2018. This might be due to an increasing population and demand for agricultural production ( Figure 3). As per the housing and population census, the population of Ziarat increased from 32,196 people in 1981 to 160,422 people in 2017. Barren land also significantly increased from 76.5% in 1988 to 81% in 2018 due to deforestation.
As shown in Figure 4 (land use change trend and forecast), the forest area decreased from 71,005 hectares in 1988 to 50,311 hectares in 2018, depicting a loss through deforestation of 20,694 hectares in the past 30 years. On the other hand, the agriculture land increased from 4848 hectares in 1988 to 11,625 hectares in 2018 representing an increase of 6777 hectares. Figure 4 also shows the linear forecast values inferred for next 30 years. The forecast analysis shows that the forest area will decrease to 29,976.5 hectares in 2048, while the agricultural and the barren land area will increase to 17,826.7 hectares and 279,259 hectares, respectively, in the coming three decades. Table 3 below shows the accuracy assessment of the land-use maps of Ziarat district. Accuracy assessment was performed in Arc Map using the confusion matrix, where 140 training points were selected for each temporal map. The confusion matrix shows producer's and user's accuracy for each land-use category of the respective yearly maps. Accuracy from the perspective of the map designer is termed as the producer's accuracy and it shows the probability of whether land cover has been classified correctly. User's accuracy tells us the correctness of the map from a user's viewpoint. It tells us the probability of how often the classified class on the map will be present on the ground [36,37]. It is quite high and above 90%. The overall accuracy of the classification for 1988, 1998, 2008 and 2018 was 0.97, 0.98, 0.98, and 0.97. Kappa analysis was undertaken if the performance of the classification did well in comparison to randomly assigning values. This ranged from −1 to 1 and a value of zero represents that the classification was no better random value assignment [36,37]. The Kappa coefficient was 0.96, 0.98, 0.97 and 0.96, respectively, for each consecutive year.  Figure 4 also shows the linear forecast values inferred for next 30 years. The forecast analysis shows that the forest area will decrease to 29,976.5 hectares in 2048, while the agricultural and the barren land area will increase to 17,826.7 hectares and 279,259 hectares, respectively, in the coming three decades.    Table 3 below shows the accuracy assessment of the land-use maps of Ziarat district. Accuracy assessment was performed in Arc Map using the confusion matrix, where 140 training points were selected for each temporal map. The confusion matrix shows producer's and user's accuracy for each land-use category of the respective yearly maps. Accuracy from the perspective of the map designer is termed as the producer's accuracy and it shows the probability of whether land cover has been classified correctly. User's accuracy tells us the correctness of the map from a user's viewpoint. It tells us the probability of how often the classified class on the map will be present on the ground [36,37].     Table 3 below shows the accuracy assessment of the land-use maps of Ziarat district. Accuracy assessment was performed in Arc Map using the confusion matrix, where 140 training points were selected for each temporal map. The confusion matrix shows producer's and user's accuracy for each land-use category of the respective yearly maps. Accuracy from the perspective of the map designer is termed as the producer's accuracy and it shows the probability of whether land cover has been classified correctly. User's accuracy tells us the correctness of the map from a user's viewpoint. It tells us the probability of how often the classified class on the map will be present on the ground [36,37].      The overall average height of all the plots was 4.92 m while the overall average diameter was 15.22 cm. The overall average basal area was 0.76 m 2 . Table 5 below shows the total biomass and the total volume of juniper forest using three different above-ground allometric equations. Total biomass calculated from all three equations were 37,281.51 kg, 24,812.12 kg and 37,518.67 kg respectively. The volume It can be noted that there was a small difference between the values of the equations of Jenkin et al. [15] and Chave et al. [33]. But there was a visible difference in the biomass and volume calculated using the equation of Ali [32].  Figure 5 represents the correlation between the parameters of tree as specified in Tables 3 and 4. There is a good correlation between the height and diameter I-e 0.89, and this does not vary among the allometric equation. This shows that the height of the tree increases with the increase in diameter of the juniper tree.  Table 6 shows an overall plot summary of carbon stock in above ground pool using the three equations as mentioned earlier. Plot 17 contains the highest amount of above ground biomass i.e., 22.48, 22.58 and 25.30 ton/ha, and the lowest above ground biomass is found in plot 21 i.e., 2.37, 1.99, and 2.01 ton/ha. The overall average living biomass of the forest of all plots using the three equations is 8.34, 7.79 and 8.4 ton/ha, respectively. Table 6. Comparison of the above ground carbon using the allometric equations from Jenkin et al. [15], Ali [32], and Chave et al. [33].

ABG (t/ha) Using Chave et al. (2005) Allometric Equation
Plot 1 Table 7 below shows the carbon estimation of the pools: soil, shrubs and litter. The respective plots where there is highest and lowest carbon stock have been highlighted. The overall average soil carbon of the juniper forest is 24.35 ton/ha. The overall carbon stored in shrubs and litter is 4.66 and 1.52 ton/ha, respectively.  The bold values are the high and the low values in all the data set.

Discussion
Ziarat juniper forest was declared a biosphere reserve in 2013. Such forests will play a pivotal role in future in furthering the cause of carbon sequestration [10]. However, they face a visible threat of deforestation and forest degradation [38]. Major drivers of deforestation and forest degradation are population expansion, agriculture intensification, fuel wood consumption, poor regeneration, illegal cutting, overgrazing, canopy dieback, mistletoe attack, periodic drought and medicinal use of the juniper tree [38,39]. For identifying deforestation, a quantitative assessment of the land-use change was performed for four land use classes namely, forest, agriculture, water and barren land. Focus of the study was reduction in forest cover, which according to our results was accounted for by the increase in agricultural area.
There is an inverse relationship between forest resources and both population as well as the amount of land used for agriculture. With population increase, more and more forests are encroached upon for harvesting firewood and timber, thereby increasing the non-forest area in Ziarat district [40]. The population increase has led to the increase in dependency on agriculture [12], thus having a synergetic effect on the forest. Besides other drivers, Archeothobium oxycedri, known as dwarf mistletoe, also damages the Juniperus species [41]. Sarangzai et al. [42] also reported the spread of the same parasitic and infectious plant in the juniper forest of Ziarat. Harvesting fuelwood is also one of the major drivers of deforestation since there is no other source of fuel to keep the local people warm in winter [43]. Even though natural gas has been provided to Ziarat, the pressure in gas pipelines is extremely low thus compelling residents to cut trees in the winter season [44]. Most of the juniper forest have low seedlings and show no adequate regeneration thus slowing down recovery time [45].
Urban areas were not classified, because most of the construction sites are either muddy houses or wooden-built sites, thus making it extremely difficult to separate its spectral information in Landsat imagery. The area also has myriad soil forms and colors resulting in the overlapping of the pixels during classification. Even the forest area was difficult to classify since the juniper forest is a sparse forest with less tree density. However, it was classified very well when observed simultaneously with Google Earth imagery. The results of this study are also comparable to the results of a study conducted by WWF and SUPPARCO (Pakistan Space & Upper Atmosphere Research Commission) in 2012 on the Juniper Forest of Ziarat using SPOT data [11]. This suggests that the Landsat data are good for forest classification and land-use change.
Ismail et al. [14] estimated the above ground carbon per hectare of Juniperus communis. The total number of juniper trees counted were 278 and the estimated average carbon was 1.96 ton/ha with a total basal area of 12.28 m 2 . This value is substantially less compared to our estimated carbon of 8.34 ton/ha and a total basal area of 16.06 m 2 using the same equation for Juniperus excelsa.
Since no specific equations exist for the Juniperus excelsa species, it is recommended that the forest department in collaboration with academia develop an allometric equation for it. If developing an equation is not possible or not within the capacity, it is recommended that the three equations used in the study are used for any future carbon stock study. The three equations used in this study gave very similar results.
The data were collected from 21 plots having 585 trees. The results show that the soil contains more carbon than trees, which may be attributed to the compact soil of the area, low temperature conditions and the less dense/sparse nature of the forest. It may also be due to the age of the forest, which is very old thus, accumulating humus for thousands of years. One of the major reasons for less biomass in the tree is small height and the very sparse nature of the forest reserve.
The juniper forest of Ziarat is an extremely rare forest and needs to be protected and sustained. The forest area has decreased from 10,025 hectares in 1960 to 53,092 hectares in 2010 and is a clear manifestation of threats faced by the juniper forest of Ziarat, Balochistan. [5,11]. This comprises multiple factors. The first one is the population of the area that has dramatically increased (about 400%) from 1981 to 2017 as per the population and housing census. Besides marble mining at few places, there is no such industrial activity in the Ziarat district. Therefore, the population must depend on natural resources for their daily subsistence. The agriculture area has also increased (2% from 1988 to 2018) as indicated in the previous section of this study. Other factors include medicinal use of berries, climate change, timber extraction, tourism and poor forest management by the forest department. All these factors are clearly posing an enhanced threat to preservation of these ancient forests of Ziarat. Keeping the current scenario and past practices in view, it is most probable that the forest area will keep on shrinking as mentioned by various studies presented in Figure 6.
Moreover, the protection regime should not be only limited to Ziarat but also to the juniper forest resources in Loralai, Harnai, Quetta and Pishin. Despite all the services the juniper forest provides, it has not been taken care of in a viable way resulting in its deforestation and degradation. If a similar trend of disregard for this precious and ancient forest continues, we might only study about this archeological heritage in archives. For this purpose, the forest department must take practical steps for the conservation of the juniper forest of Balochistan.

Conclusions
The study concluded that the ancient juniper forest of Ziarat is an important carbon sink, storing a significant amount of carbon in all its pools i.e., above-ground live tree biomass, soil, shrubs and litter. The soil of the juniper forest stores more carbon than the living biomass due to low canopy cover and scattered growth of the trees. Similarly, the litter contains more carbon stock than the shrubs since the quantity of litter found in the plots was higher than the shrubs. Land-use maps showed a tremendous change in the forest cover of Ziarat Balochistan from 1988 to 2018 with a decreasing forest trend and increasing agricultural trend. Furthermore, carbon stock of juniper has also decreased over the past three decades due to excessive deforestation and forest degradation. Clearly, the juniper forest is threatened and if the same pace of deforestation continues, it is very likely that this forest will soon be wiped out. Additionally, the carbon stock of the juniper forest was assessed using three different equations which gave similar results, so it is suggested that these equations may be used for future carbon stock studies on the juniper forest. It is also concluded that, if no appropriate steps are taken towards the conservation of this forest, we may lose the ancient world biosphere reserve in a very short time.
Supplementary Materials: The following are available online at https://www.mdpi.com/1999-490 7/12/1/51/s1, Table S1: Linear forecast model. Funding: The authors are grateful to Universiti Putra Malaysia (UPM) and NUST for providing financial support from the postgraduate RnD fund from S and T Mega Project to conduct this study.