Quantitative assessments of water-use efficiency in Temperate Eurasian Steppe along an aridity gradient

Water-use efficiency (WUE), defined as the ratio of net primary productivity (NPP) to evapotranspiration (ET), is an important indicator to represent the trade-off pattern between vegetation productivity and water consumption. Its dynamics under climate change are important to ecohydrology and ecosystem management, especially in the drylands. In this study, we modified and used a late version of Boreal Ecosystem Productivity Simulator (BEPS), to quantify the WUE in the typical dryland ecosystems, Temperate Eurasian Steppe (TES). The Aridity Index (AI) was used to specify the terrestrial water availability condition. The regional results showed that during the period of 1999–2008, the WUE has a clear decreasing trend in the spatial distribution from arid to humid areas. The highest annual average WUE was in dry and semi-humid sub-region (DSH) with 0.88 gC mm-1 and the lowest was in arid sub-region (AR) with 0.22 gC mm-1. A two-stage pattern of WUE was found in TES. That is, WUE would enhance with lower aridity stress, but decline under the humid environment. Over 65% of the region exhibited increasing WUE. This enhancement, however, could not indicate that the grasslands were getting better because the NPP even slightly decreased. It was mainly attributed to the reduction of ET over 70% of the region, which is closely related to the rainfall decrease. The results also suggested a similar negative spatial correlation between the WUE and the mean annual precipitation (MAP) at the driest and the most humid ends. This regional pattern reflected the different roles of water in regulating the terrestrial ecosystems under different aridity levels. This study could facilitate the understanding of the interactions between terrestrial carbon and water cycles, and thus contribute to a sustainable management of nature resources in the dryland ecosystems.

