Monitoring of Urban Landscape Ecology Dynamics of Islamabad Capital Territory (ICT), Pakistan, Over Four Decades (1976–2016)

In the late 1960s, the Islamic Republic of Pakistan’s capital shifted from Karachi to Islamabad, officially named Islamabad Capital Territory (ICT). In this aspect, the ICT is a young city, but undergoing rapid expansion and urbanization, especially in the last two decades. This study reports the measurement and characterization of ICT land cover change dynamics using Landsat satellite imagery for the years 1976, 1990, 2000, 2010, and 2016. Annual rate of change, landscape metrics, and urban forest fragmentation spatiotemporal analyses have been carried out, along with the calculation of the United Nations Sustainable Development Goal (SDG) indicator 11.3.1 Land Consumption Rate to the Population Growth Rate (LCRPGR). The results show consistent increase in the settlement class, with highest annual rate of 8.79% during 2000–2010. Tree cover >40% and <40% canopy decreased at an annual rate of 0.81% and 0.77% between 1976 to 2016, respectively. Forest fragmentation analysis reveals that ‘core forests of >500 acres’ class decreased from 392 km2 (65.41%) to 241 km2 (55%), and ‘patch forest’ class increased from 15 km2 (2.46%) to 20 km2 (4.54%), from 1976 to 2016. The LCRPGR ratio was 0.62 from 1976 to 2000, increasing to 1.36 from 2000 to 2016.


Introduction
The conversion of one land cover type to another is one of the most visible and rapid changes that the earth is experiencing and these conversions have profound social and environmental impacts at multiple scales. In the past, the drivers of land cover changes have been classified into two categories: proximate and distant, or indirect [1]. Recently, the phenomenon of tele-coupling between places has been recognized as key to understanding how distant drivers relate to proximate drivers and influence local landscape changes [2][3][4]. Some of the key drivers of landscape changes include economic, technological, institutional and policy, cultural, and demographic factors [5].
The growth of urban areas leads to land cover change in many parts of the world, especially in developing countries [3,6]. Intense urbanization and increase in anthropogenic activities reflect the scope, intensity, and frequency of human interference, and the changes they cause in ecological processes and systems in the urbanized areas [7]. The urban areas consist of only 1-6% of the earth's land surface, yet they have enormous impacts on the functioning and service of local and global ecosystems, by modifying local climate conditions, eliminating and fragmenting native habitats, generating anthropogenic pollutants, etc. [8]. The spatial pattern of an urban landscape is a result of the Land 2020, 9, 123 3 of 23 Clarke et al. [28] proposed a framework to combine remote sensing and spatial metrics for improved understanding and representation of urban dynamics to come up with alternative conceptions of urban spatial structure and change. Particularly with regards to urban forestry, a few studies have focused towards the assessment, mapping, and monitoring of urban forest parameters in fast-growing cities of developing countries. For example, Gong et al. [29] carried out a 30-year forest fragmentation study over Shenzhen Special Economic Zone (SEZ), a city which was established in 1979 in Southern China. Huang et al. [30] utilized satellite images of 77 metropolitan areas in Asia, US, Europe, Latin America, and Australia to calculate and analyze seven spatial metrics (area weighted mean shape, area weighted mean patch fractal dimension index, centrality, compactness index, compactness index of the largest patch, ratio of open space, and density). According to their analysis, the compactness, density, and regularity of urban areas in developing regions generally exceeded the levels reported throughout developed countries [30]. Dewan et al. [31] studied the dynamics of land use/cover changes though landscape fragmentation analysis in Dhaka Metropolitan, Bangladesh, computing and analyzing the following metrics: Number of patches, Patch density, Landscape shape index, Largest patch index, Mean patch size, Area-weighted mean fractal dimension, Interspersion and juxtaposition, Contagion, and Shannon's diversity index.
In Pakistan, like other developing countries, most urban development is haphazard, typically lacking appropriate planning strategies [10]. Pakistan's urbanization rate is the highest in South Asia, and by 2030, Pakistan will have more people in cities than in rural areas. Growing population and rapid development is causing prime agricultural land to be encroached and also causing loss of tree cover [32][33][34][35][36]. In the late 1960s, the capital of the Islamic Republic of Pakistan was shifted from Karachi to Islamabad (officially named Islamabad Capital Territory (ICT)). The masterplan of ICT was developed by the famous Greek architect and town planner C. A. Doxiadis [37]. In terms of a planned new capital, ICT is similar to planned new post-colonial capitals/relocations as in the cases of Brasilia (Brazil), Nur-Sultan (named as Astana from 1998 to 2019, and Akmola previously) in Kazakhstan, and Canberra (Australia) [38][39][40]. In the recent decades, with various ongoing development activities, ICT has been struggling with rapid urbanization and gigantic levels of pollution from industrial, residential, and transportation sources. In terms of population, ICT is considered as the most diverse city of Pakistan with a large percentage of immigrants and foreigner population [41]. Unprecedented influx of migrants and population increase has resulted in urban sprawl and conversion of fertile agricultural land and green cover into concrete-a clear deviation from the original ICT master plan [39,40]. Uncontrolled population growth in ICT due to rapid urbanization has deteriorated the living environment, and increased the adverse ecological impacts on human health, flora, and fauna [42].

