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

: 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 coefﬁcient 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 signiﬁcant 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, ﬁres, 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.


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  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 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.  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.  University of Maryland

Forest Change Footprint Information Extraction
Results from time-series analysis of Landsat images in characterizing global forest extent and change from 2000 through 2020. The data are used for an accuracy assessment. 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.

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.

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 Remote Sens. 2022, 14, 744 6 of 19 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 2 in the study area, with a forest coverage fraction of 7.47%.

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. 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 2 in the study area, with a forest coverage fraction of 7.47%.

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. 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 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. Table 2. Bands/indices used in this study.

Band/Index Name Description Calculate Method Reference
Blue, Green, Red, NIR, SWIR1, SWIR2 Original bands from the TM, ETM+, and OLI.
The SLC-off data were removed, and the harmonization function [26] was used for consistency processing along with ETM+ and OLI data to obtain the annual sequence images of six bands.
Normalized difference moisture index (NDMI) Normalized difference indices generated by TM, ETM+, and OLI sensors.
Normalized difference vegetation index (NDVI) Normalized difference indices generated by TM, ETM+, and OLI sensors.
Normalized difference snow index (NDSI) Normalized difference indices generated by TM, ETM+, and OLI sensors.

Enhanced vegetation index (EVI)
A vegetation index calculated from three bands of TM, ETM+, and OLI sensors. EVI = 2.5 * ((NIR -Red)/(NIR + 6 * Red -7.5 * Blue + 1)) [37] Tasseled cap brightness (TCB) It is derived from spectral data and the tasselled-cap transformation algorithm. The algorithm can compress spectral data into several bands of physical scene features with minimal information loss.

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].

Type Level Description Thresholds
Disturbance Serious Land use type changes, deforestation, and forest fires caused complete changes in the surface; for example, the transition from forests to agricultural land and buildings.

Moderate
Due to different reasons such as selective logging, drought, or pests and diseases, the forest has been severely disturbed.

Light
Local changes in the forest, forest disturbances can be reflected in high-resolution images, such as plenter-thinning in the management process.
The transition from non-forest land use types to forest types is mainly through afforestation, and areas with better climate conditions can also be self-regulated or community succession through forest ecosystems. <−500

Moderate
The opposite process of moderate disturbance indicates the change of forest structure, such as the change from sparse forest to dense forest.

Light
Due to afforestation or self-recovery of forests, the density of forest structure gradually becomes higher, which can be observed in high-resolution images.

−350 ≤magnitude < −200
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].

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).

•
HGFC was produced by the Global Land Analysis and Discovery Laboratory at the University of Maryland, in partnership with the Global Forest Watch, and provides annually updated global-scale forest loss data derived using Landsat time-series imagery (https://storage.googleapis.com/earthenginepartners-hansen/GFC-2020-v1.8 /download.html, accessed on 10 December 2021) [9]. This dataset was generated based on multi-source remote sensing data, such as Landsat and MODIS from 2000 to 2020, combined with bagged decision tree classification methods. In contrast to the Hansen product, we extended the time to 1990 and used a change detection algorithm based on spectral-temporal segmentation, a near-automated change detection algorithm that has the advantage of requiring less input data and having a high-detection efficiency compared to the classification method. We downloaded and synthesized HGFC data from 2000 to 2020 for validation, including the disturbance and recovery bands. • Validation samples were obtained through visual interpretation of high-resolution Google Earth historical images. 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 2 of forest was disturbed, accounting for 28.64% of the total forest area. The average disturbed area was 426.88 km 2 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 km 2 , and reached a peak of 1410.23 km 2 in 1997. After 2009, the disturbed forest area was less than 300 km 2 . 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 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 2 of forest was disturbed, accounting for 28.64% of the total forest area. The average disturbed area was 426.88 km 2 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 km 2 , and reached a peak of 1410.23 km 2 in 1997. After 2009, the disturbed forest area was less than 300 km 2 . 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 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  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 2 of forest was restored, accounting for 29.66% of the total forest area. The average annual recovery area was 442.01 km 2 , 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 km 2 in 1990, and reached the peak of 1464.34 km 2 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. km 2 of forest was restored, accounting for 29.66% of the total forest area. The average annual recovery area was 442.01 km 2 , 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 km 2 in 1990, and reached the peak of 1464.34 km 2 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. 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 2 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. 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 2 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.  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 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  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%). 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.

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.
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.
The high producer and user accuracies of the disturbance and recovery categories indicate that the LandTrendr method is robust in obtaining disturbance and recovery foot- 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.

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. 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. 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.

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%.    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 2 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 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 2 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.

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 km 2 and 40.89 km 2 , 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

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 km 2 and 40.89 km 2 , 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. In the Gilgit Baltistan of Pakistan, the forest disturbance area was 40.89 km 2 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 km 2 , 168.89 km 2 , and 117.26 km 2 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 CO 2 emission/retention trends  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.

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.

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.