Temperate Eurasian Steppe (TES), the largest grassland belt in the world, stretches from Hungary to the northeast of China, lying between the boreal forest in Russia and the vast desert area in Central Asia. It spans around one-fifth of the longitude with an area of over 1.3×10 7 km 2 [9]. This eco-zone plays a significant role in the regional and global carbon and hydrological circulations, and also offers essential ecosystem services to maintain life [10][11][12]. Much of the region is controlled by the continental climate. With large spatiotemporal variations in precipitation, the water shortage issue is severe [13,14]. Recent climate change draws extra pressure on this region. Observed evidence is that regional annual precipitation decreased during the past 30 years [15]. Such trend increases the regional drought stress due to the higher water loss from the land surface. According to the recent analysis from Mohammat, Wang [16], the possibility of summer drought is enhancing since the 1990s.
Efforts have been done to understand the carbon and water cycles in this region at different spatiotemporal scales. Precipitation has been recognized as a primary limiting factor for vegetation biomass production in the arid and semi-arid ecosystems [17,18], while the energy availability (i.e., radiation) becomes more important with the amelioration of water condition [19]. Further results indicate that the terrestrial water availability can significantly influence the ecosystem WUE. Site level study showed that ecosystem WUE stimulated with increasing precipitation in a temperate steppe of Inner Mongolia [20]. But the regional scale results of WUE exhibit a two-stage pattern in response to water availability [21][22][23]. That is, WUE first increases with lower drought stress but decreases under the humid environment. Meanwhile, studies also suggested that some vegetational indexes (e.g. LAI) could be important predictors of WUE [21,24]. However, the terrestrial aridity condition has not been explicitly divided in the previous studies, which could induce certain concerns on the regional results. For example, in the widely studied "arid/semi-arid" of Inner Mongolia, the water conditions and grassland types are diverse, ranging from the humid meadow steppe in the northeast to an extremely arid desert steppe in the southwest. The aridity level varies greatly throughout the region (Fig  1). If we simply use the entire region to study the effects of climatic change on dryland ecosystems, then the results could be affected by including those relative humid areas. Furthermore, a great imbalance of research exists in this region. Maestre, Salguero-Gomez [25] summarized that most regional studies about global change were conducted in the eastern and westernmost part of TES, including the Mongol Steppe and East-European Steppe, whereas other subregions in TES were largely under-investigated, especially the Kazakh Steppe. Therefore, it is necessary to investigate those less-researched areas to improve our understanding.
In order to better address the issues, we used a late version of Boreal Ecosystem Productivity Simulator (BEPS) to simulate the carbon and water cycles in TES. The prototype model is a multi-cycle combined terrestrial biogeochemical model and has been proved to work well in grassland research [26,27]. To further improve the model performance in TES, we modified some important parameters and algorithms in NPP simulation. Thereafter, the long-term observation based Aridity Index (AI) was used to classify the region from arid to humid. In used widely to estimate regional or global terrestrial carbon and water fluxes [28][29][30][31][32]. Detail model descriptions could be found in previous studies [32][33][34][35][36]. Only NPP and ET calculations and major model modifications are introduced here. NPP calculation. The NPP is calculated by subtracting vegetation autotrophic respiration (R a ) from gross primary production (GPP): The GPP is calculated with a "two-leaf" scheme. The canopy is separated into sunlit and sunshade. Sunlit and sunshade LAI can be achieved using total leaf area index (LAI) and clumping index (O) of specific vegetation type. The instantaneous photosynthesis rate is calculated by the Farquhar model [37]. The maximum carboxylation rate (V c,max ) is an important parameter to determine the photosynthesis assimilation. Instead of using a preset value for the entire region, we incorporated a multi-factor based model from Zhang and Zhou [38]. The model was developed by compiling a series of experimental results from previous observations. The equations were established based on the statistical functions between observed V c,max and each single environmental factors. The model has been successfully validated using multiobservation from various vegetation types [38,39].
The model considers the effects of temperature, water, atmospheric CO 2 : Where V m,0 is the maximum carboxylation rate under an optimum environmental condition, f(Temp), f(SW), f(CO 2 ) are functions of atmospheric temperature (Temp)(˚C), surface soil water content(SW)(%), CO 2 concentration(CO 2 ) (Pa) that influence the value of V c,max .
R a is separated into two parts: maintenance respiration (R m ) and growth respiration (R g ). R g is estimated as proportions of GPP for different plant component according to previous studies. For R m , we replaced Bonan [40] equation in the original model with a component-specific scheme: where B i is a biomass of plant component i, A i is a carbon allocation fraction in i, p m,i and p g,i are coefficients of maintenance respiration and growth respiration, respectively. Q is the temperature sensitivity factor of respiration calculated as function of temperature following Arora [41]: sunshade parts: Where T p,lit and LAI p,lit are transpiration and LAI of sunlit leaves, and T p,shade and LAI p,shade are transpiration and LAI of sunshade leaves, respectively.
Then, T p and E s are calculated using the Penman-Monteith equation for no snow cover areas. For calculating E s , surface resistance is replaced by soil resistance according to Monteith [42].
E p is calculated as following: Where S ci is the intercept daily solar radiation (J m -2 day -1 ), A water is the absorptivity of solar radiation for water, LH v is the latent heat of vaporization (= 2.5×10 6 J kg -1 at 0˚C) and P i is the daily intercepted rainfall by canopy, which is determined as a proportion of LAI: Where k i is defined as the rainfall interception coefficient. The value is 0.3 mm LAI -1 d -1 . PREP is the daily precipitation.
Sublimation only occurs when snow exists and S p and S s are estimated by following equations according to Liu [2003]: Where A snow is the absorptivity of solar radiation for snow, LH v is the latent heat of sublimation (= 2.8×10 6 J kg -1 at 0˚C), SNOW is the snow water equivalent (mm), S is the daily total income solar radiation and m snow is a coefficient for the fraction of solar radiation transferred to latent heat by sublimation, which is set to 0.12 according to Saunders, Munro [43].

Descriptions of sites for model validation
The NPP results were tested against the field measurement results. The observed data were collected from Eastern Kazakhstan (KS sites), Xinjiang (XJ sites) and Inner Mongolia (IM sites). Each of these sites represents a typical vegetation pattern in the local area. No specific permits were required for the described field studies because all field sites were located on public land. The biomass data were sampled from four replicates within test sites. Total biomass was defined as the sum of the herbaceous biomass, root, litter carbon pools. The peak-season living aboveground biomass and litter were measured by destructive sampling of 1.0 m 2 , roots were collected by excavating a square of 1.0 × 1.0 m to a depth of 50 cm. For each test site, productivity was calculated from five repeated observations at the sample plots. We both run the original model (i.e., BEPS) and the revised model to test if our modification could improve the model performance in the study area.
Four grassland eddy covariance (EC) sites' data were used to validate the ET outputs. Three of them are located in Northern Kazakh Steppe and the other one is located in Inner Mongolia. Daily step data of latent energy (LE) with complete meteorological information were used to evaluate the corresponding model results. The data of three Kazakh Steppe sites were collected from the Fluxnet website (http://www.fluxdata.org/) and the data of Inner Mongolia site was offered by Meng Xiangxin and Fu Congbin. We selected the daily data of the EC sites with complete meteorological information needed for model simulation. All data were gap-filled and quality checked before the application.
The sites' information is displayed in Table 1.