Literature Review-ICT Mapping and Monitoring
In the last 20 years, several studies have been conducted on the ICT, regarding land cover change, biomass estimation, water quality monitoring, and temperature increase using satellite datasets. In this section, we present a synthesis of published work regarding land cover change dynamics in the ICT.
Adeel [43] identified urban growth potential through land use for ICT zone IV (Figure 1), based on SPOT-5 2.5 m panchromatic dataset and population census data, and found that nearly 63% of zone IV carries a 'High' to 'Very High' future growth potential, which is mainly located close to Islamabad Expressway. This work uses satellite imagery and field data from one year (2007) and does not report spatiotemporal change dynamics. Butt et al. [35] studied the metropolitan development in ICT, based on growth direction and expansion trends from the city center, for the period 1972-2009 using Landsat satellite images. Using Principal Components Analysis (PCA), band ratios, and supervised classification methods, they found that the urban development had expanded by 87.31 km 2 in 38 years. Butt et al. [44] conducted a study on land cover change analysis over Simly dam watershed, ICT; the results derived from maximum likelihood supervised classification showed tree cover loss of up to 26% and 6% increase in settlements from 1992-2012, based on Landsat 5 TM and SPOT-5 imagery, respectively. Similarly, the watershed analysis of Rawal dam, ICT using Landsat 5 TM imagery, showed Land 2020, 9, 123 4 of 23 3% degradation of tree cover and 2% gain of settlement from 1992-2012 [45]. Another ICT land cover change dynamics study conducted by Hassan et al. [46] utilized 30 m Landsat 5 TM data for 1992 and 2.5 m SPOT-5 data for 2012, using the maximum likelihood algorithm for image classification. The study revealed a decrease in forest cover of approximately 49% and over 213% gain of settlement area from 1992-2012. Sohail et al. [47] conducted a study to assess the water quality index and analyze the major change in land cover types, vegetation cover, rate of urbanization and its possible impact on groundwater resources, vegetation, and barren land. They used Landsat images for the years 1993,1997,2002,2007,2013, and 2017 for the assessment and mapping of land cover dynamics; according to their findings, from 1993 to 2017, vegetation areas decreased by 101.77 km 2 , surface water was reduced by 1.10 km 2 , barren land was reduced by 2.90 km 2 , while built-up lands expanded by 105.77 km 2 .
A comparison of Beijing, China and ICT for the role of vegetation in "controlling the eco environmental conditions for sustainable urban environment" was performed by Naeem et al. [42], where they used Gaofen-1 (GF-1) and Landsat-8 Operational Land Imager (OLI) satellite imagery with 8 m and 30 m spatial resolution, respectively. They evaluated various scenarios and models for future development to predict future spatial patterns in both cities. Another study was conducted by Naeem et al. [48] to study the association between green space characteristics, analyzed through landscape metrics, and land surface temperature for sustainable urban environments comparing Beijing, China and ICT.
Khalid et al. [49] conducted a study to quantify the decline of forest reserves and associated temperature variations in a relatively unexplored biodiversity hotspot of ICT, the Margalla Hills National Park (MHNP). In this work, Landsat satellite imagery from 1992, 2000, and 2011 was used to monitor the changes in forest cover and statistical significance tests were used to determine the significance of temperature variation associated with a shift in land cover classes. The study finds that deforestation and forest degradation by local communities is an ongoing practice in MHNP; this necessitates the promotion of conservation practices to minimize ecological disturbances here [49]. Batool and Javaid [50] carried out a study on the assessment of Margalla Hills forest by using Landsat imagery for years 2000 and 2018, and report that the forest cover has decreased from 87% in 2000 to 74% in 2018, whereas built-up area has increased from 5% in 2000 to 7% in 2018, and open land in the study area increased from 2% in 2000 to 7% in 2018. Mannan et al. [51] conducted a study using Landsat imagery, Markov Chain, and Cellular Automata on Margalla Hills, focusing on the quantitative assessment of spatiotemporal land use and land cover changes during 1998, 2008, 2018, and a simulation of 2028. In addition, a forest inventory survey was conducted for biomass and carbon sink estimations. This work shows that the forest area has reduced from 409.36 km 2 to 392.31 km 2 and settlement area has increased from 14.97 km 2 to 39.66 km 2 from 1998 to 2018. The average yearly biomass and carbon losses were 50.34 Gg/ha/yr and 31.33 Gg C/ha/yr, respectively.
The ICT is a relatively new and spatially heterogeneous city surrounded by the Himalaya mountainous dense forest as compared to other fast growing and expanding cities like Dhaka, Bangladesh [31,[52][53][54][55], New Delhi, India [56], Beijing, China [42,48,57,58], Shanghai, China [59][60][61], Tokyo, Japan [27,62], etc. Based on the literature review, we observed that most studies of the ICT land cover dynamics have used different remotely sensed data, methods, definitions, and classification schemes, and have provided diverse results. Most studies which have analyzed the land change dynamics in the ICT focus on the overall analysis of land-cover and land-use change, and a detailed analysis of the landscape ecology and urban forestry characteristics is missing. There has further been very little focus in these studies towards urban landscape metrics and indicators of sustainable urban growth such as LCRPGR.

Study Objectives
In this paper, well established, proven, and articulated research methodology, satellite datasets, and definitions of features were adopted with the goal to systematically achieve the following defined objectives: Land 2020, 9, 123 5 of 23

•
Detection, measurement, and characterization of land cover features using Landsat medium resolution freely available satellite data (1976, 1990, 2000, 2010, and 2016) and determine the annual rate of change in land cover classes at 10 years interval.

•
Landscape metrics and forest fragmentation spatiotemporal analysis, to estimate and report changes in the ICT urban ecosystem over forty years.

•
Calculation of SDG indicator number 11.3.1 "Land Consumption Rate to the Population Growth Rate (LCRPGR)."

Study Area
The ICT is the capital city of Islamic Republic of Pakistan ( Figure 1) located in the Potohar plateau. It comprises an area of 906 km 2 including mountains and uneven plains exceeding 1,175 m in height above the mean sea level [63]. According to the 2017 national population census, the total population of ICT is approximately two million, which makes it the ninth largest city of Pakistan. It has a humid and sub-tropical climate with four distinct seasons: autumn, spring, summer, and winter. The temperatures vary from 13 • C in January to 38 • C in June. ICT consists of five planning zones: Zone I, II, and V are reserved for planned urban development, while the remaining two zones, III and IV, are managed as National Park and the rural fringes. Marble and chemical factories, steel mills, flour mills, oil units, pigments, paints, and pharmaceutical manufacturing plants are some of the main industries in the city [34].

Study Objectives
In this paper, well established, proven, and articulated research methodology, satellite datasets, and definitions of features were adopted with the goal to systematically achieve the following defined objectives:  Detection, measurement, and characterization of land cover features using Landsat medium resolution freely available satellite data (1976, 1990, 2000, 2010, and 2016) and determine the annual rate of change in land cover classes at 10 years interval.  Landscape metrics and forest fragmentation spatiotemporal analysis, to estimate and report changes in the ICT urban ecosystem over forty years.  Calculation of SDG indicator number 11.3.1 "Land Consumption Rate to the Population Growth Rate (LCRPGR)."

Study Area
The ICT is the capital city of Islamic Republic of Pakistan ( Figure 1) located in the Potohar plateau. It comprises an area of 906 km 2 including mountains and uneven plains exceeding 1,175 m in height above the mean sea level [63]. According to the 2017 national population census, the total population of ICT is approximately two million, which makes it the ninth largest city of Pakistan. It has a humid and sub-tropical climate with four distinct seasons: autumn, spring, summer, and winter. The temperatures vary from 13 °C in January to 38 °C in June. ICT consists of five planning zones: Zone I, II, and V are reserved for planned urban development, while the remaining two zones, III and IV, are managed as National Park and the rural fringes. Marble and chemical factories, steel mills, flour mills, oil units, pigments, paints, and pharmaceutical manufacturing plants are some of the main industries in the city [34].

Datasets
To study the land cover dynamics and spatial analysis of ICT, various datasets have been collected from primary and secondary sources. The data collected from the primary sources include Landsat medium resolution satellite imagery and field observations using Geographical Positioning System (GPS) receiver, while the secondary or ancillary data consist of population census and ground truth data from Very High Resolution Satellite (VHRS) imagery.

Landsat Satellite and Digital Elevation Model (DEM) Data
For land cover mapping, we used 30 m spatial resolution orthorectified and cloud free Landsat Multispectral Scanner System (MSS), Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), and Operational Land Image (OLI) sensor images. The TM, ETM+, and OLI sensors have the same 30 m spatial resolution, while the MSS has a spatial resolution of 57 m. The overall ICT region is covered within one Landsat scene (185 × 185 km). A total of five images for years 1976,1990,2000,2010, and 2016 were downloaded from the United States Geological Survey (USGS) Earth Resources Observation and Science (EROS) archive. All satellite data was collected for the months of October and November, to avoid cloud cover (Table 1). A 30 m spatial resolution Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) from USGS was used to understand the topography of the area and to prepare the study area maps.

