Urban Growth Patterns and Forest Carbon Dynamics in the Metropolitan Twin Cities of Islamabad and Rawalpindi, Pakistan

The unchecked and unplanned expansion of urban areas has led to the conversion of millions of green areas to gray areas. The recent urban growth patterns of Pakistan’s metropolitan twin cities, Islamabad and Rawalpindi, is a matter of concern for the surrounding green areas. The present study aimed to categorize and quantify the land-use and land-cover change (LULCC) patterns and the corresponding impacts on the forest carbon dynamics around Islamabad and Rawalpindi. Multispectral satellite images for the year 1990 (Landsat 5 TM) and 2020 (Landsat 8 OLI) were used to determine, quantify, and compare the LULCC inside and around the twin metropolitan cities. Field inventory surveys in the reserved forests of Rawalpindi and Islamabad were also conducted to determine the amount of stored carbon in these forests. Our results showed an accelerated annual urban expansion (i.e., an increase in the built-up area) of 16.49% and 26.72% in Rawalpindi and Islamabad, respectively, during the study period. Similarly, the amount of barren land and agricultural land was reduced at an annual rate of 2.08% and 2.18%, respectively, in Rawalpindi and 0.25% and 1.04% in Islamabad. A reduction in the area of barren mountains also occurred at an annual of 2.26% in Islamabad, while it increased by 4.16% in Rawalpindi. The amount of carbon stored in the reserved forests of Islamabad stood at 139.17 ± 12.15 Mg C/ha while that of Rawalpindi was 110.4 ± 13.79 Mg C/ha. In addition, total stored forest carbon was found to have decreased from 544.70 Gg C to 218.05 Gg C in Rawalpindi, while in Islamabad it increased from 2779.64 Gg C to 3548.16 Gg C. Investment in ecological urban planning, sustainable cities, and appropriate land-use planning is recommended to curb the degradation and conversion of the surrounding green areas of Rawalpindi and Islamabad.


Study Area
The present study was carried out in Rawalpindi and Islam 33.662883° Lat, 73.086373° Lon (Figure 1). With an estimated population lion inhabitants [12,33], Islamabad is the capital city of Pakistan and s area of 906.84 km 2 [34]. The minimum and maximum annual temperatur −3 °C in winter to 45 °C in summer, while mean annual rainfall ranges in winter and 1846 mm in summer [35]. On the other hand, with an esti of about 2 million inhabitants, Rawalpindi covers an area of 259 km 2 [3 is rich in flora and fauna. The dominant forest species of the area include Pinus roxburghii, Quercus incana, Acacia modesta, Acacia nilotica, Olea ferrug spp. These species act as important carbon sinks and help conserve soil. T the study area are the habitat of important wildlife species, including Macaca mulatta, Muntiacus muntjac, Canis aureus, Catreus wallichii, and Al The twin cities are also of great economic importance, given that Islam pindi contributed USD 2.86 billion and USD 3.96 billion, respectively, to tic product (GDP) of Pakistan, out of the total GDP of USD 263.7 billion

Satellite Data and Pre-Processing
Landsat 8 OLI (operational land imager) for 2020 and Landsat 5 TM per) (WRS Path 150, Row 37) were downloaded from the US Geologic website. Cloud-free images for the month of February were used, becau

