Assessing Terrestrial Ecosystem Resilience using Satellite Leaf Area Index

: Quantitative approaches to measuring and assessing terrestrial ecosystem resilience, which expresses the ability of an ecosystem to recover from disturbances without shifting to an alternative state or losing function and services, is critical and essential to forecasting how terrestrial ecosystems will respond to global change. However, global and continuous terrestrial resilience measurement is fraught with di ﬃ culty, and the corresponding attribution of resilience dynamics is lacking in the literature. In this study, we assessed global terrestrial ecosystem resilience based on the long time-series GLASS LAI product and GIMMS AVHRR LAI 3g product, and validated the results using drought and ﬁre events as the main disturbance indicators. We also analyzed the spatial and temporal variations of global terrestrial ecosystem resilience and attributed their dynamics to climate change and environmental factors. The results showed that arid and semiarid areas exhibited low resilience. We found that evergreen broadleaf forest exhibited the highest resilience (mean resilience value (from GLASS LAI): 0.6). On a global scale, the increase of mean annual precipitation had a positive impact on terrestrial resilience enhancement, while we found no consistent relationships between mean annual temperature and terrestrial resilience. For terrestrial resilience dynamics, we observed three dramatic raises of disturbance frequency in 1989, 1995, and 2001, respectively, along with three signiﬁcant drops in resilience correspondingly. Our study mapped continuous spatiotemporal variation and captured interannual variations in terrestrial ecosystem resilience. This study demonstrates that remote sensing data are e ﬀ ective for monitoring terrestrial resilience for global ecosystem assessment.


Introduction
Terrestrial ecosystems face increasing pressure from both climate change and anthropogenic activities [1][2][3]. As many environmental components are changing, it is less likely that an ecosystem will remain unaffected. Thus, explaining the trends and processes of ecosystem changes within the context of large-scale alterations of environmental conditions constitutes a research priority of global relevance [4,5]. The concept of resilience has been proposed to analyze the response of ecosystems to short-term climatic disturbances [6]. Resilience is quantified by measuring the rate of the ecosystem required to return the ecosystem to the state that existed before the stress [7], and it generally focuses on the disturbance from which the ecosystem could recover to its original state without leading to ecosystem shift to an alternative state [6,8]. Many studies have quantified terrestrial ecosystem resilience using in-situ data [9][10][11][12][13]. Field-measurement-based studies can have high accuracy and reliability. However, these methods have been applied mainly at the regional scale, and global resilience measurement is almost absent from the literature.
Remote sensing can provide information over broad spatial extents that is difficult or impossible to obtain from field studies [14]. Moreover, continuous remote sensing data, obtained over several decades, enable the terrestrial ecosystem resilience measurement of long time-series [15]. Because of the characteristics of remote sensing and the key role of vegetation in ecosystems, vegetation has generally been the main object of satellite-based terrestrial resilience measurements [16]. Until recently, satellite-based resilience measurements have been focused mainly on the recovery capability of ecosystems to a specific type of disturbance. These methods generally detect the disturbance first then monitor the response of vegetation. Ecosystem recovery capability from fire has been widely studied in many types of terrestrial ecosystems [17][18][19], and measuring ecosystem recovery capability from drought is becoming a hot area of research [20]. In addition, various types of ecosystem resilience to other disturbances have also been studied extensively [9,21,22]. Compared with field-measurement-based research, remote sensing data increase the spatial coverage of resilience measurement. However, remote-sensing-based studies that measure terrestrial resilience from one specific type of disturbance are imperfect because ecosystems typically are being disturbed continuously, and because multiple disturbances may have an interactive effect on the ecosystem. Thus, detecting a complete disturbance and capturing the corresponding ecosystem anomaly are laborious tasks to achieve, especially at a global scale.
To determine terrestrial ecosystem recovery capability to multiple disturbances, several studies have directly detected anomalies of vegetation from remote sensing data to measure terrestrial ecosystem resilience. Global engineering resilience has been mapped using the satellite Normalized Difference Vegetation Index (NDVI) and reanalysis data [23]. Harris et al. used the length of time for vegetation to recover from disturbances as an indicator of ecosystem resilience and investigated spatial and temporal patterns of resilience of vegetation across southern Africa   [24]. Moreover, recent studies have developed several satellite-based indexes, which have been demonstrated to be potential key components of resilience. Seddon et al. [25] developed the Vegetation Sensitivity Index (VSI), which identified the global terrestrial resilience of ecosystems during the period of 2000-2013. Verbesselt et al. [26] measured the resilience of tropical forests using the temporal autocorrelation of vegetation based on the critical slowing down theory. These studies enabled measurement of global resilience, but the results are multiyear average values that cannot characterize the resilience variation or underlying causes.
In this study, we measured terrestrial ecosystem resilience by monitoring vegetation anomalies caused by disturbances. Because of the characteristics of the remote sensing data, this study focused only on disturbances that could lead to changes in vegetation. Resilience is related to the maximum stress from which an ecosystem can recover and is negatively correlated with the time it takes for an ecosystem to recover from disturbance [27,28]. Thus, we characterized the maximum stress from which an ecosystem can recover (Ms) and the time required for the ecosystem to recover from this maximum stress (Rt) as a ratio to measure terrestrial resilience [29]. Because Leaf Area Index (LAI) was directly estimated from remote sensing data using physical models, and LAI exhibited high accuracy for vegetation monitoring, we selected two LAI datasets (GLASS AVHRR LAI product and GIMMS AVHRR LAI 3g product) and two disaster factors (drought and fire) in this study to monitor vegetation anomalies, measure terrestrial ecosystem resilience, and analyze spatiotemporal variation. Overall, this study had three specific objectives: (1) to map global terrestrial ecosystem resilience using an advanced remote sensing product, (2) to characterize the spatiotemporal variations of global terrestrial ecosystem resilience, and (3) to analyze impact of potential drivers on terrestrial resilience distribution and variation.