Data preparation
Aridity index classification map. The AI map from the Consultative Group for International Agriculture Research (CGIAR) was used (http://www.csi.cgiar.org). The dataset follows the standard form proposed by United Nations Environment Programme (UNEP). It defines the AI as the ratio of Mean Annual Precipitation (MAP) to Mean Annual Potential Evapotranspiration (MAE): According to the original classification, terrestrial ecosystems are classified into 5 levels: hyper-arid (AI < 0.03, HAR), arid (0.03 < AI < 0.20, AR), semi-arid (0.20 < AI < 0.50, SAR), dry and sub-humid (0.50 < AI < 0.65, DSH) and humid (AI > 0.65, HU). However, few HAR areas exist in TES, so we combined that category with AR (Fig 1).
Model input data. Land cover map: MODIS land cover product for 2001 (MOD12Q1, V004, IGBP global vegetation classification, 1km resolution, http://ladsweb.nascom.nasa.gov/ data/) was used as model input because it is of high accuracy and has wide usage. The study region was extracted in Arc GIS 10.0. We used the version with the IGBP global vegetation classification. The TES map was extracted and resampled to an 8km resolution to fit the model input.
Daily meteorological data: Daily temperature, precipitation, radiation and specific humidity are required to drive the model. We imported the dataset of Global Meteorological Forcing Dataset for Land Surface Modeling (http://rda.ucar.edu/datasets/ds314.0/). The dataset is based on global observation datasets and NCEP/NCAR reanalysis. And it was updated and improved recently by using the results from the World Meteorological Organization (WMO) Solid Precipitation Measurement Inter-comparison and Global Precipitation Climatology Project (GPCP) daily product. Then it was evaluated by the second Global Soil Wetness Project (GSWP-2) dataset. The original dataset of precipitation and temperature were superimposed on re-sampled 8km Worldclim data generated from weather stations to reduce the deviation from spatial variation of topography. The regional MAT and MAP outputs in this study are also based on the daily temperature and precipitation from this dataset. Atmosphere CO 2 data: Monthly atmospheric CO 2 data were collected from Mauna Loa Observatory (MLO), Hawaii (20˚N, 156˚W) (http://cdiac.esd.ornl.gov/ftp/trends/co2/ maunaloa.co2). The atmospheric CO 2 concentration is assumed to be uniform across the region.
Soil texture data: Soil texture data were collected from the Global Soil Dataset for use in Earth System Models (GSDE, Available at: http://globalchange.bnu.edu.cn/research/soilw). The data are displayed as the volumetric percentages of silt, clay and sand. The original spatial resolution is 30 arc-seconds (about 1 km at the equator) and was bi-linearly interpolated into 8 km resolution.
The protocol of this study is available at: dx.doi.org/10.17504/protocols.io.h3mb8k6.

Spatiotemporal distribution of meteorological conditions in TES
Mean annual temperature (MAT) shows an obvious zonal characteristic from AR to HU subregions. (Fig 3A). High temperature areas are in the desert neighboring areas such as the Kyzylkum Desert, the Karakum Desert and the Gobi Desert, which are mainly characterized as AR. Low temperature areas are in the northwestern part of Mongolia, the northeastern part of Kazakhstan, the mountainous Altai region connected to Siberia in Russia, and the Tianshan Mountains, which stretches over Kyrgyzstan, Tajikistan and China. During the study period, the MAT of all sub-regions was increasing but no significant trend was found (Fig 3C).
Mean annual precipitation (MAP) increases from AR to HU, ranging from 167.3mm to 477.1mm (Fig 3B and 3C). Deserts bordering areas of AR experience little precipitation (<100mm). High MAP areas consist of HU, DSH and scattered SAR near the northwestern border with Russia, Mongolia, Inner Mongolia in China, and the open mountain grasslands and alpine meadows in Central Asia. The MAP in HU and DSH showed greater variations than those in relatively barren areas. Decreasing MAP trends were found in all sub-regions. However, only the MAP in SAR decreased significantly during the decade (p < 0.05, R 2 = 0.516).

Validation of simulated results
Comparisons between modeled NPP and observed NPP for the three sets of field measurements are showed in  The above statistical comparisons indicate that our modifications could improve the model performances and the model results are credible and suitable for this regional study.

Spatiotemporal distributions and variations of NPP, ET and WUE
The spatiotemporal distributions and variations of the NPP, ET and WUE are showed in Fig 6  and Fig 7. There was no significant temporal trend of the NPP in each sub-region during the 1999-2008 period (Fig 6A). However, the mean values were very different across the four subregions. The annual NPP of DSH and HU were significantly larger than AR and SAR. The highest value of 212.3 gC/m 2 was in DSH and the lowest value of 39.16 gC/m 2 was in AR. The NPP shows a clear zonal distribution from near-forest to near-desert (Fig 6A). High productivity grassland mainly concentrated in these areas connected with the forested areas, such as meadow steppe in the east part of Mongol Steppe, which is joined to Great Khingan and Siberia, alpine-meadow steppe around Tianshan and Altai Mountains, and Russian steppe areas that are close to the North Caucasian Forest. The NPP of these areas are generally higher than 100 gC/m 2 per year. During the decade, The NPP reductions were found in many parts of these areas, especially the Northern part of Kazakhstan and the alpine meadows near the mountainous area in Central Asia. The NPP mainly increases in the grassland of Inner Mongolia and Russia. Grasslands with lower productivity are those areas near the deserts (mainly with AR dominated), for which the NPP are typically lower than 30 gC/m 2 per year. In these barren areas, annual NPP variations were smaller and the trends were less identifiable (Fig 6B  and 6C). According to the regional statistics, during the period, more than 60% of the total areas showed a decreasing trend of NPP (Fig 7B).
The largest mean annual ET was 335.7 mm in HU, which was significantly higher than the other sub-regions. The values in DSH, SAR and AR were 241.9mm, 215.28 mm, 175.65 mm, respectively. The ET showed decreasing trends in all four sub-regions (i.e., arid, semi-arid, dry and semi-humid and humid sub-regions), but the trend was significant (P<0.05) only in SAR (Fig 7A). Grasslands with high ET are located in the mountainous Tajikistan, Kyrgyzstan, Eastern and western part of Kazakhstan (Fig 6A); Grasslands with moderate ET are located in the areas with high precipitation but with relatively low temperatures, such as the eastern part of Inner Mongolia. The ET decreased in more than 80% areas of TES during the decade ( Fig  7B). The largest reduction was found in SAR, approximately 10% areas of the sub-region showed significant decreasing trends. They were mainly distributed in the Southwestern part of Mongolia. The largest annual reduction was found in the Central Kazakh Steppe (Fig 6B  and 6C).
During the decade, the largest mean annual WUE was found in DSH, with a value of 0.88 gC mm -1 , and the lowest value was found in AR with 0.22 gC mm -1 . The inter-annual trends of WUE were not significant in all sub-regions except AR, where the WUE showed a significant enhancement (p < 0.05) (Fig 7A). An overall increasing trend of WUE was found throughout the region. Exceptions were some areas near the northern borders in Kazakhstan and Mongolia, where the NPP decreased significantly. However, even for those areas with negative trends, their WUE reductions were very small (Fig 6C). The largest WUE increases were temperature (MAT); (c) temporal variations of MAP and MAT from 1999 to 2008 in different AI sub-regions: AR (red lines), SAR (blue lines), DSH (magenta lines) and HU (green lines). https://doi.org/10.1371/journal.pone.0179875.g003 Water-use efficiency in Temperate Eurasian Steppe found in AR, where nearly 80% areas showed increase trends and more than 10% of the areas exhibited significant (p < 0.05) or very significant trends (p < 0.01) (Fig 7B).

The response of NPP, ET and WUE to climatic factors
Spatial distribution of correlationships between the NPP, ET and WUE and the climatic factors are displayed in Fig 8. The results showed that the correlationships between the MAT and NPP were ambiguous throughout the region. The positive relationships were mainly found in HU. The ET were negatively correlated with the MAT, which could come from the Water-use efficiency in Temperate Eurasian Steppe negative relationship between temperature and precipitation [48]. Significant correlations were found in the mountainous southern part of Kazakh Steppe. For the WUE, the grasslands in HU showed higher positive correlationships with the MAT than the other subregions (Fig 9).
The NPP were positively correlated with the MAP in all four sub-regions. The areas with significant correlationships were mainly in the semi-arid Central part of Mongolia, the Northern part of and the Western part of Kazakhstan and the desert steppe in the Central part of Kazakhstan (Fig 8). The regional ET exhibited an overall significant positive correlationship with the MAP, especially in the Kazakh Steppe. The significance of the correlationships clearly increased with aridity level. The area percentage of significant positive correlationship ranged from 68.7% in AR to 20.3% in HU (Fig 9). For the WUE, more than 80% areas in AR and HU were negatively correlated with the MAP, which were mainly in the Central and South end of Kazakh Steppe. The correlationships were not very clear in SAR and DSH sub-regions.

Evaluation and comparison of the model results
In this study, our modified BEPS was tested and used to simulate the carbon and hydrological cycles in TES. The model was validated against multi-source datasets, which indicated its ability to accurately simulate the regional grassland ecosystem. Research comparison shows that our NPP results are more comparable with those models that share the same canopy photosynthesis scheme (i.e., Farquhar biochemical model), but much lower than the result from the CASA model (a light use efficiency model). However, it seems that the CASA model over-estimates NPP because the value declined significantly after an improvement by Xing, Xu [49] (Table 2). Furthermore, different uses of input data could be another source of difference. For example, for BEPS, if topography is considered to interpolate meteorological data, the production decreases significantly in grassland simulation [50]. The soil effects in sparse vegetation region cause the bias of the MODIS LAI model in arid/ semi-arid areas [51,52].
Regarding to the ET results, the results produced in this study are very similar to the results from other remote sensing based models. . The result indicates that WUE could increase with lower aridity stress in a certain range, but decline under the humid environment. This two stage pattern corresponds to the regional results in the conterminous US and East Asia [21,58].
According to the AI classification, TES mainly consists of SAR (60.48%) and AR (23.14%). The AR and SAR sub-region are typically with high MAT and low MAP according to our regional investigation (Fig 3). The continentality is more obvious in the Kazakh Steppe due to its longer distance from the sea. In this large water-deficient region, the climate continues to become drier. The temporal results indicated that during the decade, the temperature did not show statistically significant trend, but the warming tendency is obvious over all the subregions. At the longer temporal scales, the regional warming has been agreed by all the observations available. The steady significant warming trend has been recorded for more than 60 years in Central Asia [59]. According to the latest IPCC report [1], Northern Eurasia is among the land areas with the strongest warming signals. Chen, Wang [60] further indicated that the temperature increase magnitude in this region is approximately twice as the average of the entire Northern Hemisphere. Comparing to the general consensus on the regional warming, the precipitation pattern is more variable. In the study decade, the temporal trend is varied over the different sub-regions. While in the dominated sub-regions (i.e., SAR and AR), MAP tends to decrease more significantly. The decadal pattern consists with the centennial regional analysis from Ma and Fu [61] and a global investigation of drylands [25,62]. Our results also support the predictions from Rind, Goldberg [63], who indicated that the precipitation in mid-and high-latitude of North Hemisphere would decline due to the projected climate change. However, Chen, Huang [64] reported that MAP was slightly increasing during 1930-2009 in the Central Asia. Nevertheless, all these studies agreed that the variability of MAP in arid/semi-arid areas is and will be increasing. The increasing variability of MAP, coupled with the regional temperature rising triggered the warm and dry tendency.
Under this condition, the regional WUE exhibits an increasing trend. However, it could not be concluded directly that the grasslands were getting better. Based on the regional results, we even find the decreasing trends of NPP in many areas. Hence, the regional increasing of WUE was mainly attributed to the ET reduction based on our results (Fig 6). This trend is more significant in the arid areas, especially the AR and SAR areas in the Northern part of Water-use efficiency in Temperate Eurasian Steppe Kazakhstan and the Southern part of Mongolia. Considering the good positive correlationship between MAP and ET, the decreased regional MAP could, therefore, indirectly induce the WUE enhancement. Besides, the ecosystem resilience to drier environmental condition could also contribute to the WUE increase [65]. For those species adapting to the extreme environment, they have developed strong stress tolerance abilities to survive effectively, e.g., reducing autotrophic respiration, conserving production into the soil and decreasing transpiration by stomatal regulation in response to the climatic deterioration [66][67][68]. Ecosystem evidence is that the NPP in AR showed higher resilience than the other sub-regions. In fact, only some of these resistance mechanisms were captured by the model via the LAI input and description of physiological and biochemical processes (e.g. NPP allocation, vegetation autotrophic respiration, and VPD calculation). Some other proved mechanisms are beyond current model prediction. The special aboveground morphology (e.g., leaf degradation) could enhance carbon assimilation efficiency with small LAI; the root structure of desert plants facilitates water absorption ability from deeper soils [69][70][71]. These processes secure the successful survive of desert plants under harsh environment conditions. Zhang, Li [70] further suggested that oversimplification of these processes in the terrestrial models would cause certain under-estimation of carbon sequestration in Central Asia. Furthermore, the prolonged drier environmental condition could also lead to the shift of species composition to become greater resistant to the drought events [72]. The physiological properties of drought-tolerant species could generate higher rates of water and carbon dioxide exchange, which is benefit for ecosystem functioning and sustainability [73,74]. Hence, the resilience effect of the arid ecosystems should be stronger than our current prediction. Currently, the CO 2 fertilization effect is regarded as an important contributor to the increasing C accumulation in terrestrial ecosystems [75,76]. Its impact to the temperate grasslands, however, is not very clear. Based on the global satellite observations, Donohue, Roderick [77] showed that 14% increase in atmospheric CO 2 led to a 5 to 10% increase in C assimilation in warm, arid environments. While the regional studies indicated that CO 2 fertilization might only play a minor role in the temperate grassland ecosystems. Mu, Zhao [78] indicated that although CO 2 fertilization has a strong impact on the carbon assimilation, its impact on temperate grasslands is the weakest with a contribution of only 0.3% of the total NPP variations. The results from the TEM showed that CO 2 fertilization could only counteract 1.4% of a NEP decrease [79]. The study also suggested that the greatest part of C assimilation increasing attributed to a rise in CO 2 was distributed in the sub-tropical and tropical ecosystems. Considering that CO 2 fertilization has little impact on ET [80], we suggested that the CO 2 fertilization will lead to the increasing of WUE, but the extent of this impact varied spatially. Further study is largely needed to promote our understanding to this process.

Study area
The temporal period in this study avoids the significant political impacts from the huge land revolutions like the collapse of the former USSR. However, the impacts of land use and cover change (LUCC) on carbon and water cycling still widely exist over the region. The current LUCC and its spatial heterogeneity are closely related to the human activity and governmental decisions. On the one hand, inappropriate human use and activities induce grassland degradation in arid/semi-arid areas. Human dominated activities like animal husbandry are regarded as the major factors leading to the negative land conversions [25,81]. In some traditional fine pastures (e.g. Xilinhot Steppe and Hulun Buir Steppe), the quality of grass resource is much worse than the nomadic period [82,83]. According to the assessment from Yusupov [84], overgrazing is the major factor to account for the land degradation in Central Asia with a contribution of 44%. Because of the recent boom of economic development, the impacts from other human activities such as mining and oil production are increasing as well [85,86]. On the other hand, the regional conservation programs were settled to address the land degradation problems. The most well-known projects are from the Chinese government. induced a land recovery from deserts to grasslands or forests [87]. Recent evaluation pointed out that the grassland area increased by 77,993 km 2 from 2000 to 2009, with 29,432.71 Gg C. yr −1 carbon accumulation increase in Inner Mongolia [88]. In addition, the current trend of urbanization moves a large amount of the rural population to the cities and contributing to the recovery of grasslands in the rural areas [89,90]. These factors thus interact with each other and influence the regional pattern of WUE. Further assessments of the LUCC impact will be largely needed in the future studies.
The role of water to grassland ecosystems varies according to different regional features of vegetation physiology and adaptation. The plant growth in arid ecosystem is mainly limited by the vegetational limitations, such as the high root to shoot rate and the low stomatal conductance. These biophysical features are unavoidable trade-offs between the survival in the extreme environment and the higher growth rate. This type of limitation weakens with a more favorable water condition, while the biochemical limitations such as light and nutrient inputs become stronger. Interestingly, we found that in this region, the WUE shares a similar pattern of the MAP correlationship in AR and HU sub-regions. We consider that this similarity should be attributed to different underling mechanisms. For the humid areas, precipitation is sufficient for vegetation growth and transpiration, so too much water input can induce ineffective water use and higher biogeochemical requirements. For example, Austin and Vitousek [91] showed that foliar nutrient availability decreased with increasing annual precipitation, which indicates the shift of major limitation to biochemical constrains. Thus, we found that in HU, the NPP showed the weakest correlation with the MAP. For the areas in AR, which is mainly characterized with desert steppe, both the vegetation production and evapotranspiration are limited by vegetational constraints (e.g. LAI). Hence both of them increase with extra water inputs. Our results corroborate the conclusions from [92-94], who proposed that in the arid and semi-arid lands, NPP is closely correlated with precipitation and hence the rain-use efficiency is a good indicator of NPP. With great potential to lose water in AR, the ET increases almost linearly in a wide range of water input [41]. However, vegetation production cannot increase continuously to match ET due to the intrinsic limitation from vegetation physiology and increasing nutrient deficiency [95,96]. Therefore, the combination induced the negative correlationship between MAP and WUE. Likewise, the response pattern is similar in SAR subregion. However, for this sub-region, the more favorable environment allows higher vegetation growth potential, and the water loss potential is weaker than that of AR sub-region (Fig 9). As a result, the corresponding response of WUE to MAP is less significant.

Conclusion and implications
In this study, we modified and used BEPS to simulate carbon and hydrological cycles in a typical dryland ecosystem, Temperate Eurasian Steppe. The Aridity Index classification contributes to specifically distinguish areas under various terrestrial water availability conditions and reveal different responses and underling mechanisms to climatic change. Results showed that this region tends to become drier with the increasing temperature and decreasing precipitation. Although regional WUE showed an increasing trend, it could not be concluded that grassland was recovering since regional NPP was decreasing. We consider that the trend comes from a combination of ET reduction induced by lower precipitation, and ecosystem resilience to drier environmental condition. According to the prediction from the multimodel ensemble-mean [97], the regional aridity, persistent droughts and extreme climate will increase. In the long run, climatic change will lead to further vegetation degradation and deteriorative land cover conversions, e.g. desertification. The aridity aggrandizement will, thus, put more influence on multi-aspect of grassland ecosystem, e.g. carbon sequestration, soil quality and species biodiversity. Since TES mostly covers the developing countries, both the water resources and the grassland ecosystem are essential to the human-being living and social wellness [98,99]. We have already noticed the importance to protect the grassland, but how to manage and preserve its ecological services and functions efficiently should be the next task for both international scientific community and local governments.
We find the two-stage pattern of WUE and the various responses to mean annual precipitation under different aridity levels in this region. The different responses in AR and SAR should also be noticed in future studies. We believe our findings could contribute to the mechanical understanding of carbon and hydrological circulations in dryland ecosystems and offer evidences to future terrestrial model development and improvement.
Supporting information S1 File. The websites to achieve data used in this study. This work used eddy covariance data acquired by the FLUXNET networks. We acknowledge the financial support to the eddy covariance data harmonization provided by CarboEuro-peIP, FAO-GTOS-TCO, Ileaps, Max Planck Institute for Biogeochemistry, National Science Foundation, University of Tuscia, University Laval and Environment Canada and US Department of Energy and the database development and technical support from Berkeley Water Center, Lawrence Berkeley National Laboratory, Microsoft Research eScience, Oak Ridge National Laboratory, University of California-Berkeley, University of Virginia.
We also greatly acknowledge Dr. Chen Jingming from Nanjing University for offering the original model, data and instructions, Dr. Luo Yiqi from the University of Oklahoma for helping to revise the paper and give suggestions, and Dr. Fu Congbin and Dr. Meng Xiangxin from Chinese Academy of Science for providing eddy covariance data in Tongyu. Without their kind helps, this work could not be finished.