*Article* **The Forest Change Footprint of the Upper Indus Valley, from 1990 to 2020**

**Xinrong Yan 1,2 and Juanle Wang 1,2,3,4,\***


**Abstract:** The upper Indus Valley is the most important and vulnerable water tower in the South Asian subcontinent, which provides a vital water supply for 230 million people in the basin. Forests play an important role in water conservation in this region, and the security of upstream forests forms the foundation downstream water and food security. However, a big challenge is to effectively monitor the dynamics of the forest in this region. Thus, we used the LandTrendr spectral-temporal segmentation algorithm combined with 8203 scenes of multi-source remote sensing data to study the forest change footprint in the upper Indus Valley. The overall accuracy of LandTrendr extraction for forest disturbance and recovery was 86.01%, and the Kappa coefficient was 0.73. The results showed the following: (1) From 1990 to 2020, the area of forest recovery was 1.01% more than that of disturbance, 70% of disturbance occurred between 1990 and 2001, and 60% of recovery occurred between 1999 and 2012. (2) Although the overall trend of forest disturbance and recovery was balanced, there were significant differences in forest management status among the different regions. Nepal has the highest forest stability, India has the largest area of forest disturbance, and Pakistan and China have the largest areas of forest recovery. (3) India's Himachal Pradesh and Jammu and Kashmir are the two provinces with the largest disturbed areas, primarily due to grazing, fires, and commercial tree planting. Pakistan's North-West Frontier, Azad Kashmir, and China's Tibet Ali region were major contributors to the recovery, which was driven by afforestation policies in both countries. This study provides an important data base and monitoring method for planning land and forest use in Indus Valley countries, protecting fragile environments, and promoting policies for the Sustainable Development Goals.

**Keywords:** forest disturbance; forest recovery; footprint information; LandTrendr spectral-temporal segmentation algorithm; upper Indus Valley

#### **1. Introduction**

Forests are the main component of terrestrial ecosystem and the largest "carbon pool" on land, which plays an important role in regulating climate and mitigating global warming [1]. Monitoring forest disturbance and recovery has received abundant attention over the past few decades, especially to identify the important role of forests in curbing climate warming and achieving "carbon neutrality" [2]. The United Nations Sustainable Development Goal (SDG) 15 calls for the protection, restoration, and promotion of the sustainable use of terrestrial ecosystems and sustainable management of forests as well as a series of assessment indicators [3]. Forest disturbance is the main factor that affects forest growth, structure, and function. This disturbance is largely due to the degradation and disappearance of forests caused by environmental changes such as drought, hurricanes, geological disasters, and various human activities [4,5]. Forest recovery involves a series

**Citation:** Yan, X.; Wang, J. The Forest Change Footprint of the Upper Indus Valley, from 1990 to 2020. *Remote Sens.* **2022**, *14*, 744. https://doi.org/ 10.3390/rs14030744

Academic Editor: Wenxin Zhang

Received: 20 December 2021 Accepted: 3 February 2022 Published: 5 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

of favorable conditions that can promote the growth of forests, including internal and external factors. In areas with more favorable climatic conditions, forest disturbances can be recovered through succession of forest communities or internal self-regulation of ecosystems [6]. In high latitude and ecologically fragile areas, the implementation of forest protection policies plays a positive role in forest recovery [7,8]. Using the Landsat pixel scale, forest disturbance was defined as partial or complete disappearance of the forest canopy. Forest recovery is defined as the inverse ratio of forest loss, which represents an increase in forest cover of the tree canopy [9].

