Forest Resistance and Resilience to 2002 Drought in Northern China

Drought can weaken forest activity and even lead to forest mortality, and the response of different forest types to drought can be diverse. Deciduous broadleaf forest (DBF) and deciduous needleleaf forest (DNF) are two of the majority forest types in northern China. In this region, a severe drought event happened in 2002. However, due to the lack of data, the spatio-temporal characteristics of the ecosystem stability of different forest types here remain unclear. In this study, we used a machine learning approach (model tree ensemble, MTE) to drive fluxsite gross primary productivity (GPP), combined with remote sensing-based GPP and a vegetation index data (EVI), to analyze the spatial and temporal characteristics of resistance and resilience of DNF and DBF in northern China during and after the 2002 drought. The results showed that the DBFs were more acclimatized to the drought, while the resistance and resilience of DNF and DBF were diverse under different consecutive drought events. These results could be suggestive for forest conservation and vegetation modeling.


Introduction
As the largest terrestrial carbon sink, forests absorb nearly half of the carbon dioxide of the whole terrestrial ecosystem absorption [1]. Previous studies have shown that drought can affect the ecological stability of forest ecosystems by reducing forest productivity, or even lead to forest mortality [2][3][4][5]. However, spatial heterogeneity exists in the response of forests to drought. The response of different vegetation types to drought may differ [6][7][8][9]. Physiological activities of vegetation, such as photosynthesis and carbon allocation, might respond to drought differently as well [10,11]. Due to insufficient data, studies in some regions, including northern China, are still lacking.
A widely used data source for ecosystem stability is remote sensing. These products could cover the whole world. For example, a study based on MODIS-NDVI showed that drought influenced vegetation in North China seasonally [12]. However, this index is an indirect reflection of vegetation characteristics with spectral information, as most of the information came from canopies [11].
Due to the lack of direct measurements-such as from tree-rings-in the vast northern China, machine learning approaches could be employed for this circumstance. Machine learning has the ability to simulate the real world by training measured data [13], usually with a higher accuracy compared to other simulation methods. Therefore, machine learning has been gradually introduced for vegetation simulation in recent years [14][15][16]. These methods have been validated to perform well in simulating both spatial and temporal machine learning has been gradually introduced for vegetation simulation in recent years [14][15][16]. These methods have been validated to perform well in simulating both spatial and temporal dynamics of photosynthesis [17,18]. However, machine learning methods have been less involved in studying ecosystem stability.
Evidence showed that dry and hot extremes in northern China have been increasing in recent decades [19]. However, the ecosystem stability of forests to drought in this area is still unclear. In this study, we used a machine-learning simulated gross primary productivity (MTE-GPP), two remote sensing-based vegetation products (MODIS-EVI and MODSI-GPP), and the state-of-the-art drought index SPEI, to explore the ecosystem stability of forests to 2002 drought in northern China. With this study, we tried to answer three questions: (1) how did the forests respond to drought during and after the event? (2) was there any difference between the spatio-temporal characteristics of ecosystem stability of needleleaf forests and broadleaf forests in the northern China? (3) what are the similarities and differences in the performance of the three data sources: measured tree-ringwidth, machine learning-simulated GPP, and remote sensing-based vegetation data? Results from this study show diverse resistance and resilience of two forest types to drought in 2002, as well as the results from MTE-GPP and MODIS-GPP, which were the most consistent. These findings might be suggestive for forest conservation and vegetation modeling, respectively.

Study Area
Our research objects were selected from the whole northern China. The climate here transitions from a temperate continental climate to a temperate monsoon climate. Despite the climate heterogeneity of northern China, forests distribute broadly in this area, helping maintain the sustainability of the ecosystem ( Figure 1, the gray-colored background). Forest cover in northern China has been significantly altered in recent decades by afforestation, reforestation, and other human activities, which might have made the response of the forests to drought even more complicated. To disentangle the response characteristics of the typical forests in this region, we focused on forests where the landcover remained stable during our study period (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015).