Remote Sensing Data
The Global LAnd Surface Satellite (GLASS) LAI dataset with Advanced Very High-Resolution Radiometer (AVHRR) reflectance data, acquired from the Beijing Normal University, was used for this study [30][31][32] (Table 1). This long time-series data spanned from 1982 to 2016 at a spatial resolution of 0.05 • over an eight-day interval. In the GLASS LAI dataset, the physical methods are based on the inversion of canopy radiative transfer models through general regression neural networks (GRNNs). The GLASS LAI are both spatially and temporally continuous compared with other global remote-sensing LAI products. Moreover, the GLASS LAI dataset provides the longest time series, approximately four decades, and has been widely used for global vegetation monitoring and ecosystem assessment [33][34][35]. The Global Inventory Monitoring and Modeling System (GIMMS 3g) LAI dataset with an Advanced Very High-Resolution Radiometer (AVHRR), acquired from the National Oceanographic and Atmospheric Administration (NOAA), was also used for the analysis. The GIMMS AVHRR LAI 3g dataset spanned from 1982 to 2011 at a spatial resolution of 0.083 • and a 15-day interval [36]. The GIMMS AVHRR LAI 3g selected the latest version, termed the third generation NDVI dataset (GIMMS NDVI3g) as input data, which specifically aims to improve data quality in the high latitudes where the growing season is shorter than two months. It has also improved calibration that is tied to the Sea-Viewing Wide-Field-of-View Sensor, as opposed to earlier versions of GIMMS NDVI datasets that were based on inter-calibration with the Systeme Probatoire d'Observation de la Terre (SPOT) sensor. Moreover, the maximum value composite (MVC) method was adopted to remove atmosphere noise. The GIMMS AVHRR LAI 3g dataset is one of the high-quality global LAI datasets [4,32].
Fire data were acquired from the Along Track Scanning Radiometer (ATSR) World Fire Atlas (WFA). The ATSR-WFA provides a global fire monitoring service by using data acquired by the ATSR-2 and AATSR sensors from European Space Agency satellites. ATSR-WFA spanned from 1995 to 2012 at a spatial resolution of 0.05 • over a one month interval [37]. This dataset consists of two versions which were developed based on two different algorithms [38]. The version of algorithm 1 was used for disturbance detection in this study.
Land-cover data were obtained from Global Mosaics of the standard Moderate-resolution Imaging Spectroradiometer (MODIS) land cover type data product (MCD12Q1) for statistical analysis in this study [39]. The dataset boundaries are -180.0 • <= longitude <= 180.0 • ; -64.0 • <= latitude <= 84.0 • . The MCD12Q1 dataset spans 2001 to 2012 at an annual time step. The land-cover dataset has two version of spatial resolution, 5' * 5' and 0.5 • * 0.5 • , and the version of 0.5 • spatial resolution was used for analysis in this study. This dataset consists of 17 general land-cover types, which includes 11 natural vegetation classes, three developed and mosaicked land classes and three nonvegetated land classes [40,41].