Population Census Data
For the calculation of SDG Goal 11 indicator 11.3.1 LCRPGR, we used the data collected from secondary sources including demographic data (primary census abstracts for the year 1972, 1998, and 2017) from the Pakistan Bureau of Statistics. We used actual census data of 1972 and 1998 to determine the shape of the population curve, which was exponential, as expected. Using the fitted exponential regression function, we derived the growth rate through the exponential function derivative-Equation (1), and then Equation (2) was utilized to determine the approximate population for the year 1976. The population data used for the calculation of LCRPGR is given in Table 2.
Growth rate = (ab)e ac (1) where a and b are regression coefficients, c represents the population variable, POP 1976 = Estimated Population in 1976, and POP 1972 = Past Population in 1972.

Ground Truth Data
In this study, we collected ground truth data from two sources: the primary ground data was collected using a Global Positioning System (GPS) receiver in October 2016 during field survey at various locations of ICT, while for the secondary ground data we relied on sub-meter VHRS imagery from Google Earth. A total of 205 samples were gathered, which were utilized as reference data for the satellite image classification and accuracy assessment of results. Using the Google Earth temporal VHRS data, we visually identified and geo-tagged the sites witnessing land degradation, urban development, and tree loss. The geo-tagged dataset helped to validate Landsat-based land cover changes.

