# Global critical soil moisture thresholds of plant water stress

Zheng Fu, Philippe Ciais, Jean-Pierre Wigneron, Pierre Gentine, Andrew Feldman, David Makowski, Nicolas Viovy, Armen Kemanian, Daniel Goll, Paul Stoy, et al.

► To cite this version:

Zheng Fu, Philippe Ciais, Jean-Pierre Wigneron, Pierre Gentine, Andrew Feldman, et al. Global critical soil moisture thresholds of plant water stress. Nature Communications, 2024, 15 (1), pp.4826. 10.1038/s41467-024-49244-7. hal-04612805

HAL Id: hal-04612805

https://hal.science/hal-04612805v1

Submitted on 30 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Global critical soil moisture thresholds of plant water stress

Received: 7 November 2023

Accepted: 22 May 2024

Published online: 06 June 2024

<figure><img src="image_3.png" /></figure>

Check for updates

Zheng Fu<sup>1,2</sup>✉, Philippe Ciais<sup>2</sup>, Jean-Pierre Wigneron<sup>3</sup>, Pierre Gentine<sup>4</sup>, Andrew F. Feldman<sup>5,6</sup>, David Makowski<sup>7</sup>, Nicolas Viovy<sup>2</sup>, Armen R. Kemanian<sup>8</sup>, Daniel S. Goll<sup>2</sup>, Paul C. Stoy<sup>9</sup>, Iain Colin Prentice<sup>10,11</sup>, Dan Yakir<sup>12</sup>, Liyang Liu<sup>2</sup>, Hongliang Ma<sup>13</sup>, Xiaojun Li<sup>3</sup>, Yuanyuan Huang<sup>1</sup>, Kailiang Yu<sup>2</sup>, Peng Zhu<sup>14</sup>, Xing Li<sup>15</sup>, Zaichun Zhu<sup>16</sup>, Jinghui Lian<sup>2</sup> & William K. Smith<sup>17</sup>

During extensive periods without rain, known as dry-downs, decreasing soil moisture (SM) induces plant water stress at the point when it limits evapo-transpiration, defining a critical SM threshold ($\theta_{crit}$). Better quantification of $\theta_{crit}$ is needed for improving future projections of climate and water resources, food production, and ecosystem vulnerability. Here, we combine systematic satellite observations of the diurnal amplitude of land surface temperature (dLST) and SM during dry-downs, corroborated by in-situ data from flux towers, to generate the observation-based global map of $\theta_{crit}$. We find an average global $\theta_{crit}$ of 0.19 m$^3$/m$^3$, varying from 0.12 m$^3$/m$^3$ in arid ecosystems to 0.26 m$^3$/m$^3$ in humid ecosystems. $\theta_{crit}$ simulated by Earth System Models is overestimated in dry areas and underestimated in wet areas. The global observed pattern of $\theta_{crit}$ reflects plant adaptation to soil available water and atmospheric demand. Using explainable machine learning, we show that aridity index, leaf area and soil texture are the most influential drivers. Moreover, we show that the annual fraction of days with water stress, when SM stays below $\theta_{crit}$, has increased in the past four decades. Our results have important implications for understanding the inception of water stress in models and identifying SM tipping points.

The critical soil moisture threshold ($\theta_{crit}$) of plant water stress is defined as the soil moisture (SM) level at which evapotranspiration becomes SM limited in that environment¹. Below this threshold, a marginal reduction of SM reduces evapotranspiration and increases sensible heat emissions and surface temperature², making the air above the canopy warmer and drier, which in turn further reduces evapotranspiration and plant carbon dioxide uptake³⁻⁵. The control of energy partitioning regimes across $\theta_{crit}$ determines local climate through land-atmosphere coupling and can amplify warming during droughts⁶,⁷. A better knowledge of $\theta_{crit}$ is thus important for land-atmosphere interactions⁵, for climate studies⁸⁻¹⁰ and for understanding the vulnerability of ecosystems and crop yields to drought⁹.

The relationship between SM and the evaporative fraction (EF), defined as the ratio of evapotranspiration to net radiation, shows two distinct regimes<sup>2,5,11,12</sup> (see Fig. 1a). When SM is higher than $\theta_{crit}$, the system is non-water limited (energy limited) and SM does not impact evapotranspiration⁵. In contrast, when SM is lower than $\theta_{crit}$, the capacity of plants to extract soil water by roots and xylem transport becomes progressively reduced. The system becomes SM limited, and evapotranspiration decreases with decreasing SM until leaves fully close their stomata, direct evaporation at the soil surface ceases, or roots are no longer able to take up soil water (the wilting point)¹⁰. The overall EF-SM relationship (increasing below $\theta_{crit}$ and then plateauing) is conceptually well established, but a spatially explicit understanding

<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 1 | Quantifying the critical soil moisture (SM) threshold during soil dry-downs using the evaporative fraction (EF)-SM method and the land surface temperature diurnal amplitude (dLST)-SM method.** An example of estimating SM thresholds ($\theta_{crit}$) from the EF-SM method (a) and the dLST-SM method (b) using all dry-downs at Hainich beech forest site (DE-Hai, Supplementary Table 1).</figcaption></figure>

c Comparison between the SM thresholds estimated from the dLST-SM method and EF-SM method across all sites. The median and the 25th, 75th percentiles are shown for each biome. The dashed line is the 1:1 line while the red line is fitted line using least squares regression.

with accurate global maps of $\theta_{crit}$ are lacking, due to a lack of global high-frequency observations of EF<sup>11,13-15</sup>. Therefore, the factors that control the global variations in $\theta_{crit}$ are poorly known. Earth system modelers have adopted simple parametric representations of EF-SM relationship and $\theta_{crit}$ to describe soil water stress and land-atmosphere feedbacks<sup>5</sup>, leading to model biases which hinder our ability to predict drought and its ecosystem impacts<sup>8,16-18</sup>. Some model-based analyses have used the concept of critical soil water potential<sup>1,19</sup>, but current land surface models are using soil moisture rather than soil water potential, and global observation-based analyses of critical thresholds are still missing.

Satellite observations of surface SM with frequent revisit and global coverage based on microwave sensors in the L-band with stronger penetration capacity<sup>20</sup> can be combined with land surface temperature (LST) to assess the relationships between SM and the surface energy partitioning<sup>21,22</sup> (Methods). Instead of global satellite evapotranspiration products based on models with uncertain parametrizations<sup>23,24</sup>, we used here the diurnal evolution of LST as a direct observable signature of shifts in surface energy partitioning regimes<sup>25,26</sup>. Specifically, the land-surface temperature diurnal amplitude (dLST) starts to increase below $\theta_{crit}$ when ecosystems plunge into

the water-limited regime<sup>26-29</sup>. An increased dLST, for a given amount of net radiation, is directly linked to a decrease in EF and thus increased SM stress<sup>3</sup>. dLST is positively associated with sensible heating but negatively associated with EF and SM<sup>27-29</sup>. Evaporative regimes have been characterized with observed dLST-SM relationships across Africa, but not yet globally<sup>14</sup>, leaving a gap in our understanding of $\theta_{crit}$ across the globe.

To quantify the global spatial distribution of $\theta_{crit}$, we selected extensive periods without rainfall known as SM dry-downs<sup>14,30,31</sup> when the transition from energy to water limitation is likely to happen. The validity of the dLST-SM approach to determine $\theta_{crit}$ was demonstrated by comparing its results to the classical EF-SM method<sup>2,5,9,10</sup> at sites of the global network of flux tower measurements (Methods). Three global satellite SM datasets (SMAP-IB, SCA-V and SMOS-IC) and two LST datasets (Copernicus and MODIS) covering the period from April 2015 to December 2020 were then used to produce a global map of $\theta_{crit}$. Uncertainties were estimated based on an ensemble of 18 members from different pairings of SM and LST datasets, including the uncertainty on $\theta_{crit}$ from the dLST-SM relationship (Methods). Additionally, explainable machine learning models (random forest) were applied to gain insights on the climatic, biotic and edaphic factors
controlling the spatial variations of θcrit. Based on our global map of θcrit, we further calculated the fraction of days in a year when SM is below θcrit using time series of SM from satellite data and the ERA5-Land reanalysis³², to investigate the long-term trends of plant exposure to water stress over the last 40 years. Finally, we evaluated how land surface models of Earth System Models participating in the Coupled Model Inter-comparison Project Phase 6 (CMIP6) simulate the patterns of θcrit, and discussed their biases compared to our observation-based maps.

## Results and discussion

### Consistency of θcrit derived from the EF-SM and dLST-SM methods

Using soil dry-downs observed at 44 flux tower sites, we calculated daily dLST using in-situ daily maximum and minimum outgoing long-wave radiation (Methods) and compared θcrit defined as the break-point when dLST increases with decreasing SM, with the value calculated from the flux data as a tipping point of the EF-SM relationships, as previously done in refs. 2,5,9,10. An example is shown for the Hainich beech forest site in Germany (Supplementary Table 1) where the EF-SM relationship during dry-downs defines a θcrit of 0.192 ± 0.009 m³/m³ (± standard error) and the dLST-SM relationship gives a very similar estimate of 0.191 ± 0.005 m³/m³ (Fig. 1a & b). Across all the sites spanning a large range of aridity and plant functional types, the two approaches show consistent results (r = 0.87, Fig. 1c), in line with previous theoretical and observational studies⁸,¹⁴,²⁵.