Climate and Soil Data
Droughts were quantified using monthly time-series data from the Self-Calibrating Palmer Drought Severity Index (SC-PDSI) based on the Climatic Research Unit (CRU) TS3. 10.01 datasets from 2000 to 2009 at a 0.5 • resolution [42]. The SC-PDSI replaced the empirically derived climatic characteristic (K) and duration factors with values calculated automatically based on the historic climatic data of a specific location [43]. The computation of the SC-PDSI involves the classification of relative soil moisture conditions within 11 categories: extremely wet, severely wet, moderately wet, slightly wet, incipient wet spell, near normal, incipient dry spell, slightly dry, moderately dry, severely dry, and extremely dry [44].
Mean annual precipitation (MAP) and mean annual temperature (MAT) used in this study were extracted from the monthly precipitation and monthly temperature datasets in the Climate Research Unit (CRU) [45] version TS4.01, which spanned the 20th century and, therefore, covered our study period (1982-2011) with a spatial resolution of 0.5 • . MAP was calculated by summing monthly precipitation to generate annual precipitation then averaging the annual precipitation. MAT was calculated by averaging the monthly temperature to generate the annual temperature then averaging the annual temperature. The gridded datasets were obtained from more than 4000 meteorological stations and interpolated based on spatial autocorrelation functions [46,47]. These climate datasets were thoroughly evaluated and applied in global ecology and climate studies [48][49][50][51].
Soil moisture data were obtained from the United States Department of Agriculture-Natural Resources Conservation Service (USDA-NRCS) for statistical analysis in this study [52]. The global soil moisture regime data are based on an interpolation of over 20,000 climatic stations that were input into a soil water balance model to estimate soil moisture regimes. The soil moisture regime map data are rasterized on a 2-min grid cell. This soil moisture data provide one of the few high-quality soil moisture datasets on a global scale.

Data Preprocessing
To match all data we used in this study, we interpolated all datasets to a spatial resolution of 0.5 • . For fire data, we calculated the proportion of the fire area in each area of 0.5 • * 0.5 • . For the other data, we used bilinear interpolation for resampling. Averaging on a pixel level to coarser spatial resolutions is helpful for reducing data noise because the noise would be offset in the process of up-scaling the resample [53].
As our study concentrated on only the disturbance that did not lead to ecosystem regime changes, we derived the LAI grids, where the land-cover type was unchanged, from the GLASS LAI and GIMMS AVHRR LAI 3g datasets for further analysis. An unchanged pixel is defined as a pixel always being the same land-cover type during the 2001-2012 period. To diminish the impacts of disturbances caused by human activities, the grids that had urban and built-up land-cover type were excluded in further statistical analysis.
To increase the signal-to-noise ratio without greatly distorting the signal, a Savitzky-Golay filter was used pixel by pixel for smoothing long time-series data before calculating resilience index [54][55][56][57][58]. In this study, we set the half-width of the smoothing window size to 3 to get a smoother result and a higher computational efficiency.