Satellite Data and Pre-Processing
Landsat 8 OLI (operational land imager) for 2020 and Landsat 5 TM (thematic mapper) (WRS Path 150, Row 37) were downloaded from the US Geological Survey (USGS) website. Cloud-free images for the month of February were used, because this month has less Sustainability 2021, 13, 12842 4 of 14 atmospheric haze and fewer ground reflectance changes [22]. To geometrically register the image with an ariel photograph, all images were first projected to the Universal Transverse Mercator (UTM) projection system (zone 49 N), and a 30 m resolution of the pixel size was set. The images were then corrected radiometrically to convert the digital number (DN) to a numerical value by using ENVI 5.3 in Equations (1) and (2), as follows: where L λ represents radiance, QCAL are the DNs, QCALmin is the quantized calibrated minimum pixel value, Lmin λ is the spectral radiance at QCALmin, QCALmax is the maximum calibrated quantized pixel value, and Lmax λ is the spectral radiance scale at QCALmax. In Equation (2), p λ represents top of atmosphere (TOA) reflectance, d represents the Earth-Sun distance, while ESUN λ is the solar irradiance and θ is the elevation of the Sun in degrees. The TOA reflectance to surface reflectance conversion was performed using the atmospheric correction method before performing LULCC detection. By using a fast line-of-sight atmospheric analysis of hypercubes (FLAASH), the corrected image provides better results for the estimation of LULCC and C dynamics [37]. Furthermore, it also reduces the effects of water vapors, aerosols, humidity, and haze, according to [38,39].
The study area was clipped using a vector shapefile prepared from a topographic sheet (1:20,000) obtained from the Survey of Pakistan (SOP).

Extraction of LULCC Maps and Analysis
For LULCC classification, the images were enhanced with feature enhancement in order to better perceive land-use separability. In addition, the maximum likelihood method for image classification was used in ENVI 5.3, Equation (3): where gi = the land use class, x = n-dimension data, n = the no. of bands, p(w i ) = the probability that class wi appeared in the wi class, |Σi| = the co-variance matrix data for the wi class, Σ i −1 = the inverse matrix, and mi = the vector.
To analyze land-use cover, the study area was divided into six classes: forest land (FL), agricultural land (AL), barren land (BL), built-up area (BA), barren mountains (BM), and water body (WB). A description of the classification categories is given in Table 1. The image accuracy was accessed using 100 random training samples for each class. The training samples were selected from the Landsat image gallery because most of the LULCC types maintained stable characteristics. Based on this Landsat imagery, it was necessary to assess the accuracy of the classified image [40]. The combination of Band 543 in Landsat 5 TM and Band 643 in Landsat 8 provided the visually distinguished features of all the classes [41].

Accuracy Assessment
The confusion matrix in Envi 5.3 was used to measure the accuracy of the classification, in combination with Equation (4): where K is the kappa coefficient, the class number is denoted by i, where N is the total number of classified values compared to truth values, m i,i is the number of values belonging to the truth class i that have also been classified as class i (i.e., values found along the diagonal of the confusion matrix), C i is the total number of predicted values belonging to class i, and G i is the total number of truth values belonging to class i.

Carbon Dynamics
In order to calculate stored carbon, random sampling was performed across 300 circular sample plots of about 0.1 hectares (ha) in Rawalpindi and Islamabad, as described by [42]. In Rawalpindi, the samples were taken from the Loi Bher (RF), Dhamial (RF), Adiala (RF), and Takht Pari (RF) forest reserves. Meanwhile, in Islamabad, the sampling was conducted around the Margallah Hills National Park (MHNP), Tumair (RF), and Banni Gala (RF) forest reserves. All trees greater than 2.5 cm diameter at breast height (DBH) were measured in the plots. The heights of 20 trees per species were measured using an Abney level, which was sufficient enough to construct the plot-specific diameter-height function. The volume (m 3 /ha) of all the trees were calculated according to [43]. The stem biomass (SBM), the biomass expansion factor (BEF), and the root-shoor ratio R, as described by [43], were used to calculate the tree biomass (TBM), as per Equation (5): Furthermore, from the center of the 0.1 ha sample plot, the understory biomass (USV BM ), consisting of mainly herbs and shrubs, was calculated by taking a sample of all the herbs and shrubs found within the 9 m radius sample plots, which were harvested and kept in labelled bags for further analysis. The wet weight was measured at harvesting time, while the dry weight was measured after drying them in an oven at a temperature of 72 • C for 48 h. The biomass was calculated as dry weight as described by [44]. Similarly, litter and dead wood were collected for the calculation of the dead wood and litter biomass (DWLT BM ), while the ecosystem biomass (E BM ) was calculated using Equation (6): Soil samples at the depth of 0-15 cm and 15-30 cm were taken with the help of a soil auger. The volume of all soil samples was measured using a core sample of 0.0001256 m 3 , with a diameter = 4 cm and length = 10 cm. All the soil samples were weighed to measure their bulk density and organic content. Furthermore, the tree biomass was converted to carbon by using the factor 0.5, according to [45]. For carbon change, the product of carbon ha −1 and total forest cover (ha) gives the total carbon stock of the area at specific temporal periods [46][47][48]. The percentage change in carbon stock at different time periods was calculated as per Equation (7): where, ∆CC is the change in carbon content, CC 2 is the latest period, and CC 1 is the base year period.

