Detecting Mountain Forest Dynamics in the Eastern Himalayas

Abstract

Forest dynamics is critical to forested ecosystems, and considerable efforts have been devoted to monitoring long-term forest dynamics with the goals of sustainable management and conservation of forests. However, little attention has been given to mountain forests, which are more challenging to monitor due to complex topography, weather, and their distribution. We developed a 30-m resolution tree-canopy cover (TCC) and forest change dataset for the Eastern Himalayas from 1986 to 2021. The tree-canopy cover estimation was validated against estimates from the space-borne Global Ecosystem Dynamics Investigation (GEDI), demonstrating strong consistency (R-square greater than 0.81). A comprehensive assessment for the forest change dataset was performed using 448 visually interpreted points and reported high accuracy of the dataset, i.e., 97.7% and 95.9% for forest loss and gain, respectively. Higher producer and user accuracies were reported for forest loss (PA = 78.0%, UA = 60.9%) than these for forest gain (PA = 61.7%, UA = 56.7%). The results indicated that (1) the mean tree-canopy cover in the region increased by 2.76% over the past three decades, from 40.67% in 1990 to 43.43% in 2020, suggesting the forests have improved during the period; (2) forest loss was identified for a total area of 6990 km² across the study area, which is less than the 10,700 km² identified as forest gain; (3) stronger forest gains were found at elevations greater than 3000 m asl, indicating faster forest growth in high elevations likely influenced by the warming temperatures in the Eastern Himalayas.

1. Introduction

Forests comprise the largest terrestrial ecosystem [1]. They provide important ecosystem services to society, ranging from encouraging water conservation to [2] providing habitats for species [3], improving landscape function [4], regulating climate [5], and maintaining the balance of the global ecosystem via carbon sequestration [6,7]. Forests have undergone great changes in recent decades, including both forest loss and recovery [8]. The Global Forest Resources Assessments (FRA) 2020 report stated that there was a decrease of more than 80 million hectares of primary forest globally since 1990, accounting for 7.2% of the total primary forest [1,9]. Forest loss has become a substantial and complex threat to Earth, leading to negative environmental impacts, including biodiversity loss [10,11], the disruption of water cycles [12,13], increasing soil erosion [14], the destruction of local livelihoods [15], and the increase in carbon dioxide emissions to the atmosphere. The general perception is that global forest change occurs mainly in lowland areas with dense populations [16]. This statement aligns with many previous studies indicating most forest change occurs at low elevations due to human activities, with only negligible forest loss in the mountains [17,18,19]. However, recent studies have reported that mountain forests are also undergoing significant changes due to climate change and agricultural expansion [16,20,21]. Compared to other forest ecosystems, references and studies regarding mountain forests are scarce, and changes in mountain forests remain unclear. Mountain forests account for approximately 28% of the world’s forests and provide many substantial ecological, economic, hydrological, and social values [22]. They are biodiversity warehouses with 25 out of the 34 global biodiversity hotspots, and the per-unit biodiversity of mountain forests is typically higher than that of adjacent lowland forests [23]. In mountain areas, the temperature increase was larger than the global average over the past decades, and extreme climate events were more frequent in mountain regions than in other ecological systems [24,25,26,27,28]. Due to the complex and variable climate conditions, mountain forests may be more vulnerable than lowland forests to environmental changes. Therefore, it is imperative to gain a comprehensive understanding of the status and changes of mountain forests over the long term.