Model Tree Ensemble GPP
Gross primary productivity (GPP) reflects the efficient of ecosystem photosynthesis, which quantify aspects of vegetation growth. Unfortunately, on the ground, there are only sparsely in situ flux measurements available for the GPP partitioning. Instead, we used a machine-learning method-model tree ensemble (MTE) [18]-to generate regional GPP from in situ flux measurements [20]. This method has been validated and used for regional to global carbon uptake analysis [16,21]. This GPP (hereafter MTE-GPP) was simulated at a spatial resolution of 0.1 degree and monthly temporal resolution, from 1981 to 2015. We resampled the MTE-GPP into 500 m before analyzing. The variables utilized in our simulation were listed in Table A1. Details of the data-driven product is available in Yao et al.'s work published in 2018 [20].

Satellite-Derived Vegetation Growth Data
Remote sensing data could reflect long time series and large-scale spatial information of vegetation. To compare with the response generated from MTE-GPP, the satellite-derived MODIS 500 m 16-day GPP product (MOD17A2H collection 6) was used, hereafter referred to as 'GPP' directly. Another type of satellite-derived vegetation index (enhanced vegetation index, EVI) has been widely used in previous researches about the responses of forests to droughts. The EVI describes the growth condition of canopies, and has been suggested as having improved sensitivity over high biomass regions, e.g., forests [22], compared to the other popular vegetation indices such as normalized difference vegetation index (NDVI). The responses of EVI to climate factor might be different from GPP, while only a few studies compared their performances on the responses of forests to drought [10,11]. To explore the differences between GPP and EVI, the MODIS 500 m 16-day vegetation index product (MOD13A1 Version 6) was chosen as EVI data source in our study, hereafter referred to as 'EVI' directly. It contains the global EVI values since 2001 to present. Both of the two datasets in our study area were mosaiced and extracted via Google Earth Engine (GEE). After that, we regenerated monthly EVI and GPP from these 16-day data, following an official procedure from MODIS team which considered days as weights [23]. The regeneration was performed in MATLAB.