Methodology
Our methodology ( Figure 2) on land cover assessment and spatial analysis of ICT consists of (i) pre-processing of Landsat images, (ii) classification of Landsat images, (iii) accuracy assessment of land cover products, (iv) temporal land cover assessment, (v) urban landscape metrics analysis, (vi) urban forest fragmentation analysis, and (vii) LCRPGR calculation.

Ground Truth Data
In this study, we collected ground truth data from two sources: the primary ground data was collected using a Global Positioning System (GPS) receiver in October 2016 during field survey at various locations of ICT, while for the secondary ground data we relied on sub-meter VHRS imagery from Google Earth. A total of 205 samples were gathered, which were utilized as reference data for the satellite image classification and accuracy assessment of results. Using the Google Earth temporal VHRS data, we visually identified and geo-tagged the sites witnessing land degradation, urban development, and tree loss. The geo-tagged dataset helped to validate Landsat-based land cover changes.

Methodology
Our methodology (Figure 2) on land cover assessment and spatial analysis of ICT consists of (i) pre-processing of Landsat images, (ii) classification of Landsat images, (iii) accuracy assessment of land cover products, (iv) temporal land cover assessment, (v) urban landscape metrics analysis, (vi) urban forest fragmentation analysis, and (vii) LCRPGR calculation.

Pre-Preprocessing Landsat Images
Pre-processing steps for Landsat images included Top of Atmosphere (ToA) reflectance conversion, terrain illumination correction, and layer stacking of spectral bands. A reduction in between-scenes variability was accomplished through normalization for solar irradiance with a two-step process. First, we converted all pixels' Digital Numbers (DNs) to radiance values using the bias and gain values, which were scene-specific and given in the metadata file of the respective scene. Second, we converted radiance data to ToA reflectance [64]. Terrain illumination correction

Pre-Preprocessing Landsat Images
Pre-processing steps for Landsat images included Top of Atmosphere (ToA) reflectance conversion, terrain illumination correction, and layer stacking of spectral bands. A reduction in between-scenes variability was accomplished through normalization for solar irradiance with a two-step process. First, we converted all pixels' Digital Numbers (DNs) to radiance values using the bias and gain values, which were scene-specific and given in the metadata file of the respective scene. Second, we converted Land 2020, 9, 123 8 of 23 radiance data to ToA reflectance [64]. Terrain illumination correction was conducted using an empirical rotation model proposed by Tan et al. [65]. Layer stacking was performed to construct multi-spectral images at 30 m spatial resolution. The false color composite (FCC) of bands 542 from the multi-spectral images was used to visually interpret land cover features.