A Satellite-Based Terrestrial Resilience Index
In this study, LAI was selected as the proxy of vegetation to measure the terrestrial resilience because of its high accuracy on vegetation monitoring. Based on the concept of resilience which has been put forward [54,55], we characterized terrestrial resilience as a ratio of the magnitude of maximum stress from which the LAI can recover (Ms) to the length of time (Rt) required for the LAI to recover to its normal range. This index only focuses on the disturbances that lead to vegetation that is worse than the normal level, and that do not cause ecosystem regime shifts [19,54].
We defined the LAI anomalies caused by disturbances that occurred when the LAI became lower than its normal range. We used the mean and the standard deviation of LAI to scope the normal variation range. In the period where the LAI was below the normal range, we detected the lowest LAI and marked this time point as the starting point for recovery. At this time point, the difference value between the lowest LAI and the normal range was defined as Ms (Figure 1), which denotes the maximum stress from which an ecosystem can endure and recover. The time point where the LAI came back to the normal range was defined as the end point of recovery, and the time span between start point and end point was defined as Rt. These two components were combined into an overall definition of terrestrial resilience. The ratio of Ms to Rt was defined as the terrestrial ecosystem resilience of the start time point. The time-series LAI was continuously detected to find the next anomaly from the time point behind the last end point. The time and duration of anomalies are different, indicating that in some years several anomalies occur while in some years no anomalies occur. Some rules were set up for annual terrestrial resilience calculation: The resilience values for which the recovery time crossed more than one year belonged to the resilience of the start point ( Figure 2). The resilience values for years where no anomaly was found, or where the last recovery process covered the entire year, were set as equal to the resilience value of the last year. The annual resilience values of years where more than one anomaly was found were calculated by averaging the values of resilience detected in the year. The fluctuation amplitude of vegetation in different land-cover types were various when disturbances occurred. Meanwhile, we found that large numbers of anomalies caused by disturbances were neglected when we used standard deviations of LAI in different land-cover types. To include more LAI anomalies caused by disturbances and exclude the LAI anomalies caused by data noise, a parameter (β) was integrated into the equation to adjust the normal range for different land-cover types: where Ms is the maximum amplitude of LAI when it is below the normal range, denoting the maximum stress from which an ecosystem can endure and recover; LAI ij(lowest) is the LAI in month j and year i, which is the lowest LAI in the period of anomaly; MEAN j is the month j mean LAI in recent decades (1982-2016 for GLASS LAI, 1982-2011 for GIMMS AVHRR LAI 3g); β is the parameter used to adjust the normal range; STD j is the month j LAI standard deviation in recent decades (1982-2016 for GLASS LAI, 1982-2011 for GIMMS AVHRR LAI 3g); and Rt is the length of time required for the LAI value to recover to its normal range from the maximum stress. The temporal resolution of two LAI datasets is different, causing the Rt to be a multiple of 8 (for GLASS LAI) and 15 (for GIMMS AVHRR LAI), respectively. In addition, the LAI values in the same area among various LAI datasets are not exactly the same. Thus, a min-max normalization was applied to normalize the two annual terrestrial resiliences for comparison and validation, respectively.

Regression Analysis
Linear regression was used to simulate the terrestrial resilience trend for each pixel during the study period, and the least squares method was used to determine the temporal trends: where θ slope is the slope of the terrestrial resilience trend line, n is the number of simulated years, i is the ordinal number of the year, and Resilience i is the terrestrial resilience value of i years. θ slope reflects the trend of the terrestrial resilience during the study period. Generalized additive regression models (GAM) is a generalized linear model in which the linear predictor depends linearly on unknown smooth functions of some predictor variables, and interest focuses on inference about these smooth functions [56,57]. GAM were originally developed by Trevor Hastie and Robert Tibshirani to blend properties of generalized linear models with additive models [58]. The mathematical description of the GAM can be given as below: where α is the usual model intercept, the functions f i are different nonlinear functions on variable x i and are modeled by splines, x 1 through x i represent the i values for the predictor variables, and ε i is the random variable. In this study, f (x) is the terrestrial ecosystem resilience. The relationship between terrestrial ecosystem resilience and mean annual precipitation (MAP) and mean annual temperature (MAT) was assessed by GAM.

Model Determination and Validation
Terrestrial resilience generated from remote sensing data generally consists of real terrestrial resilience and sensor and data noise [59]. In this study, we applied three methods to reduce the influence of noise on terrestrial resilience measurement, namely, resampling to coarser spatial resolutions, a Savitzky-Golay filter, and β optimization. In order to find the optimal β, drought and fire were selected as the two main disturbances [13]. We selected 52 natural disturbances from 24 points globally in the period of 1982-2011 (Table 2), since the variation ranges of LAI are various within different land-cover types [30]. All points were located in national parks in areas without artificial disturbance, and the 24 points covered the all land-cover types (Figure 3). We observed whether the LAI anomalies and fire or drought were matched by varying the value of β. The optimal value of β would lead to a high proportion of anomaly with disturbance and a low proportion of anomaly without disturbance. Finally, we determined the β for each land-cover type (Table 3).   With the optimal beta (β), global resilience was calculated. We took two approaches to validate the result. Firstly, we compared the resilience calculated from GLASS LAI and GIMMS AVHRR LAI 3g ( Figure 4). Secondly, the pattern of terrestrial resilience was compared with previous researches. Ecosystem resilience was a relatively loose concept without physical meaning, and previous studies were based on similar but different theoretical hypotheses; therefore, only the pattern of resilience was compared to evaluate global terrestrial ecosystem resilience. Keersmaecker et al. quantified global vegetation resilience using the 1981-2006 GIMMS NDVI dataset and analyzed the average situation of resilience during this period [23]. Similar patterns between our ecosystem resilience and this vegetation resilience were observed ( Figure 5). These two datasets showed relatively high resilience in central and western Europe and the northern Great Lakes, and low resilience in western Australia. In addition, Harris et al. used decay times of NDVI trends as a means to characterize the potential resilience of key biomes in southern Africa, and these results showed a pattern similar to our result [24]. Overall, the methods used in this study were reasonable for mapping global terrestrial ecosystem resilience. In addition, our method of determining global terrestrial ecosystem resilience avoided the null value areas that existed in the results of previous studies, and the method improved the temporal resolution of terrestrial ecosystem resilience from a multiyear average to annum.