Drought Index
The standardized precipitation evapotranspiration index (SPEI) was used for drought detecting in this study. The SPEI is a state-of-the-art meteorological drought index. It considered not only temperature and/or precipitation as the traditional drought index (e.g., PDSI and SPI), but also considered evapotranspiration (ET), which described the average water balance of the terrestrial ecosystem [24]. In this study, we downloaded the 12-month windowed SPEI (SPEI-12) dataset from the Global SPEI Database (https://spei.csic.es/database.html, 18 January 2021),since the annual water balance has been proved to be strongly correlated with forest growth [25]. The spatial resolution of this dataset is 0.5 degree. To make the pixels location fit to the MODIS EVI and GPP for spatial statistical analyses, we resampled the SPEI pixel into 500 m following the nearest neighbor rule.

Vegetation Category
The MODIS landcover data (MCD12Q1, version 006) and the vegetation map of China (1:1,000,000) were used for extracting and categorizing in our study. The vegetation map of China from the Institute of Botany of the Chinese Academy of Sciences [26] is derived from multi-source data, including forest inventory data, field samples, and remote sensing data [27]. Hence, this map reflects the vegetation in China more precisely than present remote sensing-based landcover data. This map is also used for MTE-GPP simulation among different vegetation ecosystems (see Table A1). Unfortunately, this map is neither dynamic nor near-real-time. Given the huge landcover changes during the last two decades in the northern China, we compared the vegetation map of China with the 500 m MCD12Q1 from 2001 to 2015, then selected those pixels which did not change in MCD12Q1 during 2001 to 2015, as well as were classified as the same forest type in both the vegetation map of China and MCD12Q1 (hereafter 'stable pixels'). The vegetation map was transformed from vector to raster format with a 500 m spatial resolution. As a result, most of the stable pixels were deciduous broadleaf forest (DBF, 859,587 pixels in total) and deciduous needleleaf forest (DNF, 57,792 pixels in total). The spatial distribution of forest pixels for this study is shown in Figure 1.

Defining Levels of Water-Balance Condition and Consecutive Droughts
The water balance dynamics of different regions are spatially heterogeneous. In the meantime, the vegetation has generally adapted to certain local hydrothermal conditions. Therefore, we employed the commonly used meteorological definition, which was based on percentiles of water index, to categorize the water balance condition of each grid individually. Here, the levels of water-balance condition are defined by the thresholds of SPEI percentile in Table 1. To count the percentile of historical SPEI of each pixel, we firstly sorted the mean SPEI12 of growing season (April to October) from 1902 to 2015 for each stable forest pixel. We did not include year 1901 because the 12-month window of SPEI12 made the valid SPEI12 values were available since December 1901. Then we sorted the 114 years' SPEI12 of each pixel. We used the 11th lowest (which was approximately the 10th percentile of 114 values) SPEI12 values as the threshold for extreme drought, which corresponded to happening once per decade. In the meantime, the values between rank 12 to rank 28 (about the 25th percentile, suggested the climate events happen once per four years) corresponded to moderate drought, as well as the values between rank 29 (25th percentile) to rank 85 (75th percentile) suggested normal years. This classification was also used for detecting normal years for the ecosystem resistance and resilience calculation in this study.
The acclimation of vegetation to droughts might be differing as the frequency of perturbance changing [8,28,29]. Therefore, whether consecutive droughts happened before and after our most concerned year was also considered and discussed. Specifically, we explored whether the ecosystem stability of forests were different between droughts that were preceded either by normal or drought year. We also tested the recovery of forests by whether a subsequent drought happened after the drought year.

Ecosystem Resistance and Resilience
The ecosystem resistance stands for the ability of an ecosystem to maintain its status under disturbance, and the ecosystem resilience describes the capacity to recover from the disturbance [28]. Both of them are part of the ecosystem stability. To quantitate the resistance and resilience to drought in our study area, we utilized the two equations following Isbell et al.'s research [28]: and where Y n is the forest growth during normal years among 2001-2015, Y e is the forest growth during the event year, and Y e+i is the forest growth i years after the event year, i = 1, 2, and 4. We decided to use the changing i because legacy effect widely exists, and was suggested among one to four years for forests [29,30]. As a result, the water balance of the forests in the succeeding four years should also be considered. Hereafter, resistance would be mentioned as 'Rt', and resilience is 'Rs'. Both of the Rt and Rs values were unitless, which made it convenient to compare the stability across varies biomes.

Standardization
A routine standardization (mean/standard deviation, µ/σ) method was used to preprocess the time series of monthly EVI, GPP, and MTE-GPP. The standardized results were then comparable for our research.

Pearson Correlation Analysis
The Pearson correlation analysis was employed to analyze if there is a trade-off between the ecosystem resistance and resilience, for each forest pixel. The student-t test was used to test the significance.

Analysis of Variance (ANOVA)
To detect if there is significant difference between DNF and DBF, one-way analysis of variance (one-way ANOVA) was conducted by using the MATLAB function 'anova1'. The F-test result could describe the significance.

Spatio-Temporal Characteristics of Droughts
The time series of mean growing season SPEI 12 ( Figure 2a To classify the drought condition spatially, we retrieved the minimum SPEI12 during growing season in 2002 for each forest pixel. Unsurprisingly, most of the forests in northern China experienced drought under different levels in 2002 ( Figure 2b). There were 97.37% of the forest pixels had a SPEI12 value lower than zero, in which 53.00% were under severe droughts (−2 < SPEI12 ≤ −1.5), and 23.65% experienced extreme droughts (SPEI12 ≤ −2). The severe and extreme drought regions were clustered in the northeast and east of northern China, which were covered with dense forests. Only 18.62% of the forest pixels were under moderate drought condition (−1.5 < SPEI12 ≤ −1). The result suggested that large areas of the forests in northern China experienced severe-to-extreme droughts in 2002.
1 Figure 2. Spatial and temporal characteristics of SPEI12: (a) The time series of mean growing season (April to October) SPEI 12 during 1902-2015 (blue line) and the overall percentiles of the stable pixels (solid red lines for 10th and 90th percentile, which stands for extreme drought and extreme wet event thresholds; dashed red lines for the 25th and 75th percentile, which stands for moderate drought and moderate wet even thresholds; and the gray-colored range marked the SPEI values of normal years); (b) the spatial distribution of water-balance condition in 2002 for the stable forest pixels, the climate events were classified with individual SPEI12 thresholds calculated from the 1902-2015 records of each pixels respectively.
To compare the areas of the DNF and DBF under drought, we counted the ratio of DNF and DBF pixels under vary water condition separately. The majority of the DNFs experienced drought events, in which 31.45% encountered extreme drought, and 68.42% were under moderate drought in growing season in 2002. Only 0.13% of the DNFs were under normal to wet conditions. To the DBF pixels, 59.33% experienced extreme drought, 24.81% were under moderate drought, and 15.69% were under normal to wet conditions, respectively ( Figure 3). Since drought events occurred extensively among our study objects, pixels under drought events in 2002 were selected for analyzing.

Interannual Variation of Forest Growth
The The similarity and differences among variations of the EVI, GPP, and MTE-GPP would lead to diverse ecosystem stability results, which might be reflected in spatial distribution or the differences between DNFs and DBFs.

Spatial Distribution of Ecosystem Stability of the Forests to Drought in 2002
We first explored the spatial distribution of the ecosystem stability of the forests in our study area. Hence, we calculated the resistance during 2002, and the resilience 1-, 2-, and 4-years after 2002, for each stable forest pixel.

The Resistance during Drought in 2002
Generally, the range of RtEVI and RtGPP was similar, which was [1.16, 84.38] and [1.00, 76.58], respectively. The range of RtMTE-GPP was the smallest, between 1.25 and 45.17 ( Figure  5a, outliers were excluded). This difference was due to the higher magnitude of MTE-GPP  The similarity and differences among variations of the EVI, GPP, and MTE-GPP would lead to diverse ecosystem stability results, which might be reflected in spatial distribution or the differences between DNFs and DBFs.

Spatial Distribution of Ecosystem Stability of the Forests to Drought in 2002
We first explored the spatial distribution of the ecosystem stability of the forests in our study area. Hence, we calculated the resistance during 2002, and the resilience 1-, 2-, and 4-years after 2002, for each stable forest pixel.

The Resistance during Drought in 2002
Generally, the range of Rt EVI and Rt GPP was similar, which was [1.16, 84.38] and [1.00, 76.58], respectively. The range of Rt MTE-GPP was the smallest, between 1.25 and 45.17 (Figure 5a, outliers were excluded). This difference was due to the higher magnitude of MTE-GPP in 2001, which resulted in a relatively sharp decrease in the drought year 2002. The spatial distribution of resistance generated from EVI (Rt EVI ) and GPP (Rt GPP ) showed similar pattern (Figure 6a vs. Figure 7a), albeit the Rt EVI were generally higher than the Rt GPP . More than that, the spatial distribution of Rt MTE-GPP was distinctly different from that of Rt EVI and Rt GPP (Figure 8a). In the north area, the high values of Rt MTE-GPP was not only higher, but also more spatial-clustered than the other two spatial distributions. Meanwhile, many pixels illustrated opposite results between Rt MTE-GPP and Rt EVI with Rt GPP . The pixels with Rt EVI and Rt GPP higher than 100 often had a Rt EVI value lower than 20 in both the east and the south of the study area. The results implicated different response regimes of MODIS-EVI, MODIS-GPP, and MTE-GPP to drought. in 2001, which resulted in a relatively sharp decrease in the drought year 2002. The spatial distribution of resistance generated from EVI (RtEVI) and GPP (RtGPP) showed similar pattern (Figure 6a vs. Figure 7a), albeit the RtEVI were generally higher than the RtGPP. More than that, the spatial distribution of RtMTE-GPP was distinctly different from that of RtEVI and RtGPP (Figure 8a). In the north area, the high values of RtMTE-GPP was not only higher, but also more spatial-clustered than the other two spatial distributions. Meanwhile, many pixels illustrated opposite results between RtMTE-GPP and RtEVI with RtGPP. The pixels with RtEVI and RtGPP higher than 100 often had a RtEVI value lower than 20 in both the east and the south of the study area. The results implicated different response regimes of MODIS-EVI, MODIS-GPP, and MTE-GPP to drought.   in 2001, which resulted in a relatively sharp decrease in the drought year 2002. The spatial distribution of resistance generated from EVI (RtEVI) and GPP (RtGPP) showed similar pattern (Figure 6a vs. Figure 7a), albeit the RtEVI were generally higher than the RtGPP. More than that, the spatial distribution of RtMTE-GPP was distinctly different from that of RtEVI and RtGPP (Figure 8a). In the north area, the high values of RtMTE-GPP was not only higher, but also more spatial-clustered than the other two spatial distributions. Meanwhile, many pixels illustrated opposite results between RtMTE-GPP and RtEVI with RtGPP. The pixels with RtEVI and RtGPP higher than 100 often had a RtEVI value lower than 20 in both the east and the south of the study area. The results implicated different response regimes of MODIS-EVI, MODIS-GPP, and MTE-GPP to drought.
The spatial distribution of resilience 2-years after the 2002 drought was different from the 1-year after resilience. Rs calculated from EVI, GPP, and MTE-GPP all illustrated higher recovery, especially in the northeast of the study area (Figures 6c, 7c and 8c). These changes supported the perspective that forests generally need more than one year to recover from annual droughts.
When it came to the fourth year (2006) after the 2002 drought, the overall mean growing season SPEI12 was exactly at the 25th percentile, which was the threshold of normal year and moderate drought. The changes of the resilience were unevenly distributed as Rs EVI-4 , Rs GPP-4 , and Rs MTE-GPP-4 shown (Figures 6d, 7d and 8d). Some areas showed continuously increasing resilience, while some others' resilience decreased. The inconsistency might be due to the spatial heterogeneity of hydrothermal conditions. These results called for further comparison considering the water balance condition pixel-by-pixel.

Stability of DBF and DNF to Drought in 2002 3.4.1. Comparison of the Resistance during 2002 between DBF and DNF
The resistance of DNF pixels (blue box) and DBF pixels (pink box) during 2002 is shown in Figure 5a. The median resistance values of DNF calculated from EVI, GPP, and MTE-GPP were 17.85, 17.45, and 21.78, respectively. The median resistance values of DBF were 21.86, 18.95, and 11.26, respectively. Resistance of DNF based on MTE-GPP was higher than that based on EVI and GPP. Meanwhile, resistance of DBF was significantly lower from MTE-GPP. Although the values among EVI, GPP, and MTE-GPP were different, the significant differences are prevalent between NF and BF in the three datasets (Table A2).

Comparison of the Resilience 1-, 2-, and 4-Year after Drought between DNF and DBF
The resilience of stable DNF and DBF pixels 1-, 2-, and 4-years after the drought in 2002 was shown in Figure 8b Figure A2), the resilience in different years indicated that both the DNF and DBF would fully recover in four years. Furthermore, similar to the resistance, significant differences are prevalent between the resilience of DNF and DBF in the three datasets (Table A2).

Comparison among Interannual Variation of EVI, GPP, and MTE-GPP
MTE-GPP and GPP generally showed similar variations along the 1-3 years after the drought year. However, the lower GPP value before and during 2002 drought led to higher resistance from the MODIS-GPP. This bias may be due to the simplification of the LUE modeling approach employed by the MODIS-GPP, which could not precisely depict the vegetation photosynthesis before, during, and after disturbances. On the other hand, the subtle differences between variation of EVI and GPP might come from the different vegetation activities that EVI and GPP reflected. EVI is an index indirectly reflected the physical characteristics of the plant canopy such as phytochrome and water content [31], while GPP is a physiological indicator that is directly influenced by a range of ecological factors [16,32]. Therefore, the response of EVI to drought can be considered as the response of physical condition of leaves to drought. This information indicated that if one meteorological drought event did not damage the canopy, the EVI might not show response to this drought. Prior studies have shown that EVI only responded to drought when the canopy was damaged [33]. Hence, MTE-GPP might be more suitable for assessing the response of forest to seasonal to annual droughts than both MODIS-GPP and MODIS-EVI.

Differences between Ecosystem Stability of DNF and DBF
The spatial trade-off between resistance and resilience of all DNF and DBF in northern China was shown in this study (Figures 6-8), which corresponded to previous studies [28]. However, the differences of ecosystem stability between DNF and DBF showed complicated characteristics. None of the ecosystem stability of DNF and DBF were similar, as the oneway ANOVA results shown in Table A2. This result implies that the balance between resistance and resilience might be biome-or-species-specific. Research on certain species has drawn similar conclusions-the stability of different biomes are heterogeneous [7,34]. Besides, the resilience of the DBF 1-4 year after the drought was relatively stable, which suggested faster recovery of the DBF than the DNF. Meanwhile, the DNF recovered more than the DBF in the second year after the drought. The result indicated that temporal trade-off of the DBF and DNF might not only exist along timescale for decades, but also along shorter timescale such as interannual. Moreover, higher Rs 2-years after the drought for the DNF from GPP and MTE-GPP indicated that the GPP is more sensitive to annual drought than EVI.

Comparison of Ecosystem Stability Considering Consecutive Drought
Research has mentioned the importance of drought frequency in studying stability of forest to drought [8,29]. For this reason, the influence of consecutive droughts was considered and discussed in this section.
It is difficult to fully explore the effects of varies consecutive drought combinations for multi-years. To simplify the scenarios, we only considered four typical water balance conditions in this study. First, we divided the pixels of stable forests by whether the year 2001 was a normal year for the pixel, which might influence the result of the resistance. Second, we classified the pixels by whether the year after the drought was a normal year. For the rest of the years, we only picked pixels under normal water balance conditions. This selection criterion could simplify and better help to quantify the recovery of the forest after no more than three consecutive years of drought. As a result, four climate combinations were grouped. Figures 9-11 showed the classified results generated from EVI, GPP, and MTE-GPP, respectively. All of the results showed that the forests that experienced drought in 2001-which were all DBF pixels-had a higher resistance in 2002, or the DBF could quickly acclimate to drought, and would decrease less under two-year consecutive droughts than under a single annual drought. The results from MTE-GPP and the MODIS products were slightly different. Generally, the resilience in 2006 was lower than in 2004 without experiencing drought. The result suggested that the response mechanism of forests to water-balance might involve more driving factors than temperature and water. Additionally, for pixels in which 2003, 2004, and 2006 were normal years, the DBF illustrated increasing resilience for all the three years after drought. Hence, the DBF might have acclimated better to various hydrothermal conditions than the DNF, and could even enhance the growth when the climate condition is suitable.

Conclusions
Based on one machine-learning simulation result and two remote sensing-derived datasets of vegetation growth, we found that EVI and GPP showed the response of forests to drought diversely, because they reflected the ecosystem stability of forests from different perspectives. Although a spatial trade-off exists broadly in forests in northern China, this correlation is non-linear. Moreover, the ecosystem stability is community-, biome-, or species-specific. The resistance and resilience also illustrated temporal trade-offs at a yearly timescale, as well as that the DBF was better acclimatized to drought disturbances. These results can provide information for forest management to consult in the future for

Conclusions
Based on one machine-learning simulation result and two remote sensing-derived datasets of vegetation growth, we found that EVI and GPP showed the response of forests to drought diversely, because they reflected the ecosystem stability of forests from different perspectives. Although a spatial trade-off exists broadly in forests in northern China, this correlation is non-linear. Moreover, the ecosystem stability is community-, biome-, or species-specific. The resistance and resilience also illustrated temporal trade-offs at a yearly timescale, as well as that the DBF was better acclimatized to drought disturbances. These results can provide information for forest management to consult in the future for a higher-frequency extreme events context. The comparison among MODIS-EVI, MODIS-GPP, and MTE-GPP suggested that MTE-GPP could be used for evaluating the response of forests to drought. Last but not least, both the MTE and satellite-based approaches of simulating GPP can be further improved by considering the variables and processes in more details.

Data Availability Statement:
The data presented in this study are available on request from the corresponding websites or authors.

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