In recent decades, the rapid development of remote sensing technologies and the emergence of new satellite sensors and analytical methods have contributed to long-term forest monitoring [29]. In addition, remote sensing can provide periodic and seamless observations of Earth, fostering the development of large-scale forest monitoring. Numerous studies have reported forest changes at both global and local scales [9,19,30,31,32]. However, they focused primarily on short-term forest loss, neglecting forest recovery as well as long-term monitoring. In addition, many previous forests change studies were based on the examination of vegetation indices (e.g., NDVI). These indices not only represent the status of forests but are also influenced by other vegetation. To improve the accuracy of forest-change monitoring, an approach based on tree-canopy cover (TCC) has been proposed in recent years. Tree-canopy cover represents the fraction of tree canopy coverage within a pixel area as a “continuous field” [33]. The Tree Cover layer in the Moderate-resolution Imaging Spectroradiometer (MODIS) Vegetation Continuous Fields (VCF) is an annual continuous global tree-canopy cover dataset at a 250-m resolution that has been widely used in global- or local-scale assessments [34]. Because many forest changes occur in patches smaller than a 250-m footprint, one need for forest change monitoring is repeated observations at finer resolution such as those from Landsat [29]. Sexton et al. [35,36] proposed a method to produce a tree-canopy cover and forest change dataset based on long-term Landsat observations. Feng et al. [37] assessed the dataset via visual interpretation, and the accuracy is confirmed. While the tree-canopy-cover-based method has been proposed, there are few studies focusing on the long-term changes of mountain forests, where complex topography and weather present challenges to forest-change monitoring.

The Himalayas experienced a rapid change in climate, which was estimated to be far greater than the global average, along with the significant increase in extreme climate events in recent decades [26,27,38,39,40,41,42]. In this study, we choose the Eastern Himalayas as an example to investigate its forest status and change. The Eastern Himalayas could be an excellent test area for monitoring mountain forest dynamics because it is both the highest mountain in the world and the most sensitive to temperature and climate variations [43,44]. This study aims to derive Landsat-based tree-canopy cover in the Eastern Himalayas from 1986 to 2021 and monitor forest loss and forest gain. In addition, we provide a comprehensive assessment for both TCC and forest change datasets in the mountain area. Monitoring forest dynamics can provide long-term spatiotemporal consistent information that is valuable for understanding the change of mountain forests and ensuring their sustainability.

2. Study Area and Datasets

2.1. Study Area

The Eastern Himalayas (longitude: 88°22′–95°26′E, latitude: 26°41′–29°54′N) is located at the edge of the Tibetan Plateau, extending from the Yadong River in the west to Brahmaputra in the east (Figure 1). The Eastern Himalayas are a part of the Alpine–Himalayan mountain chain with high elevations [45]. The total area covers 156,699 km², and elevation ranges from 42 m asl to 7527 m asl. According to data from the Climatic Research Unit (CRU) in 2020, the monthly average temperature within the study area ranged from −3.2 °C to 13.5 °C and monthly precipitation ranged from 1 mm to 340 mm.

The complex topography and the steep south-to-north elevation gradient have resulted in many diverse physiographic landscapes and regional climates [46]. At high elevations, temperatures are very low, even remaining below 0 °C year-round. Precipitation decreases from the southern foothills to the northern foothills across the Eastern Himalayas. The northern portion of the Eastern Himalayas is characterized by less precipitation and a drier climate because high-elevation mountains block water vapor from the south. On the other hand, the southern portion of the Eastern Himalayas is predominantly influenced by the Indian monsoon [26]. The Eastern Himalayas are some of the most ecologically sensitive areas in the world, with many vegetation zones, including, from southeast to northwest, montane forests, shrubs, alpine meadows, alpine steppes, and alpine deserts [41]. Spruce is the dominant forest species widely distributed in the Eastern Himalayas [41].

2.2. Datasets

2.2.1. Landsat Data