Spatial Distribution and Temporal Variation of Satellite-Based Terrestrial Resilience
The annual mean terrestrial resilience was divided into five categories globally: low resilience (values of 0.0-0.2), relatively low resilience (0.2-0.4), medium resilience (0.4-0.6), relatively high resilience (0.6-0.8), and high resilience (0.8-1). Half of the resulting pixels showed low and relatively low resilience, and these were located in arid, semiarid, and temperate regions. Regions with medium resilience were found in Europe and tropical regions, more specifically, in the Amazon, central Africa, southeastern Australia, the Caribbean, and Southeast Asia. The pattern of relatively high resilience areas was similar to that of medium resilience areas, but with some subtle differences. For example, the area of relatively high resilience was only one-third as large as the area of medium resilience. Very few groups of terrestrial pixels, except those located in the Amazon and Southeast Asia, showed high resilience according to the multiyear average terrestrial resilience. Consequently, our model did not provide information regarding the resilience of regions consisting mostly of desert or of extremely sparse vegetation.
To compare the terrestrial resilience in different ecosystems, the mean resilience values of different land-cover types were determined ( Figure 6). Most forest ecosystems (e.g., broadleaf forest) showed relatively high resilience; however, deciduous needleleaf forest showed relatively low resilience. The gap in resilience between deciduous needleleaf forest and deciduous broadleaf forest was huge; in resilience from GIMMS AVHRR LAI 3g, the resilience of deciduous needleleaf forest was 58% lower than that of deciduous broadleaf forest. The resilience of mixed forest (mean resilience value (from GLASS LAI): 0.52) was slightly lower than that of deciduous broadleaf forest. Among the shrubland, open shrubland had a slightly lower value of resilience than that of close shrubland. Woody savannas and savannas showed large uncertainty on resilience value between the results from GLASS LAI and GIMMS AVHRR LAI 3g. In resilience from GLASS LAI, the resilience value of woody savannas was almost the same as the resilience value of savannas. However, in resilience from GIMMS AVHRR LAI 3g, resilience of savannas was only half of the woody savannas' resilience. Resilience of grassland and cropland from GLASS LAI were lower than that from GIMMS AVHRR LAI 3g. The resilience value of barren or sparsely vegetated was relatively low. This study adopted a linear regression model to characterize the variation of global terrestrial resilience. Figure 7 shows the results of the per-pixel linear trend analysis based on the results of two terrestrial resilience datasets. The blue and red colors indicate decreasing and increasing trends, respectively. Resilience decrease mainly appeared in the upper northern latitudes of Europe and, most notably, in Russia, western Finland, Sweden, Romania, eastern France, western Germany, eastern Poland, and Bulgaria. A trend of decline was also found in eastern North America, the Caribbean, southeastern Latin America, and southeastern China-regions that had a rather strong rate of decline of more than 2%/year (terrestrial ecosystem resilience from GLASS LAI). A trend of increase in resilience was found in central India, southern Africa, southern Latin America, and western America-regions that had a rate of increase of more than 2%/year (terrestrial ecosystem resilience from GLASS LAI). A relatively large discrepancy between resilience calculated from GLASS LAI and resilience calculated from GIMMS AVHRR LAI 3g was observed in most of Mexico and northern Australia, where the resilience trend calculated from GLASS LAI was higher than that calculated from GIMMS AVHRR LAI 3g. By comparison, the resilience trend calculated from GLASS LAI and resilience trend calculated from GIMMS AVHRR LAI 3g achieved good consistency at the middle and high latitudes of the Northern and Southern Hemispheres. This is probably because the time period and quality of LAI products were different [32]. In order to compare the variation of terrestrial resilience in different ecosystems, the mean resilience variation of different land-cover types was also calculated using the results from GLASS LAI (Figure 8). Among these ecosystems, resilience in evergreen broadleaf forest, deciduous broadleaf forest, and closed shrubland showed strong fluctuations. The fluctuation of resilience in evergreen broadleaf forest and deciduous needleleaf forest started to weaken in the 21st century, while the fluctuation of resilience in closed shrublands became strong after 2003. Resilience in woody savannas and savannas showed a continuous relatively strong fluctuation in recent decades. Resilience in evergreen needleleaf forest, mixed forest, cropland, and cropland/natural vegetation mosaic showed similar variation patterns, which decreased during the 20th century but started to increase in the 21st century. Resilience in deciduous needleleaf forest, open shrublands, and barren or sparsely vegetated land have remained steady in recent decades. Resilience in grassland and permanent wetlands has exhibited fluctuations, but has not shown obvious changes or trends in recent decades.