Classification of Landsat Satellite Images
The image classification procedure is to instinctively categorize all pixels in an image into land cover classes [66]. In this study, the maximum likelihood classification method was used for land cover-mapping from Landsat images. The maximum likelihood classification method is a supervised classification technique, which works on the basis of multivariate normal probability density function of categories, and utilizes both the variance and co-variance of the spectral response of each unknown pixel to assign it to a particular category [67,68]. In this study, five land cover classes were specified: Tree cover >40% canopy, Tree cover <40% canopy, Settlement, Soil, and Water. Almost 35 training samples were collected for each land cover class to classify the various Landsat images. The classified objects with an area smaller than the Minimum Mapping Unit (MMU) (i.e., 1 ha~3 × 3 pixels) were fused with the neighboring land cover classes [69]. We adopted on-scene digitization technique for land cover change detection. First, we overlaid 2016 base land cover map on the multi-spectral Landsat image of 2010. We traced patches where land changes had occurred, leaving unchanged patches unmodified for consistency [70]. We followed the similar approach to detect land cover change between 2000 and 2010, using 2010 land cover map as a base layer. This process was repeated to generate land cover change maps between 1976-1990, 1990-2000, 2000-2010, and 2010-2016.

Accuracy Assessment of Land Cover Products
Assessment of classification accuracy of 1976,1990,2000,2010, and 2016 land cover maps was carried out to determine the quality of information derived from the data. Random sampling method was adopted to assess the accuracy of satellite derived land cover products. Accuracy was assessed using 125 reference points, based on ground truth data and satellite visual interpretation. The comparison of classification results and reference data was carried out statistically using error matrices. In addition, the non-parametric kappa statistic was computed for each classified map to measure the accuracy of results, as it not only accounts for diagonal elements but also for all elements in the confusion matrix [71].

Temporal Land Cover Assessment
In this study, within the 906 km 2 total area of ICT, the land cover maps (1976, 1990, 2000, 2010, and 2016) were compared in terms of the area. According to Puyravaud [72], the annual rate of change is based on the compound interest law, considering non-linear change across the timeline to estimate the percentage change per year. In this study, the annual rate of change was calculated using Equation (3) proposed by Puyravaud [72].
where r is the annual rate of change in percentage, and A 1 and A 2 is the area at earlier time t 1 and later time t 2 , respectively, and ln denotes the natural logarithm function.

Urban Landscape Metrics Analysis
Landscape metrics are used in this study to quantify the spatial patterns of land cover categories. The landscape metrics can be defined as quantitative and comprehensive measurements showing spatial diversity at a specific scale and resolution [73,74]. In this study, landscape metrics analysis was performed at land cover class level to quantitatively analyze spatial structures and patterns of the topographically and biophysically heterogeneous and diverse landscape of the ICT [74].
Under the three habitat categories Patch Density and Size, Shape and Edge, and Proximity/Isolation, a total of nine landscape metrics were used to quantify structural changes:  Table 3). The landscape metrics were calculated using the FRAGSTATS v4 software tool. The output of the landscape metric analysis depends upon the spatial resolution of the data [75]; in this particular study, 30 m spatial resolution data has been chosen from 1990 to 2015 and 60 m spatial resolution for 1976. Measures total length of edge per unit area. It explains complexity of patch shape.

Mean Shape Index (MSI)
Measures complexity of patch shape compared to a standard shape of the same size. Its value increases with complexity of shape.
Mean Perimeter to Area Ratio (MPAR) Patch shape complexity, based on perimeter length to patch area. It explains shape complexity without standardization to a standard Euclidean shape (square).

Proximity/Isolation
Mean Proximity Index (MPI) Measure of connectedness of a habitat class. Considers size and proximity of all patches with the same habitat type inside a specified search radius.

Mean Euclidean Nearest
Neighbor Distance (MED)

Meter (m)
Measures minimum edge-to-edge distance to the nearest neighboring patch of the same type. It explains connectedness or isolation in landscape or habitat class.

Urban Forest Fragmentation Analysis
Forest fragmentation is the splitting up of large contiguous forest fields into smaller or less contiguous areas. A number of events or activities can lead to forest fragmentation including road formations, woodcutting, forest conversion to agriculture, forest fires, and human conflict over forest patches [25]. To assess forest fragmentation in ICT, the land cover map is divided into two major categories: forest and non-forest. Forest class consists of tree cover greater and less than 40% canopy cover while non-forest class comprises of settlement, soil, and water land cover classes. The outcome of forest fragmentation analysis was represented into six categories: Patch, Edge, Perforated, Core (<250 acres), Core (250-500 acres), and Core (>500 acres). These categories are signs of forest ecosystem quality and can be used to estimate the amount of fragmentation present in a landscape and the potential habitat impacts [25,76] temporal and spatial scales [12,77]. If the LCRPGR ratio value lies between 0 ≤ LCRPGR ≤ 1, it shows the simultaneous increase of population growth rate (PGR) and land consumption rate (LCR), but the land consumption rate is much slower than the population growth rate. On the other hand, if LCRPGR > 1, it reflects the simultaneous increase of PGR and LCR, with a faster LCR than the PGR. To estimate the LCRPGR, satellite data derived land cover maps were used for the years 1976, 1998/2000 and 2016/2017, and census data was used for the specific years ( Table 2). The indicator 11.3.1 was assessed at the local level for ICT using the mathematical expressions currently proposed by UN-Habitat, given below in Equations (4)

