Forest disturbance decreased in China from 1986 to 2020 despite regional variations

Human activities have altered disturbance patterns in many parts of world, but there is no quantitative information on patterns and trends of forest disturbance regimes in China. We applied a spectral-temporal segmentation approach over all available Landsat data to map individual disturbance patches and characterize the patterns and trends in disturbance rate, size, frequency, and severity across China’s forests. From 1986 to 2020, about 39.7% of China’s forests were disturbed with an annual rate of 1.16 ± 0.41% yr−1. The disturbance decreased at a rate of −390 ± 142 km2 yr−1, primarily driven by the effective implementation of forest protection policy since 2000s. The rate, frequency, and size of disturbance generally intensified in Southeast, but weakened in Northeast China. Our high-quality, spatially explicit disturbance map provides an essential data layer to understand the landscape-scale drivers of forest dynamics and functions for important but less understood pan-temperate forest regions. The rate, frequency and extent of forest disturbances decreased in China over 1986–2020 probably due to forest protection policies, but while they intensified in the southeast, they weakened in the northeast, according to 30 m resolution Landsat satellite data

D isturbances play a critical role in shaping forest structure and dynamics and therefore have profound impacts on ecosystem functions and biodiversity 1,2 . Satellite observations have shown that about 5.7% of global forest (ca. 2.3 million km 2 ) was lost to stand-replacing disturbances from 2000 to 2012 3 , and the drivers of forest loss varied considerately among regions 4 . Meanwhile, low severity, non-stand replacing and diffuse disturbances may be more prevalent, affecting much larger forested areas. Multiple lines of evidence have demonstrated that climate change and human activities increased the extent and magnitude of forest disturbances in many parts of the world 5 . Intensifying disturbance regimes, together with anthropogenicdriven exacerbation of chronic environmental drivers (e.g., rising temperature and water stress), will likely increase tree mortality and force forests toward shorter-statured and younger stands, and thus have cascading consequences for biodiversity and climate mitigation 6 . Therefore, quantification of spatial patterns and trends of disturbance regime is key to understanding the processes underlying the shifts in forest dynamics and functions.
Disturbances are spatially contiguous landscape processes that are best characterized by earth observations at fine to moderate spatial resolutions (e.g., 10-50 m) 1 . Recently, there is a growing interest in using Landsat observations to map long-term changes in forest disturbance regimes over large regions [7][8][9][10] . For example, about 17% of European forests were disturbed during 1986-2016 11 with increasing disturbance rates due to forest operations, including timber harvest 12 , even though Landsatbased detection changes may bias the trends of harvested areas 13,14 and natural disturbances (e.g., fires) 15 . In Russian boreal forests, climate warming has significantly increased fire extent and severity in recent decades 16 . In North America, patterns, trends, and severity of forest disturbances were also well characterized by recent studies [17][18][19] . China ranks the 5th around the world in terms of forested area, together with Europe and North America, comprises most of world's temperate forests. Recent studies showed that disturbances and regrowth in temperate forests contributed to over half of the carbon sink in the world's forests 15,[20][21][22] . To date, there is little quantitative information available on disturbance regimes and their changes over time in China's forests, which represents a significant spatial gap in understanding the underlying drivers for structure and function in temperate forests.
China's forests are heavily disturbed and primarily secondgrowth forests recovering from extensive logging since the early 20th century. Excessive timber extraction and poor management have led to forest loss and degradation, with cascading effects on regional environmental and ecological conditions. Since the late 1970s, the Chinese government has implemented several forestrelegated ecological restoration projects to reverse this trend 23,24 . Notably, the Natural Forest Protection Program (NFPP) covers 27.5% of China's land area launching in 1999 and represents a paradigm shift in forest management policy 25,26 . Recent studies have recognized the contribution of these ecological restoration projects in preserving and enhancing forest covers, quality and ecosystem services, such as carbon sequestration [27][28][29][30] . However, there is no knowledge on how changes in forest restoration policy affected the spatial and temporal patterns of disturbance regimes in China's forests. Recent global-scale disturbance data is often limited by the temporal extent (e.g., post-2000 3,4 ), to a particular disturbance type (e.g., wildfire 31 ), or at spatial resolutions that are often too coarse to capture landscape-scale disturbance dynamics 32 . Therefore, the main objective of this study is to provide a high quality, spatially explicit disturbance dataset to characterize the disturbance regimes in China's forests.
Here, we used a spectral-temporal segmentation method to map and characterize the disturbance regimes of China's forest at 30-m spatial resolution for the period 1986-2020 utilizing all available Landsat time series images (i.e., TM, ETM + and OLI) over the growing season (June-September) (# of Landsat images = 56634, Supplementary Fig. 1) and 31,225 reference plots ( Supplementary Fig. 1). We followed the definition of disturbance by 33 : "any relatively discrete event that disrupts the structure of an ecosystem, community, or population, and changes resource availability or the physical environment". We detected disturbances within forests by tracking the abrupt changes in timeseries spectral indices, which were indicative of removal or change in tree cover associated with forest conversion to other land use types, management activities at various intensities (e.g., logging), or forest fires; and mapped reforestation/afforestation disturbances outside of forests by tracking forest expansion. Our overarching goal is to understand the patterns, trends and drivers of disturbance regimes in China's forests during 1986-2020. Specifically, (1) What are the spatial patterns of disturbance regimes across China's forests? (2) How have the forest disturbance regimes changed during 1986-2020? (3) What are the roles of restoration policy in shaping the disturbance regimes across China's forests? Here we demonstrate that about 39.7% of China's forests were disturbed with an annual rate of 1.16 ± 0.41% yr −1 , and the annual disturbed forested decreased at a rate of −390 ± 142 km 2 yr −1 from 1986 to 2020 primarily driven by the effective implementation of forest protection policy since 2000s. The disturbance regimes generally intensified in Southeast and weakened in Northeast China, reflecting a regional shift in restoration and conservation policy.