Driving Factors of Satellite-Based Terrestrial Resilience
On a global scale, mean annual precipitation (MAP) had a significant influence on terrestrial resilience. The latter manifests high sensitivity to precipitation for MAP values below 2000 mm. Within that high sensitivity range, we found that terrestrial resilience on a global scale was negative and decreased markedly as MAP dropped below around 500 mm (Figure 9). Above that value, the resilience transitions to positive. The impact of mean annual precipitation on terrestrial resilience remained steady and positive, with mean annual precipitation in the range of 2000 to 3000 mm. We found a huge uncertainty of MAP impact on terrestrial resilience between resilience from GLASS LAI and resilience from GIMMS AVHRR LAI 3g when MAP rose above 3000 mm.
The sensitivity of the terrestrial resilience to mean annual temperature (MAT) was not as strong as to mean annual precipitation (MAP). We found contrasting results of mean annual temperature's impact on terrestrial resilience between resilience from GLASS LAI and resilience from GIMMS AVHRR LAI 3g when mean annual temperature fell below around -5°C. As mean annual temperature rose, the effect of mean annual temperature on terrestrial resilience started to become vanishingly small. The negative effect of mean annual temperature was observed on terrestrial resilience when mean annual temperature rose up over 25°C on a global scale.
We found that soil condition exhibits a potential influence on terrestrial resilience. In most udic regimes, terrestrial resilience was high, while terrestrial resilience in most aridic regimes was low (Figures 7 and 10). For permafrost and inter forest regimes, the value of terrestrial resilience had large uncertainty and was not the main focus in this study, since effects of snow on remote sensing data were relatively large.
In addition to the effects of precipitation and temperature, we also found heavy fluctuations of terrestrial resilience were related to disturbance frequency. In recent decades, we observed three heavy fluctuations of natural disturbance frequency in 1989, 1995, and 2001, respectively. In each year of dramatically increasing disturbance frequency, the terrestrial resilience always significantly went down ( Figure 11). This relationship from GLASS LAI was more obvious than that from GIMMS AVHRR LAI 3g. Moreover, the correlation between terrestrial resilience and disturbance frequency was undefined when all the data of past three decades were taken into account.