Accuracy Assessment of Land Cover Products
According to the confusion matrix analysis, 90% overall accuracy and kappa coefficient value of 0. 85 (Table A1 in Appendix A).

Temporal Land Cover Assessment
The overall settlement in the ICT increased while tree cover classes (greater and less than 40% tree canopies) decreased in forty years (1976 to 2016). The estimated land cover area for each image is summarized in Table 4 and spatial distributions are presented in Figure 3.   (Table 4).
Based on the temporal land cover data in the ICT, urban landscape settlement increased at an annual rate of 4.34% since 1976, with the highest annual rate of 8.79% during 2000-2010. The tree cover >40% canopy decreased at an annual rate of 0.81% between 1976 to 2016 and tree cover >40% canopy declined at an annual rate of 0.77% in forty years (1976-2016). The tree cover >40% canopy witnessed the highest loss rate at 2.38% per annum between 2000-2010, and tree cover <40% canopy experienced tree cover loss at 2.72% during 2010-2016 (Table 5).  (Table 4).
Based on the temporal land cover data in the ICT, urban landscape settlement increased at an annual rate of 4.34% since 1976, with the highest annual rate of 8.79% during 2000-2010. The tree cover >40% canopy decreased at an annual rate of 0.81% between 1976 to 2016 and tree cover >40% canopy declined at an annual rate of 0.77% in forty years . The tree cover >40% canopy witnessed the highest loss rate at 2.38% per annum between 2000-2010, and tree cover <40% canopy experienced tree cover loss at 2.72% during 2010-2016 (Table 5).

Urban Landscape Matrix Analysis
Landscape matrix analysis was performed at class level, and the results for each class are reported in Figure 4 and Table A2 in Appendix A. The parameters Number of patches (NP), Largest Patch Index

Urban Landscape Matrix Analysis
Landscape matrix analysis was performed at class level, and the results for each class are reported in Figure 4 and Table A2 in Appendix A. The parameters Number of patches (NP), Largest Patch Index (LPI), Edge Density (ED), and Mean Proximity Index (MPI) for the settlement class overall increased from 1976 to 2016. Largest Patch Index (LPI) increased to 16% in 1990 from 13% in 1976, and then reduced to less than 10% in 2000, 2010, and 2016. For tree cover greater than and less than 40% canopy, Number of patches (NP) and Edge Density (ED) initially increased from 1976 to 2000, and then declined from 2000-2016. For tree cover greater than and less than 40% canopy, soil, and settlement classes, the Mean Patch Size (MPS) and Mean Euclidean Nearest Neighbor Distance (MED) sharply declined from 1976 to 1990 and remained lower. In the Mean Euclidean Nearest Neighbor Distance (MED), tree cover greater than and less than 40% canopy, soil and settlement classes remained less than 250 m throughout from 1976 to 2016. From 1976 to 2010, the Mean Radius of Gyration (MRG) decreased for tree cover greater than and less than 40% canopy, and then slightly increased during the years 2010 to 2016. For tree cover >40%, the canopy Mean Proximity Index (MPI) slightly increased between 1976 to 2016 while for tree cover >40%, the canopy MPI decreased from 1976 to 1990, grew between 1990 to 2010, and then again sharply reduced from 2010 to 2016. In all the land cover classes, the Mean Shape Index (MSI) value remained greater than 1 with slight fluctuations. In the Mean Perimeter to Area Ratio (MPAR), the settlement class increased from 640

Urban Forest Fragmentation Analysis
The forest fragmentation analysis results are presented in Figures 5 and 6, and Table A3 in Appendix A. The analysis shows that the overall core forests of >500 acres decreased from 391. Appendix A. The analysis shows that the overall core forests of >500 acres decreased from 391.98 km 2 (65.41%) in 1976 to 241.44 km 2 (40.29%) in 2016. In the forty years' time span, the patch class increased to 20 km 2 (4.54%) in 2016 from 14.74 km 2 (2.46%) in 1976. The perforated forest fragmentation class was 35.21 km 2 (5.88%) in 1976, which increased 41.79 km 2 (7.23%) in 2000, and then again declined to 22.74 km 2 (5.18%) in 2016. For the edge class, a minor increase was observed from 1976 to 2000 and then a drop in 2016.