The upper Indus Valley is one of the most critical water towers in the world [10], and ecosystem security in this region is directly linked to the welfare of 230 million individuals downstream. The Karakoram, Hindu Kush, Ladakh, and Himalayas in the upper reaches form diverse forest ecosystems, which are important carbon sinks and water conservation areas. The security of forest ecosystems is related to the security of water and agricultural resources in the Indus Valley basin. However, forest disturbances in this region have been a concern in recent decades. Studies have shown that the annual deforestation rate in this area reaches 2.2% [11]. A previous study by Rashid et al. (2017) showed that the evergreen broadleaved forest in Kashmir was degraded, and aboveground biomass and carbon storage were lost, thereby adversely affecting carbon fixation in this region [12]. There is evidence that forests in the upper Indus Valley were destroyed, and severe forest disturbances devastated the stability of the original forest ecosystem, leading to catastrophic floods in 1992 and 2010. In response to these problems, the governments of relevant countries have issued several commercial logging bans, but this has not prevented the loss of upstream forest areas. The rate of deforestation and loss of total forest area in this region are still higher than the global average [11,13]. Thus, it is urgent to obtain information on the disturbance and recovery of forests in recent decades through various monitoring methods to provide a reference for the formulation of forest policies in this region.

Forest disturbance and recovery information can be obtained through many technical means such as regular field surveys of forest resources. China, Canada, and Finland have developed relatively complete forest resource survey techniques, but these require higher human and economic costs. The upper Indus Valley includes five countries, and it is difficult to organize forest surveys and obtain forest disturbance and recovery data for all countries. Remote sensing technology is currently considered the most effective and lowest cost method, providing spatiotemporal awareness of forest change on multiple scales [14,15]. The principle of remote sensing technology for monitoring disturbance and recovery is to track the change information of the vegetation canopy spectrum in a time-series image. When the band or index representing the vegetation canopy suddenly or gradually exceeds a certain threshold, the forest is considered disturbed or recovered [16]. Time-series images provide important, consistent, and continuous data sources for the acquisition of forest footprint information in forest-covered areas. Early forest change monitoring based on remote sensing mainly came from the comparative analysis of the interpretation results of two or more temporal remote sensing images. However, because of the limitation of the timeliness of interpretation of samples and other reasons, such methods generally span many years, such as 5 or 10 years, ignoring the short-interval footprint information in the process of forest change, and it is difficult to achieve continuous monitoring of forest change information from a large area.

In the upper Indus Valley, Joshi et al. (2014) interpreted images of the central Himalaya in 1979, 1999, and 2009 by means of artificial visual interpretation, in which forest stands are reduced and forest areas are affected by degradation and isolation [17]. To obtain information on forest loss in the Himalayas, Qamer et al. (2016) used remote sensing images during two time periods from 1990 to 2000 and 2000 to 2010 [18]. Owing to the limitation of computing power or samples, these methods did not describe the forest loss footprint information in detail, and the images of the three periods could not dynamically restore the process of forest change over 21 years. However, most of these were short-term and partial studies, making it difficult to create an overall assessment of the forest status of the entire upper Indus Valley. With the accumulation of remote sensing images and the enhancement of data computing capability, forest change monitoring has gradually developed into a time-series image combined with a change detection algorithm. Since the release of Earth big data computing platforms, such as Google Earth Engine (GEE), increased data computing power has enabled forest change monitoring to use higher spatial resolution data. [19]. In recent years, the big data cloud computing platform for Earth observation has increasingly integrated interference and recovery monitoring algorithms. This includes the vegetation change tracker [20], breaks for additive seasonal and trend [21], vegetation regeneration and disturbance estimates over time [22], continuous change detection and classification [23], and continuous monitoring of land disturbance [24]. These algorithms are widely used in forest event monitoring, such as forest fires, deforestation, and degradation, and are important for recognizing long-term changes in forests.