Landsat satellite data constitute the longest record of Earth observation with a relatively high spatial resolution of approximately 30 m [47]. Therefore, Landsat data have been extensively applied to investigating how forests have been changing for decades. Landsat data are publicly available on the United States Geological Survey (USGS) website (https://earthexplorer.usgs.gov/, (accessed on 9 April 2022)). We downloaded Landsat-5 Thematic Mapper (TM), Landsat-7 Enhanced Thematic Mapper Plus (ETM+), and Landsat-8 Operational Land Imager (OLI) images during the growing season of 1986 to 2021. The study area consists of 15 tiles of the Landsat World Reference System 2 (WRS-2), including 1983 images in total.

2.2.2. MODIS VCF

The MODIS Vegetation Continuous Fields (VCF) MOD44B Collection 6 data product is a sub-pixel-level representation of surface vegetation cover. The MOD44B product produces at a 250-m spatial resolution from monthly composites of the MODIS surface reflectance data from 2000 [34]. The MOD44B dataset includes three cover layers, i.e., percent tree cover layer, percent non-tree cover, and percent non-vegetated layer. In this study, the MODIS VCF tree cover data were selected to provide tree cover training data. The MOD44B dataset was downloaded from the USGS website (https://earthexplorer.usgs.gov/, (accessed on 9 April 2022)).

2.2.3. GEDI

The Global Ecosystem Dynamics Investigation (GEDI) measurements observe the Earth’s surface between a longitude of 51.6°N and 51.6°S both day and night [48], covering temperate and tropical forests. The returned lidar waveform measured by the GEDI laser can be used to derive products such as the canopy cover fraction, canopy top height, and topographic surface elevation. The GEDI dataset is publicly available on the EARTHDATA website (https://earthdata.nasa.gov/, (accessed on 17 April 2022)). We used the GEDI L2B canopy cover data from 2019 and 2020 to validate our Landsat-derived tree-canopy cover.

2.2.4. Ancillary Data

The Advanced Land Observing Satellite (ALOS) World 3D-30 m (AW3D30) is a 30-m resolution global digital surface model (DSM) dataset produced by the Japan Aerospace Exploration Agency (JAXA). The AW3D30 is produced using data archived from the Remote-sensing Instrument for Stereo Mapping (PRISM) onboard the ALOS satellite [49]. We used this DSM to analyze tree-canopy cover and forest change at different elevations. The dataset was downloaded from (https://www.eorc.jaxa.jp/ALOS/, (accessed on 3 May 2022)). We obtained the temperature and precipitation data from the Climate Research Unit (CRU) (http://www.cru.uea.ac.uk/, (accessed on 11 May 2022)) at a resolution of 0.5° × 0.5°. In addition, we used high-resolution satellite imagery from Google Earth Pro and Normalized Difference Vegetation Index (NDVI) time-series data derived from Landsat to aid visual interpretation. The NDVI profile was derived from Landsat NDVI using the Google Earth Engine platform.

3. Methodology

3.1. Generation of Tree-Canopy Cover

TCC can be estimated using a piecewise linear function of surface reflectance and temperature parameters as was proposed by Sexton et al. in 2013 [35]:

$$C_{i,t}=f\left(X_{i,t}\right)+\epsilon$$

(1)

where $X$ is a vector of temperature and surface reflectance, $\epsilon$ is the error produced by $f$ applied to $X$ in the estimates, the subscript $i$ represents a pixel’s location, and the subscript $t$ represents the time indexed by year.

The TCC regression tree model was fitted using CatBoost [50], an implementation of gradient boosting, which is a type of machine learning boosting technique that relies on an ensemble of sequentially built models to minimize the overall prediction error [51]. The MODIS MOD44B tree cover estimates were used as training data. Landsat surface reflectance data were aggregated to 250 m using a spatial mean estimator. In this way, a joint sample of tree cover and surface reflectance variables was prepared to create a training dataset for the Eastern Himalayas each year. Then, the fitted model was then used in the original Landsat data with a resolution of approximately 30 m to estimate TCC at a finer scale from 1986 to 2021. In this study, we only analyzed forest dynamics in the study area from 1990 to 2020 due to a large amount of missing data in earlier years.

3.2. Detection of Forest Change

Forest change detection was based on the bi-temporal class-probabilities model proposed by Sexton et al. in 2015 [36]. The model was used to detect forest change in each pixel using TCC and its corresponding error as the input raster. While the model is applicable to land cover and any bi-temporal biophysical attribute (e.g., biomass, tree-canopy height, diversities of land covers), it is most pertinent to detecting forest change using remotely sensed characteristics. The detection of forest change was performed in two steps.

Step one: Define “forest” as a kind of land cover category wherein TCC exceeds a predefined threshold. The probability of a pixel belonging to “forest” is calculated as the function:

$$p\left(F\right)\stackrel{\scriptscriptstyle\mathrm{def}}{=}p\left(c>c^{*}\right)={{\displaystyle \int }}_{c^{*}}^{100}p\left(c\right)dc$$

(2)

where $p\left(F\right)$ is the probability of a pixel being identified as “forest”, $p\left(c\right)$ is the probability density function of TCC, $c$ is the value of TCC, and $c^{*}$ is the predefined threshold of TCC. The probability of a pixel belonging to “non-forest” is then calculated as 1 − $p\left(F\right)$.

Step two: Calculate the probability of forest dynamics. Based on the probability of forest in a pixel for each of two times, each pixel was assigned one of four forest dynamics categories: Stable forest (FF), stable non-forest (NN), forest loss (FN), and forest gain (NF). These four probabilities were calculated as follows:

$$p{\left(FF\right)}_i=p\left(F_{i,t_1},F_{i,t_2}\right)=p\left(F_{i,t_1}\right)\times p\left(F_{i,t_2}\right)$$

(3)

$$p{\left(NN\right)}_i=p\left(N_{i,t_1},N_{i,t_2}\right)=\left(1-p\left(F_{i,t_1}\right)\right)\times \left(1-p\left(F_{i,t_2}\right)\right)$$

(4)

$$p{\left(FN\right)}_i=p\left(F_{i,t_1},N_{i,t_2}\right)=p\left(F_{i,t_1}\right)\times \left(1-p\left(F_{i,t_2}\right)\right)$$

(5)

$$p{\left(NF\right)}_i=p\left(N_{i,t_1},F_{i,t_2}\right)=\left(1-p\left(F_{i,t_1}\right)\right)\times p\left(F_{i,t_2}\right)$$

(6)

where $p{\left(FF\right)}_i$, $p{\left(NN\right)}_i$, $p{\left(FN\right)}_i$, and $p{\left(NF\right)}_i$ represent the four forest dynamics from time $t_1$ to time $t_2$. Each pixel was assigned to the most probable class.

To better interpret the different long-term forest changes, the following four categories were identified: Forest loss without gain, forest gain without loss (forest newborn), forest gain after loss (forest recovery), and forest loss after gain.

3.3. Assessment of Tree-Canopy Cover

The GEDI lidar altimeters have been collecting data on forest structure since 2019. In this study, the GEDI L2B canopy cover data, including 2769 observation points in 2019 and 3731 observation points in 2020 (Figure 2), were used to validate Landsat-derived tree-canopy cover for the same years. The selected data points were distributed both in areas with dense forests and areas with sparse trees. For a point from GEDI, if it meets one of the following criteria, we set canopy cover to zero: (a) The point was classified as snow/ice, (b) the point was located in water, or (c) the canopy height of the point was 0 m.

3.4. Assessment of Forest Change

3.4.1. Sampling Design

A stratified sampling design was developed to assess the accuracy of the forest change results. The changed area in the Eastern Himalayas was much smaller than the unchanged area in the region according to our forest change results. To increase the sampling efficiency and improve the representation of forest change samples, we used a preliminary version of the forest change layer to divide the Eastern Himalayas into two strata (i.e., changed and unchanged).

3.4.2. Point Selection

A set of 225 points was designed to be randomly collected for each of the changed and unchanged strata. The sampling probabilities of strata were calculated as follows:

$$p\left(i|K\right)=\frac{n_k}{N_k}$$

(7)

where $p\left(i|K\right)$ represents the probability of a pixel being sampled, $n_k$ is the number of desired pixels, and $N_k$ is the total number of pixels in this stratum. $n_k$ was set to 225 for each stratum. Each pixel was given a random number $p^{*}$, and pixels that met the criteria $p^{*}$ < $p\left(i|K\right)$ were selected. A total of 448 sample points were selected across the Eastern Himalayas. The sampling probability for the changed stratum was 0.00126%, and that for the unchanged stratum was 0.00017%.

3.4.3. Visual Interpretation

Forest dynamics were visually interpreted by experts by examining the historical high-resolution satellite imagery in Google Earth Pro. To improve its reliability, phenology based on Landsat’s NDVI was also used. The NDVI profile of each sample point was derived from Landsat data using the Google Earth Engine. For example, a sample point (93.8182°E, 27.7992°N) is located in the eastern part of the Eastern Himalayas (Figure 3). High-resolution satellite imagery from Google Earth Pro shows that large areas of trees were cleared in 2014. There was a concomitant sharp decline in NDVI (Figure 4). The high-resolution imagery from 2019 showed that trees had recovered after several years, and NDVI showed a visible recovery. Therefore, it is reasonable to assume that this place experienced forest loss in 2014 with a forest recovery process that started after the loss event.

3.4.4. Validation Metrics

The labeled reference points were used to assess the accuracy of forest changes by calculating the confusion matrix and the accuracy metrics, i.e., overall accuracy (OA), producer’s accuracy (PA), and user’s accuracy (UA) [37,52]. Weights were applied to points to remove the bias of the disproportional sampling between changed and unchanged pixels. Each point’s weight for the changed ($w_i$), and unchanged ($w_j$) strata were calculated as follows:

$$w_i=\frac{p\left(j|K_j\right)}{p\left(i|K_i\right)+p\left(j|K_j\right)}$$

(8)

$$w_j=\frac{p\left(i|K_i\right)}{p\left(i|K_i\right)+p\left(j|K_j\right)}$$

(9)

where $p\left(i|K_i\right)$ and $p\left(j|K_j\right)$ represent the inclusion probability of a pixel being sampled in the changed and unchanged stratum, respectively. In this example, the point’s weight for the changed stratum was 0.1163, while that for the unchanged stratum was 0.8837.

Accuracy assessment of forest change was calculated using a confusion matrix [53]. For each validation point ($i$), the agreement between the estimate (y_i) and reference ($y$) was defined as follows:

$$y_i=\left\{\begin{array}{c}1 if {\widehat{c}}_i=c_i\\ 0 if {\widehat{c}}_i\ne c_i\end{array}\right.$$

(10)

The OA indicates the probability that a pixel was classified correctly, and was defined as the ratio of diagonal elements of the confusion matrix to the total weighted number of sample points, $n_a$:

$$OA={{\displaystyle \sum }}_{i=1}^{n_a}{y}_i\times w_i / {{\displaystyle \sum }}_{i=1}^{n_a}{w}_i$$

(11)

The PA and UA were calculated as follows:

$$P{A}_C={{\displaystyle \sum }}_{i=1}^{n_c}{y}_i\times w_i / {{\displaystyle \sum }}_{i=1}^{n_c}{w}_i$$

(12)

$$U{A}_C={{\displaystyle \sum }}_{i=1}^{n_{\widehat{c}}}{y}_i\times w_i / {{\displaystyle \sum }}_{i=1}^{n_{\widehat{c}}}{w}_i$$

(13)

where $n_c$ represents points labeled as type c (i.e., forest loss or gain) by the reference, and $n_{\widehat{c}}$ represents points identified as type c.

4. Results

4.1. Spatial Distribution of Tree-Canopy Cover

The tree canopy is densest in the southeast portion of the Eastern Himalayas, thinning from the southeast to the northwest (Figure 5). Trees are mainly distributed on the southern slopes of the Eastern Himalayas. Areas with zero tree-canopy cover accounted for 36% of the whole study area, mainly found in the northern slopes of the Eastern Himalayas.

The areas associated with different tree-canopy cover show a bimodal pattern with peaks at approximately 0% to 5% and 85% to 96% (Figure 6a). Tree canopy was found primarily above sea level to below 5000 m asl, with the tree-canopy cover reaching the densest values in elevations between 500 m asl and 3000 m asl and becoming nearly absent in areas above 5000 m asl (Figure 6b). Tree-canopy cover also increases gradually from west to east with increasing longitude; the densest distribution is found between 93°E and 95.4°E (Figure 6c). South of latitude 29.2°N, tree-canopy cover in the Eastern Himalayas decreases, with the densest distribution between 26.6°N and 27.6°N. In contrast, tree canopy increases in density as latitude increases above 29.2°N (Figure 6d).

4.2. Spatio-Temporal Changes of Tree-Canopy Cover

Forests in the Eastern Himalayas increased in density from 1990 to 2020 (Figure 7). The mean tree-canopy cover in the region increased by 2.76% over 31 years, from 40.67% in 1990 to 43.43% in 2020, with an annual increase of 0.11% tree-canopy cover per year (i.e., percent cover = 0.11x × year − 175.46, R² = 0.65, p < 0.001). Trees covered 56,650 km² of the study area in 1990 and 60,460 km² in 2020, increasing by 3810 km² over the 31-year period.

Tree canopy gains occurred more frequently than tree canopy losses in the Eastern Himalayas from 1990 to 2020, with the strongest increases occurring in the high elevations (i.e., roughly 3000~4500 m asl) (Figure 8). In contrast, tree canopy losses in the study area were smaller in magnitude and occurred primarily at lower elevations where human activities are more intense.

4.3. Validation of Tree-Canopy Cover

A comparison of our Landsat-derived TCC values with GEDI-derived values displays strong agreement, with R-squared values of 0.87 and 0.81 for GEDI observations in 2019 and 2020, respectively (Figure 9a,b), indicating Landsat-derived tree-canopy cover estimates are reliable in the Eastern Himalayas.

4.4. Forest Loss and Gain

The area of forest gain in the Eastern Himalayas was much greater than that of forest loss from 1990 to 2020 (Figure 10). In addition, forest gain exceeded forest loss in every elevation zone except below 500 m asl. More than 10,700 km² in the Eastern Himalayas demonstrated forest gain, with an average increase of 350 km² per year over the 31-year period. Approximately 6990 km² of forest were lost in the Eastern Himalayas in the past three decades, with an average loss of 230 km² per year over the 31-year period. Forest gain and forest loss mainly occurred between elevations of 3000 m asl and 4500 m asl. Barely any forest gain and loss occurred at elevations greater than 5000 m asl since trees were nearly absent there.

The spatial distribution of forest changes was analyzed across the entire study area. Significant forest losses were concentrated in the lowlands of the Eastern Himalayas (Figure 11a), which have been subjected to deforestation for decades. In contrast, forest gains were mainly found at higher elevations. The mountain forests exhibited a net gain in each longitude and latitude zone of the Eastern Himalayas (Figure 11b,c). Over the 31-year period, more than 3870 km² of forest were recovered after loss. Approximately 1420 km² of the forest experienced loss without recovery detected, and the area was the smallest among the four categories and mainly distributed across lower elevations (Figure 12). Forest gain without loss detected was identified in a 5135 km² area, which has the largest area among the four categories and exceeded other categories in every elevation zone above 1000 m asl.

4.5. Validation of Forest Change

We performed an accuracy assessment of the forest change dataset using a total of 448 visually identified sample points (Figure 13). The results of the accuracy are shown in

Table 1. Accuracy of forest change detection.

Accuracy Metrics	Forest Loss (%)		Forest Gain (%)
OA	97.7		95.9
	Loss	No loss	Gain	No gain
PA	78.0	98.3	61.7	97.6
UA	60.9	99.3	56.7	98.1

. Forest loss had higher accuracy than forest gain with an OA of 97.6%, PA of 78.0%, and UA of 60.9%. The OA of forest gain was 95.9%, and the PA and UA were 61.7% and 56.7, respectively. There are four kinds of forest dynamics categories, including change status (i.e., loss, gain) and unchanged status (i.e., no loss, no gain). The producer and user accuracies of the unchanged pixels were much larger than that of the changed pixels. It was likely that many unchanged sample points were located on bare lands at high elevations. The overall accuracy of forest change was high (>95%) in this study, indicating the reliability of Landsat scenes in identifying forest changes.

5. Discussion

5.1. Mountain Forest Changes

This analysis reveals that the mountain forests of the Eastern Himalayas exhibited a net gain of forest over the three decades. This is likely the result of human activities and climate change [19]. The United Nations Framework Convention on Climate Change defines forests as having a minimum tree cover greater than 10–30% [54]. Using this definition, the Eastern Himalayas contained between 69,385 km² and 76,083 km² of forest in 1990 and increased to between 71,176 km² and 76,620 km² of forest in 2020. Forest gains were distributed in the southern slopes of the Eastern Himalayas, and there were no gains in the northern slopes due to the area’s high elevation and associated dry, cold climate. More specifically, the forest gains were mainly distributed at high elevations greater than 3000 m asl, which agrees with the results from Zheng et al. [26]. Forest loss occurred mainly at lower elevations, likely due to human activities being more intense at lower elevations. Similar patterns were also reported in other regions worldwide [17,18,55]. For example, Curran et al. [18] point out the lowland forests in Indonesian Borneo decreased by 56% from 1985 to 2001. The study from Aide et al. [17] in the Andes showed that forest loss mainly accrued at elevations below 1500 m asl due to agriculture expansion, while above 1500 m asl, it was dominated by forest gain. In addition, our findings were consistent with the results from Song et al. [19] who researched global land change from 1982 to 2016. Song et al. [19] pointed out that tree cover in mountain systems of Earth has increased from 1982 to 2016. Curtis et al., estimated the drivers of forest loss [56]. Although the data are of a coarser resolution and only span 2000 to 2015, they can still be considered a reference for understanding the causes of forest loss in the Eastern Himalayas (Figure 14), and forestry (i.e., large-scale clearcutting within managed tree plantations and forests) was the primary cause of forest loss in the Eastern Himalayas, accounting for 60% of the forest loss. Commodity-driven deforestation leading to permanent land use change was the secondary driver of forest loss in the Eastern Himalayas, accounting for 28%. The felling of trees for paper pulp, timber, and plantation expansion is the major commercial activity that adversely affected the forests. Forest loss for the remaining areas was attributed to wildfires (7%) and shifting agriculture (5%).

5.2. Forest Dynamics Monitoring

The long-term Landsat-derived tree-canopy cover provided a consistent dataset for capturing complex mountain forest dynamics at a high resolution and large scale. Compared to our estimates of forest loss, the Global Forest Change product produced by Hansen et al. [9] seriously underestimated forest loss in the Eastern Himalayas. Our findings on the underestimation issue agreed with those reported by Milodowski et al. [57] in the Amazon. From 2001 to 2020, more than 5630 km² of the Eastern Himalayas experienced forest loss, while only 1440 km² of co-occurring forest loss were detected by the Hansen Global Forest Change product. This difference is likely due to the fact that the Hansen Global Forest Change product only considered forest loss under conditions of a stand-replacing disturbance or complete canopy removal [9], while we were able to include weaker disturbances such as selected logging and forest degradation when calculating forest loss. For example, this sample point (91.7642°E, 26.9048°N) is located in the southern part of the Eastern Himalayas (Figure 15). Satellite imagery from Google Earth Pro shows that large areas of trees were disturbed due to the construction of roads. Forest losses were relatively accurately identified by our study, while much forest loss information in this area was lost in the Global Forest Change product (Figure 16).

5.3. Corelation between TCC and Climate Factors

Since climate change directly influenced the phonology of vegetation [42,58,59,60], we investigated the correlation between Landsat-derived tree-canopy cover and key climate variables, e.g., temperature and precipitation. According to the CRU dataset, the Eastern Himalayas experienced rapid warming with an increase of 0.44 °C over the past three decades, from 9.47 °C in 1990 to 9.91 °C in 2020, with an annual increase of 0.025 °C per year (i.e., temperature = 0.025x × year + 9.4575, r = 0.53, p = 0.002) (Figure 17a). This result aligned with previous studies that the Himalayas experienced significant warming during the recent decades [38,42,61]. The annual precipitation in the study area decreased by 94 mm with significant fluctuation, from 1678 mm in 1990 to 1584 mm in 2020 (Figure 17b). The fluctuation of precipitation agrees with the data reported by Shrestha et al. [62], which pointed out that precipitation from 78 station records showed significant fluctuation both on decadal and annual scales. The annual mean tree-canopy cover from 1990 to 2020 showed a significant positive correlation with the annual mean temperature in the Eastern Himalayas (r = 0.47, p = 0.008), and no statistical correlation with annual precipitation (r = −0.10, p = 0.593) (Figure 17c,d). A similar conclusion was made by Li et al. [43] in the Eastern Himalayas. Zheng et al. [26] also pointed out that tree growth in the Eastern Himalayas, especially the southeastern Tibetan Plateau, showed an overall increase in growth with a rapid increase after 1982, which was likely due to continuous warming. In addition, our findings also aligned with previous studies reported by Mohammat et al. [63] in Inner Asia, which demonstrated that the increase in NDVI was more strongly correlated with temperature than with precipitation in high-altitude regions. Therefore, the temperature was likely the dominant driver of the increase in tree-canopy cover in the Eastern Himalayas.

5.4. Validation and Uncertainties

Both Landsat-derived tree-canopy cover estimates and forest change data in the Eastern Himalayas were assessed. The Landsat-derived tree-canopy cover estimates were found to be strongly consistent with GEDI observations in 2019 and 2020, indicating that the tree-canopy cover estimates are reliable in the Eastern Himalayas. The lower agreement at lower- and higher tree-canopy cover was likely due to the fact that GEDI’s data were somewhat overestimated at the higher end of the tree cover range and underestimated at the lower end of the tree cover range [64]. In addition, the capture of within-canopy gaps in GEDI’s measurements [65] could also contribute to the differences in the Landsat-derived tree-canopy cover estimates. Our validation results indicate that forest change derived from tree-canopy cover is able to capture forest dynamics (loss and gain). However, forest gain was found to be less accurate than forest loss, which agrees with results published by Hansen et al. [9]. This finding is likely due to the long and slow process of canopy recovery and the high uncertainty associated with visual identification of forest gain, which was also reported by Feng et al. [37].

6. Conclusions

In this study, we validated Landsat-derived tree-canopy cover, as well as forest loss and forest gain in the Eastern Himalayas. Even though the method and accuracy of forest detection had been improved in recent years, forest change detection remains a challenge, especially for mountain forests. Understanding the accuracy and uncertainties is critical for researchers and governments to use this data. In particular, the space-based laser altimeters, such as ICESat-2 and GEDI, provide the ability to more accurately detect forest structure for forest change analysis. Given the high consistency between our Landsat-derived tree-canopy cover estimates and GEDI-based estimates, our methods show promise for use in mountain forest monitoring, as well as other forest systems, and further research should be conducted to ensure that this technique provides accurate monitoring for forest change in more areas.

During the past few decades, with the vast number of newly available satellite and airborne sensors [66], more remote sensing data can be integrated and fused with Landsat data to improve the spatial and temporal scale and consistency of forest change detection. Additional data from these sensors would dramatically increase our ability to monitor forest structure and change dynamics. Tree-canopy cover could be a key variable for incorporating these multi-source satellite data and promoting robust forest change detection and real-time monitoring.