Results
Disturbance regimes. Satellite observations showed that the total forested area was 209.8 million ha in 2000, which were consistent with the national forest inventory data, global forest resource assessment, and other estimates from satellite observations [34][35][36] . The forested area was mainly distributed in Northeast (NE, 20%), Southwest (SW, 30%) and Southeast (SE, 37%) China. Our results showed that about 83.3 million ha forest was disturbed in the period 1986-2020, accounting for 39.7% of China's forested area, with an annual rate of 1.16 ± 0.41% (mean ± sd, hereafter) ( Fig. 1). Nationally, forest fires constituted only a small percentage of disturbances (~5%), with regional hotspots (>30%) mostly located in the deciduous needle leaf forests of the Great Xing'an Mountains in NE region, and the mountainous coniferous forests of the southeast Tibet in SW region ( Supplementary Fig. 3a). However, climate warming increased the percentage of forest fires in Great Xing'an Mountains in NE region ( Supplementary  Fig. 3b).
We used the user's and producer's accuracy as the main accuracy assessment metrics. Comparisons between the visually derived reference data and Random Forest (RF) classified disturbance maps showed that the user's and producer's accuracy was 0.76 and 0.96, respectively, with an overall accuracy of 0.96 (CI 95% = (0.955, 0.963), Kappa score = 0.93) (Supplementary Table 1). The temporal accuracy (i.e., year of disturbance) was within 3 years, consistent with other studies 11 . Our spatially explicit disturbance maps also showed significant positive correlations with the Landsat-derived global forest cover changes between 2001 and 2019 (r 2 = 0.44, Supplementary Figs. 4, 3) and the national forest fire records between 2003 and 2009 (r 2 = 0.5, Supplementary Fig. 5). These results collectively suggested that the spectral-temporal segmentation method used in this study generally captured the spatiotemporal characteristics of disturbance in China's forests.
There were substantial spatial variations in the disturbance regimes across China's forests (Fig. 2, Supplementary Table 2). Among 0.5°grid cells, the annually disturbed forest area averaged 10.17 ± 11.78 km 2 yr −1 , with a highly left skewed distribution ( Supplementary Fig. 6), in which 75% of disturbance patches had an annual disturbed area <17.95 km 2 yr −1 (Supplementary Table 2). There were a few regions (e.g., mountainous coniferous forests in SW) that experienced relatively larger percent of annually disturbed forest area due to higher fire activity levels ( Supplementary Fig. 3). The percentage of annually disturbed forest, normalized by forest area in each grid cell, roughly showed the similar patterns with the annually disturbed forest area ( Fig. 2a and b). The annual disturbance frequency was 1.70 ± 1.67 patches yr −1 per km 2 of forest area with higher frequency in SE region, where forest management practices (e.g., plantation) provided the majority of timber in China since 2000s. The patch size of disturbance was 1.35 ± 0.45 ha, with larger patches located both in SE region due to intensive plantations and in NE region due to frequent, large forest fires. The patch size distribution was also left skewed, in which 75% of grid cells were 1.56 ha or smaller. The disturbance severity, calculated as the changes  of spectral vegetation greenness index 1 year before and 1 year after disturbance (i.e., ΔNDVI = NDVI t-1 − NDVI t+1 ), was 0.36 ± 0.51, with higher severity of disturbances generally located in the north of the Yangtze River probably due to higher percentage of fire occurrences.
Besides mapping of disturbances within forests, we also mapped the pattern and rate of forest expansion via reforestation/afforestation outside forests across China (See Methods). Consistent with recent studies 29,30 , our results showed that forest covers increased substantially in the ecologically sensitive regions, such as the Loess plateau, the upper reach of Yangtze River, and the northern Beijing ( Supplementary Fig. 7), which also experienced a relatively low and declining rate of forest disturbance (Figs. 2 and 3). The spatial synchrony between forest cover increases and disturbance decreases implied that forest restoration and conservation policy significantly affected the extent and dynamics of forest cover and disturbance in China. We also observed that the provinces with largest increases in forest area were located in SW region while simultaneously experiencing increased disturbance rates ( Supplementary Fig. 7, Fig. 3), suggesting the regional balance between forest conservation and timber production.
Trends in disturbance regimes. Disturbance regimes changed profoundly between 1986 and 2020, but their trends of change differed among regions (Fig. 3, Supplementary Table 3). At the national scale, the annual disturbance area decreased at rate of 390 ± 142 km 2 yr −1 (p < 0.01) from 1986 to 2020 (Fig. 1b), which was −1.6 ± 0.58% relative to the annually disturbed forested area (24489 ± 8707 km 2 yr −1 ). The national trend in disturbance rate was mainly driven by the decreasing trends in NE region ( Fig. 1c and Fig. 3a) while there was no significant trend in overall disturbance rate in SW and SE regions ( Fig. 1). Trends in disturbance regimes also exhibited a north-south gradient (Fig. 3). Generally, the disturbance regime (rate, frequency, and size) intensified in SE region and weakened in NE region. The disturbance severity generally showed a decreasing trend across China's forests with only a few areas in SE region that showed an increasing trend. The different trends of disturbance regime between northern and southern China broadly reflected the regional restoration strategies to promote ecological functions of forests in the slower growing regions (e.g., NE region) and sustain timber products in the faster growing regions (e.g., SE, SW region).
Role of forest restoration policy. To understand the underlying drivers of changing disturbance regimes, we used the generalized linear model (GLM) to investigate the effects of social-economic factors, climate changes, vegetation types, soil properties, and variability in Landsat time series on the percent of annually disturbed forest area at 0.5°grid cells. The GLM analysis explained about 42% of variation in annual disturbed area, in which forest restoration had a greater impact compared with other factors under investigation. Although the non-stationarity of variability in Landsat images may affect the detection of trends in disturbance, we found that its influence was relatively small (Fig. 4). Vapor pressure deficit, Normalized difference vegetation index (NDVI), and soil organic content had negative impacts on annually disturbed forest area, suggesting forests with favorable growth conditions experienced lower disturbance rates. On the contrary, regions with higher elevation and population density, and shorter distance from roads were more likely to experience higher disturbance rates (Fig. 4).
Of six major ecological restoration projects, the RSFP (Yangtze River and Zhujiang River Shelter Forest Projects) and NFPP were selected as the major drivers of forest management policy. However, different influences on annually disturbed forest area were observed by these two projects. There is decreasing trend within the NFPP region but increasing trend within the RSFP region ( Supplementary Fig. 8), indicating trends within in the NFPP region governed the overall decreased disturbance rate in China's forests. Specifically, annual disturbance rate was reduced from 1.14 ± 0.47% during 1985-1999 to 0.81 ± 0.57% during 2000-2020 at the national level. However, after the NFPP was implemented in 1999, the annual disturbance rate was significantly reduced (p < 0.05) by almost half from 1.21 ± 0.48% during 1986-1999 to 0.59 ± 0.42% during 2000-2020 within the NFPP regions (Fig. 5). The most significant decreases in disturbance rate occurred in NE region, the Loess Plateau, and the upper Yangtze River region, where timber harvest was significantly reduced, but reforestation/afforestation substantially increased forest covers ( Supplementary Fig. 7) to enhance ecological functions of forests. In contrast, the annual disturbance rate was only slightly reduced (p > 0.05) from 1.09 ± 0.46% to 0.98 ± 0.61% outside the NFPP regions. There was also a sharp contrast in the spatial patterns of trends for annually disturbed forest between 1985-1999 and 2000-2020, indicating significant changes in disturbance regime at the turn of the 21st century (supplementary Fig. 9). Therefore, our findings suggested the NFPP played a substantial role in changing disturbance regimes across China's forests since its implementation.