Normalized Difference Vegetation Index (NDVI)
The NDVI is an important tool that helps to determine the spatiotemporal change in forest cover [49]. NDVI is the most reliable means of determining the intensity of forest degradation and vegetation loss [50]. The NDVI was calculated for the two periods of 1990 and 2020 for Islamabad and Rawalpindi using Equation (8), in ENVI 5.3: where NIR is the near-infrared band and RED is the red band.

Urban Growth
We used GIS and remote sensing data and techniques to categorize and quantify the LULCC across the Rawalpindi and Islamabad areas for a thirty year period (1990 to 2020). Our results showed that the twin cities have experienced an enormous expansion of urban growth during the past three decades. For instance, in 1990, the built-up areas in Islamabad and Rawalpindi were 3018.42 ha and 2561.02 ha, respectively, and have expanded by 26.72% and 16.49 a year, respectively, for the two cities during the study period (Tables 2 and 3). In 1990, the built-up area in Islamabad was 3.33% (3018.42 ha) of its total area, which expanded to 30.01% by 2020 (27,218.06 ha) (Table 3). Similarly, for Rawalpindi, it was 9.88% in 1990, and has increased to 58.08% by 2020 (Figures 2 and 3). For the same study period, barren land has decreased in Rawalpindi at a rate of 2.08% annually as compared to 0.25% for Islamabad (Table 3). Decreasing trends were observed in agricultural land at a rate of 2.18% and 1.04% per annum for Rawalpindi and Islamabad, respectively. Furthermore, the share of barren mountain area decreased from 31.88% to 10.16% in Islamabad, but increased from 4.98% to 11.22% in Rawalpindi ( Figure 3). Interestingly, an increase of 0.92% per annum in the forest area has been identified for the Islamabad area during the study period. Forest land in Rawalpindi has been reduced from 4930.55 ha in 1990 to 1973.82 ha in 2020 (Table 2). which expanded to 30.01% by 2020 (27,218.06 ha) (Table 3). Similarly, for Rawalpindi, it was 9.88% in 1990, and has increased to 58.08% by 2020 (Figures 2 and 3). For the same study period, barren land has decreased in Rawalpindi at a rate of 2.08% annually as compared to 0.25% for Islamabad (Table 3). Decreasing trends were observed in agricultural land at a rate of 2.18% and 1.04% per annum for Rawalpindi and Islamabad, respectively. Furthermore, the share of barren mountain area decreased from 31.88% to 10.16% in Islamabad, but increased from 4.98% to 11.22% in Rawalpindi (Figure 3). Interestingly, an increase of 0.92% per annum in the forest area has been identified for the Islamabad area during the study period. Forest land in Rawalpindi has been reduced from 4930.55 ha in 1990 to 1973.82 ha in 2020 ( Table 2).
The accuracy of these figures has been assessed, and the Kappa coefficient for Islamabad was found to be 0.92 and 0.87 for the years 2020 and 1990, respectively. The Kappa coefficient for Rawalpindi was 0.90 and 0.89 for the years 2020 and 1990, respectively.   Decreasing trends were observed in agricultural land at a rate of 2.18% and 1.04% per annum for Rawalpindi and Islamabad, respectively. Furthermore, the share of barren mountain area decreased from 31.88% to 10.16% in Islamabad, but increased from 4.98% to 11.22% in Rawalpindi (Figure 3). Interestingly, an increase of 0.92% per annum in the forest area has been identified for the Islamabad area during the study period. Forest land in Rawalpindi has been reduced from 4930.55 ha in 1990 to 1973.82 ha in 2020 ( Table 2).
The accuracy of these figures has been assessed, and the Kappa coefficient for Islamabad was found to be 0.92 and 0.87 for the years 2020 and 1990, respectively. The Kappa coefficient for Rawalpindi was 0.90 and 0.89 for the years 2020 and 1990, respectively.  The accuracy of these figures has been assessed, and the Kappa coefficient for Islamabad was found to be 0.92 and 0.87 for the years 2020 and 1990, respectively. The Kappa coefficient for Rawalpindi was 0.90 and 0.89 for the years 2020 and 1990, respectively.

