Global critical soil moisture thresholds of plant water stress

Domain: earth_science
Published: 2024-06-06
Editorial outcome: minor_revision

Abstract
During extensive periods without rain, known as dry-downs, decreasing soil moisture (SM) induces plant water stress at the point when it limits evapotranspiration, defining a critical SM threshold (θ crit). Better quantification of θ 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 θ crit. We find an average global θ crit of 0.19 m3/m3, varying from 0.12 m3/m3 in arid ecosystems to 0.26 m3/m3 in humid ecosystems. θ crit simulated by Earth System Models is overestimated in dry areas and underestimated in wet areas. The global observed pattern of θ 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 θ 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.

Introduction
The critical soil moisture threshold (θ crit) of plant water stress is defined as the soil moisture (SM) level at which evapotranspiration becomes SM limited in that environment1. Below this threshold, a marginal reduction of SM reduces evapotranspiration and increases sensible heat emissions and surface temperature2, making the air above the canopy warmer and drier, which in turn further reduces evapotranspiration and plant carbon dioxide uptake3–5. The control of energy partitioning regimes across θ crit determines local climate through land‐atmosphere coupling and can amplify warming during droughts6,7. A better knowledge of θ crit is thus important for land-atmosphere interactions5, for climate studies8–10 and for understanding the vulnerability of ecosystems and crop yields to drought9.

The relationship between SM and the evaporative fraction (EF), defined as the ratio of evapotranspiration to net radiation, show two distinct regimes2,5,11,12 (see Fig. 1a). When SM is higher than θ crit, the system is non-water limited (energy limited) and SM does not impact evapotranspiration5. In contrast, when SM is lower than θ 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)10. The overall EF–SM relationship (increasing below θ crit and then plateauing) is conceptually well established, but a spatially explicit understanding with accurate global maps of θ crit are lacking, due to a lack of global high-frequency observations of EF11,13–15. Therefore, the factors that control the global variations in θ crit are poorly known. Earth system modelers have adopted simple parametric representations of EF–SM relationship and θ crit to describe soil water stress and land-atmosphere feedbacks5, leading to model biases which hinder our ability to predict drought and its ecosystem impacts8,16–18. Some model-based analyses have used the concept of critical soil water potential1,19, but current land surface models are using soil moisture rather than soil 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 capacity20 can be combined with land surface temperature (LST) to assess the relationships between SM and the surface energy partitioning21,22 (Methods). Instead of global satellite evapotranspiration products based on models with uncertain parametrizations23,24, we used here the diurnal evolution of LST as a direct observable signature of shifts in surface energy partitioning regimes25,26. Specifically, the land-surface temperature diurnal amplitude (dLST) starts to increase below θ crit when ecosystems plunge into the water-limited regime26–29. An increased dLST, for a given amount of net radiation, is directly linked to a decrease in EF and thus increased SM stress3. dLST is positively associated with sensible heating but negatively associated with EF and SM27–29. Evaporative regimes have been characterized with observed dLST–SM relationships across Africa, but not yet globally14, leaving a gap in our understanding of θ crit across the globe.

To quantify the global spatial distribution of θ crit, we selected extensive periods without rainfall known as SM dry-downs14,30,31 when the transition from energy to water limitation is likely to happen. The validity of the dLST–SM approach to determine θ crit was demonstrated by comparing its results to the classical EF–SM method2,5,9,10 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 θ crit. Uncertainties were estimated based on an ensemble of 18 members from different pairings of SM and LST datasets, including the uncertainty on θ 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 reanalysis32, 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
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 longwave radiation (Methods) and compared θ crit defined as the breakpoint 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 m3/m3 (± standard error) and the dLST–SM relationship gives a very similar estimate of 0.191 ± 0.005 m3/m3 (Fig. 1a & b). Across all the sites spanning a large range of aridity and plant functionality types, the two approaches show consistent results (r = 0.87, Fig. 1c), in line with previous theoretical and observational studies8,14,25.

Global distribution of θ crit
Global high frequency LST and SM observations from multiple satellites during dry downs 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 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 m3/m3. Even though satellites only probe surface SM, whereas plants may be sensitive to stress from surface and root zone moisture deficits, surface SM has been shown to be equally skillful for identifying evapotranspiration regime changes as deeper soil layers or root zone SM measurements33,34. 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, Akbar33 and Fu, Ciais35.

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 evaporative demand often exceeds 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 potentials36. 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 toward 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 m3/m3) than maize, wheat and potato (0.20 m3/m3, 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 America37. 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 environments1,14,19,33,38–40. Land surface models often have a lower θ crit model parameter in arid biomes38,40,41. The map of ecosystem-scale isohydricity from remotely sensed observations showed that the anisohydric behavior is more common in arid ecosystems39. By quantifying the soil water potential threshold, Bassiouni, Good1 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 index38, the Land Surface Water Index42,43, and SM anomalies40. 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 allow 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.

Drivers of global variation in θ crit
To evaluate the possible mechanisms controlling the spatial variation of θ 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 θ 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, Good1, 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 θ crit and used a comprehensive set of environmental variables to identify the main drivers of global θ crit variations. Partial dependence analysis further showed that θ 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 θ 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 m2/m2, the θ crit decreases; above that, further increases in LAI are less important (Fig. 3c). This suggests that low θ 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 θ crit (Fig. 3c) further supports the interpretation that wetter ecosystems can sustain more leaves without compromising transpiration, given that SM rarely drops below θ 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 stress35. On the other hand, the negative response of θ crit to the sand fraction is consistent with the fact that sandy soils have lower SM wilting points44. Indeed, sand fraction regulates the dependence of water potential to SM and water potential is the primary driver of plant water stress41,45. Sandy soils have a lower soil water content for the same critical soil water potential for plant stress46,47, which explains the negative dependence of θ crit to sand fraction. Our results also show that a higher leaf nitrogen content is associated with a lower θ crit, consistent with the fact that plants with a higher leaf nitrogen content have a larger resistance to drought48. We also find that θ crit shows a positive dependence on precipitation frequency but a negative dependence on shortwave radiation (Fig. 3f). More frequent precipitation events49–51 and lower shortwave radiation help reduce water stress, and thus appear to favor an adaptation toward higher θ crit.

Global distribution of the fraction of stressed days and its trend over 1979–2020
The global map of θ 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 θ 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 θ crit from the ERA5-Land reanalysis during 2016–2020 confirm the high yearly FSD in dryland ecosystems and gives a similar spatial pattern, but a lower mean FSD (Fig. 4b). This is because the θ 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 drivers52, this increasing trend of FSD may also be attributed to increased frequency of drought and heatwaves, resulting into an overall decline of SM53,54. This result is in line with recent findings from Jiao, Wang55 and Denissen, Teuling38 based on independent data, suggesting a regime shift from energy to water limitation in relation to an overall decline of SM.

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

Comparison with Earth System models
Finally we diagnosed θ 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 θ crit than in the observation-based map (Fig. 5a, Supplementary Figs. 11–12, Fig. 2g) and significantly underestimated θ 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 CO2 uptake. Conversely, models significantly overestimated θ crit in dry regions and failed to capture the observed very low θ crit values in arid areas (Fig. 5b, Supplementary Fig. 11), which could partly explain why ESMs underestimate both gross and net CO2 fluxes in dryland ecosystems57,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 transpiration59–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 θ 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 θ crit but also on their simulation of water uptake and transport when SM is lower than θ crit. For example, models reduce gas exchange at different rates when ecosystem becomes water-limited62,63. This leads to differences in their water and carbon simulations64 and better θ 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 θ 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 values65,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 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 θ crit being a key variable characterizing the coupling between soil-plant continuum and the atmosphere. Yet our ability to characterize θ 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 θ crit across the globe. We also showed the usefulness of hourly LST data from geostationary satellites to understand ecosystem water stress67,68. By uncovering the relationships between θ crit and environmental factors, including climatic, biotic and edaphic variables, we further added mechanistic understanding of the environmental factors driving the global variation in θ 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)69, AmeriFlux70,71 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 filling72,73. Data were processed following a consistent and uniform processing pipeline72. 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-down detected 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 law74 (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 = ε × σ × LST4 (1)

dLST = ((LWmax/σ)1/4 - (LWmin/σ)1/4) / ε1/4 (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)20, single channel vertical polarization (SCA-V, SMAP_L3_SM_P)75 and Soil Moisture and Ocean Salinity inversion IC (SMOS-IC)76. 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, Jackson77, applied to the SMAP mono-angular dual-polarized brightness temperature18. 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 temperature78. 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)79. There-by, 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 202076. Based on the recent study from Li, Wigneron20, 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 database80,81. Across all ISMN in-situ measurements, Li, Wigneron20 found that the biases of SMAP-IB, SCA-V and SMOS-IC are 0.002, 0.008 and −0.054 m3/m3, 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)82,83. 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 studies30,31,84–86. 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 highlight 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, Corti5, Seneviratne, Lüthi87, Koster, Dirmeyer88, Koster, Suarez89, Teuling, Seneviratne90, Feldman, Short Gianotti14). 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 or dLST = a + b(SM - θ crit) if SM < θ crit