Satellite-Based Terrestrial Resilience
Many different theories and potential methodologies of resilience have been developed since the concept of ecosystem resilience was put forward [6,[60][61][62][63][64]. Because ecosystem resilience is a broad concept, previous researches has developed many proxy indicators that have been used to present the recovery capability of ecosystem (e.g., resilience, resistance, return interval) [8,13,59,[65][66][67][68][69]. Every proxy indicator can reflect the recovery capability of an ecosystem based on different perspectives, although the similar calculation methods of these indicators. In our study, we measured the terrestrial resilience from the site at which the LAI began recovery, indicating that the resilience we measured did not include a capability of resistance. This can explain how deciduous needle-leaf forest had a strong capability to resist natural disasters but, in our study, had a low resilience value.
Despite many theories and methodologies for resilience measurement, long time-series terrestrial ecosystem resilience monitoring has faced several obstacles [70,71]. One of the main obstacles is choosing the appropriate data to represent the state of the ecosystem. The development of remote sensing has provided long time-series data of earth's surface [31,72]. Among advanced remote sensing products, LAI exhibits very high precision and long time series [25,[72][73][74]. Compared with Gross Primary Productivity (GPP) and biomass, LAI is a measurable index in the field, and LAI is directly estimated from remote sensing data using physical models, while GPP and biomass are usually generated from LAI, NDVI or other vegetation indexes and based on statistical modeling [75][76][77]. For this reason, LAI was selected to measure terrestrial resilience in this study.

Distribution and Variation of Satellite-Based Terrestrial Resilience
In this study, we measured terrestrial resilience in a broad sense to monitor the ability of an ecosystem to return to equilibrium or a steady state following disturbances [29,78]. Previous studies have suggested that terrestrial resilience can be affected by several component, such as climate, hydrology, radiation, and soil condition [9,19,79]. Climate, in particular, has a significant effect on terrestrial resilience [80,81]. Díaz-Delgado found that rainfall had a positive effect on ecosystem resilience response to fire in Mediterranean forest communities [19]. In tropical areas, mean annual precipitation has been found to have a positive effect on vegetation resilience [26,82]. Our study found that mean annual precipitation has more significant impact on terrestrial resilience than mean annual temperature globally in the most recent four decades. Moreover, huge uncertainty of MAP impact on terrestrial resilience between resilience from GLASS LAI and resilience from GIMMS AVHRR LAI 3g was observed when MAP rose above 3000 mm. This is probably due to the influence of clouds and haze in areas with more than 3000 mm MAP, increasing the uncertainty of the LAI datasets [9,26]. Patterns of relationships between terrestrial resilience and climate condition varied among different regions ( Figure 12). In deciduous needleleaf forest and savannas, mean annual precipitation did not show a significant impact on terrestrial resilience. The growth of deciduous needleleaf forest is limited by temperature because of its location. For savannas, the increase of precipitation has a contrasted impact on its growth between the wet season and dry season [83,84]. It is worth noting that mean annual precipitation not always has positive effect on all terrestrial ecosystems. Our study found in deciduous broadleaf forest, mixed forest, and closed shrublands that the effect of mean annual precipitation began to convert to negative when mean annual precipitation was excessive. In our study, no consistent relationships were found between mean annual temperature and terrestrial resilience on a global scale. This was also indirectly confirmed by a previous study that analyzed the relationship between temperature and tropical vegetation resilience [26,85]. Precipitation and temperature both have a huge influence on the state of vegetation, while the impact of mean annual temperature is insignificant on terrestrial resilience that was generated from the vegetation index. There are probably two reasons for this: firstly, the impact of mean annual temperature was offset when it was calculated on a global scale because the impacts of mean annual temperature were contrast among different ecosystems; secondly, the range of temperature variation was narrower and the process of temperature influence slow and long, which is not reflected in statistical results based on nearly four decades of data.
Previous researches have suggested that soil condition has a significant influence on terrestrial resilience. Soil can provide essential nutrition to vegetation for recovery from disturbances [52,86]. In this study, we found that resilience was higher in most udic moisture regimes and lower in most aridic moisture regimes ( Figure 10). Moreover, the results of the global mean resilience showed that savannas had low resilience. This agrees with the results of Ponce-Campos, who studied ecosystem resilience under large-scale alteration of the hydroclimatic conditions [79]. They found that savanna water-use efficiency decreased with increasing aridity, which indicates a low resilience of savannas to drought. This could be considered evidence of a relationship between terrestrial resilience and soil moisture regimes, because savanna ecosystems are located largely in aridic moisture regimes at a global scale.
Several studies have found a positive correlation between the time interval between wildfires and terrestrial ecosystem resilience, which indicates that the occurrence of frequent wildfires leads to low resilience [19]. Wildfire can alter ecosystem heterogeneity and also lead to degradation of ecosystem function. Although spatial heterogeneity can enhance the terrestrial resilience of an ecosystem [87], a degraded ecosystem may have lower resilience [80]. For this reason, the relationship between disturbances and resilience remains unclear. In this study, to explore the relationship between disturbance frequency and terrestrial ecosystem resilience, we calculated the number of annual anomalies in recent decades ( Figure 11). We observed three dramatic rises of disturbance frequency and three significant drops of resilience correspondingly. The reason for this is probably the dramatically increasing frequency of disturbance changes soil characteristics, reducing nutrient availability, and increasing soil susceptibility, which might prolong the time required for recovery from a disturbance. Furthermore, the short-term spatial legacies of disturbance can persist for years to decades [2]. High-frequency disturbances have an interaction effect and create a novel disturbance with potential impacts on the ecosystem [88]. Furthermore, the correlation between terrestrial resilience and disturbance frequency was undefined when all the data of past three decades were taken into account. This is probably because the impacts of disturbance frequency's slight fluctuation on terrestrial resilience can be offset by other environmental factors.
For the resilience variation among various ecosystems, we found the fluctuation of resilience in evergreen broadleaf forest and deciduous needleleaf forest have started to weaken in 21st century, while the fluctuation of resilience in closed shrublands became strong after 2003. Resilience in evergreen needleleaf forest, mixed forest, cropland, and cropland/natural vegetation mosaic decreased in the 20st century, then started to increase in the 21st century. Obvious change was not found from precipitation and soil conditions in these ecosystems [89]. This is probably because of the effects of El Niño, which was strongly monitored in 1997 and 1998 and led to disturbance frequency change [90]. Other causes need to be analyzed by future works.
In our study, we did not measure terrestrial resilience in certain highly managed ecosystems such as cropland/natural vegetation mosaic or urban and built-up areas, because they are not natural ecosystems. They might be well organized and highly productive, but they are subject to collapse without management due to their lack of a capacity to recover from disturbance [29,91].

Limitations and Recommendations for Future Research
Although remote sensing has advantages for large-scale and long-term environmental monitoring and leads to outstanding global terrestrial resilience measurement performance compared with field studies, it has three known limitations. Firstly, imperfect removal of cloud and noise in remote sensing products can affect terrestrial ecosystem resilience measurement, although these products promote high accuracy [26]. Secondly, remote sensing products cover relatively limited aspects, focusing mainly on vegetation, compared with the complexity of ecosystems [92]. Finally, the limitation of the number of field studies constrains extensive validation. In addition, a lack of uniformity of terrestrial resilience leads to a low availability of extant field data.
To make large-scale resilience monitoring more applicable globally, it is urgent to collect more resilience-related information for ecosystem resilience measurement. Most previous studies have used mainly temporal indicators to measure ecosystem resilience (e.g., temporal autocorrelation, temporal variance), but recent evidence has shown that spatial pattern dynamics could lead the change of resilience [8]. Several spatial indicators have been developed, and this facilitates more comprehensive methods of resilience measurement based on multi-indicators. Future research will consider both temporal information and spatial information when developing resilience indexes to improve the accuracy of global resilience measurement.

Conclusions
Terrestrial ecosystem resilience is an important component of ecosystem health and ecosystem function. The identification of large-scale metrics to quantify terrestrial ecosystem resilience and its dynamics remains a vital strategy for global ecosystem assessment. In this study, we mapped global terrestrial ecosystem resilience based on the long time-series GLASS LAI and GIMMS AVHRR LAI 3g products. The result showed that, globally, evergreen broadleaf forest exhibited the highest resilience and deciduous needleleaf forest exhibited the lowest resilience. On a global scale, the increase of mean annual precipitation has a positive effect on terrestrial resilience enhancement, while we found no consistent relationships between mean annual temperature and terrestrial resilience. For terrestrial resilience dynamics, we observed three dramatic raises of disturbance frequency in 1989, 1995, and 2001, respectively, and three significant drops of resilience correspondingly. On a regional scale, the distribution and dynamic of resilience have also been analyzed based on different land cover. Our study avoided the problems of some areas having null values encountered in previous studies, and this method improved the time scale of terrestrial ecosystem resilience to the interannual scale. This study demonstrated that satellites are effective for monitoring terrestrial resilience for global ecosystem assessment. With the availability of longer records of remotely sensed products, future research will consider both temporal information and spatial information when developing resilience indexes, which will improve the accuracy of resilience measurement.
Author Contributions: J.W. and S.L. conceived and designed the experiments; J.W. performed the experiments; J.W. analyzed the data; J.W. and S.L. contributed to the analysis of the results; all of the aforementioned contributed towards writing the manuscript. All authors have read and agreed to the published version of the manuscript.