Forest Carbon Dynamics
In Islamabad, stored carbon in the forest's upper story, lower story, and in dead and dry matter was found to be 78.86 ± 9.94 Mg/ha, 2.1 ± 1.1 Mg/ha, and 1.9 ± 0.7 Mg/ha, respectively. The stored carbon in the ecosystem was 139.17 ± 10.48 Mg/ha for Islamabad. For Rawalpindi, the carbon stored in the upper story, lower story, and in dead and dry wood was 65.16 ± 10.39 Mg/ha, 1.5 ± 1.12 Mg/ha, and 1.7 ± 0.72 Mg/ha, respectively. The stored carbon in the ecosystem for the Rawalpindi forest area was 110.47 ± 11.07 Mg/ha ( Figure 4). In 1990, the total carbon stored in Islamabad was 2779.64 Gg, which has increased to 354,816 Gg in 2020, whereas in Rawalpindi, a decreasing trend was observed in the stored carbon, which has dropped to 218.5 Gg C in 2020 from 544.70 Gg C in 1990.
The mean value of the NDVI for Islamabad was 0.07 in 1990, and this has increased to 0.09 in 2020 ( Figure 5). The mean NDVI value calculated for Rawalpindi was 0.03 in 1990, which has decreased to 0.01 in 2020 ( Figure 6).
For Rawalpindi, the carbon stored in the upper story, lower story, and in dead and dry wood was 65.16 ± 10.39 Mg/ha, 1.5 ± 1.12 Mg/ha, and 1.7 ± 0.72 Mg/ha, respectively. The stored carbon in the ecosystem for the Rawalpindi forest area was 110.47 ± 11.07 Mg/ha ( Figure 4). In 1990, the total carbon stored in Islamabad was 2779.64 Gg, which has increased to 354,816 Gg in 2020, whereas in Rawalpindi, a decreasing trend was observed in the stored carbon, which has dropped to 218.5 Gg C in 2020 from 544.70 Gg C in 1990.

OR PEER REVIEW
The mean value of the NDVI for Islamabad was 0.07 in 1990, and this has in to 0.09 in 2020 ( Figure 5). The mean NDVI value calculated for Rawalpindi wa 1990, which has decreased to 0.01 in 2020 ( Figure 6).

R PEER REVIEW
The mean value of the NDVI for Islamabad was 0.07 in 1990, and this has in to 0.09 in 2020 ( Figure 5). The mean NDVI value calculated for Rawalpindi was 1990, which has decreased to 0.01 in 2020 ( Figure 6).