Discussion
Our results showed that about 1.16 ± 0.41% yr −1 , equivalent to 39.7% of total forested area was disturbed in China from 1986 to 2020. This figure was higher than the rate in European forests (0.55% yr −1 ) 11 , Canadian boreal forests (0.37~0.49 % yr −1 ) 17 , and contiguous United States (0.9~1.09% yr −1 ) 18,37,38 . However, there was a significant decline in disturbance rate across China's forests due to the effective implementation of protective forest policy for ecological restoration. On the contrary, increases in disturbance rate and severity that many parts of world were experiencing, were caused by either climate-driven wildfires in boreal forests and western U.S. 39 , management and droughts in Europe's forests 12,40 , or deforestation in tropical forests 3,41 . Changes in disturbance regime in China's forests were primarily governed by forest operation activities, ranging from more intensive timber harvest and plantation to less intensive thinning or understory treatment, similar to findings by 4 , suggesting restoration practices were key factors in regulating the disturbance regimes, and thus the structure and functions in China's forest ecosystems.
Our study provides high-quality quantitative and spatially explicit characterization of disturbance regime and its changes. Our findings highlight the role of forest restoration and conservation in driving the spatial and temporal variability in disturbance regime across China's forests. The NE region dominated by cold temperate forests is experiencing a decreasing disturbance rate. Historically, natural forests in NE region were one of the major timber providers in China. However, the NFPP first set strict limits on area, tree species, and age for selective logging to maintain a high-quality forest in 1999, followed by even stricter policy implemented in 2015 to prohibit all commercial logging in this region. These policies closely followed a shifting management objective of regional forest resources from timber production to sustainable supply of ecological services at the turn of the 21st century. In addition, wildfire historically was the dominant natural forest disturbance in NE region with a return interval between 30 and 100 years 42 . The strict fire exclusion policy, especially after the "black dragon" fire in 1987, has greatly lengthened the fire return interval in the region 43 . Since fires only accounted for <10% of disturbed area in NE region, the decreased disturbance regime was thus mainly caused by the reduced timber extraction. However, climate change will likely intensify the fire disturbance, compensating for the reduction in timber harvest in NE region (Supplementary Fig. 3b). In contrast, there was a significant increase in plantation in SE China to meet rising demand for forest products. As a result, more intensive forest management practices (e.g., timber harvest) led to higher disturbance rate and frequency, and larger individual patches in the region. Such trends are expected to continue in the near future. Forest fires in SE region accounted for around 84% of total fire occurrences in China, which were largely of anthropogenic origin with high frequency and relatively small patches 44 . However, there was an overall decreasing trend in forest area disturbed by fires from 1986 to 2020 ( Supplementary Fig. 3), suggestive of increasing intensity of forest operations in Southern China.
Our spatially explicit forest disturbance dataset provides a baseline to assess the impacts of changing disturbance on forest dynamics, biodiversity, and climate forcing. Forest management and wildfire are major disturbance types in temperate and boreal forests, but with large regional variations 4 . Our study confirmed that natural disturbance was a minor disturbance type in China's forests; Instead, forest management activities were the dominant drivers of disturbance regimes and had stronger effects on the extent, composition, and stand structure across China's forests. Over the past few decades, afforestation/reforestation projects have led to a continuous increase of forest cover which was also observed in the recent national forest inventory statistics 45,46 . Recent studies have found that increases in forest cover contributed greatly to increasing leaf area index or land greening 32,[47][48][49] , and the enhancement of China's terrestrial carbon sink over the past three decades 29,[50][51][52] . On the other hand, forest conservation projects, such as the NFPP, have effectively decreased the disturbance rate in China's forests and to what extent such decreasing disturbances contributed to leaf area index and carbon sequestration remains to be quantified. Our spatially explicit forest disturbance dataset provides a baseline to disentangle the contribution of increasing forest cover and changing disturbance regime to forest dynamics, carbon cycle, and climate feedback [53][54][55] . However, it is also worthwhile to note that although intensified afforest/reforestation has significantly increased forest covers in China, there is an increasing need to enhance forest structure and quality via effective management 26 .
Although the Landsat represents the longest and one of the best satellite observations to detect forest disturbances over large spatial and temporal scales 7 , the spatial, spectral, and temporal resolution offered by Landsat data are still not equipped to separate subtle structural and compositional changes caused by different silvicultural practices and chronic, low intensity disturbances. Thus, we lumped different types of disturbances together including disturbance events that caused abrupt changes in canopy cover or biomass through permanent forest losses, or different silvicultural practices (e.g., timber harvest and other silvicultural options), though we tried to separate forest fires from other types of disturbance. We also only mapped reforestation/ afforestation by tracking forest expansion without detecting postdisturbance establishment practices (e.g., planting and assisted regeneration) within forests that increased forest covers or biomass because of lack of high-quality reference data to distinguish forest recovery and these establishment practices. Future work should thus aim for the improved attribution algorithms that more explicitly consider the complex interactions between human and natural disturbance processes in China's forest ecosystems 56 . We also only mapped the largest disturbances for each 30 × 30m 2 pixel, in which there was only one disturbance event recorded for the whole 36-year period. For the short-rotation forest systems (e.g., fast-growing plantations in subtropical regions) or sequential disturbances (e.g., post-fire salvage logging), our approach may underestimate the disturbance rate if two or more disturbance events occurred in the past three decades. Finally, we did not further consider changes in disturbance regimes between natural and planted forests, because of a lack of detailed spatial and temporal distribution of these forests. Despite these limitations, remote sensing offers the only feasible way to map disturbance at large spatiotemporal scales and our quantitative analysis of patterns and trends in forest disturbances provides a first step towards better understanding of the drivers and ongoing changes in China's forest ecosystems.
In conclusion, we demonstrated that 1.16 ± 0.41% of China's forest was disturbed annually from 1986 to 2020. However, there was a decline in disturbance rate across China's forests mainly due to the effective implementation of protective forest policy. Generally, we found disturbance regime weakened in Northeast and increased in Southeast China. Our high-quality spatially explicit characterization of disturbance map provides an essential data layer for vegetation models to improve their capacity to represent disturbances and their roles in regulating ecosystem services provided by China's forests.