a if SM ≥ θ crit

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’91 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 Gianotti14, 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 Criterion92, 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 regimes25,26. 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 stress3. dLST is positively associated with sensible heating but negatively associated with EF and SM27–29. Evaporative regimes and θ crit estimating have been characterized previously with observed dLST–SM relationships across some regions, such as Africa14 and site level28, 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 θ 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, Akbar33 and Fu, Ciais35 revealed that these thresholds also provide information deeper into the subsurface, and proved that surface and root zone SM are often similarly skillful for identifying evapotranspiration regime changes based on in-situ observations. Feldman, Short Gianotti34 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 root zones. 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 Program93 (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, Ramankutty94. The Global Map of Irrigation Areas (Version 5) was downloaded from the website of The Food and Agriculture Organization95. 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, Turubanova37. 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 difference96. 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 θ crit9. Predictor variables with low predictive power were removed from the random forest model to avoid overfitting. Following Green, Ballantyne97, 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)98. 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 precipitation99), 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 variables97.

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 time 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, Hansen100 based on daily satellite observations acquired by the Advanced Very High Resolution Radiometer. Song, Hansen100 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, Ciais35 have analyzed the temporal dynamics of θ crit at five flux tower sites with at least 15 year 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 values101, 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.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors present an analysis of the ‘critical soil moisture threshold’ (theta_crit) based on remotely sensed and model data and validated with tower data. The theta_crit value is where soil moisture becomes limiting and evapotranspiration (ET) fluxes being to decrease. Understanding the value of theta_crit can help in creating and parameterizing models of vegetation stress and its effects on ET. There have be numerous recent studies of theta_crit, including a similar study the authors have previously conducted with global sets of tower-based observation (see ref 6 in particular). The main addition here appears to be doing a similar thing with satellite data as well as comparing with theta_crit derived from model outputs and evaluating trends in theta_crit through time.\nOverall, the study is well written and the subject matter is of interest to a wide audience. I suggest a few interconnected issues which should be addressed.\n1) Given the number of prior studies in this subject area (e.g. refs 1, 6, 7, 16, and others cited) the authors do not do a sufficient job of placing the finding here within the context of these earlier works. These prior studies also examined theta_crit, or similar quantities, and so how is what is learned here different, similar, and building upon this prior body of work. We’re similar patterns observed with respect to aridity, LAI, veg type, etc.? This is lacking in the manuscript.\n2) In a similar vein, little is present about what hydrologic theory would suggest these values are. For instance plant wilting points and other soil moisture thresholds are often defined with respect to the soil water potential. Wilting points are typically specified as being the same wilting point across textures, with the volumetric water content different from this based on soil retention curves. While I see that the authors have soil textures, they should be able to translate their observations into water potentials which would represent more physical properties driving the movement of water.\n3) Similarly ecological theory about plant stress in various environments has a long history. What does this body of literature say we should expect in more (or less) arid environments or different ecosystems. Which types of ecosystems are expected to shut down quicker or later? This has been approached in a variety of ways including LSM model parameterization, remote sensing of isohydricity, and other studies on vegetation stress thresholds, and so on. What new information about these trends/patterns/relationships by the study the authors present.\nAll the above points are interconnected and represent a lack of placement of this studies results within the context of prior work. Given that these theta_crit values are not directly to be used in models (if they are then the authors need to compare against what values models are currently using) a better effort needs to be made to demonstrate the value of the work done here. This can be achieved by properly building off of prior understanding.", "reviewer_2": "This paper presents a global assessment of the critical soil moisture thresholds using satellite-based soil moisture and land surface temperature (LST) estimates, unveiling global patterns of plant hydraulic strategies adapted to atmospheric dryness and soil water availability. A global map of critical soil moisture thresholds was derived by analyzing the relationship between diurnal LST range and soil moisture during dry-down periods (days after rainfall). The approach is validated with in-situ data from global flux towers (44 sites). The study showed that Earth System Models underestimate the spatial variability of the critical soil moisture threshold, indicating potential model deficiencies. Lastly, the paper revealed an increasing trend in the fraction of days with soil moisture below the critical threshold during the past four decades based on ERA5-Land reanalysis data.\nThis study identifies the critical soil water threshold, its global patterns, and the environmental factors, which provide important insights into ecosystem water-stress regimes and soil moisture constraints on global water cycles. The paper is well-written, and the results are clearly and effectively presented. However, I have a few questions and comments regarding the potential confounding factors for critical soil moisture threshold determination and the analysis of ERA5-Land data, that I hope the authors could clarify or consider.\nFirst, the analysis could be enhanced by a thorough evaluation of the confounding variables affecting the dLST-SM and EF-SM relationships, thereby, the temporal variation of critical soil moisture threshold. In this study, a single soil moisture threshold was derived for each pixel/site, yet it is possible that such a threshold present temporal variation linked to atmosphere and vegetation changes. As suggested by Feldman et al. (2019, ref 16), VPD consistently mediates the threshold inferred from the dLST-SM relationship in African grasslands. Therefore, it would be beneficial to include confounding variables in the statistical analysis and understand whether and how they may affect spatial patterns of critical soil moisture threshold.\nSecondly, the extensive data gaps in the soil moisture threshold derived from ERA5-Land data presented, as seen in Fig. 4 and Fig. 9, require a more detailed explanation. The extent of missing data, which is much higher than satellite-based analysis, is quite surprising since ERA5-Land is a reanalysis product. I’m interested in understanding the main cause of data gaps: lack of dry-down periods, limited range of SM, or missing data. If the main cause is the failure of the breakpoint model, it may suggest that there are major inconsistencies in SM or LST data between ERA5-Land and satellite estimates. The large spread of the scatters between ERA5-Land derived and satellite-derived theta_crit in Supplementary Fig. 9b partly indicates such a discrepancy. An explanation and discussion of these inconsistencies would be helpful to strengthen the analysis.\nAnother minor suggestion is to include a simple plot of theta_crit between satellite-derived v.s. ESM-model estimates, perhaps summarized by biome or climate types. This could provide a more straightforward illustration of how ESMs underestimate the threshold in wet regions and overestimate the threshold in water-limited regions than the maps shown in Fig. 5.\nMinor comments:\nMain text\nL164-165: This is an interesting hypothesis. A more detailed elaboration on how recently expanded cropland is associated with lower critical soil moisture than established ones would be helpful.\nFig. 4: More context on the ERA5_Land analysis would be helpful for clarity. Did the ERA5-Land analysis use the theta_crit derived from ERA5-Land? Which relationship was used, dLST-SM or ET-SM?\nL263: While the ESMs’ theta_crit estimates show a much smaller spatial variability, the relative patterns are quite similar to that of the satellite-derived values, at least in some models.\nMethod\nL295: Does daily minimum LW occur during nighttime?\nL302: To avoid confusion, consider using the official product name for “SCA-V” (e.g., SMAP_L3_XX). I was initially confused that “SCA-V” is soil moisture retrieved from another satellite, before reading the method second in detail.\nL316: It would be helpful to provide the temporal availability of the SMOS-IC dataset.\nL389: Could you please double-check the reference here? After a brief read of Reference 7, I didn’t find relevant results on dLST-SM relationship, but maybe I missed it. In any case, a brief summary of previous work on dLST-SM for theta_crit estimating would be helpful.\nL411: I am curious if there are any potential explanations about the consistently higher estimates from satellite versus flux tower in the low theta_crit range.\nL417: Could you please elaborate on the third point about how the high temporal variations could lead to uncertainties in the estimates?\nL467: It would be helpful to have the final selected hyperparameters as a reference for other and future studies.", "reviewer_3": "Fu et al. introduce a novel approach for estimating the critical soil moisture threshold and the fraction of water-limited days, deviating from methods previously limited to reanalysis and climate models. The global map generated in this study is of significant value, establishing a foundational reference for future comparisons. Given my own focus on soil moisture regime and critical values estimation, I find the figures in this study particularly exciting. The writing is logical, and the methods are clearly articulated. While I have only a few comments and suggestions, I recommend acceptance of the paper after minor revisions. Hsin Hsu (Wish not to remain anonymous)\n1.The paper shows that the Aridity Index is a primary determinant of spatial variation in θcrit. This also suggests that annual variability in Aridity Index can lead to annual variability in θcrit. However, the annual FSD is calculated based on climatological θcrit rather than year-specific values. This approach potentially compromises the fidelity of the FSD trend estimation. Although the authors reference a study (Line 507) indicating limited changes in θcrit, another research (Hsu and Dirmeyer, 2023) suggests significant changes in θcrit over many locations in one century. The FSD examination of this study (40 year- reanalysis-global) is between the two mentioned extreme cases (15 year-observation-few sites vs. 100 year with +1%CO2 per year in climate model). Given that the period used in all of these studies are within similar increasing rate of C02, the evidence of 15-year observation analysis might not effectively support authors’ statement. The authors could consider a decadal θcrit estimation from ERA5-land as supplementary information and/or moving the argument of temporal variation of θcrit from method to main text with a more comprehensive discussion.\nHsu and Dirmeyer, 2023. Uncertainty in Projected Critical Soil Moisture Values in CMIP6 Affects the Interpretation of a More Moisture-Limited World. Doi: 10.1029/2023EF003511\n2.Maybe I miss this somewhere: The paper lacks clarity regarding the sampling method for drydown used to estimate θcrit for climate models. It is unclear whether the author selects identical months as observations or the model's peak growing seasons. Discussing the implications of this choice is essential, as differences in background climate (and the aridity index) between observations and model outputs may lead to variations in estimated θcrit.\n3.Regarding the comment #3 on divergence of hydroclimate among data products. how does author think the inherent divergence of simulated soil moisture distribution (Koster et al. 2009) affects the interpretation here. Does this extends beyond model-specific vegetation/soil dynamics as discussed in lines 231-236?\nKoster et al. 2009: On the Nature of Soil Moisture in Land Surface Models. Doi: 10.1175/2009JCLI2832.1\n4.While the method using dLST appears straightforward, the main text could benefit from a physical description of the connection between using Evaporative Fraction (EF) and dLST, explaining why this approach is effective (Line 84-85). This addition would enhance the audience's understanding of the methodology.\nLine 77. typo: atmosphere"}, "author_response": "We greatly appreciate the opportunity to revise our manuscript, and would like to thank all reviewers for their valuable expertise and constructive comments. We have carefully revised the manuscript following the reviewers’ suggestions. Consequently, our manuscript has been considerably improved. Please see below our point-by-point responses in blue text following reviewer comments. The line numbers referred to are for the clean version of the revised manuscript (non-track-change version).\n---------------------------------------------------------------------------------------------------------------------\nReviewer #1 (Remarks to the Author):\n1.1 The authors present an analysis of the ‘critical soil moisture threshold’ (theta_crit) based on remotely sensed and model data and validated with tower data. The theta_crit value is where soil moisture becomes limiting and evapotranspiration (ET) fluxes being to decrease. Understanding the value of theta_crit can help in creating and parameterizing models of vegetation stress and its effects on ET.\nThere have be numerous recent studies of theta_crit, including a similar study the authors have previously conducted with global sets of tower-based observation (see ref 6 in particular). The main addition here appears to be doing a similar thing with satellite data as well as comparing with theta_crit derived from model outputs and evaluating trends in theta_crit through time.\nOverall, the study is well written and the subject matter is of interest to a wide audience. I suggest a few interconnected issues which should be addressed.\nThank you very much for supporting this paper and recognizing the significance of our study for a wide audience. We have carefully revised the manuscript following your suggestions. Consequently, our manuscript has been considerably improved. Please see below our detailed responses.\n1.2 1) Given the number of prior studies in this subject area (e.g. refs 1, 6, 7, 16, and others cited) the authors do not do a sufficient job of placing the finding here within the context of these earlier works. These prior studies also examined theta_crit, or similar quantities, and so how is what is learned here different, similar, and building upon this prior body of work. We’re similar patterns observed with respect to aridity, LAI, veg type, etc.? This is lacking in the manuscript.\nThank you very much for this constructive comment. We agree that prior studies investigated this subject area. Bassiouni, Good 1 examined the critical soil water potential threshold based on a soil water balance model and an inverse modeling framework. Several previous observational studies also reported θ at flux tower sites2 or regional level3. Although these previous works contributed to improve our knowledge on θ at sites level and highlighted the importance of scaling observed θ globally, they did not provide global maps of θ , due to a lack of global high-frequency observations of EF, leaving a gap in our understanding of θ at the global scale.\nPlease also note that current land surface models are using soil moisture rather than soil water potential.\nGoing beyond previous studies, we expended an approach initially developed over Africa to generate the first observation-based global map of θ by combining systematic satellite observations and in-situ data from flux towers. The approach to map θ is based on the diurnal amplitude of the land surface temperature and daily soil moisture measurements during periods of consecutive days without rain, using in-situ flux tower networks and satellite observations from different sensors, including geostationary satellites to measure the diurnal amplitude of land surface temperature. This observation-based map of θ allowed us to evaluate the possible mechanisms controlling the global variation of θ , and to diagnose simulated θ from current Earth System Models. Our study thus evaluated the global distribution of observation-based θ and its drivers. We also show that models differ greatly from observations, and that models have spatially too smooth gradients of θ , underestimating θ in drylands and overestimating it in wet biomes. The novelties of this study include:\n1. A global comparison of two approaches to estimate θ values: Evaporative Fraction decreases during soil moisture dry-downs, which only applicable at in situ flux tower sites, and land surface temperature diurnal amplitude increases during soil moisture dry-downs, applicable globally from satellite observations;\n2. The first observation-based global map of θ , derived by combining satellite observations and in-situ data from flux towers.\n3. Exploration of the drivers of global θ variations using a very comprehensive set of drivers and explainable machine learning;\n4. Uncovering the distribution of the fraction of stressed days when soil moisture stays below θ and its trend over 1979-2020;\n5. Applied the framework to CMIP6 models and diagnosed the θ from current Earth System Models simulations;\n6. Brought together multiple lines of evidence to better understand θ and its implications.\nFollowing your suggestions, we have added more discussion of previous works and highlighted our findings and novelties in the revised manuscript.\nLine 79-82:\n“Some model-based analyses have used the concept of critical soil water potential1, 4, 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.”\nLine 202-205:\n“This result is consistent with Bassiouni, Good 1, 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 θ and used a comprehensive set of environmental variables to identify the main drivers of global θ variations.”\nLine 215-217:\n“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 stress2.”\nLine 92-93:\n“Evaporative regimes have been characterized with observed dLST–SM relationships across Africa, but not yet globally3, leaving a gap in our understanding of θ across the globe.”\nLine 150-155:\n“We further use SM data with different soil layers from ERA5-Land (Methods) and compare the θ 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 θ is highly correlated with θ derived from deep soil layers (Supplementary Fig. 6), showing that θ obtained from surface SM can provide information deeper into the subsurface, consistent with the results of flux tower observations reported by both Dong, Akbar 5 and Fu, Ciais 2.”\nLine 179-192:\n“The spatial distribution of θ in this study aligns with previous findings in ecological theory regarding plant stress across various environments 1, 3, 4, 5, 9, 10, 11. Land surface models often have a lower θ model parameter in arid biomes8, 9, 11. The map of ecosystem-scale isohydricity from remotely sensed observations showed that the anisohydric behavior is more common in arid ecosystems10. By quantifying the soil water potential threshold, Bassiouni, Good 1 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 index9, the Land Surface Water Index12, 13, and SM anomalies11. However, these indicators are not direct measures of water stress. The θ values quantified in our study reflect the long-term adaptation of ecosystems to aridity regimes. θ is simple to define and is a direct measure of water stress, but θ remains not observed and our study allows to compare it across biomes. θ can also be used to quantify the time spent below θ and understand how recent climate trends have affected the exposure of ecosystems to water stress.”\nLine 284-304:\n“Based on the dLST–SM relationship from multiple satellite observations, this study provides the geographical distribution and assessment of the variations of θ across the globe. We also showed the usefulness of hourly LST data from geostationary satellites to understand ecosystem water stress6, 7. By uncovering the relationships between θ and environmental factors, including climatic, biotic and edaphic variables, we further added mechanistic understanding of the environmental factors driving the global variation in θ . It sheds light on potential tipping points of water stress impairing ecosystem functioning, and future research will aim to use the new understanding built based on the map of θ 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. The ten state-of-the-art ESMs that we tested failed to accurately reproduce the magnitude of θ , 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 θ 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.”\n1.3 2) In a similar vein, little is present about what hydrologic theory would suggest these values are. For instance plant wilting points and other soil moisture thresholds are often defined with respect to the soil water potential. Wilting points are typically specified as being the same wilting point across textures, with the volumetric water content different from this based on soil retention curves. While I see that the authors have soil textures, they should be able to translate their observations into water potentials which would represent more physical properties driving the movement of water.\nThank you for this critical comment. We fully agree that soil water potential represents better the physical properties in driving the movement of water than soil moisture. Following your suggestions, we tried to converted our map of θ into a map of critical soil water potential (y ) using soil pedotransfer functions, which showed the map of y is not realistic, with more negative y values in the driest regions and less negative ones in wet ecosystems (Figure R1).\nThe bias is likely due to the uncertainty in soil pedotransfer functions and parameters. It is well known that the parameter distributions in pedotransfer functions are poorly constrained and prevent confident transformation of soil moisture to soil water potential. For example, even relatively small variations in a single parameter of the van Genuchten model cause soil water potential to vary by an order of magnitude over a wide range of soil moisture8. Please look at the ref 8 for more details about the large uncertainty in converting soil moisture to soil water potential using soil pedotransfer functions and parameters. Given the large uncertainties in soil pedotransfer functions and parameters, we choose not to add the map of critical soil water potential into the manuscript.\nFigure R1. Estimated soil matric potential threshold by converting the θ based on soil pedotransfer function and parameters.\n1.4 3) Similarly ecological theory about plant stress in various environments has a long history. What does this body of literature say we should expect in more (or less) arid environments or different ecosystems. Which types of ecosystems are expected to shut down quicker or later? This has been approached in a variety of ways including LSM model parameterization, remote sensing of isohydricity, and other studies on vegetation stress thresholds, and so on. What new information about these trends/patterns/relationships by the study the authors present.\nThank you for this valuable suggestion. The spatial distribution of θ in this study is in agreement with previous studies in ecological theory about plant stress in various environments1, 3, 4, 5, 9, 10, 11. Land surface models often have a lower θ model parameter in arid biomes8, 9, 11. The map of ecosystem-scale isohydricity derived from remotely sensed observations showed that the anisohydric behavior is more common in arid ecosystems10. Note that this plant trait is not a direct indicator of water stress. By quantifying the soil water potential threshold, Bassiouni, Good 1 showed that water uptake strategies in arid locations are generally more drought resistant. But 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 index9, the Land Surface Water Index12, 13, and SM anomalies11. However, these indicators are not direct measures of water stress.\nQuantifying the inception of water stress – the θ , as done in our study, is a prerequisite for understanding the response rates of gas exchanges to SM stress. θ is simple to define and is a direct measure of water stress (Evaporative Fraction drops), but θ remains not observed and our study allows to compare it across biomes. The θ values quantified in our study reflect the long-term adaptation of ecosystems to aridity regimes. Ecosystems with low θ often have greater adaptive capacity to deal with water stress because low θ reflects a resistance to soil dryness. In addition, explainable machine learning models (random forest) were applied in our study to gain insights on the climatic, biotic and edaphic factors controlling the spatial variations of θ . Based on our map of θ , we can also calculate the fraction of stressed days (FSD, Fig. 4 in the manuscript) to explore how many days in a year that ecosystem are water-limited.\nObservation based models of evapotranspiration and gross primary productivity (e.g., light use efficiency models) typically assume fixed plant functional type values14, 15 to define SM stress thresholds, that are used across regions and climate. Our study provides spatially explicit parameterizations of plant water stress as a function of envirometal drivers that could be incorporated in future model iterations to improve the representation of plant water stress and its spatial variations.\nIn the revised manuscript, we now added this discussion as suggested.\nLine 179-192:\n“The spatial distribution of θ in this study aligns with previous findings in ecological theory regarding plant stress across various environments 1, 3, 4, 5, 9, 10, 11. Land surface models often have a lower θ model parameter in arid biomes8, 9, 11. The map of ecosystem-scale isohydricity from remotely sensed observations showed that the anisohydric behavior is more common in arid ecosystems10. By quantifying the soil water potential threshold, Bassiouni, Good 1 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 index9, the Land Surface Water Index12, 13, and SM anomalies11. However, these indicators are not direct measures of water stress. The θ values quantified in our study reflect the long-term adaptation of ecosystems to aridity regimes. θ is simple to define and is a direct measure of water stress, but θ remains not observed and our study allows to compare it across biomes. θ can also be used to quantify the time spent below θ and understand how recent climate trends have affected the exposure of ecosystems to water stress.”\n1.5 All the above points are interconnected and represent a lack of placement of this studies results within the context of prior work. Given that these theta_crit values are not directly to be used in models (if they are then the authors need to compare against what values models are currently using) a better effort needs to be made to demonstrate the value of the work done here. This can be achieved by properly building off of prior understanding.\nThank you very much for these constructive comments. In this study, we are trying our best and making the most of current knowledges to generate the first observation-based global map of θ by combining systematic satellite observations and in-situ data from flux towers. We have followed your suggestions and carefully revised the manuscript as above. As a consequence, our manuscript has been considerably improved. Please see our detailed responses above. We hope that our responses and revision of the manuscript are satisfactory to you.\n2.1 This paper presents a global assessment of the critical soil moisture thresholds using satellite-based soil moisture and land surface temperature (LST) estimates, unveiling global patterns of plant hydraulic strategies adapted to atmospheric dryness and soil water availability. A global map of critical soil moisture thresholds was derived by analyzing the relationship between diurnal LST range and soil moisture during dry-down periods (days after rainfall). The approach is validated with in-situ data from global flux towers (44 sites). The study showed that Earth System Models underestimate the spatial variability of the critical soil moisture threshold, indicating potential model deficiencies. Lastly, the paper revealed an increasing trend in the fraction of days with soil moisture below the critical threshold during the past four decades based on ERA5-Land reanalysis data.\nThis study identifies the critical soil water threshold, its global patterns, and the environmental factors, which provide important insights into ecosystem water-stress regimes and soil moisture constraints on global water cycles. The paper is well-written, and the results are clearly and effectively presented. However, I have a few questions and comments regarding the potential confounding factors for critical soil moisture threshold determination and the analysis of ERA5-Land data, that I hope the authors could clarify or consider.\nWe greatly appreciate the reviewer’s positive comments and support of the manuscript. We have thoroughly revised the manuscript following your suggestions. Consequently, our manuscript has been considerably improved. Please see below our detailed responses.\n2.2 First, the analysis could be enhanced by a thorough evaluation of the confounding variables affecting the dLST-SM and EF-SM relationships, thereby, the temporal variation of critical soil moisture threshold. In this study, a single soil moisture threshold was derived for each pixel/site, yet it is possible that such a threshold present temporal variation linked to atmosphere and vegetation changes. As suggested by Feldman et al. (2019, ref 16), VPD consistently mediates the threshold inferred from the dLST-SM relationship in African grasslands. Therefore, it would be beneficial to include confounding variables in the statistical analysis and understand whether and how they may affect spatial patterns of critical soil moisture threshold.\nThank you for this critical comment. While other factors limit evapotranspiration besides SM and the linear dependency is a simple approximation, many previous studies have showed that the EF–SM framework provides a good first-order representation of regimes of land–atmosphere coupling (Fig. 5 from Seneviratne, Corti 16), both in models and observations (e.g., Seneviratne, Corti 16, Seneviratne, Lüthi 17, Koster, Dirmeyer 18, Koster, Suarez 19, Teuling, Seneviratne 20, Feldman, Short Gianotti 3). Here we provided a global analysis based on this first-order theoretical and empirically verified framework.\n7\n[Redacted]\nFig. 5 from Seneviratne, Corti 16. Definition of soil moisture regimes and corresponding evapotranspiration regimes. This paper has been cited over 4000 times by the field since 2010.\nWe fully agree that it is possible that the critical soil moisture threshold present temporal variation linked to atmospheric and vegetation changes. But as Feldman et al. 2019 concluded that, the actual impacts of these confounding factors (e.g. VPD, net radiation) on θ are minor, and that soil moisture remains to be the dominant factor in partitioning surface energy in the transition from non-water limited to water limited regimes. It is also important to note that Feldman et al. 2019’s analysis on confounding factors quantified how short‐term changes of confounding factors influenced dLST–SM relationships. As stated in Feldman et al. 2019, caveats exist in the analysis of confounding factors, including high covariation between regressors, assumption of linear interactions, and inability to account for auto-regressive terms due to irreconcilable time gaps between drydowns. Given these uncertainties, we choose to not conduct such an analysis. Instead, we used random forest models and uncovered the relationships between θ and environmental factors, including climatic, biotic and edaphic factors, to evaluate the possible drivers in controlling the spatial variation of θ . We believe this is much more important than considering the minor effects of confounding factors on θ based on short-term temporal variations.\nIn the revised manuscript, we now added this discussion according to your suggestion.\nLine 392-397:\n“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 16, Seneviratne, Lüthi 17, Koster, Dirmeyer 18, Koster, Suarez 19, Teuling, Seneviratne 20, Feldman, Short Gianotti 3). Here we provided a global analysis based on this first-order theoretical and empirically verified framework.”\n2.3 Secondly, the extensive data gaps in the soil moisture threshold derived from ERA5-Land data presented, as seen in Fig. 4 and Fig. 9, require a more detailed explanation. The extent of missing data, which is much higher than satellite-based analysis, is quite surprising since ERA5-Land is a reanalysis product. I’m interested in understanding the main cause of data gaps: lack of dry-down periods, limited range of SM, or missing data. If the main cause is the failure of the breakpoint model, it may suggest that there are major inconsistencies in SM or LST data between ERA5-Land and satellite estimates. The large spread of the scatters between ERA5-Land derived and satellite-derived theta_crit in Supplementary Fig. 9b partly indicates such a discrepancy. An explanation and discussion of these inconsistencies would be helpful to strengthen the analysis.\nThank you for this constructive comment. We compared the daily ERA5-Land SM and dLST with those from satellites for a day in 2020 (Supplementary Fig. 10). We found that the primary biases between ERA5-Land and satellites lie in SM data, rather than dLST data, emphasizing the importance of global soil moisture measurements. Following your suggestions, we have added the explanation and discussion in the revised manuscript.\nSupplementary Fig. 10. Comparison of daily ERA5-Land SM and dLST with those from satellites for a day in 2020.\nLine 537-543:\n“We noted that there are some data gaps in θ 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 θ (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).”\n2.4 Another minor suggestion is to include a simple plot of theta_crit between satellite-derived v.s. ESM-model estimates, perhaps summarized by biome or climate types. This could provide a more straightforward illustration of how ESMs underestimate the threshold in wet regions and overestimate the threshold in water-limited regions than the maps shown in Fig. 5.\n0.5\n0.4\n0.3\n0.2\n0.1\nSupplementary Fig. 11. Comparison between multi-model mean θ and observation-based θ grouped by climate types based on the aridity classification.\nLine 257-264:\n“We found that the models showed less spatial variability of θ than in the observation-based map (Fig. 5a, Supplementary Figs. 11-12, Fig. 2g) and significantly underestimated θ 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 uptake. Conversely, models significantly overestimated θ in dry regions and failed to capture the observed very low θ values in arid areas (Fig. 5b, Supplementary Fig. 11), which could partly explain why ESMs underestimate both gross and net CO fluxes in dryland ecosystems21, 22.”\nMain text\n2.5 L164-165: This is an interesting hypothesis. A more detailed elaboration on how recently expanded cropland is associated with lower critical soil moisture than established ones would be helpful.\nThank you for this comment. Following your suggestion, we have added more explanations in the revised manuscript.\nLine 476-478:\n“We also tested the hypothesis that the areas of recent cropland expansion over drier marginal lands should be associated with a decrease of θ , as more crop species adapted to dry environments would be selected.”\n2.6 Fig. 4: More context on the ERA5_Land analysis would be helpful for clarity. Did the ERA5-Land analysis use the theta_crit derived from ERA5-Land? Which relationship was used, dLST-SM or ET-SM?\nYes, the ERA5-Land analysis used the θ derived from ERA5-Land SM and dLST. We now stated it clearly in the revised manuscript as suggested.\nLine 528-529:\n“Following satellite datasets analysis, the θ for ERA5-Land was estimated using ERA5-Land SM and dLST.”\nLine 712-713:\n“The global distribution of θ using ERA5-Land surface SM and dLST.”\n2.7 L263: While the ESMs’ theta_crit estimates show a much smaller spatial variability, the relative patterns are quite similar to that of the satellite-derived values, at least in some models.\nThank you for this suggestion. We have modified the sentence to read “Although the relative patterns of θ 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 θ , suggesting the uncertain projection of current and future response of carbon uptake and evapotranspiration to droughts.”\nMethod\n2.8 L295: Does daily minimum LW occur during nighttime?\nThe daily minimum LW occurs either in the morning (e.g., 5 am) or at night (e.g., 19 pm), and this timing varies across different sites and days.\n2.9 L302: To avoid confusion, consider using the official product name for “SCA-V” (e.g., SMAP_L3_XX). I was initially confused that “SCA-V” is soil moisture retrieved from another satellite, before reading the method second in detail.\nThank you for this suggestion. We have added the official product name for “SCA-V” (SMAP_L3_SM_P) as suggested.\nLine 334-336:\n“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)23, single channel vertical polarization (SCA-V, SMAP_L3_SM_P)24 and Soil Moisture and Ocean Salinity in version IC (SMOS-IC)25.”\n2.10 L316: It would be helpful to provide the temporal availability of the SMOS-IC dataset.\nChanged as the reviewer suggested.\nLine 348-349:\n“SMOS-IC has 25 km resolution and two to four-day revisit from 1 January 2011 to 31 December 202025.”\n2.11 L389: Could you please double-check the reference here? After a brief read of Reference 7, I didn’t find relevant results on dLST-SM relationship, but maybe I missed it. In any case, a brief summary of previous work on dLST-SM for theta_crit estimating would be helpful.\nThank you for your careful reading. We have deleted the reference 7 here and added a brief summary of previous work per your suggestions.\nLine 420-423:\n“Evaporative regimes and θ estimating have been characterized previously with observed dLST–SM relationships across some regions, such as Africa3 and site level26, showing that the dLST–SM relationship is an effective method to estimate θ . Here we applied this method to the global scale using multiple satellite observations.”\nLine 88-92:\n“Specifically, the land-surface temperature diurnal amplitude (dLST) starts to increase below θ when ecosystems plunge into the water‐limited regime27, 26, 28, 29. An increased dLST, for a given amount of net radiation, is directly linked to a decrease in EF and thus increased SM stress30. dLST is positively associated with sensible heating but negatively associated with EF and SM26, 28, 29.”\n2.12 L411: I am curious if there are any potential explanations about the consistently higher estimates from satellite versus flux tower in the low theta_crit range.\nThank you for this comment. We have added the possible explanations in the revised manuscript per your suggestion.\nLine 452-455:\n“We noted that the q estimated from satellites is a bit higher than that of flux towers in the low q 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.”\n2.13 L417: Could you please elaborate on the third point about how the high temporal variations could lead to uncertainties in the estimates?\nThank you for this suggestion. We have elaborated on the third point in the revised manuscript as suggested.\nLine 449-452:\n“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 θ values.”\n2.14 L467: It would be helpful to have the final selected hyperparameters as a reference for other and future studies.\nChanged as the reviewer suggested.\nLine 504-507:\n“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.”\n3.1 Fu et al. introduce a novel approach for estimating the critical soil moisture threshold and the fraction of water-limited days, deviating from methods previously limited to reanalysis and climate models. The global map generated in this study is of significant value, establishing a foundational reference for future comparisons. Given my own focus on soil moisture regime and critical values estimation, I find the figures in this study particularly exciting. The writing is logical, and the methods are clearly articulated. While I have only a few comments and suggestions, I recommend acceptance of the paper after minor revisions. Hsin Hsu (Wish not to remain anonymous)\nWe appreciate Dr. Hsu’s positive comments and very helpful suggestions. We have thoroughly revised the manuscript following your suggestions. Please see below our detailed responses.\n3.2 1.The paper shows that the Aridity Index is a primary determinant of spatial variation in θcrit. This also suggests that annual variability in Aridity Index can lead to annual variability in θcrit. However, the annual FSD is calculated based on climatological θcrit rather than year-specific values. This approach potentially compromises the fidelity of the FSD trend estimation. Although the authors reference a study (Line 507) indicating limited changes in θcrit, another research (Hsu and Dirmeyer, 2023) suggests significant changes in θcrit over many locations in one century. The FSD examination of this study (40 year- reanalysis-global) is between the two mentioned extreme cases (15 year-observation-few sites vs. 100 year with +1%CO2 per year in climate model). Given that the period used in all of these studies are within similar increasing rate of CO2, the evidence of 15-year observation analysis might not effectively support authors’ statement. The authors could consider a decadal θcrit estimation from ERA5-land as supplementary information and/or moving the argument of temporal variation of θcrit from method to main text with a more comprehensive discussion.\nHsu and Dirmeyer, 2023. Uncertainty in Projected Critical Soil Moisture Values in CMIP6 Affects the Interpretation of a More Moisture-Limited World. Doi: 10.1029/2023EF003511\nThank you for this constructive comment. Following your suggestion, we have moved the argument of temporal variation of θ from method to main text with a more comprehensive discussion by referring to Hsu and Dirmeyer 31.\nLine 247-253:\n“We acknowledge that θ may change over time. Based on model outputs analysis, Hsu and Dirmeyer 31 found significant temporal variations in θ across many locations spanning 100 years. Conversely, another study analyzed the temporal dynamics of θ at five flux tower sites with at least 15 years of measurements and found no significant trend over time2. This underscores the need for future research to gain a better understanding of the temporal dynamics of θ through longer observations. We considered here that the temporal dynamics of θ should not hamper our trend analysis, given that even if θ changes, its magnitude over 40 years is minimal.”\n3.3 2.Maybe I miss this somewhere: The paper lacks clarity regarding the sampling method for drydown used to estimate θcrit for climate models. It is unclear whether the author selects identical months as observations or the model's peak growing seasons. Discussing the implications of this choice is essential, as differences in background climate (and the aridity index) between observations and model outputs may lead to variations in estimated θcrit.\nThank you for this suggestion. The full year data of satellite observations were used and results were found to be similar when only growing season data were used. Following satellite datasets analysis, the full year data of model outputs were also used. This is now clearly stated in the revised manuscript.\nLine 560-564:\n“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 θ pixel-by-pixel by means of the EF–SM relationship.”\n3.4 3.Regarding the comment #3 on divergence of hydroclimate among data products. how does author think the inherent divergence of simulated soil moisture distribution (Koster et al. 2009) affects the interpretation here. Does this extends beyond model-specific vegetation/soil dynamics as discussed in lines 231-236?\nKoster et al. 2009: On the Nature of Soil Moisture in Land Surface Models. Doi: 10.1175/2009JCLI2832.1\nThank you for this critical comment. We agree that different models led to different simulated soil moisture values, which could contribute to the differences between observation-based θ and models-based θ values. But we found that all models consistently showed less spatial variability of θ than in the observation-based map (Supplementary Fig. 12), suggesting that our result did not depend on the inherent divergence of simulated soil moisture distribution. Please also note that we only used the surface soil moisture for all models so that the effects of the inherent divergence of simulated soil moisture from different models are muted.\nFollowing your suggestion, we have discussed the effects of inherent divergence of simulated soil moisture distribution from different models by referring to Koster, Guo 32.\nLine 566-570:\n“We noted that different models led to different simulated SM values32, and this inherent divergence of simulated SM distribution could also contribute to the differences between observation-based θ and models-based θ values. But we found that all models consistently showed less spatial variability of θ than in the observation-based map, suggesting that our result did not depend on the inherent divergence of simulated SM distribution.”\n3.5 4.While the method using dLST appears straightforward, the main text could benefit from a physical description of the connection between using Evaporative Fraction (EF) and dLST, explaining why this approach is effective (Line 84-85). This addition would enhance the audience's understanding of the methodology.\nLine 88-92:\n“Specifically, the land-surface temperature diurnal amplitude (dLST) starts to increase below θ when ecosystems plunge into the water‐limited regime27, 26, 28, 29. An increased dLST, for a given amount of net radiation, is directly linked to a decrease in EF and thus increased SM stress30. dLST is positively associated with sensible heating but negatively associated with EF and SM26, 28, 29.”\n3.6 Line 77. typo: atmosphere\nChanged as suggested.\n1. Bassiouni M, Good SP, Still CJ, Higgins CW. Plant water uptake thresholds inferred from satellite soil moisture. Geophysical Research Letters 47, e2020GL087077 (2020).\n2. Fu Z, et al. Critical soil moisture thresholds of plant water stress in terrestrial ecosystems. Science Advances 8, eabq7827 (2022).\n3. Feldman AF, Short Gianotti DJ, Trigo IF, Salvucci GD, Entekhabi D. Satellite‐based assessment of land surface energy partitioning–soil moisture relationships and effects of confounding variables. Water Resources Research 55, 10657-10677 (2019).\n4. Bassiouni M, Higgins CW, Still CJ, Good SP. Probabilistic inference of ecohydrological parameters using observations from point to satellite scales. Hydrology and Earth System Sciences 22, 3229-3243 (2018).\n5. Dong J, Akbar R, Short Gianotti DJ, Feldman AF, Crow WT, Entekhabi D. Can Surface Soil Moisture Information Identify Evapotranspiration Regime Transitions? Geophysical Research Letters, e2021GL097697 (2022).\n6. Xiao J, Fisher JB, Hashimoto H, Ichii K, Parazoo NC. Emerging satellite observations for diurnal cycling of ecosystem processes. Nat Plants 7, 877-887 (2021).\n7. Wen J, Fisher JB, Parazoo NC, Hu L, Litvak ME, Sun Y. Resolve the Clear‐Sky Continuous Diurnal Cycle of High‐Resolution ECOSTRESS Evapotranspiration and Land Surface Temperature. Water Resources Research 58, e2022WR032227 (2022).\n8. Novick KA, et al. Confronting the water potential information gap. Nat Geosci 15, 158-164 (2022).\n9. Denissen J, et al. Widespread shift from ecosystem energy to water limitation with climate change. Nature Climate Change 12, 677-684 (2022).\n10. Konings AG, Gentine P. Global variations in ecosystem-scale isohydricity. Global Change Biology 23, 891-905 (2017).\n11. Li X, et al. Global variations in critical drought thresholds that impact vegetation. National Science Review 10, nwad049 (2023).\n12. Du J, et al. Synergistic satellite assessment of global vegetation health in relation to ENSO‐induced droughts and pluvials. Journal of Geophysical Research: Biogeosciences 126, e2020JG006006 (2021).\n13. 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).\n14. Kolby Smith W, et al. Large divergence of satellite and Earth system model estimates of global terrestrial CO2 fertilization. Nature climate change 6, 306-310 (2016).\n15. Pei Y, et al. Evolution of light use efficiency models: Improvement, uncertainties, and implications. Agricultural and Forest Meteorology 317, 108905 (2022).\n16. Seneviratne SI, et al. Investigating soil moisture–climate interactions in a changing climate: A review. Earth-Science Reviews 99, 125-161 (2010).\n17\n17. Seneviratne SI, Lüthi D, Litschi M, Schär C. Land–atmosphere coupling and climate change in Europe. Nature 443, 205 (2006).\n18. Koster RD, et al. Regions of strong coupling between soil moisture and precipitation. Science 305, 1138-1140 (2004).\n19. Koster RD, et al. Realistic initialization of land surface states: Impacts on subseasonal forecast skill. Journal of Hydrometeorology 5, 1049-1063 (2004).\n20. Teuling A, Seneviratne SI, Williams C, Troch P. Observed timescales of evapotranspiration response to soil moisture. Geophysical Research Letters 33, (2006).\n21. MacBean N, et al. Dynamic global vegetation models underestimate net CO2 flux mean and inter-annual variability in dryland ecosystems. Environmental Research Letters 16, 094023 (2021).\n22. Wang L, et al. Dryland productivity under a changing climate. Nature Climate Change, 1-14 (2022).\n23. Li X, et al. A new SMAP soil moisture and vegetation optical depth product (SMAP-IB): Algorithm, assessment and inter-comparison. Remote Sensing of Environment 271, 112921 (2022).\n24. Jackson TJ. III. Measuring surface soil moisture using passive microwave remote sensing. Hydrological processes 7, 139-152 (1993).\n25. Wigneron J-P, et al. SMOS-IC data record of soil moisture and L-VOD: Historical development, applications and perspectives. Remote Sensing of Environment 254, 112238 (2021).\n26. Amano E, Salvucci GD. Detection and use of three signatures of soil-limited evaporation. Remote sensing of environment 67, 108-122 (1999).\n27. Panwar A, Kleidon A, Renner M. Do surface and air temperatures contain similar imprints of evaporative conditions? Geophysical Research Letters 46, 3802-3809 (2019).\n28. Betts AK, Desjardins R, Worth D, Beckage B. Climate coupling between temperature, humidity, precipitation, and cloud cover over the Canadian Prairies. Journal of Geophysical Research: Atmospheres 119, 13,305-313,326 (2014).\n29. Thakur G, Schymanski SJ, 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. Scientific reports 12, 1-14 (2022).\n30. Gentine P, Chhang A, Rigden A, Salvucci G. Evaporation estimates using weather station data and boundary layer theory. Geophysical Research Letters 43, 11,661-611,670 (2016).\n31. Hsu H, Dirmeyer PA. Uncertainty in Projected Critical Soil Moisture Values in CMIP6 Affects the Interpretation of a More Moisture‐Limited World. Earth's Future 11, e2023EF003511 (2023).\n32. Koster RD, Guo Z, Yang R, Dirmeyer PA, Mitchell K, Puma MJ. On the nature of soil moisture in land surface models. Journal of Climate 22, 4322-4335 (2009).\nReviewer #1 (Remarks to the Author):\nThe authors have addressed my comments in their revision. Nice work.\nReviewer #2 (Remarks to the Author):\nDear authors,\nThank you very much for thoroughly addressing my previous comments! All my concerns have been comprehensively resolved and the paper has been significantly improved in clarity. I just have a few minor suggestions for the authors to consider incorporating.\nLine 536-537 (track-change manuscript): Thanks for the clarification. Could you please also specific the temporal range for this analysis? I am curious to know whether it covers 2016-2020 (consistent with satellite-based analysis) or the entire time period of ERA5-Land.\nSupplementary Fig. 9, Fig.10: The scatter plots might benefit from using a point density plot (or adjusting the transparency of the points) to mitigate the visual bias caused by overlapping. This overlap seems to visually amplify the impact of outliers and reduce the actual significance of the correlations.\nThanks again for your efforts in revising the manuscript and addressing my comments!\nCongratulations on completing an exciting paper!\nReviewer #2 (Remarks on code availability):\nThe authors have provided a ReadMe file that provides detialed information in running the code to reproduce the results. Yet, since the raw data is public and not included in the repository, I suggest providing a more detailed description of the raw data format and preprocessing processes (if any) to faciliate reproducing. For example, the first step, \"SM dry-down identification.m\" reads in \"*.mat\" files, it would be benefitial to briefly outlinet the content of these files as comment or in the ReadMe.\nReviewer #3 (Remarks to the Author):\nI appreciate the authors for addressing the comments and modifying their manuscript accordingly. The revised manuscript extensively discusses the relevant literature and uncertainty of their analysis, and provides great detail on their methods. I close my paper-review work and anticipate the paper's publication.\nMinor comment:\nThe information provided in Line 310: \"have undergone a standardized set of quality control and gap filling\" is highly redundant to Line 320: \"Data were quality controlled so that only measured...\"\n-Hsin Hsu\nGlobal critical soil moisture thresholds of plant water stress\n#NCOMMS-23-53911A\nResponse to Reviewers\nReviewer #1 (Remarks to the Author):\n1.1 The authors have addressed my comments in their revision. Nice work.\nThank you very much for your constructive comments in the process.\nReviewer #2 (Remarks to the Author):\n2.1 Thank you very much for thoroughly addressing my previous comments! All my concerns have been comprehensively resolved and the paper has been significantly improved in clarity. I just have a few minor suggestions for the authors to consider incorporating.\nWe greatly appreciate the reviewer’s positive comments and support of the manuscript. We have thoroughly revised the manuscript following your suggestions. Please see below our detailed responses.\n2.2 Line 536-537 (track-change manuscript): Thanks for the clarification. Could you please also specific the temporal range for this analysis? I am curious to know whether it covers 2016-2020 (consistent with satellite-based analysis) or the entire time period of ERA5-Land.\nThank you very much for your suggestions. As per your suggestions, we have specified the temporal range in Line 525-526: “We first compared the FSD from ERA5-Land reanalysis SM during 2016-2020 with that of satellite observations.”\n2.3 Supplementary Fig. 9, Fig.10: The scatter plots might benefit from using a point density plot (or adjusting the transparency of the points) to mitigate the visual bias caused by overlapping. This overlap seems to visually amplify the impact of outliers and reduce the actual significance of the correlations.\nThank you for this constructive comment. Following your suggestion, we have adjusted the transparency of the points in Supplementary Figs. 9-10.\n2.4 Thanks again for your efforts in revising the manuscript and addressing my comments!\nCongratulations on completing an exciting paper!\nThank you for your kind words and valuable suggestions.\n2.5 Reviewer #2 (Remarks on code availability):\nThe authors have provided a ReadMe file that provides detialed information in running the code to reproduce the results. Yet, since the raw data is public and not included in the repository, I suggest providing a more detailed description of the raw data format and preprocessing processes (if any) to faciliate reproducing. For example, the first step, \"SM dry-down identification.m\" reads in \"*.mat\" files, it would be benefitial to briefly outlinet the content of these files as comment or in the ReadMe.\nThank you for this suggestion. Following your suggestion, we have added the detailed description of the raw data format and preprocessing processes in the ReadMe file. As an example, the ‘SM_SCA-V_2020.mat’ and ‘Delta_LST_2020.mat’ files have been included in the repository.\nReviewer #3 (Remarks to the Author):\n3.1 I appreciate the authors for addressing the comments and modifying their manuscript accordingly. The revised manuscript extensively discusses the relevant literature and uncertainty of their analysis, and provides great detail on their methods. I close my paper-review work and anticipate the paper's publication.\nMinor comment:\nThe information provided in Line 310: \"have undergone a standardized set of quality control and gap filling\" is highly redundant to Line 320: \"Data were quality controlled so that only measured...\"\n-Hsin Hsu\nThank you for your positive comments and helpful suggestions. As per your suggestions, we have removed Line 320: \"Data were quality controlled so that only measured...\"."}]

Decision letter