Based on this background, this research intends to combine Landsat time-series imagery and Earth observation big data cloud computing platform GEE to carry out forest disturbance and recovery monitoring research in the upper Indus Valley. The main objectives include the following: (1) combining the time-spectral segmentation algorithm LandTrendr to obtain the forest change footprint in the upper Indus Valley from 1990 to 2020; and (2) analysis of the results of forest disturbance and recovery to obtain the spatiotemporal cognition of the forest change footprint in the upper Indus Valley.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The Indus River is a major river in Central Asia and an important source of agricultural irrigation that feeds more than 230 million people. The Indus Valley, ranging from 24◦ N to 37◦ N and from 66◦ E to 82◦ E (Figure 1), has a diverse geographical environment, covering an area of 618,580.35 km2. It borders the Karakoram Mountains and Himalayas to the northeast, the Thar Desert in India to the southeast, the Hindu Kush Mountains in Afghanistan to the northwest, the Baluchistan Plateau to the southwest, and the Arabian Gulf to the south. The Indus Valley has a subtropical climate, characterized by a distinct monsoon climate. However, owing to the influence of the high mountains in the northeast, the climate is usually between dry and semi-dry, tropical, and subtropical. The year is divided into four seasons.

The upper Indus Valley is the intersection of many mountain ranges, and the complex geographical environment has led to the development of a variety of forest ecosystems. The natural vegetation mainly includes coniferous forests (subalpine forests, arid temperate forests, humid temperate forests, and subtropical pine forests), shrub forests (arid subtropical broad-leaved forests and arid tropical thorn forests), economic forests, and shelter forests. The forests are widely distributed between 500 and 5500 m elevation, and 90% of the forests are located above 1500 m.

According to a report from the Worldwide Fund for Nature, the population boom coupled with poverty and a lack of awareness have led to illegal and unsustainable logging, excessive logging for fuel and charcoal, and increased small-scale agriculture, which continue to reduce the amount of forest cover in the upper Indus Valley. In addition, forest fires, natural disasters, climate change, pests, and diseases further contribute to the degradation and disappearance of forests [25].

**Figure 1.** The study area is in the upper Valley of the Indus River. To highlight the forests, we used a Sentinel-2 image synthesized with false colors (the darker red represents the lusher vegetation).

#### *2.2. Data Preprocessing*

Table 1 shows the sources and descriptions of multi-source remote sensing data used in this study. Landsat time-series data were used from the United States Geological Survey (USGS), including thematic mapper (TM), enhanced thematic mapper (ETM), and operational land imager (OLI) sensors. The data have the advantages of having high resolution, a short revisit cycle, and easy access and processing. GEE was used as the main data processing platform, which provided access to the atmospheric correction surface reflectance data of all Landsat sensors. We connected all available remote sensing images of vegetation growth season on the GEE platform from July 15 to October 1, which can effectively avoid images caused by seasonal snow in high-altitude mountainous areas. Landsat 7 ETM+ scan line corrector failure (SLC-off) images with data gaps were removed. A total of 8203 remote sensing images were used in this study. In order to reduce the influence of the difference in reflection wavelength between ETM+ and OLI sensors on the results, the harmonic function from OLI to ETM+ was used to process the data consistently [26]. In the high mountains of South Asia, clouds, snow, ice, and mountain shadows have a serious impact on the results of forest change monitoring. We used the image quality assessment band (QA) and the CFMask algorithm to generate masks to remove these four types of elements, and the annual image was generated by the median synthesis method [27]. According to the statistics obtained from annual synthetic images, cloud-free observation has been achieved in 90% of the forest for more than 20 years, and the average cloud-free observation time for each pixel was 25.5 years from 1990 to 2010. The annual synthetic images meet the requirements for LandTrendr algorithm fitting.


**Table 1.** Datasets used in this study.

#### *2.3. Forest Change Footprint Information Extraction*

Figure 2 shows the workflow for obtaining forest change footprint information. This process mainly includes data collection and processing, mask construction, algorithm segmentation, different levels of disturbance and recovery classification, and accuracy assessment.

**Figure 2.** Flowchart of forest change footprint monitoring in the upper Indus Valley.

2.3.1. Built Mask for Forest/Non-Forest Areas

The forest mask can effectively avoid the influence of cultivated land and dense grassland on LandTrendr segmentation results. In this study, the forest area needed to be a maximum union of the forest distribution from 1990 to 2020, which cannot be achieved by using remote sensing images or datasets of a single time phase.