Methods
Disturbance detection. We used a well-established spectral-temporal segmentation method, Landsat-based Detection of Trends in Disturbance and Recovery (LandTrendr), to detect disturbances within the Google Earth Engine (GEE) cloudcomputing platform 57,58 . The core of the LandTrendr is to extract a set of disturbance-related metrics by breaking pixel-level annual time-series spectral trajectories into linear features using Landsat observations. The LandTrendr has been widely used for change detection in various forest settings, and details about the algorithms can be found in previous publications 57 . Here we briefly described the key steps in generating the year and type of disturbances in China's forests using the LandTrendr within the GEE platform. The overall analytic flow can be found in Supplementary Fig. 10.
First, we generated annual spectrally consistent time-series data by using all available, good quality (cloud cover ≤ 20) Tier 1 Landsat 5 (Thematic Mapper), Landsat 7 (Enhanced Thematic Mapper Plus), and Landsat 8 (Operational Land Imager) images acquired during the peak growing seasons (June 1-September 30) from 1986 to 2020. The peak growing seasons were selected to exclude compounding influences from ice, snow, and soil, and to maximize the spectral changes after forest disturbances. To tackle the spectral inconsistency among c-e showed the histogram distribution of percent of disturbed forest annually (% yr −1 ) before (blue) and after (red) NFPP at national (c), within the NFPP region (d), and outside the NFPP region (e). Blue lines in a and b indicated the NFPP boundary (also see Supplementary  Fig. 8a).
Landsat sensors, we harmonized spectral values via linear transformations according to band-respective coefficients presented in 59 . Clouds, cloud shadows, snow, and water were masked out using the Fmask algorithm 60 . The annual band composites at 30-meter spatial resolution during 1986-2020 were computed using the Medoid method 61 .
Secondly, we ran the LandTrendr using five spectral indices, including two spectral bands (shortwave infrared I and II that were B5 and B7), tasseled cap wetness (TCW), normalized burn ratio (NBR), and normalized difference vegetation index. These five indices were effective indictors to represent vegetation greenness and structures, and were commonly used for detecting changes in forest disturbance and recovery 62 . For each spectral index, the LandTrendr produced a set of parameters to describe a possible disturbance event at the pixel level, including spectral values at pre-disturbance level (preval), magnitude of change (mag), duration (dur) and rate of change (rate), and the signal-to-noise ratio (dsnr) (n = 5). Using these five spectral indices, we generated a stack of disturbancerelated parameter layers (n = 25, 5 spectral indices × 5 parameters), which were later used to detect and classify disturbances using machine learning models derived from reference data (described below).