### Global distribution of θcrit

Global high frequency LST and SM observations from multiple satellites during drydowns are then used to quantify the spatial distribution of θcrit (Methods). To calculate daily dLST, we used the MODIS Terra and Aqua satellites and the Copernicus dataset based on a constellation of geostationary satellites. Even though MODIS only passes over the earth four times a day while the geostationary satellites data (Copernicus) have 24 observations per day, allowing us to define more accurately the diurnal amplitude of LST, both observations match well with each other (Supplementary Figs. 1–2). Over Siberia and India where no geostationary data are available, we used only MODIS. For quasi-daily SM, we used satellite all-weather data from SMAP-IB, SCA-V and SMOS-IC, which show a similar pattern of the number of dry-downs per year over each point of the Earth (Supplementary Figs. 3–4). Here, the SM drydowns were defined as periods with at least five (SMAP-IB and SCA-V) or four (SMOS-IC) consecutive overpass measurements over intervals longer than 10 days during which SM is persistently decreasing (Methods). The areas with the largest number of dry-downs are in central America, Argentina, central Europe, eastern Europe and eastern Australia (Supplementary Fig. 3). On the other hand, only few dry-downs could be used to infer θcrit in wet regions such as the Amazon, central Africa and southern China. Grid points with no clear dry-downs to calculate θcrit were masked.

The global maps of θcrit obtained with the three SM and the two dLST satellite datasets show consistent patterns (Fig. 2a–f). Different maps from paired SM and dLST observations with three θcrit estimates (mean, and plus or minus one standard error, see Methods) provide an ensemble of 18 members. The median and standard error across all members of the ensemble shown in Fig. 2g represent our best estimate of the global distribution θcrit and its uncertainty. The relative uncertainty of θcrit, defined as the ratio of standard error to the median value of the 18 ensemble members, is less than 10% over most areas (Fig. 2h). Moreover, despite a mismatch in spatial scales, the value of θcrit from satellites at a global resolution of 25 by 25 km, is significantly correlated to the local estimate calculated at point-scale flux tower measurements (Supplementary Fig. 5). The median value of θcrit over the global vegetated areas is 0.19 m³/m³. Even though satellites only probe surface SM, whereas plants may be sensitive to stress from surface and

rootzone moisture deficits, surface SM has been shown to be equally skillful for identifying evapotranspiration regime changes as deeper soil layers or rootzone SM measurements³³,³⁴. We further use SM data with different soil layers from ERA5-Land (Methods) and compare the θcrit values derived from ERA5-Land SM layer 1 (0–7 cm depth), layer 2 (7–28 cm) and layer 3 (28–100 cm). We found that surface θcrit is highly correlated with θcrit derived from deep soil layers (Supplementary Fig. 6), showing that θcrit obtained from surface SM can provide information deeper into the subsurface, consistent with the results of flux tower observations reported by both Dong, Akbar³³ and Fu, Ciais³⁵.

The lowest θcrit values were observed in dryland ecosystems over the western United States, western Argentina, eastern Brazil, South Africa, northwestern China and Australia (Fig. 2g). In those dryland regions, plant hydraulic features adapted to conditions when evapotative demand often exceed soil water supply, are likely to minimize θcrit through mechanisms of sustained SM extraction by roots and transport by xylem, even at low soil water potentials³⁶. Conversely, the highest θcrit values were found in humid ecosystems such as Indonesia, south-eastern China, south-eastern United States, and Uruguay (Fig. 2g). Differences of θcrit between biomes were found to be significant, with increasing θcrit from dry shrublands, grasslands, and savannas towards temperate, boreal and tropical forests (Supplementary Fig 7a). Similar patterns were found across climate types, with increasing θcrit from hyper-arid, arid, and semi-arid ecosystems towards humid ecosystems (Supplementary Fig. 8).

We performed a more detailed analysis of the θcrit differences between cropland types, based on the expectation that θcrit should be affected by the choice of cultivars and by management practices such as irrigation. We found that θcrit varied among different crop species (Supplementary Fig. 7b), with rice (mostly irrigated) having significantly higher values (0.28 m³/m³) than maize, wheat and potato (0.20 m³/m³, *p* < 0.05). Moreover, θcrit tended to increase with increasing irrigation (Supplementary Fig. 7c). We also tested the hypothesis that the areas of recent cropland expansion over drier marginal lands should be associated with a decrease of θcrit (Methods). Most new cropland expansion occurred over drier areas, such as in southern Sahel, central Highlands and Zambia over Africa, and the Cerrado and Chaco plains in South America³⁷. Over the ‘new’ cropland areas that were cultivated after 2003 according to the high resolution map of ref. 37, we verified that θcrit was lower on average than on established cropland areas (Supplementary Fig. 7d). More insights on regional patterns of θcrit and its impacts for yields could be gained based on regional management and cultivars information, which is beyond the scope of this study.

The spatial distribution of θcrit in this study aligns with previous findings in ecological theory regarding plant stress across various environments¹,¹⁴,¹⁹,³³,³⁸⁻⁴⁰. Land surface models often have a lower θcrit model parameter in arid biomes³⁸,⁴⁰,⁴¹. The map of ecosystem-scale isohydricity from remotely sensed observations showed that the anisohydric behavior is more common in arid ecosystems³⁹. By quantifying the soil water potential threshold, Bassiouni, Good¹ showed that water uptake strategies in arid locations are generally more drought resistant. Note that soil water potential is rarely measured in situ, and land surface models are using soil moisture rather than soil water potential. Different vegetation water stress in arid and humid ecosystems have also been recognized in many other studies, based on the ecosystem limitation index³⁸, the Land Surface Water Index⁴²,⁴³, and SM anomalies⁴⁰. However, these indicators are not direct measures of water stress. The θcrit values quantified in our study reflect the long-term adaptation of ecosystems to aridity regimes. θcrit is simple to define and is a direct measure of water stress, but θcrit remains not observed and our study allows to compare it across biomes. θcrit can also be used to quantify the time spent below θcrit and understand how recent climate trends have affected the exposure of ecosystems to water stress.
<figure><img src="image_6.png" /><figcaption>**Fig. 2** | The global distribution of critical soil moisture threshold ($\theta_{crit}$). The global distribution of estimated $\theta_{crit}$ using Copernicus land surface temperature diurnal amplitude (dLST) and soil moisture (SM) from SMAP-IB (a), SCA-V (b) or</figcaption></figure>

SMOS-IC (c). The global distribution of estimated $\theta_{crit}$ using MODIS dLST and SM from SMAP-IB (d), SCA-V (e) or SMOS-IC (f). The median $\theta_{crit}$ (g) and its relative uncertainty (h) across 18 ensemble members of $\theta_{crit}$ (Methods).

## Drivers of global variation in $\theta_{crit}$

To evaluate the possible mechanisms controlling the spatial variation of $\theta_{crit}$ (Fig. 2g), we excluded croplands and used random forest models (Methods) with 35 candidate factors, including soil properties, vegetation structure, plant hydraulic traits and climatic variables (Supplementary Table 2). Based on a recursive feature elimination algorithm (Methods), a subset of 11 most influential predictors were selected in the final ‘best’ model, which explains 74% of the global spatial variation in $\theta_{crit}$ (Fig. 3a). The aridity index, defined as the ratio of mean annual potential evapotranspiration to precipitation, was identified as the most important factor; followed by leaf area index (LAI) and the sand fraction of soil texture (Fig. 3a). This result is consistent with Bassiouni, Good¹, who evaluated the relation between critical soil water potential and aridity index based on a soil water balance model and an inverse modeling analysis. But our study rather focused on observation-based $\theta_{crit}$ and used a comprehensive set of environmental variables to identify the main drivers of global $\theta_{crit}$ variations. Partial dependence analysis further showed that $\theta_{crit}$ decreases with a higher aridity index (Fig. 3b) and sand fractions (Fig. 3d) but becomes higher at higher LAI (Fig. 3c). A lower aridity index reflects wetter climates where a higher $\theta_{crit}$ can be interpreted as an adaptation trait in view of the low risk of plants to be exposed to a water limited regime. We noted that below a leaf area index of about 2.5 m²/m², the $\theta_{crit}$ decreases; above that, further increases in LAI are less important (Fig. 3c). This suggests that low $\theta_{crit}$ in arid areas are also