We obtained all available remote sensing images for July 1 to October 1 from 1990 to 2020, and based on the forest extract rule, a normalized difference vegetation index (NDVI) > 0.6 and short-waved infrared radiation (SWIR) < 0.1 generated the annual forest distribution area. Zhu et al. (2012) showed that the NDVI of the forest growth season is always above 0.6, but sometimes other vegetation types, such as crops, grass, and shrubs, could also have an overall forest NDVI value above 0.6. Since forests are generally dark in the SWIR band compared to other vegetation types, a surface reflectance threshold of 0.1 in the overall band 7 excludes other vegetation types that may have high NDVI values [14,28]. The maximum forest union from 1990 to 2020 showed that there was a total forest area of 46,192.25 km<sup>2</sup> in the study area, with a forest coverage fraction of 7.47%.

#### 2.3.2. Detection Methods for Forest Disturbance and Recovery

LandTrendr is a set of spectral-temporal segmentation algorithms based on remote sensing image pixels that are useful for change detection in a time series of moderateresolution satellite imagery. The time segmentation algorithm is considered an effective method for detecting forest disturbance and recovery. The trajectory of the generated spectral time-series data had almost no interannual signal noise [29]. The algorithm uses a time segmentation strategy based on regression and a point-to-point fitting spectral index as a time function, allowing for the capture of slow-evolving processes such as recovery and unexpected events [30]. Interactive Data Language (IDL) initially implemented LandTrendr, and later Google engineers ported LandTrendr to the GEE platform [19,31].

The GEE framework nearly eliminates the onerous data management and image preprocessing aspects of the IDL implementation. LandTrendr combined with GEE also simplifies tedious data management and image preprocessing by directly accessing geospatial datasets. Figure 3 shows the process of the LandTrendr algorithm for extracting forest disturbances. The discrete original value of the spectral index or band was divided into a series of straight lines and breakpoints. From the segmentation results, the start year, end year, duration, and magnitude of the spectral index or band of forest disturbance can be easily obtained.

**Figure 3.** LandTrendr pixel time-series segmentation. Image data were reduced to a single band or spectral index and then divided into a series of straight-line segments by breakpoint (red points) identification. For example, this figure shows the segmentation process of the NDVI, AB, and CD, which represent the steady state of the forest, where BC is fitted as a disturbance event, and CD is an inverse process of BC and is identified as a recovery process.

The LandTrendr algorithm requires a series of surface reflectance bands and spectral indices as inputs when performing spectral-time segmentation. The first band or index was used for the main segmentation, and the other bands were used for fitting to supplement the best results. Cohen et al. (2018) showed that using the 13 bands and indices in Table 2 can effectively reduce the error rate of the LandTrendr segmentation [32]. Table 2 shows the bands and indices used in this study. After obtaining the segmentation results from LandTrendr, we filtered the segmentation results using a forest mask to remove the effects of dense grassland and cropland.



2.3.3. Classification of Forest Disturbance and Recovery Levels

Three levels of disturbance and recovery were classified to further study the extent of forest disturbance and recovery (Table 3). In the segmentation results of the LandTrendr algorithm, the "magnitude" band records the variation amplitude of disturbance and recovery in different years. Different levels of disturbance and recovery were obtained by reclassifying this band. Referring to previous study [41], we obtained the thresholds of disturbance and recovery using the visual interaction method between classification results and high-resolution remote sensing images from Google Earth.


**Table 3.** Different levels of disturbance and recovery [41].

Serious disturbance and strong recovery categories indicate distinct processes of forest change, such as clear-cutting or barren afforestation. Moderate disturbance, moderate recovery, light disturbance, and light recovery may be changes in forest growth trends, such as varying levels of forest degradation or the growth process from young forests to mature forests [41].

#### 2.3.4. Accuracy Assessment

In this study, we validated the forest disturbance and recovery results using two data sources: HGFC datasets and Google Earth images. The validation data were classified into four categories: disturbance, recovery, and both disturbance and recovery (disturbance + recovery).


