Response of Ecosystem Water Use Efficiency to Drought over China during 1982–2015: Spatiotemporal Variability and Resilience

Ecosystem water use efficiency (WUE describes carbon-water flux coupling in terrestrial ecosystems. Understanding response and resilience of WUE to drought are essential for sustainable water resource and ecosystem management under increasing drought risks over China due to climate warming. Here we analyzed the response of ecosystem WUE to drought (spatiotemporal variability and resilience) over China during 1982–2015 based on an evapotranspiration (ET) dataset based on the model tree ensemble (MTE) algorithm using flux-tower ET measurements and satellite-retrieved GPP data. The results showed that the multiyear average WUE was 1.55 g C kg−1 H2O over China. WUE increased in 77.1% of Chinese territory during the past 34 years. During drought periods, the ecosystem WUE increased mainly in the northeast of Inner Mongolia, Northeast China and some regions in southern China with abundant forests but decreased in northwestern and central China. An apparent lagging effect of drought on ecosystem WUE was observed in the east of Inner Mongolia and Northeast China, the west and east regions of North China and the central part of Tibetan Plateau. Some ecosystems (e.g., deciduous needle-leaf forests, deciduous broadleaf forests, evergreen broadleaf forests and evergreen needle-leaf forests) in Central China, Northeast and Southwest China exhibited relatively greater resilience to drought than others by improving their WUE. Our findings would provide useful information for Chinese government to adopt a reasonable approach for maintaining the structure and functions of ecosystems under drought disturbance in future.


Introduction
Drought is an intermittent disturbance when water supply does not satisfy the demand for a long time [1,2], which can bring profound effects on carbon-water fluxes in terrestrial ecosystems at global and regional scales [3,4]. Drought-induced water stress could constrain vegetation growth, trigger wildfire and lead to crop yield reduction, grassland degradation and biological invasion [5,6]. Both frequency and intensity of global drought would increase in the context of climate warming until the year 2100 [7][8][9]. Knowledge about how terrestrial ecosystems adapt to and recover from drought stress is obviously critical for governmental agents to implement effective measures for sustainable ecosystem development and water resources management under future drought disturbances [10,11].
Water use efficiency (WUE) is a widely accepted parameter in describing carbon-water flux coupling between terrestrial ecosystems and the atmosphere and linking biological (e.g., vegetation transpiration and photosynthesis) and physical (e.g., soil evaporation) processes [12,13]. It reflects how much water is consumed by an ecosystem while it obtains carbon unit, which can be quantified as the ratio of carbon uptake (measured by GPP-gross primary productivity, NPP-net primary productivity or NEP-net ecosystem productivity) and water consumption (measured by ET-evapotranspiration or precipitation) [14][15][16]. Many previous studies have investigated the spatiotemporal changes of WUE. For example, Zhou et al. (2017) showed that underlying WUE (GPP × VPD 0.5 /ET) increased slowly during 1901-1975 but increased rapidly from 1976 to 2010 over the world based on terrestrial biosphere models [17]. Tang et al. (2014) found that the WUE had a similar latitudinal trend on earth, increasing from the subtropics to about 51 • N and then decreased in higher latitudes [18]. Ponce-Campos et al. (2013) indicated that the WUE varied with water availability at the global scale, with an increasing trend from moderate to extreme drought and a decreasing trend in wet conditions [19].
The ecosystem WUE is vulnerable to climate change because various abiotic factors (e.g., radiation, wind, precipitation and temperature) can have a huge impact on it [11,20,21]. Drought disturbances (less water available than in normal conditions) can influence plant transpiration and soil evaporation directly, and ultimately affect ecosystem WUE and productivity [2,22]. In the last three decades, the spatiotemporal responses of WUE to drought have been widely investigated. For example,  calculated ecosystem WUE in China using the Boreal Ecosystem Productivity Simulator (BEPS) model, which indicated that WUE increased in Northeast China and central Inner Mongolia but WUE was reduced in Central China when drought happened from 2000 to 2011 [23]. Yu et al. (2008) indicated that WUE estimated from site measurements enhanced in drier regions (relative to precipitation-abundant regions [24]) in eastern China because plants decreased stomatal conductance to adapt to water-limited conditions [25]. Moreover, Yang et al. (2016) investigated the relationship between global WUE and wetness index (WI) and found that WUE rose in arid ecosystems but declined in semi-arid/sub-humid ecosystems during drought periods mainly because of distinct sensitivities of different ecosystems to changes in hydro-climatic conditions [26]. However, Huang et al. (2017) suggested that the global WUE negatively (both positively and negatively) responded to drought in arid (humid) ecosystems during 2000-2014 based on the spearman correlation between the WUE and Standardized Precipitation-Evapotranspiration Index (SPEI) [20]. The inconsistency may arise from distinct data source, time periods, biomes and the way WUE is defined in different studies. Although previous studies explored how drought affected the WUE at the global and point scales, the responses of regional WUE across different biomes and eco-climatic zones to drought were still not well understood. Holling (1973) proposed the concept of system resilience to disturbances [27], which has been widely applied to terrestrial ecosystems in the past few decades [19,28]. The terrestrial ecosystem resilience represents the ability of an ecosystem to absorb disturbances from hydro-climate (e.g., drought) and to sustain or recover its structure and functions with ever-changing environments [28]. A resilient ecosystem could sustain or increase WUE to e nsure its productivity under water shortage conditions due to drought [29]. Ponce Campos et al. (2013) explored the relationship between WUE and drought during the drought in the early 21st century and found that Australian and North American ecosystems had strong biome-scale resilience according to increased WUE in the dry year [19]. Sharma et al. (2017) examined ecosystem resilience at the river basin scale in India during 2000-2014 based on the MODIS (Moderate Resolution Imaging Spectroradiometer) dataset and indicated that most basins in India were not resilient to disturbances [29].
Recently, previous studies have explored the spatiotemporal changes in WUE and the resilience in China [30,31], but how drought affected the WUE and the ecosystem resilience to drought in different biomes and nine sub-regions over China are not understood. In this study, the objectives are to (1) explore the spatiotemporal responses of WUE to drought over China during the period from 1982 to 2015 based on model tree ensemble-based evapotranspiration (MTE ET) and satellite-retrieved Global Land Surface Satellite (GLASS) Gross Primary Productivity (GPP) product; (2) assess the ecosystem resilience to drought in China at different scales (e.g., vegetation type and eco-climate regions). It will provide critical advice to the government for making reasonable ecosystem and water resources management strategies under natural and anthropogenic disturbances (e.g., drought).

Data Sources
Annual ET data used in this study were obtained from Li et al. (2018) [32], which estimated through machine learning approach (model tree ensemble, MTE) integrating eddy covariance ET and remote-sensing data from 1982 to 2015 at spatial resolution (hereafter MTE-ET, spatial resolution: 0.1 • ) [33,34]. Previous studies have demonstrated the capacity of MTE-ET to predict ET, even in regions where training data are not enough [32,35]. We used annual GLASS GPP product from 1982 to 2015 (hereafter GLASS-GPP, spatial resolution: 0.05 • ) produced by National Earth System Science Data Sharing Infrastructure, National Science and Technology Infrastructure of China (http://www.geodata.cn). It was estimated based on the Bayesian algorithm ensemble of eight light use efficiency models integrating the MODIS (Moderate Resolution Imaging Spectroradiometer) dataset, AVHRR (the Advanced Very High Resolution Radiometer) dataset and meteorological data [36]. The GLASS GPP was widely validated through 155 global eddy covariance measurements which contained nine terrestrial ecosystem types and applied in previous ecosystem carbon cycles studies at global and regional scales [36,37].
The widely used Standardized Precipitation-Evapotranspiration Index (SPEI, spatial resolution: 0.5 • ) developed by Vicente-Serrano et al. (2010) was applied which combined both flexible time scales and the sensitivity to changes in evaporative demand compared with other drought indexes (e.g., SPI-Standardized Precipitation Index and PDSI-Palmer Drought Severity Index) [38][39][40]. The 12-month SPEIs were selected to define annual drought severity during the period of 1982-2015 which can be directly acquired from the SPEI base version 2.5 products (http://digital.csic.es/handle/ 10261/153475). Annual gridded precipitation data (spatial resolution: 0.1 • ) was downloaded from the National Meteorological Information Center of the China Meteorological Administration (http://data.cma.cn/) at 613 meteorological stations in China, which was interpolated to raster data in order to be consistent with the resolution of ET and GPP through the Inverse Distance Weight (IDW) method which was weighted average interpolation based on the distance between points. Land cover map of 2015 (spatial resolution: 300 m ) from the European Space Agency Climate Change Initiative (http://maps.elie.ucl.ac.be/CCI/) was also aggregated to 0.1 • by using the nearest neighbor method and also used to classify vegetation types (we reclassified the 38 land cover types into nine categories in this study, namely croplands-CRO, grasslands-GRA, evergreen broadleaf forests-EBF, deciduous broadleaf forests-DBF, evergreen needle-leaf forests-ENF, deciduous needle-leaf forests-DNF, mixed forests-MF, shrub lands-SHR and non-vegetation regions-NOV, Figure 1a) over China [41].
In addition, the study was carried out at nine sub-regions according to the Chinese climatic and administrative boundaries [42,43] and four eco-climatic zones according to multiyear average precipitation in China, as shown in Figure 1b. Moreover, GLASS-GPP and SPEI datasets were also uniformly aggregated to 0.1 • using the nearest neighbor resampling method.  [20,26,45]. In order to investigate the lagging effect, we established linear regression models between 159 WUE and current-year SPEI as well as multivariable linear regression models among WUE, current-

168
The ecosystem resilience was considered as the ability of an ecosystem which could resist 169 external disturbances (such as drought) and sustain the same structure and functions under changing 170 conditions [19,49]. In this study, we used a dimensionless index (Rd) defined by Sharma et al. [29,50] 171 to quantify terrestrial ecosystem resilience to drought on a per-pixel basis: where Rd was ecosystem resilience index, WUEd was the WUE in driest year (the driest year can be 173 determined according to the minimum of annual precipitation on each pixel) and WUEm was the 174 mean annual WUE.

175
The value of Rd can be classified into four categories [29]: if Rd was greater than or equal to 1, the 176 grid cell was considered to have resilience and it meant that the ecosystem can sustain its productivity

Trend and Correlation Analysis
In our study, we assessed the trends of annual WUE using the least-square regression model during 1982-2015 on a per-pixel basis [44]. A positive value of slope indicated an increasing trend, and vice versa. We applied Pearson correlation coefficient to examine the relationship between WUE and SPEI in China. Previous studies found the lagging effect of drought on ecosystem WUE and the lagging effect meant that previous-year drought could affect current-year WUE in the ecosystem [20,26,45]. In order to investigate the lagging effect, we established linear regression models between WUE and current-year SPEI as well as multivariable linear regression models among WUE, current-year, previous-year SPEI on each pixel as follows [20,26]: where WUE current and SPEI current represented current-year WUE and drought, SPEI previous represented previous-year SPEI. k (including k 1 and k 2 ) and b were slope and intercept of linear model, respectively. Akaike Information Criterion (AIC) was applied to assess the degree of fitting optimization and to determine whether the inclusion of previous-year SPEI improved the regression model between WUE and SPEI [46]. If the difference of AIC values between Equations (1) and (2) was higher than 2, the Equation (2) is optimized [47,48].

Resilience Analysis
The ecosystem resilience was considered as the ability of an ecosystem which could resist external disturbances (such as drought) and sustain the same structure and functions under changing conditions [19,49]. In this study, we used a dimensionless index (R d ) defined by Sharma et al. [29,50] to quantify terrestrial ecosystem resilience to drought on a per-pixel basis: where R d was ecosystem resilience index, WUE d was the WUE in driest year (the driest year can be determined according to the minimum of annual precipitation on each pixel) and WUE m was the mean annual WUE.
The value of R d can be classified into four categories [29]: if R d was greater than or equal to 1, the grid cell was considered to have resilience and it meant that the ecosystem can sustain its productivity by increasing or maintaining WUE during drought periods. The ecosystem was considered slightly non-resilient if 0.9 < Rd < 1, moderately non-resilient if Rd lied between 0.8 and 0.9, and severely non-resilient if Rd < 0.8.
Moreover, all analyses were performed in MATLAB R2017a (MathWorks, Natick, MA, USA), and maps were drawn using ArcMap 10.2 (Environment System Research Institute, Redlands, CA, USA).

Spatial Distribution and Spatiotemporal Trends in Annual WUE over China
The spatial distribution of mean annual WUE showed large regional differences during the period 1982-2015 ( Figure 2)  The annual WUE showed a significantly increasing trend (slope: 8. among the nine sub-regions. For different vegetation types, the highest WUE was found for 193 deciduous broadleaf forests (2.48 g C kg −1 H2O), followed by deciduous needle-leaf forests (2.35 g C

Effects of Drought on WUE in China during 1982-2015
The correlation between WUE and SPEI showed a large spatial heterogeneity in China from 1982 to 2015 (Figure 4). The WUE positively correlated to SPEI in approximately 65.4% of Chinese territory (significant in 13.7% of Chinese territory at the 0.05 level) whereas negatively correlated to SPEI in approximately 34.6% of Chinese territory (significant in 4.5% of Chinese territory, mainly in the northeast of Inner Mongolia and Northeast China and for DBF, DNF and SHR).

240
In order to explore the drought lagging effects on WUE in China, in our study we also calculated 241 Pearson correlation coefficient between WUE and previous-year SPEI. As shown in Figure   In order to explore the drought lagging effects on WUE in China, in our study we also calculated Pearson correlation coefficient between WUE and previous-year SPEI. As shown in Figure 5, the spatial distribution of correlations mostly coincides with that shown Figure 4, which indicated that the WUE responded to current-year drought and previous-year drought in the same direction in most parts of China. However, the WUE was positively correlated with current-year SPEI but negatively correlated with previous-year SPEI in the central of Inner Mongolia and some areas of South China, Central China and Southeast China. WUE was negatively correlated with current-year SPEI but positively correlated with previous-year SPEI in the central of Southwest China, north and east parts of Northeast China. In the DBF, DNF and SHR vegetation types, the WUE showed positive correlations with previous-year drought but negative correlations with current-year drought, which led to even greater effects on carbon and water exchange processes during long-lasting droughts. In our study, we also used Akaike Information Criterion (AIC) determine whether the inclusion of previous-year SPEI improved performances of the regression model between WUE and SPEI, which could also reveal the lagging effects of previous-year drought on WUE. Figure 5 showed that the difference in AIC of the two equations was less than −2 mainly in the east parts in Inner Mongolia and Northeast China, west and east parts in North China and central parts in Tibetan Plateau, which indicated that previous-year drought lagging effects on WUE in these regions were relatively larger than other regions.

283
At vegetation type scale (Figure 6c), DNF, DBF, EBF and ENF exhibited relatively higher 284 ecosystem resilience to drought (Rd ≥ 1) during the past 34 years. The croplands, shrub lands and  Figure 6a shows the spatial distribution of ecosystem resilience to drought in China from 1982 to 2015, which evaluated the ability of an ecosystem to absorb and recover from the drought disturbances. On the whole, the terrestrial ecosystems were mostly resilient to drought in southeastern part of China and non-resilient to drought in northwestern China. Ecosystems in approximately 36.0% of Chinese territory were resilient to drought, which were mainly distributed in central parts

296
In our study, we estimated multiyear-averaged WUE (1.55 g C kg −1 H2O) and found large spatial

299
(2019) [31,32]. The highest annual WUE in Northeast China is probably due to its large areas of 300 deciduous forests (DBF and DNF) [23,30] while the higher annual WUE in Central China, North

301
China, Southeast China and Southwest China may be because of their suitable temperature, sufficient 302 water supply [23] or large amount areas of crop lands [51,52]. In previous studies, the WUE of forests 303 was usually higher than croplands, and grassland was usually the lowest [29,53,54] in those studies could be partly attributed to the different environmental conditions, plant species, At vegetation type scale (Figure 6c), DNF, DBF, EBF and ENF exhibited relatively higher ecosystem resilience to drought (R d ≥ 1) during the past 34 years. The croplands, shrub lands and mixed forests showed slightly non-resilient to drought with R d values varying from 0.9 to 1. Grasslands had the lowest R d (R d = 0.8148), which meant that grassland was relatively more vulnerable to drought among all vegetation types.

Discussion
In our study, we estimated multiyear-averaged WUE (1.55 g C kg −1 H 2 O) and found large spatial differences of mean annual WUE in nine sub-regions due to the heterogeneities in climate, soils and plant physiological characteristics [20,23], which agreed mostly with Sun et al. (2018) and   [31,32]. The highest annual WUE in Northeast China is probably due to its large areas of deciduous forests (DBF and DNF) [23,30] while the higher annual WUE in Central China, North China, Southeast China and Southwest China may be because of their suitable temperature, sufficient water supply [23] or large amount areas of crop lands [51,52]. In previous studies, the WUE of forests was usually higher than croplands, and grassland was usually the lowest [29,53,54]. However, the WUE in this study showed the following order: deciduous broadleaf forests > deciduous needle leaf forests > evergreen broadleaf forests > croplands > shrub lands > evergreen needle leaf forests > mixed forests > grasslands, which were similar to   [23] except for croplands. The discrepancy in those studies could be partly attributed to the different environmental conditions, plant species, planting areas and cultivation methods. Moreover, the significant increasing trend (8.5 × 10 −3 g C kg −1 H 2 O per year) of annual WUE found in this study which may be due to the increased atmospheric CO 2 concentration [17,55,56]. However, the WUE trend was much higher than Sun et al. (2018) (7.32 × 10 −4 g C kg −1 H 2 O per year) [30] which may induced by their different data sources. The significant increasing WUE in some regions of North China, Northeast China, Northwest China and Inner Mongolia may, to some extent, be owing to the implementation of ecological restoration projects (e.g., the 'Three North' Shelterbelt Development Program, the Grain-for-Green Project and the Beijing-Tianjin Sand Source Control Program) [57][58][59][60].
The relationship between WUE and SPEI in our study was positive in semi-arid which was consistent with Yang et al. (2016) [26]. The WUE in the Northeast China responded negatively with an increase in SPEI but Yang et al. (2016) showed positive correlations between WUE and drought in this region [26]. Moreover, our findings indicated that the relationship between WUE and SPEI in DBF was negative whereas was positive in Yang et al. (2016) [26] and Huang et al. (2017) [20]. These divergences could be due to differences in drought index, study areas and data sources. The different sensitivities of GPP and ET to droughts controlled the WUE-drought relationships across various sub-regions and biomasses [23,26]. For example, water supply played an important role for vegetation productivity [61] in arid and semi-arid regions (e.g., most parts in Inner Mongolia and Tibetan Plateau, Northwest China) with sparse vegetation (e.g., grasslands and croplands), and increase in precipitation (wetter) would mitigate water limitations and thus improve plant productivity in these regions [21,62], precipitation maybe was also consumed by soil evaporation, leading to an increase in ET [26,63]. However, the increased GPP was higher than ET, resulting in an increase in WUE. On the other hand, decrease in precipitation (drier) would trigger plants' water-saving adaptive strategy (partly close their stomata) and further reduce photosynthesis. The ecosystem resilience would thus decrease under drought stresses due to decreased WUE induced by reduced GPP and almost unchanged ET in such regions [11,64]. The croplands in humid and semi-humid regions (e.g., Dongbei Plains, Huabei Plains, Yangtze River basin and some regions of South China) also showed relatively lower resilience to drought (the WUE for this region showed positive response to drought) due to their vulnerability to water stresses. The crops would wilt or even die under water deficiency conditions, thus the WUE will decrease with increased ET and declined ecosystem GPP [65,66]. In some humid and semi-humid regions (Northeast China, the west of Inner Mongolia and some many areas of southern China) with abundant vegetation (i.e., DBF and DNF), energy was the main factor affecting vegetation growth. The incoming shortwave solar radiation was reduced when SPEI increased, hindering the absorption of carbon [23]. And soil evaporation and vegetation transpiration enhanced owing to additional water input, so the WUE decreased. When drought event happened, the ecosystems of those regions could absorb the disturbances from drought and sustain its GPP by increasing WUE because of their strong ecosystem resilience [26]. The above findings are generally consistent with  except for the central region of Inner Mongolia, which may be induced by the differences of data sources and time periods used [23]. Moreover, the anthropogenic factors (e.g., afforestation and deforestation) could also affect the ecosystem resilience in the west parts of in North China and Northwest China, the central of Inner Mongolia as well as some regions of Tibetan Plateau.
In our study, different vegetation types showed different WUE and their WUE showed different responses to drought. It may because GPP and ET in different vegetation types showed different sensitivities to climate change leading to the difference in carbon uptake and water consumption [65,66]. For example, grassland showed lowest WUE (Figure 2c), but its positive correlation between SPEI and WUE was the highest (Figure 4b). The growth and activity of grassland was highly dependent on water supply [67] and ecosystem productivity in grassland was closely related to precipitation [19,67]. Therefore, grassland had lowest GPP and lower ET [65,66] resulting in lowest WUE due to water limitations in arid and semi-arid regions. The growth rate and activity of grassland would decline when SPEI decreased (drought) [68] and the decreased GPP was higher than ET since GPP was more sensitive to drought than ET [26], leading to a decreased in WUE in the dry years. Moreover, DNF showed high WUE (Figure 2c), but its negative correlation between SPEI and WUE was the highest (Figure 4b). DNF was mainly distributed in humid and semi-humid regions, where energy supply was the main driver affecting vegetation functions [23,69]. Drought usually resulted in more incoming solar radiation (less precipitation and less cloud cover) and therefore in these radiation-limited, humid regions, drier, but also sunnier, conditions would imply a larger increase in GPP compared to ET, resulting in higher WUE [26,69].
Climate factors such as precipitation and temperature usually have lagging effects on vegetation growth [70,71]. In this study, we found that droughts also had lagging effects on ecosystem WUE in some regions (e.g., the east of Inner Mongolia and Northeast China, the west and east regions in North China and central parts in Tibetan Plateau). Yang et al. (2016) suggested that the relationship between WUE and SPEI in different years (previous-year and current-year) were in the opposite direction in all vegetation types over the world [26]. However, Huang et al. (2017) found that the impacts of drought in current-year and previous-year on WUE were in the same direction [20]. Their findings were generally not common in our investigation, which may be induced by the differences of data sources, methods and time periods used.
The main limitation of our study was the uncertainty of GPP and ET datasets. The ET data used in this study probably showed higher uncertainty in grasslands because few grassland flux-tower sites were applied in gridded ET estimation, which may result in slight underestimation of WUE [32]. Moreover, the interpolation of datasets with different spatial resolutions (e.g., the original spatial resolutions of ET, GPP and SPEI were 0.1 • , 0.05 • and 0.5 • ) may also cause uncertainties. Especially, high resolution drought monitoring datasets would undoubtedly be helpful to improve the robustness of related studies. Nevertheless, our findings revealed that the non-resilient ecosystems would encounter more challenges in ecological protection, agricultural production and carbon sequestration in the future due to increasing drought disturbances [29,72]. Policy makers could make reasonable water use strategies and efficient ecosystem management pathways to minimize the effects of drought under the warming climate.

Conclusions
In this study, we analyzed spatiotemporal variations in WUE and its response to droughts over China from 1982 to 2015. The ecosystem resilience to drought was then estimated using a dimensionless index at each grid cell. The multiyear-averaged WUE was 1.55 g C kg −1 H 2 O in China with relatively higher values in Northeast China and for deciduous broadleaf forests. The annual WUE increased (significantly increased) over 77.1% (56.2%) of Chinese territory. The annual WUE showed negative response to SPEI in forest regions in the northeast of Inner Mongolia, Northeast China and some regions of southern China as well as in deciduous broadleaf forests, deciduous needle leaf forests and shrub lands biomasses. An apparent lagging effect of drought on annual WUE was found in east parts in Northeast China, Inner Mongolia, west and east parts in North China and central parts in Tibetan Plateau. The terrestrial ecosystems showed resilience to drought in approximately 36.0% of Chinese territory mainly including several sub-regions (e.g., Central China, Northeast China, Southwest China) and vegetation types (e.g., deciduous needle-leaf forests, deciduous broadleaf forests, evergreen broadleaf forests and evergreen needle leaf forests). The results obtained in this study can improve our understanding on how ecosystem WUE responds to drought stresses and can inform local governments to take appropriate measures and strategies for sustainable water resources and ecosystem management in China.
Author Contributions: L.G. and W.L. designed and constructed this study framework; L.G. undertook data analysis and drafted the manuscript; W.L., Y.Z., F.S. and H.W. (Hong Wang). contributed to the manuscript; W.L., J.Z., H.C., B.D. and H.W. (Hongquan Wang). contributed significantly to our results section and helped to process the dataset.