related to an increasing fraction of soil exposure, highlighting the role of soil evaporation in arid areas. Thus, further evaluation and measurements of soil evaporation are needed in the future to better quantify the significance of soil evaporation in arid areas. A higher LAI being positively associated with $\theta_{crit}$ (Fig. 3c) further supports the interpretation that wetter ecosystems can sustain more leaves without compromising transpiration, given that SM rarely drops below $\theta_{crit}$ during the year. Recent studies have also shown that ecosystems with higher leaf area index have a more gradual stomatal closure in response to a SM decrease, which sustains photosynthesis in periods of low to moderate water stress³⁵. On the other hand, the negative response of $\theta_{crit}$ to the sand fraction is consistent with the fact that sandy soils have lower SM wilting points⁴⁴. Indeed, sand fraction regulates the dependence of water potential to SM and water potential is the primary driver of plant water stress⁴¹,⁴⁵. Sandy soils have a lower soil water content for the same critical soil water potential for plant stress⁴⁶,⁴⁷, which explains the negative dependence of $\theta_{crit}$ to sand fraction. Our results also show that a higher leaf nitrogen content is associated with a lower $\theta_{crit}$, consistent with the fact that plants with a higher leaf nitrogen content have a larger resistance to drought⁴⁸. We also find that $\theta_{crit}$ shows a positive dependence on precipitation frequency but a negative dependence on shortwave radiation (Fig. 3f). More frequent precipitation events⁴⁹⁻⁵¹ and lower shortwave radiation help reduce water stress, and thus appear to favor an adaptation towards higher $\theta_{crit}$.
<figure><img src="image_7.png" /><figcaption>**Fig. 3 | The factors influencing global variation of critical soil moisture threshold ($\theta_{crit}$).** a The importance of climatic, biotic and edaphic variables in controlling $\theta_{crit}$. Aridity index is defined as the ratio of mean annual potential evapotranspiration to precipitation. b–g Partial dependence plots of the top six</figcaption></figure>

predictors. The y-axis is SHAP value for corresponding predictor (X-axis). The partial dependence plots indicate the effects of individual variables on the response, without the influence of the other variables (Methods).

## Global distribution of the fraction of stressed days and its trend over 1979–2020

The global map of $\theta_{crit}$ (Fig. 2g) can also be used to understand how recent climate trends have affected the exposure of ecosystems to water stress. We calculated the fraction of stressed days (FSD) each year, selecting the days when SM is below $\theta_{crit}$ over each location of the globe (Methods). Combining SMAP-IB, SCA-V and SMOS-IC time series of SM during 2016–2020, we find that dryland ecosystems have a yearly FSD higher than 70% (Fig. 4a). The same analysis performed with daily surface SM and $\theta_{crit}$ from the ERA5-Land reanalysis during 2016–2020 confirms the high yearly FSD in dryland ecosystems and gives a similar spatial pattern, but a lower mean FSD (Fig. 4b). This is because the $\theta_{crit}$ estimated from ERA5-Land data is larger than that the satellite observations (Supplementary Fig. 9).

After removing the pixels with large land cover changes (>10%) to limit the impacts of land changes on trend analysis (Methods), using daily ERA5-Land SM since 1979, we found that the FSD has been increasing significantly over the last 40 years (Fig. 4c, d), implying that terrestrial ecosystems became exposed to more extensive periods of water stress. Over the past four decades, FSD increased globally, on average, by about one day per year (Fig. 4d). In addition to increased evaporative demand and atmospheric drivers$^{52}$, this increasing trend of FSD may also be attributed to increased frequency of drought and heatwaves, resulting into an overall decline of SM$^{53,54}$. This result is in line with recent findings from Jiao, Wang$^{55}$ and Denissen, Teuling$^{38}$ based on independent data, suggesting a regime shift from energy to water limitation in relation to an overall decline of SM.

We acknowledge that $\theta_{crit}$ may change over time. Based on model outputs analysis, Hsu and Dirmeyer$^{56}$ found significant temporal variations in $\theta_{crit}$ across many locations spanning 100 years. Conversely, another study analyzed the temporal dynamics of $\theta_{crit}$ at five flux tower sites with at least 15 years of measurements and found no significant trend over time$^{35}$. This underscores the need for future research to gain a better understanding of the temporal dynamics of $\theta_{crit}$ through longer observations. We considered here that the temporal dynamics of $\theta_{crit}$ should not hamper our trend analysis, given that even if $\theta_{crit}$ changes, its magnitude over 40 years is minimal.

## Comparison with Earth System models

Finally we diagnosed $\theta_{crit}$ by using daily EF and surface SM simulations from Earth System Models (Supplementary Table 3, Methods). We found that the models showed less spatial variability of $\theta_{crit}$ than in the observation-based map (Fig. 5a, Supplementary Figs. 11–12, Fig. 2g) and significantly underestimated $\theta_{crit}$ in wet regions (Fig. 5b, Supplementary Fig. 11), suggesting that they may underestimate the soil moisture point of inception of plant water stress in wet regions. Such a bias may lead to overly optimistic projections of the future increase of plant CO$_2$ uptake. Conversely, models significantly overestimated $\theta_{crit}$ in dry regions and failed to capture the observed very low $\theta_{crit}$ values in arid areas (Fig. 5b, Supplementary Fig. 11), which could partly explain why ESMs underestimate both gross and net CO$_2$ fluxes in dryland ecosystems$^{57,58}$. Plants growing in arid areas have evolved many adaptation strategies to survive drought, for example by reducing leaf area index, reducing plant hydraulic and stomatal conductance, and using water stored in vegetation for transpiration$^{59–61}$. These mechanisms have not been properly parameterized or fully integrated in models. Thus, our results can help to guide the research directions that can improve the simulation of SM stress.

We noted that the biases in $\theta_{crit}$ should not be directly equated with model accuracy in simulating water stress because the ability of models to simulate water stress not only depend on the value of $\theta_{crit}$ but also on their simulation of water uptake and transport when SM is lower than $\theta_{crit}$. For example, models reduce gas exchange at different rates when ecosystem becomes water-limited$^{62,63}$. This leads to differences in their water and carbon simulations$^{64}$ and better $\theta_{crit}$ estimates will not resolve these differences that drive much of the water stress impacts on gross primary productivity and evapotranspiration. However, quantifying the inception of water stress – the $\theta_{crit}$, as done here, is a prerequisite for understanding the response rates of gas exchanges to SM stress. In addition, observation based models of evapotranspiration and gross primary productivity (e.g., light use efficiency models) typically assume fixed plant functional type values$^{65,66}$ to define SM stress thresholds, that are used across regions and climate. This study provides spatially explicit parameterizations of plant water stress as a function of envirometal drivers that could be
<figure><img src="image_8.png" /></figure>