Ratio of Land Consumption Rate to the Population Growth Rate (LCRPGR)
For the years 1976-2000, the LCRPGR ratio (SDG indicator 11.3.1) was obtained to be 0.62, which highlights the simultaneous increase of Population Growth Rate (PGR) and Land Consumption Rate (LCR), but the LCR is much slower than the PGR. For the years 2000-2016, the LCRPGR ratio was obtained to be 1.36, which reflects the simultaneous increase of PGR and LCR, with a faster LCR than PGR.

Discussion
The present study was conducted in order to quantitatively analyze the landscape change dynamics in the ICT from 1976 to 2016 using medium spatial resolution Landsat images. In this study, we observed an increase in settlements over the last forty years. The dramatic land cover change in settlements is exerting severe pressure on other land cover classes, particularly tree cover and soil. Already existing urban areas can be seen expanding through rapid construction of residential blocks in the form of housing societies, industrial blocks and road expansions, leading to horizontal and vertical developments in the city and development of lavish farm houses in the vicinity areas (Figure 7). This rapid increase in urbanization is linked also with migrations. Urban growth may have positive or negative impacts on the environment but unplanned growth of urban areas has negative effects. For the economic development of the country, necessary planning is required to make urbanization helpful, as social, health, and environmental issues often accompany the process of urbanization.
In the time span of forty years, an overall decrease has been observed in the area of tree cover classes (i.e., greater than and less than 40% tree canopy) in the ICT, and most of the tree loss has occurred after the year 2000 with a corresponding increase in built-up areas. Apart from tree cutting, another important reason behind rapid tree loss or degradation is forest fires. In the Margalla Hills, forest fires usually take place during dry hot climate conditions when there is no rain for months and temperature goes up to 45 • C. According to Khalid and Ahmad [78], a total of 320 forest fires were recorded from 2002 to 2012 and approximately 8 km 2 area got burnt as a result. In the ICT, due to uncontrolled urbanization and lack of awareness, huge tree loss has been observed in the last sixteen years, i.e., 2000-2016. Another factor for forest degradation is uncontrolled grazing of livestock [78]. As such, there are no adequate plans or manageable methods to stop grazing activities. For ecotourism and public awareness drives, a number of jogging and hiking trails have been formed in Margalla Hills National Park, and a large number of visitors has severely affected the Margalla Hills National Park by dumping waste. These illegal and unmonitored activities in the forested area cause threats to the forest ecosystem [50,79].
Massive migrations have occurred over the past few years from rural to urban areas, mostly due to low cultivated land output, landlessness, sub-division of land, poor economy, and better educational and health opportunities in urban areas. The rapid increase in population has contributed towards natural resource depletion and rapid deforestation close to settlements [80].
The LCRPGR parameter is an indicator of urban sustainable development, whether urban expansion is in balance with population growth or not. According to literature review, limited scientific peer reviewed studies have been reported on the monitoring and mapping of LCRPGR. Under each SDG, a number of targets and indicators have been defined, which countries have to quantify, but most of the developing countries do not have comprehensive databases through which they can compute, quantify, and report the SDGs indicators. To the best of the authors' knowledge, only Nicolau et al. [77] and Wang et al. [81] have computed and reported scientific results of SDG proposed LCRPGR over urban areas. Nicolau et al. [77] based their study over the mainland of Portugal, while Wang et al. [81] carried out their study over mainland China, using earth observation and population census data, and reported the increase of LCRPGR value from 1.69 in 1990-2000 to 1.78 in 2000-2010. The LCRPGR related research findings from these studies show that in most cities, both horizontal and vertical urban expansions are carried out in an unplanned manner which has already effected the equilibrium of land consumption versus population increase to attain effective development goals by 2030 [77,81]. In this study over the ICT, the LCRPGR ratio was 0.62 from 1976 to 2000, which increased to 1.36 from 2000 to 2016. Based on studies conducted on the global scale, in the most of the cases LCR is higher than or equal to PGR due to high demand of luxurious occupancies in the urban areas [26,30]. Based on studies conducted on the global scale, in the most of the cases LCR is higher than or equal to PGR due to high demand of luxurious occupancies in the urban areas [26,30]. The Shenzhen SEZ and ICT cities were both developed in late 1970s. In terms of urban forest cover change, Shenzhen SEZ had been restored to ~85% (1973-2005) [29] while in this study, we detected urban forest loss of ~27%  in the ICT. In both cities (SEZ and ICT), urban forest fragmentation results revealed losses in forest patches. In South Asia, the landscape of the ICT The Shenzhen SEZ and ICT cities were both developed in late 1970s. In terms of urban forest cover change, Shenzhen SEZ had been restored to~85%   [29] while in this study, we detected urban forest loss of~27%  in the ICT. In both cities (SEZ and ICT), urban forest fragmentation results revealed losses in forest patches. In South Asia, the landscape of the ICT resembles strongly with the cities of Kathmandu (the capital of Nepal) and Thimphu (the capital of Bhutan), as they are all cities topographically surrounded by the pine trees. As compared to ICT, Kathmandu and Thimphu cities are older but highly migrant resident populated, vastly encroached, and over grazed. Although many studies have investigated the land cover dynamics of Kathmandu and Thimphu cities using temporal satellite data [82][83][84][85], detailed analysis using parameters such as landscape metrics, forest fragmentation, and LCRPGR has not been performed. In this aspect, this study serves as a methodological framework for application and analysis over other similar cities in developing countries.
Due to rapid tree loss and urbanization, the ICT has observed rapid spells of dust storms and soil erosion [86]. Soil erosion, dry temperatures, and consequent dust storms have a negative impact in the form of air pollution [87] and land degradation. The soil in ICT and the surrounding areas is shallow and has a clay composition [63]. The alluvial lands and terraces in the area tend to have low agricultural productivity and in the southern and western parts of the Potohar plateau, the soil is thin and infertile [88]. Streams and ravines cut the loose plain and cause erosion and steep slopes. This land is generally unsuitable for cultivation. However, large patches of deep, fertile soil are found in the depressions and sheltered parts of the plateau and these support small forests and agriculture [89]. Butt et al. [45] described the most significant land degradation issues in the Rawal watershed as soil erosion and loss of soil nutrients. They further showed that, between 1992 to 2012, the majority of land that was previously vegetation, bare soil, or water bodies was converted to agriculture and settlements, suggesting increased pressure on natural resources in the Rawal watershed [45].
The results of landscape analysis in this study reveal that the ICT urban landscape has become more heterogeneous, disproportional and diverse, and tree patches have declined. Alarmingly, core forests of >500 acres have declined almost 15% in forty years. Although at the individual level, the residents of ICT, civil society, and the local government are trying to recover tree cover loss by planting trees, but these initiatives should be continuous and ongoing on a regular basis to monitor the growth of trees without damaging the existing mature standing trees. The temporal forest fragmentation analysis shows that due to tree cover loss, the three categories of core forest fragmentations (i.e., <250 acres, 250-500 acres, and >500 acres) have decreased in forty years . The loss in forest fragmentation negatively influences the habitat and terrestrial biodiversity of ICT [79,90,91]. Based on landscape metrics analysis over the Dhaka metropolitan, a similar Asian developing city, Dewan et al. [31] revealed that cultivated areas and vegetated lands became highly fragmented with increasing anthropogenic disturbances and urban built up category became aggregated and convoluted.
While our analysis of ICT land cover change dynamics in this study is important and unique, there are a few caveats which are worth mentioning. First of all, for the 1976 land cover map, we have relied on approximately 57 m spatial resolution Landsat 3 Multispectral Scanner System (MSS) sensor data which is relatively coarser spatially as compared to Landsat 5 and 8 (i.e., 30 m), which may affect the spatial heterogeneity and accuracy of developed land cover maps. Second, as Landsat is a pioneer Earth Observation (EO) satellite program initiated in 1972, so acquiring remote sensing imagery of the years before that is not possible, and thus we cannot derive land cover maps before the 1960s, when the Islamic Republic of Pakistan's capital shifted from Karachi to Islamabad. Third, from Landsat 30 m medium spatial resolution satellite data, we cannot detect and delineate the boundaries of built-up areas as well as can be done from sub-meter VHRS images available from the 2000s onward. However, the VHRS datasets come with a high cost, especially when acquisitions at multiples times have to be acquired, and may not be therefore feasible for study sites in developing countries. Of course, research in this domain is ongoing with regards to utilization of VHRS data for study of urban areas, using methods like object-based image analysis, machine learning, etc. Fourth, there may be some level of uncertainty for the field measurement data, as there is always a potential for human error especially when ground truth is collected over larger areas. Fifth, in this study, we utilized only temporal optical satellite data, which may cause optical signal saturation in closed canopy forests, and atmospheric effects (i.e., cloud coverage, haze, and smog).