#### **3. Results**

#### *3.1. Spatial and Temporal Patterns of Forest Disturbance and Recovery*

Figure 4 shows the spatial distribution of forest disturbances (including serious, moderate, and light disturbances) in the upper Indus Valley from 1990 to 2020 using transitional colors. Forest disturbances are widely distributed in the middle of the study area, mainly

in the southwestern Himalayas and the southern Hindu Kush. Large areas of forest disturbance occurred in the upper reaches of the three rivers of Jhelum, Chenab, and Ravi and on both sides of the Kashmir Valley, largely distributed around cities and farmlands, between 1995 and 2015. In the past 31 years, a total of 13,233.55 km<sup>2</sup> of forest was disturbed, accounting for 28.64% of the total forest area. The average disturbed area was 426.88 km2 per year, indicating that 0.92% of forests were disturbed at different levels per year on average. The annual disturbed area showed a fluctuating downward trend. In 1990, the disturbed forest area was 1080.17 km2, and reached a peak of 1410.23 km2 in 1997. After 2009, the disturbed forest area was less than 300 km2. The disturbance area from 1990 to 2001 accounted for 70% of the total disturbed area in 31 years, indicating that the disturbance mostly occurred in the first 10 years of the study period.

**Figure 4.** Spatial and temporal footprint distribution map of forest disturbance. Four selected subareas in the figure represent disturbance result examples in different areas, including (**a**) plain area, (**b**) plain and alpine transition area, (**c**) alpine area, and (**d**) river valley area.

Figure 5 shows the spatiotemporal information of forest recovery (including strong, moderate, and light recovery) in the upper Indus Valley from 1990 to 2020, and the recovery area is widely distributed in the middle region as well as the disturbed area. Similar to the distribution of disturbances, the restored areas were also concentrated in the southwestern Himalayas and the southern Hindu Kush. In the past 31 years, a total of 13,702.55 km<sup>2</sup> of forest was restored, accounting for 29.66% of the total forest area. The average annual recovery area was 442.01 km2, and 0.96% of the forests were restored on average. The annual recovery area showed a hump distribution trend. At the beginning, the recovery forest area was 924 km2 in 1990, and reached the peak of 1464.34 km2 in 2000. The area recovered from 1999 to 2010 accounted for 54.78% of the total area of 31 years, indicating that the forest recovery was mainly concentrated in the middle part of the study period.

**Figure 5.** Spatial and temporal footprint distribution map of forest recovery. Four selected sub-areas in the figure represent recovery result examples in different areas, including (**a**) plain area, (**b**) plain and alpine transition area, (**c**) alpine area, and (**d**) river valley area.

From the results of the long-term analysis, it can be concluded that forest recovery and disturbance areas tend to be balanced in the upper Indus Valley. The recovery area was only 469 km<sup>2</sup> larger than the disturbance area, accounting for 1.01% of the forest area in the upper Indus Valley. Figure 6 shows the interannual variation characteristics of disturbance and recovery. In 1990–2000, 70% of forest disturbance occurred, followed by 60% forest recovery in 1999–2012. In terms of spatial distribution, 70% of disturbances in 1990–2000 and 60% of recovery in 1999–2012 were mostly spatially coincident. This indicates that the disturbed forests from 1990 to 2000 were restored in the following 10 to 13 years, whether natural or man-made restorations. For the entire upper Indus Valley, 2012 is an important time node. The cumulative disturbance area since 1990 is equal to the cumulative recovery area this year, reaching a point of equilibrium. After 2012, the cumulative recovery area was slightly greater than the cumulative disturbance area.

#### *3.2. Temporal and Spatial Characteristics of Different Levels of Disturbance and Recovery*

Figure 7 shows the spatiotemporal information distribution of 16 types of disturbance and recovery processes, including three different levels of combinations of disturbance and recovery. The results showed that 17.93% of the forests in the study area underwent two processes of disturbance and recovery from 1990 to 2020, especially in the areas bordering forests and farmland in the southern Himalayas and the sides of the Kashmir Valley.