Disturbance classification
Reference data. High-quality consistent reference data is key to train and classify disturbance types. To do so, we generated a total of 31225 reference points using a hierarchical approach. We first generated a large number of potential disturbance points using forest loss data from 2001 to 2020 3 . Then we separated fire disturbances from non-fire disturbances by overlaying MODIS burned area (BA) with potential disturbance points following the procedure used by 63 . Specifically, fire disturbances were determined if the MODIS BA data coincided with the Landsatderived forest loss for the fire year and 2 years postfire (i.e., t + 0, t + 1, t + 2) to account for delayed post-fire tree mortality. Following this step, we derived points as potential disturbances that consisted of fires and non-fire disturbances (including forest conversion to other land use types and silvicultural practices at various intensities). We also generated roughly the same number of points that experienced no disturbances (e.g., persistent forests), which were determined by selecting pixels with very few changes in spectral indices. These reference points, including fire, non-fire disturbances, and persistent forests, were then used to sample the timeseries spectral data from 1986 to 2020. Finally, time-series spectral data from each reference point were visually checked to make sure they accurately represented disturbance events. This process resulted into a total of 31225 reference data points, including 2356 fire disturbance points, 13,242 non-fire disturbance points, and 15,627 no disturbance points (persistent forests) (Supplementary Fig. 2).
Random forest classification. We used machine learning modeling to classify each pixel into fire disturbance, non-fire disturbance, or no disturbance. The reference data points were used to sample the LandTrendr-derived disturbance-related parameter layers described above, which resulted into a dataset consisting of disturbance types. We divided the dataset into 70% of training data, and 30% as validation data. Using the training data, a Random Forest (RF) model was trained to classify each reference point into fire, non-fire disturbance, or no disturbance. Our RF approach showed that short-wave infrared (SWIR)-based moisture indices (e.g., B7, TCW) were strong predictors for detecting forest disturbances ( Supplementary Fig. 11) likely because of their sensitivity to vegetation water content and canopy structure 64 . Finally, we applied the trained RF model to the full classification stack to consistently map the disturbance types from 1986 to 2020 across China's forests, assuming that the spectral trajectories derived from reference data period 2001-2020 can be extrapolated to the whole mapping period 1986-2020. However, note that our approach was meant to detect relatively acuate and discrete disturbances that caused canopy opening, rather than subtle changes of forest structure or composition resulted from low intensive silvicultural practices and chronic disturbances.
Year of disturbance. We used the LandTrendr to determine the year of disturbance as the onset of magnitude of spectral change. Since we ran LandTrendr on five spectral indices, there were five possible years of disturbance for each pixel. Thus, we determined the year of disturbance using the median value from at least three different indices (i.e., NDVI, NBR, TCW, B5, B7). In this way, we only kept pixels that were detected as disturbances using at least three indices, thus reducing commission errors.
The year with the greatest spectral changes generated by the LandTrendr often had an accuracy within 3 years 11 . A confidence level was also assigned to each disturbed pixel based on numbers of indices which showed possible disturbance events. Specially, low, medium, and high confidence were assigned if the disturbance was detected by three, four, or five spectral indices, respectively.
Validations. We validated the disturbance map at the pixel and national levels. At the pixel level, we validated the final map using the validation sub-sample described in the previous section. We derived a confusion matrix to report user's and producer's accuracy (Supplementary Table 1) as the main accuracy assessment metrics. At the national level, we compared forest disturbance detected in this study to available existing dataset. Specifically, we compared the area of forest fire disturbance between our study and the national fire records during 2003-2009 (Supplementary Fig. 5). We compared the disturbance rates between our study and Landsat-derived global forest cover changes from 2001 to 2019 3 (Supplementary Fig. 4).
Post-processing. We applied a series of spatial filters to minimize the unrealistic outliers from two potential sources of uncertainty, including speckle in time-series spectral trajectories or misregistration among images. This may lead to individual pixel or small patches including only a few pixels, which were (a) detected as disturbances, thus increasing the commission errors, or (b) not detected as disturbances, while their surrounding pixels were mostly disturbed, thereby increasing the omission errors. To address the issue (a), we removed all single-pixel disturbance patches through setting the minimum mapping unit as two 30 × 30 m 2 pixels (0.18 ha). To address the issue (b), we applied a 3 by 3 moving window to fill holes through assigning the year of disturbance based on the years in the surrounding pixels. Finally, we smoothed the year of disturbance by assigning the center pixel using majority rules from surrounding pixels within the 3 by 3 windows, thus accounting for artefacts associated with uncertainties in the correct identification of the disturbance year.
Characterizing disturbance regimes and their trends. We characterized the disturbance regime using five indicators within each 0.5°grid cell (n = 1946) across China's forests based on annual forest disturbance maps generated from the previous step. Within each grid cell, we calculated (1) total annually disturbed forest area (km 2 yr −1 ), (2) percentage of forest disturbed annually (% yr −1 ), as annual disturbed forest area divided by the total forested area, (3) disturbance size (ha), as the number of disturbed pixels for each individual patch using an eight-neighbor rule, (4) disturbance frequency (# of patches per 1000 km 2 forested area each year), as the number of disturbance patches per year divided by the total forested area, (5) disturbance severity (ΔNDVI = NDVI t−1 − NDVI t+1 ), as magnitude of NDVI change 1 year before and 1 year after disturbance, obtained from the LandTrendr analysis. We used (1) and (2) to characterize the disturbance rate, and (3)-(5) to describe the patch characteristics. The (2) and (4) were normalized by forest area within each grid cell, thus making them comparable among grid cells. For (3)-(5), we only calculated the patch size >0.45 ha (five 30 × 30-m 2 pixels), because patches <0.45 ha only contributed <10% of disturbance area, while greatly increasing the computation demand. Larger patches were also likely to be detected more accurately, having greater effects on forest dynamics and functions. In our analysis, the area was calculated by numbers of pixels multiplied by pixel size (0.09 ha).
We reported the five disturbance regime indicators at each grid cell, and aggregated patch-based metrics (disturbance size, frequency, and severity) within each grid cell. We visualized the spatial patterns of forest disturbance regime using arithmetic mean and reported the statistics using both arithmetic mean and median. Trends in five disturbance regime indicators at each grid cell were quantified using the Mann Kendall trend test, which is a non-parametric measure of monotonic trends in time series insensitive to outliers.
We also assessed how variations in satellite image availability affected the trend of disturbance. There is an increasing availability of Landsat images per year as Landsat mission evolves ( Supplementary Fig. 1). Although it provides us more opportunities to detect changes on the ground, this could potentially cause a false positive trend 13,14 . The number of available Landsat images was much lower from 1986 to 1999 but increased with time. As a result, there has been great interannual variability in disturbance area, with a non-significant decreasing trend. However, the number of available Landsat images has doubled since 2000, except that the number of available Landsat images was especially low in 2000 and 2012 due to the decommission of TM and ETM + sensors, respectively ( Supplementary Fig. 1). This resulted in a relatively small interannual variation in disturbance area and a decreasing trend. Although more detailed analysis (see below) indicated positive trends of Landsat observations contributed to the increasing disturbance rate during 1986-2020, our results showed that its influence was relatively small compare the other factors (Fig. 4). Furthermore, the declining rate of disturbance would be higher if the density of Landsat observation was consistent during 1986-1999, because the denser Landsat observations in the later observation period could create false positive trends.
Underlying drivers of changing patterns of forest disturbance. To understand the underlying drivers of changing patterns of forest disturbance, we used the generalized linear model (GLM) to explore the effects of climatic factors, biophysical factors (topography, vegetation, soil), social-economic factors (e.g., management, pullulation density, road density, and distance to the nearest road), and variability in Landsat observations on annually disturbed forest area at each 0.5degree grid cell. Climatic factors were used to investigate if background climates had an influence on disturbance patterns. We summarized annual mean air temperature (MAT), mean annual precipitation (MAP), and vapor pressure deficit (VPD) from 1985 to 2020 using 4-kilometer, monthly gridded TerraClimate data 65 . VPD was then only included in the analysis due to its collinearity with MAT and MAT. Forest management variables included percent of forest management area and the status of ecological projects. The percent of forest management area was calculated as percent of pixels with forest management activities divided by all forest pixels at 0.5°grid cell, using a newly derived forest management dataset for 2015 66 . The boundary of ecological project was used to determine whether disturbance patches were located in the ecological project regions (i.e., 1-within the ecological project region, 0-outside ecological project region, detailed descriptions of ecological projects were provided below). Population density (persons/km 2 ), road density (km/km 2 ), and distance to the nearest road (km 2 ) were summarized at 0.5°grid cell. Trends for Landsat observations were calculated using annual availability of all Landsat observations (Supplementary Fig. 1). Biophysical factors included topography (e.g., elevation, slope), vegetation (e.g., NDVI as a proxy for rate of growth), and soil (e.g., soil organic content, soil bulk density).
The GLM was fitted using annually disturbed forest area as dependent variables, and climatic factors, biophysical factors, management, and availability of Landsat observations as independent variables. The GLM was fitted using the glm function of the stats package in R. A parsimonious model was selected using the stepwise algorithm based on Akaike Information Criterion (AIC). The standardized regression coefficient from the GLM can be interpreted as the relative influences of independent variables. Partial plots were used to indicate direction of independent variables (Fig. 4).
Forest restoration policy and its effects on disturbances. Forest management policy has undergone substantial changes since the establishment of P.R. China in 1949. Before 1970s, excessive timber harvest decreased forest covers, degraded forest conditions, and caused serious ecological and environmental consequence, such as soil erosion, sandstorm, and flood. To restore the degraded forest ecosystems, Chinese government has implemented several national key ecological restoration projects since 1987, including Three-North Shelter Forest Program (TNSFP, also known as Green Great Wall), Yangtze River and Zhujiang River Shelter Forest Projects (RSFP), Natural Forest Protection Project (NFPP), The Grain for Green Program (GGP, also known as China's Sloping Lands Conversion Project). The TNSFP, initiated in the late 1970s, is an afforestation program, and spans half of northern China. It aims to improve environmental conditions, reduce the soil erosion, and promote the forest/grass-related products in the Three-North regions. The RSFP, initiated in 1989, is a forest restoration program, and covers the Yangtze River and Zhujiang River watersheds. It aims to reduce the soil erosion, rocky desertification, and occurrence of flood in south China. The NFPP, initiated in 1999, is a forest protection/conservation program to conserve China's natural forests, with the objectives of conserving biodiversity, protecting the water quality, preventing soil erosion and desertification, and reducing the likelihood of floods and other natural disasters associated with deforestation. Within the NFPP regions, the major restoration strategies include decreasing and adjusting the logging intensity, and accelerating afforestation in ecologically sensitive regions, especially along the Yangtze River and the Yellow River basin. The GCP, initiated in 2000, is a nation-wide ecological compensation program, which encourages farmers to convert croplands in hilly areas (>20 degree) to forests. Other than the forestrelated programs mentioned above, there are also a few other ecological projects, including the Beijing-Tianjin Sand Source Control Project and Returning Grazing Land to Grassland Project. However, these two projects focus on dryland ecosystems and thus were not included in our analysis 24,46 .
To further assess the effects of forest management policy on disturbance regimes, we compared the fraction of positive and negative trends of annually disturbed forest area at 0.05 significance level (i.e., p < 0.05) within and outside the ecological restoration projects boundary ( Supplementary Fig. 8). A higher proportion of negative trends within the ecological restoration projects boundary indicated a decreasing trend of disturbance rate caused by the effective implementation of the project. The boundaries of these projects were digitized at the province level from maps obtained from the administrative agencies.
Mapping forest expansion. Besides the mapping of disturbances within forests, we also mapped the pattern and rate of forest expansion via reforestation/afforestation outside of forests to fully capture disturbance regimes that regulated forest changes across China. Since the LT-GEE cannot map forest extent dynamics, we derived the forest extent (defined as canopy cover >20%) by merging the 30 m annual land cover map from 1990-2019 30 and an average tree canopy cover map from AVHRR (1985-2016) 32 , MODIS (2000-2020) 67 and Landsat (2000-2020) 3 . The annual forest expansion was defined as newly gained forests on previously non-forested land use types and remained forested for the rest of study period. The rate of forest expansion (km 2 yr −1 ) was mapped with each 0.5°grid cell (n = 1946) across China's forests (Supplementary Fig. 7).
Forest mask. Forest mask in 1986 was derived using the following procedures. First, a composite tree cover map in 1986 (TC compostite1986 ) was obtained by merging the tree cover map in 1986 (TC 1986 ) using observations by Advanced Very High-Resolution Radiometer (AVHRR) 32 and the tree cover map in 2000 (TC 2000 ) using observations by Landsat sensor 3 . The percent of tree cover in TC compostite1986 was determined as the maximum tree cover in TC 1986 or TC 2000 (that was TC compostite1986 = max(TC 1986, TC 2000 )). In this way, we could capture the possible forest disturbances that occurred before 2000 (e.g., TC 1986 > TC 2000 ), and the expansion of forested area from 1986 to 2000 (e.g., TC 1986 < TC 2000 ). Then the forest mask in 1986 was defined as TC compostite1986 > 20% following Liu et al., (2019). We should note that our study area did not include the newly afforested area after 2000. All analyses were performed within the forest mask, thus excluding the potential confounding factors from other land cover types. The description of TC 1986 and TC 2000 can be found in 3,32 .