Discussion
RS and GIS tools are the most effective methods around today to determine the spatiotemporal LULCC both in urban and rural areas [51]. LULCC has become an important factor in determining the socioeconomic and environmental development of an area [52]. Forests act as a carbon sink. Therefore, a change in forest cover directly affects the climatic and environmental conditions of a region [53]. Urban forests and green areas are more prone to LULCC due to different developmental activities, and so therefore categorizing and quantifying LULCC in urban areas is very important [54]. Our results indicate that between 1990 and 2020, the amount of settlement and built-up area has increased considerably. This increase has exerted pressure on other land-use classes, such as BL, AL, BM, etc., in the study area. In 1990, the built-up area around Islamabad was 3.3%, which has increased to 30.01% in 2020. In Rawalpindi, the built-up area was 9.8% of the total area in 1990, which has expanded to about 58.80% in 2020. An annual urban expansion rate of 0.88% for Islamabad and 1.63% for Rawalpindi has therefore been identified during the study period. This expansion could be attributed to an increase in real estate agents and housing companies operating within the study area. The areas have experienced a huge human influx, as people from the rural areas migrated to the capital for different purposes, including education, jobs, and other business reasons. According to [12], the accelerated urban expansion can be attributed to governmental policies which favor urban development and real estate businesses. Similarly, these governmental policies concomitantly encourage and support industrial growth and infrastructure development within and around the twin cities [55]. Industrial growth and infrastructure development are found to be important factors in the urbanization of the twin cities [12]. Population size directly affects socioeconomic development [56,57], as it determines the available of resources in relation to the social, economic, and environmental stability of any locality [58]. During the study period, an increase of 36,870.23 ha in the built-up area has been recorded in the twin cities ( Table 2). The rapid increase in the built-up area could be attributed to the increased human population in the country, now at 207.8 million people [59], with the population size of the twin cities increasing at a rapid rate. During the study period, in Islamabad, the population increased from 0.342 million to 1.129 million, while in Rawalpindi, the population has increased from 1.087 million to 2.236 million during the study period (Figure 7). Sustainability 2021, 13, x FOR PEER REVIEW relation to the social, economic, and environmental stability of any locality [58]. the study period, an increase of 36,870.23 ha in the built-up area has been recorde twin cities ( Table 2). The rapid increase in the built-up area could be attributed increased human population in the country, now at 207.8 million people [59], w population size of the twin cities increasing at a rapid rate. During the study pe Islamabad, the population increased from 0.342 million to 1.129 million, while in pindi, the population has increased from 1.087 million to 2.236 million during th period (Figure 7).  The master plan of Islamabad was made by Constantinos A. Doxiadis and Doxiadis Associates in 1950, based on the hierarchical community system, with an urban roads spacing network and sustainable conservation [60]. On the other hand, Rawalpindi has a semi-planned urbanization pattern [34,61,62]. The management of the Islamabad is under the supervision of the Capital Development Authority (CDA), which has all the protocols and guidelines for the management and expansion of the city [63], while Rawalpindi is managed by the Rawalpindi Development Authority (RDA) [64]. In addition, the weak enforcement of the law in Rawalpindi has stimulated the encroachment of government land and the development of illegal housing societies, causing haphazard urbanization [65]. In 2017, the Capital Development Authority of Islamabad (CDA) demolished 50 commercial buildings, including marriage halls and shops on illegal land worth PKR 9 billion (USD 51,660,765), as well as 1670 ha of land costing PKR 300 billion (USD 1,722,025,500), which was retrieved from the land grabbers [66]. No such action has been taken by the RDA. However, the main focus of the CDA was to maintain the standard infrastructure and beauty of the city, safeguarding areas such as Rawal Lake, the Rose and Jasmine Garden, Saidpur Village, etc. [67]. The housing sectors in Islamabad, i.e., E, F, G, H, and I, are subdivided into smaller division of 2 km 2 , i.e., G6, G7, G8, G9, etc. All these subdivisions have one community center at the midpoint, with highways at the border as well as link roads [68]. The drinking water, solid waste management, and sanitation facilities are controlled by the management who are stationed at the community center.
From 1990 to 2020, barren and agricultural land decreased at the rate of 0.25% and 1.24% annually in Islamabad, respectively, and 2.08% and 2.18% annually in Rawalpindi. This may be related to the fact that small-scale agriculture had been the main stay and source of livelihood for local communities [12]. People also used to keep land adjacent to their homesteads for animal rearing and gardening. In later years, governmental policies began to change in favor of urban development and ushered in a rapid increase in land value. These open lands around the homesteads, along with other less economically valuable agricultural and barren lands, were all converted into urban sprawls and settlements.
Our results show that, compared to Rawalpindi, forest land in Islamabad not only remained intact, but has increased at an annual rate of 0.92% during the study period, with a positive increase in the NDVI value. This positive change in the NDVI for Islamabad could be attributed to strict management and law enforcement by the capital's forest management department. Strict forest protection law enforcement in the urban area can lead to the development of protected areas, including green belts, research areas, and botanical gardens [69]. The establishment of the Margallah Hills National Park (MHNP) in Islamabad is an example of urban forest protection. Furthermore, green belts, roadside plantations, and botanical gardens increased in Islamabad, supporting the green forest cover [70]. This increase in forest area is also related to the planned forest policies of the federal government in Islamabad, whose main objective is increasing the forest area and ensuring its conservation. On the other hand, the protected and reserved forests of Rawalpindi are challenged by illegal wood harvesting, artisanal logging, and forest land encroachment [71]. For instance, illegal encroachment has been observed in the Loi Bher Reserve Forest, located in the center of Rawalpindi. This could be one the reasons for the decrease in the NDVI value for Rawalpindi during the study period. Similarly, unplanned urbanization and forest land degradation in Rawalpindi also resulted in the decrease in the amount of stored carbon in its forests. Islamabad, on the other hand, witnessed an increase in stored carbon in its forests between 1990 and 2020.
The forests of Rawalpindi, including both protected and reserved forests, are very much prone to encroachment from illegal building by housing societies, i.e., Bahira Town, Media Town, the Police Foundation, etc. Encroachment in Loi Bher Reserve Forest, located in the center of Rawalpindi, was noted from 1990 to 2000, and as a result, the mean value of the NDVI decreased from 0.678 to 0.552 in the observed period (see Figure 8).
The forests of Rawalpindi, including both protected and reserved forests, much prone to encroachment from illegal building by housing societies, i.e., Bah Media Town, the Police Foundation, etc. Encroachment in Loi Bher Reserve Fores in the center of Rawalpindi, was noted from 1990 to 2000, and as a result, the me of the NDVI decreased from 0.678 to 0.552 in the observed period (see Figure 8) Due to unplanned, uncontrolled, and unchecked urbanization in Rawalpind est carbon decreased from 544.70 Gg C to 218.05 Gg C, where in Islamabad, t carbon in the forest was 2779.64 Gg C in 1990 and increased to 3548.16 Gg C in 2 could be attributed to the increase of forest land in Islamabad and the decrease land in Rawalpindi for the same period of time. Deforestation reduces the ability to act as carbon sinks, as it is directly related to the biomass of forest. Due to defo carbon loss was observed to be 31.33 GgC/ha annually in the subtropical forest stan [72]. However, contradicting results were observed in China, in which th sink was increased from 36 MgC/ha to 41 MgC/ha from 1981 to 2000, due to fore vation policies [73]. Although the urbanization in Islamabad was more rapid com Rawalpindi, it does not affect the forest land, due to the active implication of env tal policies. Due to unplanned, uncontrolled, and unchecked urbanization in Rawalpindi, the forest carbon decreased from 544.70 Gg C to 218.05 Gg C, where in Islamabad, the stored carbon in the forest was 2779.64 Gg C in 1990 and increased to 3548.16 Gg C in 2020. This could be attributed to the increase of forest land in Islamabad and the decrease of forest land in Rawalpindi for the same period of time. Deforestation reduces the ability of forests to act as carbon sinks, as it is directly related to the biomass of forest. Due to deforestation, carbon loss was observed to be 31.33 GgC/ha annually in the subtropical forests of Pakistan [72]. However, contradicting results were observed in China, in which the carbon sink was increased from 36 MgC/ha to 41 MgC/ha from 1981 to 2000, due to forest conservation policies [73]. Although the urbanization in Islamabad was more rapid compared to Rawalpindi, it does not affect the forest land, due to the active implication of environmental policies.

Conclusions
The findings of this study show that the twin cities of Rawalpindi and Islamabad have experienced rapid urbanization during the study period. Despite expansion in the urban area, an increase in the vegetation cover was observed in Islamabad. In Rawalpindi, a decrease in the forest area was recorded. The relative increase in Islamabad's vegetation cover could be attributed to better management and implementation of urban forest protection policies. Our study signifies the importance of forest conservation in the cosmopolitan cities. Our results would assist policy makers in designing conservation and management policies that could potentially support sustainable urban ecosystems. Based on our findings, we recommend that researchers conduct a comparative study focusing on the status and dynamics of land surface temperature for the twin cities.