**Fig. 4** | The global distribution and trend in the fraction of time when soil moisture is below the critical soil moisture threshold ($\theta_{crit}$). a. b The global distribution of the fraction of time when soil moisture (SM) is below $\theta_{crit}$ using satellite observed SM (median values from SMAP-IB, SCA-V and SMOS-IC) and ERA5-Land reanalysis SM during 2016–2020. c Annual time series of the fractions of time when SM is below $\theta_{crit}$ in regions with different fraction bins over 1979–2020. The trend (Sen's slope) and its 95% confidence interval are detected using the

nonparametric trend test technique (Mann–Kendall test; *p* < 0.05). The solid line shows the median value while the shading bounds the interval of the 25th to 75th percentiles. d Spatial patterns of the temporal trend in the fraction of time when SM is below $\theta_{crit}$ with white indicating those areas with no significant changes (Mann–Kendall test; *p* > 0.05) or >10% land cover changes during 1982–2016¹⁰⁰, and colored pixels indicating areas with significant trends (*p* < 0.05).

<figure><img src="image_9.png" /></figure>

**Fig. 5** | The global distribution of critical soil moisture threshold ($\theta_{crit}$) from Earth System Models. a The multi-model mean $\theta_{crit}$ using ten Earth System Models. b The differences (multi-model mean $\theta_{crit}$ minus observation-based $\theta_{crit}$) between multi-model mean $\theta_{crit}$ and observation-based $\theta_{crit}$.

incorporated in future model iterations to improve the representation of plant water stress and its spatial variations.

Vegetation regulates the terrestrial water and carbon cycles, as it controls and adapts to changing SM availability, with $\theta_{crit}$ being a key variable characterizing the coupling between soil-plant continuum and the atmosphere. Yet our ability to characterize $\theta_{crit}$ at the global scale has been limited to date. Based on the dLST-SM relationship from multiple satellite observations, this study provides the geographical

distribution and assessment of the variations of $\theta_{crit}$ across the globe. We also showed the usefulness of hourly LST data from geostationary satellites to understand ecosystem water stress⁶⁷,⁶⁸. By uncovering the relationships between $\theta_{crit}$ and environmental factors, including climatic, biotic and edaphic variables, we further added mechanistic understanding of the environmental factors driving the global variation in $\theta_{crit}$. It sheds light on potential tipping points of water stress impairing ecosystem functioning, and future research will aim to use
the understanding built based on the map of θcrit and its environmental drivers to improve land-surface model representation of SM constraints on water and carbon cycles. In addition, we showed that the terrestrial ecosystems experienced more frequent water-stress regimes through the past four decades, with important implications on the current land carbon sink capacity. Although the relative patterns of θcrit in some models are similar to those observed, the ten state-of-the-art ESMs that we tested failed to accurately reproduce the magnitude and spatial variability of θcrit, suggesting the uncertain projection of current and future response of carbon uptake and evapotranspiration to droughts. These ESMs need to improve the simulation of SM and related processes, conforming to measurements, to provide more reliable projections of the response of terrestrial ecosystems to climate change and feedbacks between land and atmosphere. Together, these results demonstrated the global distribution of θcrit and its drivers, applications and models' performance, with important implications for understanding the inception of water stress in models and identifying tipping points of water stress that could result in widespread impairment of ecosystem functioning and loss of ecosystem services with continued climate warming.

## Methods

### Eddy covariance measurements

We used half-hourly SM, latent heat flux, sensible heat flux, and outgoing longwave radiation from the recently released ICOS (Integrated Carbon Observation System)⁶⁹, AmeriFlux⁷⁰,⁷¹ and FLUXNET2015 datasets of energy, water, and carbon fluxes and meteorological data, all of which have undergone a standardized set of quality control and gap filling⁷²,⁷³. Data were processed following a consistent and uniform processing pipeline⁷². There were 279 flux tower sites in total by combining ICOS, AmeriFlux and FLUXNET2015 datasets. We first removed 130 sites without SM or outgoing longwave radiation measurements; then dropped all wetland sites because they have a perched water table and infrequently show SM limitations. Since for some sites, there is no dry-downdetected during the peak growing season across all available years; these sites were also excluded (81 sites remaining). The evaporative fraction (EF)-SM and land surface temperature diurnal amplitude (dLST)-SM relationships in these 81 sites were evaluated site-by-site, respectively, to detect the θcrit for each site (see below). There were 44 sites with the θcrit estimates for both EF-SM and dLST-SM methods. We only used the surface SM observations because surface SM (0-10 cm, varying across sites) was measured at all sites.

### Derivation of dLST from eddy covariance measurements

At each flux tower site, we derived daily dLST using measured daily maximum and minimum outgoing longwave radiation. The outgoing longwave radiation (LW) is emitted by the surface and depends on radiometric surface temperature (LST), the Stefan-Boltzmann constant (σ) and emissivity (ε) according to the Stefan-Boltzmann law⁴ (Eq. 1). Therefore, the dLST can be calculated as Eq. 2, where LWmax and LWmin are the daily maximum and minimum outgoing longwave radiation, respectively; ε is considered as constant at the same site and same day because we are deriving dLST (not LST).

$$LW = \varepsilon \times \sigma \times LST^4 \quad (1)$$

$$dLST = ((LW_{max}/\sigma)^{1/4} - (LW_{min}/\sigma)^{1/4})/\varepsilon^{1/4} \quad (2)$$

### SM and dLST from satellite observations

We used three L-band passive daily surface SM (to a depth of 5 cm) products: Soil Moisture Active Passive (SMAP)-INRAE-BORDEAUX (SMAP-IB)²⁰, single channel vertical polarization (SCA-V, SMAP_L3_SM_P)⁷⁵ and Soil Moisture and Ocean Salinity in version IC

(SMOS-IC)⁷⁶. Both SMAP-IB (version 1) and SCA-V (version 7, L3 products) have 36 km resolution and one to three-day revisit from 1 April 2015 to 31 December 2020. The SMAP-IB algorithm is based on the two-parameter inversion of the L-MEB model, as defined in Wigneron, Jackson⁷⁷, applied to the SMAP mono-angular dual-polarized brightness temperature¹⁸. SCA-V is not independent from SMAP-IB, but their retrieval algorithms, vegetation correction and surface roughness correction are different. SCA-V is adopted as the operational baseline algorithm to estimate SM from SMAP brightness temperature⁷⁸. In the SCA-V, vegetation is accounted for by the τ-ω model as in L-MEB. However, optical depth at nadir (τNAD) is not retrieved as for SMOS-IC and SMAP-IB. Instead it is estimated from the linear relation τNAD = b × VWC between τNAD and vegetation water content (VWC)⁷⁹. Thereby, values of the b-parameter are assumed polarization independent and will be provided from a land cover look up table, and the VWC is estimated from values of the NDVI Index. SMOS-IC is derived from the two-parameter L-MEB inversion applied to the SMOS multi-angular and dual-polarized brightness temperatures. SMOS-IC has 25 km resolution and two to four-day revisit from 1 January 2011 to 31 December 2020⁷⁶. Based on the recent study from Li, Wigneron²⁰, the biases of these three SM datasets were corrected using ISMN in-situ measurements, an international cooperation to construct and maintain a global in-situ SM database⁸⁰,⁸¹. Across all ISMN in-situ measurements, Li, Wigneron²⁰ found that the biases of SMAP-IB, SCA-V and SMOS-IC are 0.002, 0.008 and -0.054 m³/m³, respectively. We thus corrected the biases of these three SM datasets by subtracting the corresponding bias.

Two land surface temperature datasets from the Copernicus Global Land Operations and Moderate Resolution Imaging Spectroradiometer (MODIS) were used. The Copernicus LST (version 2) datasets are obtained from a constellation of geostationary satellite missions: Meteosat Second Generation (MSG) and Indian Ocean Data Coverage (IODC) missions, Geostationary Operational Environmental Satellite (GOES) and Himawari (and its predecessor Multi-Function Transport Satellite - MTSAT)⁸²,⁸³. The Copernicus LST provides hourly data at a spatial resolution of 5 km covering most of the globe's land surface, but there is no geostationary coverage in parts of northern and eastern Europe, Central Asia, and the Indian subcontinent as well as parts of eastern Siberia and northern North America (Supplementary Fig. 1). The second LST datasets are from the Terra (MOD11C1) and Aqua (MYD11C1) MODIS Version 6.1 Land Surface Temperature, providing four observed LST per day (10:30 AM/PM, 1:30 AM/PM) at a 0.05-degree resolution. We calculated the daily dLST as the difference between daily maximum and minimum LST using hourly LST from Copernicus or four observed LST every day from MODIS Terra and Aqua. The bilinear interpolation algorithm was applied to resample all data into the grid resolution of 0.25 degree.

### Soil moisture dry-down identification

Dry-downs following rainfall are episodes with no rain for several consecutive days during which SM shows a short term 'pulse' rise after rain and then decays until the next rain event. At each flux tower site, a dry-down is retained for our analysis when SM decreases consecutively for at least 10 days after rainfall following previous studies³⁰,³¹,⁸⁴⁻⁸⁶. The results were similar after requiring the soil dry-down to be at least 9 or 11 days. To ensure the reliability of latent heat flux measurements (high signal-to-noise ratio), we focused on the soil dry-downs during the peak growing season for all available site-years, defined as three-month periods with the maximum mean gross primary productivity across the available years. For satellite observations, soil drydowns were defined as at least 5 (for SMAP-IB and SCA-V, one to three-day revisit) or 4 (for SMOS-IC, two to four-day revisit) consecutive overpasses (over ≥ 10 days) of decreasing SM. The full year data of satellite observations were used and results were found to be similar when only growing season data were used.
Soil dry-down periods provided the unique and consistent opportunity for us to detect the θcrit because the transition from energy limitation to soil water limitation is likely to happen during the dry-down period. However, less dry-downs are available in wet regions than in dry regions. As a results, some wet regions including a small number of dry-downs were masked as it was not possible to detect the θcrit in such situations. Note that, if some wet regions have experienced more frequent droughts recently, they will thus have more dry-downs and have been included in our analysis. This approach partly reduces the spatial coverage of our global map of θcrit, but highlights that some wet regions should be further investigated in the future with longer satellite observations.

## θcrit estimation using EF-SM and dLST-SM methods

While other factors limit evapotranspiration besides SM and the linear dependency is a simple approximation, many previous studies have showed that the EF-SM and dLST-SM framework provides a good first-order representation of regimes of land-atmosphere coupling, both in models and observations (e.g., Seneviratne, Corti⁵, Seneviratne, Lüthi⁸⁷, Koster, Dirmeyer⁸⁸, Koster, Suarez⁸⁹, Teuling, Seneviratne⁹⁰, Feldman, Short Gianotti¹⁴). Here we provided a global analysis based on this first-order theoretical and empirically verified framework. We calculated the daily EF as the ratio of observed latent heat flux to the sum of latent and sensible heat fluxes. Then, we characterized the EF-SM and dLST-SM relationship at each site, respectively, using all available soil dry-downs, from a regression between these two variables with a linear-plus-plateau model:

$$EF \text{ or } dLST = \begin{cases} a + b(SM - \theta_{crit}) & \text{if } SM < \theta_{crit} \\ a & \text{if } SM \ge \theta_{crit} \end{cases} \quad (3)$$

where a is the maximum (or minimum) value of EF (or dLST) in the absence of SM stress (energy-limited stage), b represents the slope of the linear phase (water-limited stage) between EF (or dLST) and SM, and θcrit is the critical SM threshold. θcrit represents the breakpoint until which EF (or dLST) increases (or decreases) linearly as a function of SM. The θcrit and its standard error were simultaneously estimated by least squares fit with the R software package 'segmented'⁹¹ for each site. An example to estimate the θcrit using EF-SM and dLST-SM methods is shown in Fig. 1a, b. Following Feldman, Short Gianotti¹⁴, we considered three models: SM varying only within a water-limited regime (linear model) or energy-limited regime (linear model), and SM varying within a transitional regime (linear-plus-plateau model). According to the lowest Akaike Information Criterion⁹², we selected the “best” model pixel by pixel, and the θcrit is detected when the linear-plus-plateau model is selected. Some pixels or sites did not have a defined θcrit value if there were either no dry-downs or if SM varied only within a water- or energy-limited regime, thus rendering the breakpoint analysis of dLST-SM or EF-SM impossible. Based on the EF-SM and dLST-SM relationships, there were 44 sites (Supplementary Table 1) with the θcrit estimates for both EF-SM and dLST-SM methods. The Pearson correlation and its associated statistical test were used to compare the θcrit values from the dLST-SM method with that of EF-SM method. Increasing dLST is a direct observable signature of shifts in the surface energy partitioning regimes²⁵,²⁶. An increased diurnal temperature range, for a given amount of net radiation, is directly linked to a decrease in EF and thus increased soil moisture stress³. dLST is positively associated with sensible heating but negatively associated with EF and SM²⁷⁻²⁹. Evaporative regimes and θcrit estimating have been characterized previously with observed dLST-SM relationships across some regions, such as Africa¹⁴ and site level²⁸, showing that the dLST-SM relationship is an effective method to estimate θcrit. Here we applied this method to the global scale using multiple satellite observations.

For global satellite observations, we quantified θcrit for each pixel based on dLST-SM method using all dry-downs from 1 April 2015 to 31

December 2020. There were 18 maps of θcrit in total by considering all possible combinations of three SM datasets (SMAP-IB, SCA-V and SMOS-IC) and two dLST datasets (Copernicus and MODIS) and the uncertainty of θcrit estimates (θcrit and θcrit ± standard error, from the linear-plus-plateau model), resulting from different data sources and estimating variants. The median θcrit and its relative uncertainty at each pixel were calculated across 18 ensemble members. The relative uncertainty was defined as the ratio of standard error to the median value of these ensemble members. The θcrit estimated from satellite ensembles was compared with θcrit estimated from flux tower sites using the dLST-SM method or the EF-SM method, respectively. For each site, we extracted and calculated the median θcrit values within a 3 × 3 pixel window around the site from satellites-derived θcrit map. The Pearson correlation and its associated statistical test were used to compare θcrit based on satellite observations and flux towers across 26 sites. For the remaining 18 sites, the satellite data could not be used to derive θcrit because there were either no dry-downs, or SM varied only within a water- or energy-limited regime, or the number of samples were too low, thus rendering the breakpoint analysis of dLST-SM unreliable. Note that both the measurement time periods and frequency of flux tower sites differ from those of satellites. The results show that satellite-derived θcrit is significantly correlated to that estimated independently from eddy covariance measurements (Supplementary Fig. 5). The correlations between satellite-derived θcrit and tower-derived θcrit using both the dLST-SM method and the EF-SM method are strongly significant (p < 0.01), but their correlation coefficients are not very high (r = 0.57 for the dLST-SM method and r = 0.55 for the EF-SM method), which may be due to several factors. First, the footprint size ranges from a few meters to dozens of meters for the flux tower measurements but reaches 25 kilometers for satellite observations (0.25 degree). This mismatch is expected to lead to difference of θcrit values estimated from flux towers and satellites. Second, the soil depths of measured SM from flux towers and the quality of SM measurements varied among different sites, and the depths from flux towers are also different from those of satellite observations. This could also contribute to the differences found between θcrit values. Third, daily data from both flux towers and satellites were used, and high variability and measurement errors affect the data at this short time scale. Moreover, there are 48 measurements per day for flux towers but only a few revisits per week for satellite SM. These factors could introduce biases when comparing their θcrit values. We noted that the θcrit estimated from satellites is a bit higher than that of flux towers in the low θcrit range (Supplementary Fig. 5), which may be attributed to higher SM values from satellite data compared to measurements from flux towers in arid regions because of different sampling depths between flux tower measurements and satellite observations.

While satellite-based θcrit used surface SM, both Dong, Akbar³³ and Fu, Ciais³⁵ revealed that these thresholds also provide information deeper into the subsurface, and proved that surface and rootzone SM are often similarly skillful for identifying evapotranspiration regime changes based on in-situ observations. Feldman, Short Gianotti³⁴ recently also reported that remotely sensed surface SM can capture deep water dynamics relevant to plant water uptake so that L-band satellite SM data used here are relevant to vegetation rootzones. To further test whether the θcrit obtained from surface SM also provide information deeper into the subsurface at global scale, we used dLST and SM data with different soil layers from ERA5-Land and compared the θcrit values derived from ERA5-Land SM layer 1 (0–7 cm depth) with the layers 2 (7–28 cm) or 3 (28–100 cm).

## θcrit among different biomes and the impacts of farm management on θcrit in croplands

To compare θcrit values among different biomes, the International Geosphere-Biosphere Program (IGBP) classification from MCD12C1
and Köppen climate classification map were used (Supplementary Fig. 7). We also used the aridity classification by the United Nations Environment Program⁹³ (Supplementary Table 2), and the global landmass is classified into five categories, namely, (i) hyperarid, (ii) arid, (iii) semi-arid, (iv) dry sub-humid, and (v) humid (Supplementary Fig. 8). We performed a more detailed analysis of the $θ_{crit}$ differences between cropland types, based on the expectation that $θ_{crit}$ should be affected by the choice of cultivars and by management practices such as irrigation. The geographic distribution of main staple crops was from Monfreda, Ramankutty⁹⁴. The Global Map of Irrigation Areas (Version 5) was downloaded from the website of The Food and Agriculture Organization⁹⁵. This map showed the amount of area equipped for irrigation in percentage of the total area on a raster. We also tested the hypothesis that the areas of recent cropland expansion over drier marginal lands should be associated with a decrease of $θ_{crit}$, as more crop species adapted to dry environments would be selected. For cropland expansion, we used the map of percent of cropland net gain per pixel during 2003–2019 from Potapov, Turubanova³⁷. Differences in $θ_{crit}$ between groups (different biomes or climate types or farm managements in croplands) were analyzed using the Kruskal–Wallis test, a nonparametric test of difference⁹⁶. A *p* < 0.05 was used to identify significant differences between groups.

## Drivers of global variation in $θ_{crit}$

A random forest analysis was used to identify the factors (soil property, vegetation structure, plant hydraulic traits and climate – 35 factors in total) that contribute the most to the geographic variation in $θ_{crit}$. These variables were chosen due to their relevance to soil and vegetation dynamics based on field studies, sites observations and their availability at the global scale. $θ_{crit}$ is a composite attribute, reflecting soil and vegetation attributes. We investigated the importance of 35 factors (Supplementary Table 2) that reflect the dual role of soils and vegetation in determining $θ_{crit}$⁹. Predictor variables with low predictive power were removed from the random forest models to avoid overfitting. Following Green, Ballantyne⁹⁷, we first ran a random forest model with all predictor variables included, and the predictor variables were ranked according to their permutation importance. The model was then rerun with the least important variable removed from the model, a process called recursive feature elimination (RFE)⁹⁸. Importance values were then recalculated and stored, and this process was repeated until the three most important predictor variables remained. From here, the R-squared value was tabulated based on the out-of-bag observations (-one-third of the observations), and then, the model was rerun with the next most important variable added back in (based on the importance rankings stored during RFE). The R-squared value of this model based on the out-of-bag observations was then retabulated, and should the R-squared value increase by at least 0.005, the predictor variable remained in the model (otherwise, it was removed) and the next most important variable was then added back into the model and was rerun with a new R-squared value tabulated. This process was repeated until all predictor variables had been added back into the model, and the variable combination with the highest out-of-bag R-squared value was selected for the final model. Additionally, for each model, the number of variables used at each node split (between 2 and the number of predictor variables, with a final selection of 4) and the number of trees used in the model (between 50 and 5000, with a final selection of 1500) were optimized to maximize out-of-bag R-square value. In this way, the best quality model could be developed by only including the most informative inputs. The final set of predictors included the following 11 predictor variables: sand fraction and volumetric fraction of coarse fragments to describe soil properties; leaf area index, leaf nitrogen concentration, woody density and specific leaf area to describe vegetation structures; plant hydraulic resistance to describe plant hydraulic traits; aridity index (defined as the ratio of

mean annual potential evapotranspiration to precipitation⁹⁹), mean annual precipitation frequency, shortwave radiation and vapor pressure deficit to describe climatic factors.

The selected random forest model was used to calculate Shapley values (SHAP), and thus analyze the sensitivity of the output to the input variables, and improve upon feature importance. Shapley values, based on game theory, assess every combination of predictors to determine each predictor's impact. For example, focusing on the aridity index feature, the approach tests the accuracy of every combination of features not including aridity index and then tests how adding aridity index improves the accuracy on each combination. Thus, the results from the partial SHAP dependency analysis can be used to determine the effects of individual variables on the response, without the influence of other variables⁹⁷.

## Calculating the fraction of stressed days and its trend

To explore how many days in a year that ecosystem are water-limited, we calculated the fraction of stressed days (FSD), defined as the ratio of the number of days with SM < $θ_{crit}$ to the total observed daily SM in a year for each pixel. The fraction of time when SM is below $θ_{crit}$ was computed for each year during 2016–2020 using SMAP-IB, SCA-V and SMOS-IC, respectively. Then the median value across these three SM datasets was calculated. The same analysis was also performed for the daily ERA5-Land reanalysis SM dataset because it has longer time series (1979–2020). Following satellite datasets analysis, the $θ_{crit}$ for ERA5-Land was estimated using ERA5-Land SM and dLST. We first compared the FSD from ERA5-Land reanalysis SM during 2016–2020 with that of satellite observations. Then we calculated the FSD from ERA5-Land reanalysis SM for each year from 1979 to 2020. The overall trends (Sen's slope) of the FSD in regions with mean fractions of times spent below $θ_{crit}$ within 10% to 30%, 30% to 50%, 50% to 70% and 70% to 90%, were detected, respectively, using the nonparametric trend test technique (Mann–Kendall test). To avoid the impacts of extreme values, we did not include the regions with mean fractions below 10% or above 90%. Additionally, we also used the Mann–Kendall test to evaluate trend (Sen's slope) in the fraction of time when SM is below $θ_{crit}$ at each pixel and map its trend globally. A *p* < 0.05 was used to identify statistically significant trends. We noted that there are some data gaps in $θ_{crit}$ derived from ERA5-Land datasets (Supplementary Fig. 9a) because of the failure to fit a breakpoint model, suggesting that there are some inconsistencies in SM or dLST data between ERA5-Land and satellites. The large spread of the scatters between ERA5-Land derived and satellite-derived $θ_{crit}$ (Supplementary Fig. 9b) partly indicates such a discrepancy. We thus compared the daily ERA5-Land SM and dLST with those from satellites for a day in 2020, and found that the primary biases between ERA5-Land and satellites lie in SM data, rather than dLST data (Supplementary Fig. 10).

To remove the impacts of land cover changes on the trend analysis, we masked the pixels with >10% land cover changes during 1982–2016 according to the global land changes data from Song, Hansen¹⁰⁰ based on daily satellite observations acquired by the Advanced Very High Resolution Radiometer. Song, Hansen¹⁰⁰ quantified the global land changes during the period 1982–2016 and developed an annual vegetation continuous fields (VCF, representing the land surface as a fractional combination of vegetation functional types) product consisting of tree canopy cover, short vegetation cover and bare ground cover and characterized land change over the past 35 years (0.05-degree resolution). The bilinear interpolation algorithm was applied to resample this data into the grid resolution of 0.25 degree. We also considered that the temporal dynamics of $θ_{crit}$ should not hamper the trend analysis because Fu, Ciais³⁵ have analyzed the temporal dynamics of $θ_{crit}$ at five flux tower sites with at least 15 years of measurements, and found no significant trend with time in $θ_{crit}$.
## CMIP6 ESM simulations

Ten ESMs (ACCESS-ESM1-5, BCC-ESM1, Can-ESM5, CMCC-CM2, INM-CM5, IPSL-CM6A, MIROC6, MPI-ESM1-2-HR, MRI-ESM2 and NorESM2-MM) in CMIP6 provided daily surface SM, latent and sensible heat fluxes outputs (Supplementary Table 3). Daily SM and calculated EF (from latent and sensible heat fluxes) from historical runs (2009–2014) were used for each model. Following the observational analysis, the same analysis was carried out for the ten CMIP6 models. For each model, we first selected all soil dry-downs from the full-year dataset of model outputs, defined as at least 10 consecutive days of decreasing SM, then quantified θcrit pixel-by-pixel by means of the EF–SM relationship. Multi-model mean θcrit was calculated for each pixel by averaging the θcrit across these ten models. To evaluate the θcrit performance in ESMs, we also calculated the difference between multi-model mean θcrit and observation-based θcrit. We noted that different models led to different simulated SM values¹⁰, and this inherent divergence of simulated SM distribution could also contribute to the differences between observation-based θcrit and models-based θcrit values. But we found that all models consistently showed less spatial variability of θcrit than in the observation-based map, suggesting that our result did not depend on the inherent divergence of simulated SM distribution.

## Data availability

The global critical soil moisture thresholds of plant water stress are available at https://zenodo.org/records/11183719. The eddy covariance measurements are downloaded from the ICOS (https://doi.org/10.18160/2G60-ZHAK), AmeriFlux (https://ameriflux.lbl.gov/) and FLUXNET2015 datasets (https://fluxnet.fluxdata.org/data/fluxnet2015-dataset/). SMAP-IB and SMOS-IC SM are obtained from https://ibremote-sensing.inrae.fr/. SCA-V are available on National Snow and Ice Data Center (https://smap.jpl.nasa.gov/data/). Copernicus LST are downloaded from https://land.copernicus.eu. MODIS LST are from https://lpdaac.usgs.gov/. The Global Map of Irrigation Areas (Version 5) was downloaded from the website of The Food and Agriculture Organization (https://www.fao.org/aquastat/en/geospatial-information/global-maps-irrigated-areas/latest-version). ERA5-Land reanalysis data are from https://cds.climate.copernicus.eu/. The CMIP6 data are downloaded from https://esgf-data.dkrz.de/search/cmip6-dkrz/.

## Code availability

The primary code used to generate the results is publicly available at https://zenodo.org/records/11183719.

## References

1. Bassiouni, M., Good, S. P., Still, C. J. & Higgins, C. W. Plant water uptake thresholds inferred from satellite soil moisture. *Geophys. Res. Lett.* **47**, e2020GL087077 (2020).

2. Rodriguez-Iturbe, I. Ecosystems: a hydrologic perspective of climate-soil-vegetation dynamics. *Water Resour. Res.* **36**, 3–9 (2000).

3. Gentine, P., Chhang, A., Rigden, A. & Salvucci, G. Evaporation estimates using weather station data and boundary layer theory. *Geophys. Res. Lett.* **43**, 661–611,670 (2016).

4. Gentine, P. et al. Coupling between the terrestrial carbon and water cycles—a review. *Environ. Res. Lett.* **14**, 083003 (2019).

5. Seneviratne, S. I. et al. Investigating soil moisture–climate interactions in a changing climate: a review. *Earth-Sci. Rev.* **99**, 125–161 (2010).

6. Santanello, J. A. Jr. et al. Land–atmosphere interactions: the LoCo perspective. *Bull. Am. Meteorol. Soc.* **99**, 1253–1272 (2018).

7. Zhang, P. et al. Abrupt shift to hotter and drier climate over inner East Asia beyond the tipping point. *Science* **370**, 1095–1099 (2020).

8. Schwinghackl, C., Hirschi, M. & Seneviratne, S. I. Quantifying spatiotemporal variations of soil moisture control on surface energy balance and near-surface air temperature. *J. Clim.* **30**, 7105–7124 (2017).

9. Rodríguez-Iturbe I. & Porporato A. Ecosystems: soil moisture and plant dynamics (Cambridge University Press, 2007).

10. Laio, F., Porporato, A., Ridolfi, L. & Rodriguez-Iturbe, I. Plants in water-controlled ecosystems: active role in hydrologic processes and response to water stress: II. Probabilistic soil moisture dynamics. *Adv. water Resour.* **24**, 707–723 (2001).

11. Budyko, M. I. *Climate and life* (Academic press, 1974).

12. Eagleson, P. S. Climate, soil, and vegetation: 4. The expected value of annual evapotranspiration. *Water Resour. Res.* **14**, 731–739 (1978).

13. Baldocchi, D. D., Xu, L. & Kiang, N. How plant functional-type, weather, seasonal drought, and soil physical properties alter water and energy fluxes of an oak–grass savanna and an annual grassland. *Agric. For. Meteorol.* **123**, 13–39 (2004).

14. Feldman, A. F., Short Gianotti, D. J., Trigo, I. F., Salvucci, G. D. & Entekhabi, D. Satellite-based assessment of land surface energy partitioning–soil moisture relationships and effects of confounding variables. *Water Resour. Res.* **55**, 10657–10677 (2019).

15. Koster, R., Schubert, S. & Suarez, M. Analyzing the concurrence of meteorological droughts and warm periods, with implications for the determination of evaporative regime. *J. Clim.* **22**, 3331–3341 (2009).

16. Dirmeyer, P. A., Koster, R. D. & Guo, Z. Do global models properly represent the feedback between land and atmosphere? *J. Hydrometeorol.* **7**, 1177–1198 (2006).

17. Green, J. K. et al. Large influence of soil moisture on long-term terrestrial carbon uptake. *Nature* **565**, 476–479 (2019).

18. Short Gianotti, D. J., Rigden, A. J., Salvucci, G. D. & Entekhabi, D. Satellite and station observations demonstrate water availability’s effect on continental-scale evaporative and photosynthetic land surface dynamics. *Water Resour. Res.* **55**, 540–554 (2019).

19. Bassiouni, M., Higgins, C. W., Still, C. J. & Good, S. P. Probabilistic inference of ecohydrological parameters using observations from point to satellite scales. *Hydrol. Earth Syst. Sci.* **22**, 3229–3243 (2018).

20. Li, X. et al. A new SMAP soil moisture and vegetation optical depth product (SMAP-IB): Algorithm, assessment and inter-comparison. *Remote Sens. Environ.* **271**, 112921 (2022).

21. Dong, J. & Crow, W. T. The added value of assimilating remotely sensed soil moisture for estimating summertime soil moisture-air temperature coupling strength. *Water Resour. Res.* **54**, 6072–6084 (2018).

22. Dong, J. & Crow, W. T. Land-band remote-sensing increases sampled levels of global soil moisture-air temperature coupling strength. *Remote Sens. Environ.* **220**, 51–58 (2019).

23. Badgley, G., Fisher, J. B., Jiménez, C., Tu, K. P. & Vinukollu, R. On uncertainty in global terrestrial evapotranspiration estimates from choice of input forcing datasets. *J. Hydrometeorol.* **16**, 1449–1455 (2015).

24. Bai, P. & Liu, X. Intercomparison and evaluation of three global high-resolution evapotranspiration products across China. *J. Hydrol.* **566**, 743–755 (2018).

25. Bateni, S. & Entekhabi, D. Relative efficiency of land surface energy balance components. *Water Resour. Res.* **48**, 1–8 (2012).

26. Panwar, A., Kleidon, A. & Renner, M. Do surface and air temperatures contain similar imprints of evaporative conditions? *Geophys. Res. Lett.* **46**, 3802–3809 (2019).

27. Betts, A. K., Desjardins, R., Worth, D. & Beckage, B. Climate coupling between temperature, humidity, precipitation, and cloud cover over the Canadian Prairies. *J. Geophys. Res.* **119**, 305–313,326 (2014).

28. Amano, E. & Salvucci, G. D. Detection and use of three signatures of soil-limited evaporation. *Remote Sens. Environ.* **67**, 108–122 (1999).
29. Thakur, G., Schymanski, S. J., Mallick, K., Trebs, I. & Sulis, M. Downwelling longwave radiation and sensible heat flux observations are critical for surface temperature and emissivity estimation from flux tower data. *Sci. Rep.* **12**, 1–14 (2022).

30. Akbar, R. et al. Estimation of landscape soil water losses from satellite observations of soil moisture. *J. Hydrometeorol.* **19**, 871–889 (2018).

31. Feldman, A. F. et al. Moisture pulse-reserve in the soil-plant continuum observed across biomes. *Nat. plants* **4**, 1026–1033 (2018).

32. Hersbach, H. et al. The ERA5 global reanalysis. *Q. J. R. Meteorol. Soc.* **146**, 1999–2049 (2020).

33. Dong, J., et al. Can surface soil moisture information identify evapotranspiration regime transitions? *Geophys. Res. Lett.* **49**, e2021GL097697 (2022).

34. Feldman, A. F. et al. Remotely sensed soil moisture can capture dynamics relevant to plant water uptake. *Water Resour. Res.* **59**, e2022WR033814 (2023).

35. Fu, Z. et al. Critical soil moisture thresholds of plant water stress in terrestrial ecosystems. *Sci. Adv.* **8**, eabq7827 (2022).

36. Grünzweig, J. M. et al. Dryland mechanisms could widely control ecosystem functioning in a drier and warmer world. *Nat. Ecol. Evol.* **6**, 1064–1076 (2022).

37. Potapov, P. et al. Global maps of cropland extent and change show accelerated cropland expansion in the twenty-first century. *Nat. Food* **3**, 19–28 (2022).

38. Denissen, J. et al. Widespread shift from ecosystem energy to water limitation with climate change. *Nat. Clim. Change* **12**, 677–684 (2022).

39. Konings, A. G. & Gentine, P. Global variations in ecosystem-scale isohydricity. *Glob. Change Biol.* **23**, 891–905 (2017).

40. Li, X. et al. Global variations in critical drought thresholds that impact vegetation. *Natl Sci. Rev.* **10**, nwad049 (2023).

41. Novick, K. A. et al. Confronting the water potential information gap. *Nat. Geosci.* **15**, 158–164 (2022).

42. Chandrasekar, K., Sesha Sai, M., Roy, P. & Dwevedi, R. Land Surface Water Index (LSWI) response to rainfall and NDVI using the MODIS Vegetation Index product. *Int. J. Remote Sens* **31**, 3987–4005 (2010).

43. Du, J. et al. Synergistic satellite assessment of global vegetation health in relation to ENSO-induced droughts and pluvials. *J. Geophys. Res. Biogeosci.* **126**, e2020JG006006 (2021).

44. Bonan, G. *Ecological climatology: concepts and applications* (Cambridge University Press, 2015).

45. Shao, Y. & Irannejad, P. On the choice of soil hydraulic models in land-surface schemes. *Bound. Layer. Meteorol.* **90**, 83–115 (1999).

46. Van Genuchten, M. T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. *Soil Sci. Soc. Am. J.* **44**, 892–898 (1980).

47. Mualem, Y. A new model for predicting the hydraulic conductivity of unsaturated porous media. *Water Resour. Res.* **12**, 513–522 (1976).

48. Prentice, I. C., Dong, N., Gleason, S. M., Maire, V. & Wright, I. J. Balancing the costs of carbon gain and water transport: testing a new theoretical framework for plant functional ecology. *Ecol. Lett.* **17**, 82–91 (2014).

49. Knapp, A. K. et al. Rainfall variability, carbon cycling, and plant species diversity in a mesic grassland. *Science* **298**, 2202–2205 (2002).

50. Zhang, W. et al. Ecosystem structural changes controlled by altered rainfall climatology in tropical savannas. *Nat. Commun.* **10**, 1–7 (2019).

51. Sloat, L. L. et al. Increasing importance of precipitation variability on global livestock grazing lands. *Nat. Clim. Change* **8**, 214–218 (2018).

52. Yuan, W. et al. Increased atmospheric vapor pressure deficit reduces global vegetation growth. *Sci. Adv.* **5**, eaax1396 (2019).

53. Li, W. et al. Widespread increasing vegetation sensitivity to soil moisture. *Nat. Commun.* **13**, 1–9 (2022).

54. Williams, A. P. et al. Large contribution from anthropogenic warming to an emerging North American megadrought. *Science* **368**, 314–318 (2020).

55. Jiao, W. et al. Observed increasing water constraint on vegetation growth over the last three decades. *Nat. Commun.* **12**, 1–9 (2021).

56. Hsu, H. & Dirmeyer, P. A. Uncertainty in projected critical soil moisture values in CMIP6 affects the interpretation of a more moisture-limited world. *Earth's Future* **11**, e2023EF003511 (2023).

57. MacBean, N. et al. Dynamic global vegetation models underestimate net CO2 flux mean and inter-annual variability in dryland ecosystems. *Environ. Res. Lett.* **16**, 094023 (2021).

58. Wang, L. et al. Dryland productivity under a changing climate. *Nat. Clim. Change* **12**, 981–994 (2022).

59. Smith, W. K. et al. Remote sensing of dryland ecosystem structure and function: Progress, challenges, and opportunities. *Remote Sens. Environ.* **233**, 111401 (2019).

60. Dahlin, K., Fisher, R. & Lawrence, P. Environmental drivers of drought deciduous phenology in the Community Land Model. *Biogeosciences* **12**, 5061–5074 (2015).

61. Zhou, S., Duursma, R. A., Medlyn, B. E., Kelly, J. W. & Prentice, I. C. How should we model plant responses to drought? An analysis of stomatal and non-stomatal responses to water stress. *Agric. For. Meteorol.* **182**, 204–214 (2013).

62. Medlyn, B. E. et al. Using models to guide field experiments: a priori predictions for the CO 2 response of a nutrient-and water-limited native Eucalypt woodland. *Glob. Change Biol.* **22**, 2834–2851 (2016).

63. Paschalis, A. et al. Rainfall manipulation experiments as simulated by terrestrial biosphere models: Where do we stand? *Glob. Change Biol.* **26**, 3336–3355 (2020).

64. Trugman, A., Medvigy, D., Mankin, J. & Anderegg, W. Soil moisture stress as a major driver of carbon cycle uncertainty. *Geophys. Res. Lett.* **45**, 6495–6503 (2018).

65. Kolby Smith, W. et al. Large divergence of satellite and Earth system model estimates of global terrestrial CO2 fertilization. *Nat. Clim. Change* **6**, 306–310 (2016).

66. Pei, Y. et al. Evolution of light use efficiency models: improvement, uncertainties, and implications. *Agric. For. Meteorol.* **317**, 108905 (2022).

67. Xiao, J., Fisher, J. B., Hashimoto, H., Ichii, K. & Parazoo, N. C. Emerging satellite observations for diurnal cycling of ecosystem processes. *Nat. Plants* **7**, 877–887 (2021).

68. Wen, J. et al. Resolve the clear-sky continuous diurnal cycle of high-resolution ECOSTRESS evapotranspiration and land surface temperature. *Water Resour. Res.* **58**, e2022WR032227 (2022).

69. Warm Winter 2020 Team, & ICOS Ecosystem Thematic Centre. Warm Winter 2020 ecosystem eddy covariance flux product for 73 stations in FLUXNET-Archive format—release 2022-1 (version 1.0). ICOS Carbon Portal https://doi.org/10.18160/2G60-ZHAK. (2022).

70. Novick, K. A. et al. The AmeriFlux network: a coalition of the willing. *Agric. For. Meteorol.* **249**, 444–456 (2018).

71. Lu, X. & Keenan, T. F. No evidence for a negative effect of growing season photosynthesis on leaf senescence timing. *Glob. Change Biol.* **28**, 3083–3093 (2022).

72. Pastorello, G. et al. The FLUXNET2015 dataset and the ONEflux processing pipeline for eddy covariance data. *Sci. Data* **7**, 1–27 (2020).

73. Baldocchi, D. et al. FLUXNET: a new tool to study the temporal and spatial variability of ecosystem-scale carbon dioxide, water vapor, and energy flux densities. *Bull. Am. Meteorol. Soc.* **82**, 2415–2434 (2001).

74. Luyssaert, S. et al. Land management and land-cover change have impacts of similar magnitude on surface temperature. *Nat. Clim. Change* **4**, 389–393 (2014).

75. Jackson, T. J. III Measuring surface soil moisture using passive microwave remote sensing. *Hydrol. Process.* **7**, 139–152 (1993).
76. Wigneron, J.-P. et al. SMOS-IC data record of soil moisture and L-VOD: historical development, applications and perspectives. *Remote Sens. Environ.* **254**, 112238 (2021).

77. Wigneron, J.-P. et al. Modelling the passive microwave signature from land surfaces: a review of recent results and application to the L-band SMOS & SMAP soil moisture retrieval algorithms. *Remote Sens. Environ.* **192**, 238–262 (2017).

78. Zeng, J., Chen, K.-S., Bi, H. & Chen, Q. A preliminary evaluation of the SMAP radiometer soil moisture product over United States and Europe using ground-based measurements. *IEEE Trans. Geosci. Remote Sens.* **54**, 4929–4940 (2016).

79. Jackson, T. & Schmugge, T. Vegetation effects on the microwave emission of soils. *Remote Sens. Environ.* **36**, 203–212 (1991).

80. Dorigo, W. et al. The International Soil Moisture Network: serving Earth system science for over a decade. *Hydrol. Earth Syst. Sci.* **25**, 5749–5804 (2021).

81. Dorigo, W. et al. Global automated quality control of in situ soil moisture data from the International Soil Moisture Network. *Vadose Zone J.* **12**, vzj2012-0097 (2013).

82. Freitas, S. C. et al. Land surface temperature from multiple geostationary satellites. *Int. J. Remote Sens.* **34**, 3051–3068 (2013).

83. Wan, Z. & Dozier, J. A generalized split-window algorithm for retrieving land-surface temperature from space. *IEEE Trans. Geosci. Remote Sens.* **34**, 892–905 (1996).

84. Fu, Z. et al. Uncovering the critical soil moisture thresholds of plant water stress for European ecosystems. *Glob. Change Biol.* **28**, 2111–2123 (2022).

85. McColl, K. A. et al. Global characterization of surface soil moisture drydowns. *Geophys. Res. Lett.* **44**, 3682–3690 (2017).

86. Shellito, P. J., Small, E. E. & Livneh, B. Controls on surface soil drying rates observed by SMAP and simulated by the Noah land surface model. *Hydrol. Earth Syst. Sci.* **22**, 1649–1663 (2018).

87. Seneviratne, S. I., Lüthi, D., Litschi, M. & Schär, C. Land–atmosphere coupling and climate change in Europe. *Nature* **443**, 205 (2006).

88. Koster, R. D. et al. Regions of strong coupling between soil moisture and precipitation. *Science* **305**, 1138–1140 (2004).

89. Koster, R. D. et al. Realistic initialization of land surface states: Impacts on subseasonal forecast skill. *J. Hydrometeorol.* **5**, 1049–1063 (2004).

90. Teuling, A., Seneviratne, S. I., Williams, C. & Troch, P. Observed timescales of evapotranspiration response to soil moisture. *Geo-phys. Res. Lett.* **33**, 1–5 (2006).

91. Muggeo, V. M. Segmented: an R package to fit regression models with broken-line relationships. *R. News* **8**, 20–25 (2008).

92. Akaike, H. A new look at the statistical model identification. *IEEE Trans. Autom. Control* **19**, 716–723 (1974).

93. Middleton, N. & Thomas, D. *World atlas of desertification*. ed. 2. Arnold, Hodder Headline, PLC (1997).

94. Monfreda, C., Ramankutty, N. & Foley, J. A. Farming the planet: 2. Geographic distribution of crop areas, yields, physiological types, and net primary production in the year 2000. *Glob. Biogeochem. Cycles* **22**, 1–10 (2008).

95. Siebert, S., Henrich, V., Frenken, K. & Burke, J. Global map of irrigation areas version 5. *Rheinische Friedrich-Wilhelms-Universität, Bonn, Germany/Food and Agriculture Organization of the United Nations*, Rome, Italy **2**, 1299–1327 (2013).

96. McKight, P. E. & Najab, J. Kruskal-wallis test. *The corsini encyclo-pedia of psychology*, 1-1 (2010).

97. Green, J. K. et al. Surface temperatures reveal patterns of vegetation water stress and their environmental drivers across the tropical Americas. *Glob. Change Biol.* (2022).

98. Guyon, I., Weston, J., Barnhill, S. & Vapnik, V. Gene selection for cancer classification using support vector machines. *Mach. Learn.* **46**, 389–422 (2002).

99. Zomer, R. J., Xu, J. & Trabucco, A. Version 3 of the global aridity index and potential evapotranspiration database. *Sci. Data* **9**, 1–15 (2022).

100. Song, X.-P. et al. Global land change from 1982 to 2016. *Nature* **560**, 639–643 (2018).

101. Koster, R. D. et al. On the nature of soil moisture in land surface models. *J. Clim.* **22**, 4322–4335 (2009).

## Acknowledgements

This work was financially supported by the IGSNRR (E3V30050), National Natural Science Foundation of China (31988102), CNES (5100019800) and the ANR CLAND Convergence Institute. I.C.P. acknowledges support by European Research Council funding under the European Union's Horizon 2020 research and innovation program (787203 REALM). P.G. and I.C.P. acknowledge support by the LEMONTREE (Land Ecosystem Models based On New Theory, observation and Experiments) project, funded through the generosity of Eric and Wendy Schmidt by recommendation of the Schmidt Futures program. A.F.F. was supported by an appointment to the NASA Postdoctoral Program at the NASA Goddard Space Flight Center, administered by Oak Ridge Associated Universities under contract with NASA. W.K.S. was supported by the NASA Carbon Cycle and Ecosystems Program under grant 8ONSSC23K0109. D.M. was supported by the INRAE metaprogrammes CLIMAE and XRISQUES. A.R.K. was supported by the College of Agricultural Sciences at Penn State University via USA Hatch Appropriations under Project PEN04710 and Accession 1020049. We would like to thank the ICOS Infrastructure for support in collecting and curating the eddy covariance data. This work used global eddy covariance data acquired and shared by the FLUXNET community, including these networks: AmeriFlux, AfriFlux, AsiaFlux, CarboAfrica, CarboEuropeIP, Carboltaly, CarboMont, ChinaFlux, Fluxnet-Canada, GreenGrass, ICOS, KoFlux, LBA, NECC, OzFlux-TERN, TCOS-Siberia and USCCC. The ERA-Interim reanalysis data are provided by ECMWF and processed by LSCE. The FLUXNET eddy covariance data processing and harmonization were carried out by the European Fluxes Database Cluster, AmeriFlux Management Project and Fluxdata project of FLUXNET, with the support of CDIAC and ICOS Ecosystem Thematic Center and the OzFlux, ChinaFlux and AsiaFlux offices.

## Author contributions

Z.F. and P.C. designed the study. Z.F. performed the analysis. Z.F. and P.C. wrote the paper with the inputs from all co-authors. J.P.W., P.G., A.F.F., D.M., N.V., A.R.K., D.S.G., I.C.P., P.C.S., D.Y., L.Y.L., H.L.M., X.J.L., Y.Y.H., K.L.Y., P.Z., X.L., Z.C.Z., J.H.L. and W.K.S. provided methodological suggestions and contributed to the interpretation of the results.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-024-49244-7.

**Correspondence** and requests for materials should be addressed to Zheng Fu.

**Peer review information** *Nature Communications* thanks Stephen Good, Yanghui Kang and Hsin Hsu for their contribution to the peer review of this work. A peer review file is available.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2024

¹Key Laboratory of Ecosystem Network Observation and Modeling, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing 100101, China. ²Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, Gif-sur-Yvette 91191, France. ³ISPA, INRAE, Université de Bordeaux, Bordeaux Sciences Agro, F-33140 Villenave d'Ornon, France. ⁴Department of Earth and Environmental Engineering, Columbia University, New York, NY 10027, USA. ⁵NASA Goddard Space Flight Center, Biospheric Sciences Laboratory, Greenbelt, MD 20771, USA. ⁶Earth System Science Interdisciplinary Center, University of Maryland, College Park, MD, USA. ⁷Unit Applied Mathematics and Computer Science (UMR MIA-PS) INRAE AgroParisTech Université Paris-Saclay, Palaiseau 91120, France. ⁸Department of Plant Science, The Pennsylvania State University, 116 Agricultural Science and Industries Building, University Park, PA 16802, USA. ⁹Department of Biological Systems Engineering, University of Wisconsin–Madison, Madison, USA. ¹⁰Georgina Mace Centre for the Living Planet, Department of Life Sciences, Imperial College London, Silwood Park Campus, Buckhurst Road, Ascot SL5 7PY, UK. ¹¹Ministry of Education Key Laboratory for Earth System Modeling, Department of Earth System Science, Tsinghua University, Beijing 100084, China. ¹²Earth and Planetary Sciences, Weizmann Institute of Science, Rehovot 7610001, Israel. ¹³INRAE, Avignon Université, UMR 1114 EMMAH, UMT CAPTE, F-84000 Avignon, France. ¹⁴Department of Geography, The University of Hong Kong, Hong Kong, SAR, China.

¹⁵Research Institute of Agriculture and Life Sciences, Seoul National University, Seoul, South Korea. ¹⁶Peking University Shenzhen Graduate School, Peking University, Shenzhen 518055 Guangdong, China. ¹⁷School of Natural Resources and the Environment, University of Arizona, Tucson, AZ, USA.

<figure><img src="image_10.png" /></figure>

e-mail: fuzheng@igsnrr.ac.cn