**Figure 6.** Interannual variation characteristics of disturbance and recovery. The red cross represents the cumulative area of the disturbed region equal to the cumulative area of the recovery region in the upper Indus Valley.

**Figure 7.** Spatial distribution of forest disturbance and recovery combinations with different levels. Four selected sub-areas in the figure represent different levels of disturbance and recovery examples in different areas, including (**a**) plain area, (**b**) plain and alpine transition area, (**c**) alpine area, and (**d**) river valley area.

Figure 8 shows the proportion of different combinations in the total forest area. A total of 59.71% of the forests remained stable in 31 years without disturbance or recovery. In the period from 1990 to 2020, forests that only underwent disturbances accounted for 10.72%, of which 86.51% were light, 11.62% were moderate, and only 1.86% were serious disturbances. The forests that only underwent recovery accounted for 11.64% of the total forest area, which was 0.92% more than the disturbance area. The results of different recovery level were as follows: light recovery (91.74%) > moderate recovery (7.63%) > strong recovery (0.61%).

**Figure 8.** Proportion of disturbance and recovery of different levels in forest area.

The area where both disturbance and recovery occurred accounted for 17.93% of the total forest area, indicating that this part of the forest underwent a process of disturbance and recovery conversion. The areas of light disturbance and light recovery were the largest, accounting for 56.26%, followed by moderate disturbance + light recovery (14.42%), light disturbance + moderate recovery (10.38%), and moderate disturbance + moderate recovery (9.4%). The proportions of the other five combinations were all less than 3%. The results showed that forest disturbance and recovery in the study area were mainly light disturbances, and serious disturbances and strong recovery accounted for very few (<2%).

LandTrendr not only obtained spatiotemporal information of severely disturbed and obviously restored forest areas, but also captured the natural growth trend of forests or forest degradation caused by climate change and drought. Although there is no effective way to classify these areas in detail, obtaining varying degrees of disturbance and recovery is crucial for understanding the processes of forest change.

#### *3.3. Accuracy Assessment of LandTrendr Results*

The accuracy evaluation of the HGFC data and Google Earth images on LandTrendr segmentation results shows that the algorithm can effectively monitor forest disturbance and recover the footprint (Table 4). The different classes demonstrated high producer and user accuracies.


**Table 4.** The accuracy assessment result was based on LandTrendr segmentation results and HGFC datasets.

Analysis of the HGFC data showed that the highest user accuracy disturbance class was 91.69%, whereas the user accuracy for identifying both disturbance and recovery was 65.16%. There was a 26.15% confusion between the disturbance + recovery and disturbance classes, indicating that only one disturbance process was identified in the disturbance + recovery class. There was little difference among the three classes of producer accuracy: the highest accuracy of the recovery category was 87.22%, and the lowest accuracy of the disturbance category was 85.18%. In Google image-based evaluation, the highest user accuracy of 90.05% was achieved in the disturbance class (Table 5). In producer accuracy assessment, the highest accuracy was 91.84% for the disturbance class and the lowest was 61.86% for the recovery class.


**Table 5.** The accuracy assessment result was based on LandTrendr segmentation results and samples from Google Earth images.

The high producer and user accuracies of the disturbance and recovery categories indicate that the LandTrendr method is robust in obtaining disturbance and recovery footprint information.

#### **4. Discussion**

#### *4.1. Forest Disturbance and Recovery in Different Regions*

The results showed that, although the overall trend of forest disturbance and recovery was balanced, there were significant differences in forest disturbance and recovery in different regions due to the geographical environment and management policies. The forests in the upper Indus Valley are managed by five countries. To explore the changes in forests in different countries, we created statistics and analyzed results from the perspective of forest management.

Table 6 shows the area of stability, disturbance, recovery, and both disturbance and recovery have occurred in India, Pakistan, Afghanistan, China, and Nepal. The disturbed and recovered forests of the five countries showed significantly different characteristics. From 1990 to 2020, the forest stability in Nepal was the highest (71.03%), Pakistan and China had little difference (66.97% and 64.59%, respectively), and India had the least stability

(52.98%). Figure 9 shows the proportion of disturbance and recovery areas in the forest area of the upper Indus Valley in different countries. The highest proportion of disturbance was 13.18% in India, and the least were 5.95% and 5.33% in China and Nepal, respectively. The area of forest recovery was highest in China (17.49%), followed by Pakistan (11.81%), Nepal (11.74%), India (11.41%), and Afghanistan (8.81%). In the region, where both disturbance and recovery have occurred, India accounted for the highest proportion at 22.42%, and Nepal for the lowest proportion at 11.89%.


**Table 6.** Forest disturbance and recovery statistics in the upper Indus Valley (unit: km2).

Disturbance + Recovery means both disturbance and recovery have occurred

**Figure 9.** The proportion of forest disturbance and recovery area in the upper Indus Valley.

Figure 10 shows the interannual variation characteristics of forest disturbance and recovery in the upper Indus Valley in different forest management regions. In India, the disturbance forest area has been larger than the recovery forest area for 31 years since 1990. By 2020, the disturbance and recovery forest area showed a trend towards balance, but there was still a 387.86 km<sup>2</sup> gap. From 1990 to 1994, the accumulative recovery area of Afghanistan was larger than the accumulative recovery area. After 1994, the accumulative disturbance area continued to be larger than the accumulative recovery area. Although Afghanistan's forest area in the upper Indus Valley is relatively small, the disturbance and recovery of forests did not reach equilibrium until 2020. There were two inflection points in the forest change trend in Pakistan. In 1993, the cumulative disturbance area in Pakistan exceeded the cumulative recovery area and lasted for 16 years. After reaching equilibrium in 2009, the forest change trend in Pakistan tended to recover. The interannual variation of forest in China and Nepal was similar, and the cumulative recovery area was always larger than the cumulative disturbance area after 1993. Many studies have pointed to forest disturbances in Nepal, but this study only focused on a very small portion of Nepal's forests located in the upper Indus Valley, and the results represent only that portion of forests.

**Figure 10.** Interannual variation characteristics of disturbance and recovery in in the upper Indus Valley. (**a**) India (**b**) Pakistan (**c**) Afghanistan (**d**) China (**e**) Nepal.

#### *4.2. Analysis of the Causes of Disturbance and Recovery*

There is a close relationship between the climate and vegetation in the Indus Valley. In the upper Sindh and Punjab, overgrazing and deforestation have led to the destruction of many natural vegetation types. In addition, humans have long interfered with natural water systems, and in Shivalik, deforestation has led to a marked deterioration of groundwater and vegetation cover. From 1990 to 2020, 59.71% of the forests in the upper Indus Valley remained stable, and the restored area was 0.92% greater than the disturbed area. Although the overall trend of forest change was balanced, there was still an imbalance among the five countries, especially in different provinces.

Table 7 shows the data of forest change in different administrative regions. Himachal Pradesh and Jammu and Kashmir are the two regions with the most disturbance area exceeding the recovery area, covering 362.77 km2 and 40.89 km2, respectively. Pathania et al. (2012) noted that the forest area in this region decreased with the passage of time, and grazing caused an obvious loss of plantation forest [42]. Other studies have shown that the large-scale introduction of horticultural cash crops in the region between 1998 and 2010 led to some significant changes in forest composition, with commercial cultivation leading to the decline of some important plant species [43]. We also observed that some reports indicate that the forest coverage in this area increased from 1991 to 2015, which can

be explained by commercial tree planting [44]. Forest fires, landslides, and other natural disasters are also significant causes of forest disturbance in this region [45], but it is difficult to quantify the specific area of disturbance currently caused by these reasons. These regions should be aware that changes in the composition of forest species and ecological imbalances caused by deforestation and forest degradation may cause irreversible damage to unstable and fragile mountainous areas.


**Table 7.** Forest disturbance and recovery area of 28 administrative regions in the upper Indus Valley (km2).

In the Gilgit Baltistan of Pakistan, the forest disturbance area was 40.89 km2 more than the recovery area, and the disturbance was mainly distributed in the Western border area with Afghanistan. Qamer et al. (2016) showed that the area was felled or severely degraded from 1990 to 2010 [18]. The results of this study and those of Qamer et al. (2016) show consistent temporal and spatial characteristics in the Gilgit Baltistan. The forests of Afghanistan are mainly distributed in the Eastern region, and the disturbance and recovery of the provinces are balanced, except in the Paktia, Paktika, and Nangarhar districts. A previous study on forest change in Northern Pakistan from 1990 to 2010 similarly showed that some assessment reports have grossly overstated deforestation rates, which is consistent with the findings of our study [18].

The main contribution of forest recovery came from the Khyber Pakhtunkhwa, Azad Kashmir, and Tibet, and the recovery areas were 537.05 km2, 168.89 km2, and 117.26 km2 more than the disturbance area, respectively. Forest recovery in the Khyber Pakhtunkhwa Province was primarily driven by the Khyber Pakhtunkhwa Billion Tree Forestation Project, which aims to plan, design, initiate, and implement the Green Growth Initiative in the forestry sector [46]. Khan et al. (2020) studied the temporal changes in forest cover, carbon storage, and corresponding CO2 emission/retention trends (1989–2018) in Azad Kashmir and showed that both forest cover and carbon storage increased significantly, and the research results were consistent with this study [47]. The Tibet Ali region, carried out

from 2002 biological sand prevention engineering, implemented a five-phase afforestation project, which is consistent with the results of this study.

#### *4.3. Method Limitations and Its Application*

The LandTrendr spectral-temporal segmentation algorithm results may vary greatly according to different parameters and spectral indices [48]. The fixed parameters cannot obtain the optimal result, which needs to be further determined by combining the UI program provided by LandTrendr with real ground samples. For forest disturbance and recovery, the information extraction method used in this study can obtain spatiotemporal cognition of the forest change footprint in the upper Indus Valley. However, this approach does not distinguish between the possible causes of disturbance, such as deforestation, climate change, and geological hazards. The rate band in the results provided by LandTrendr provides a possible way to distinguish between these possible causes, which require further analysis.

#### **5. Conclusions**

To describe the disturbance and recovery of forests quantitatively and objectively in the upper Indus Valley, we used multi-source remote sensing data, the LandTrendr spectraltemporal segmentation algorithm, and a remote sensing big data computing platform to complete the monitoring of forest change footprint in the upper Indus Valley. The main conclusions are as follows:

(1) The LandTrendr algorithm combined with multi-source remote sensing data and the GEE big data platform completed forest change footprint tracking of forests in the upper Indus Valley for 31 years, and the algorithm showed stable robustness and portability.

(2) Forest disturbance and recovery were widespread in 1990–2020, most of the disturbance occurred in 1990–2000, the recovery was in 2000–2010, and equilibrium was widely attained in 2012.

(3) Forest disturbance and recovery in different forest management regions showed significant differences, and forest disturbance in India and Afghanistan did not reach equilibrium by 2020. Pakistan reached equilibrium in 2009, while Nepal and China showed relatively stable and continuous trends in forest recovery.

This study can further contribute to more effective forest management policy development by identifying the spatial and temporal patterns of disturbance and recovery for quantitative assessment.

**Author Contributions:** Conceptualization, methodology, and resources, X.Y. and J.W.; software, validation, data curation, and writing —original draft preparation, X.Y.; writing—review and editing, J.W. and X.Y. All authors have read and agreed to the published version of the manuscript.

**Funding:** Special acknowledgement should be expressed to the China–Pakistan Joint Research Center on Earth Sciences and the Construction Project of the China Knowledge Center for Engineering Sciences and Technology (Grant number CKCEST-2020-2-4), which supported the implementation of this study.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

