Critical soil moisture detection and water–energy limit shift attribution using satellite-based water and carbon fluxes over China

Critical soil moisture (CSM), a tipping point of soil moisture (SM) at which surface fluxes shift from the energy-limited regime to the water-limited regime, is essential for the vegetation state and the corresponding land–atmosphere coupling. However, detecting CSM and attributing water–energy limit shifts to climate and ecosystem variables are challenging as in situ observations of water, carbon fluxes, and soil moisture (SM) are sparse. In this study, CSM was assessed over China using two satellite-based methods: (i) the difference between the correlation between SM and evapotranspiration (ET) and the correlation between vapor pressure deficit (VPD) and ET and (ii) the covariance between VPD and gross primary production (GPP). ET and GPP products were based on the Penman–Monteith–Leuning (PML) ET and GPP, Global LAnd Surface Satellite (GLASS) ET and GPP, Collocation-Analyzed Multi-source Ensembled Land Evapotranspiration (CAMELE) ET, Surface Energy Balance Algorithm for Land (SEBAL) ET, two-leaf light use efficiency (TL-LUE) GPP, and solar-induced chlorophyll fluorescence (SIF)-based (GOSIF) GPP. At flux sites, ET and GPP products were evaluated by eddy-covariance-based measurements; CSM values using two satellite-based methods were assessed using the soil moisture–evaporative fraction method. Their consistency at site scales demonstrated reliable results and applicability to regional scales. Through intercomparison, the spatial pattern of CSM from multi-source ET and GPP datasets was consistent and robust in eastern and southern China. Generally, CSM decreased from south to north. The Pearl River basin and Southeastern River basin displayed a relatively high CSM for clay-rich soils (e.g., 0.39 m³ m⁻³ using PML ET and 10 cm depth SM) and forests (e.g., 0.35 m³ m⁻³ using PML ET and 10 cm depth SM). Since CSM values were higher than the average SM at four soil layers, grassland and clay were water-limited. Thus, with increased water demand, western grasslands were more susceptible to water stress. These findings highlight the variability in CSM and the primary determinants of water–energy limit shifts, offering valuable insights into the potential water limitations on ecosystems under comparable SM circumstances.

1 Introduction

Critical soil moisture (CSM) serves as an indicator of shifts in the relationship between water and energy availability (Schwingshackl et al., 2017; Denissen et al., 2020) and is essential in shaping regional climates. Plants adjust their stomatal resistance in response to changes in soil moisture (SM) and vapor pressure deficit (VPD) (Grossiord et al., 2020; F. Li et al., 2023). Above CSM, there is no alteration in water stress with SM increases (Rodriguez-Iturbe, 2000; Seneviratne et al., 2010; Akbar et al., 2018); plants are primarily controlled by VPD. Warm and dry air above the canopy (Grossiord et al., 2020; X. Li et al., 2023) leads to a decrease in both the transpiration process as the largest part of evapotranspiration (ET) (Good et al., 2015) and gross primary production (GPP) coupled with the ET process via plant leaf stomata (Gentine et al., 2019; Liu et al., 2020). A decrease in ET, in turn, results in elevated surface temperature and VPD (Gentine et al., 2019) and leads to increased atmospheric aridity on a large spatial scale, thereby intensifying SM depletion. Below CSM, surface fluxes are primarily influenced by SM availabilities in conditions of restricted water supply. During this period, a decrease in SM results in a reduction in latent heat flux (LE) and an increase in sensible heat flux (H) (Rodriguez-Iturbe, 2000); SM and leaf conductance follow a positive linear relationship (Laio et al., 2001; Porporato et al., 2002). Previous studies have examined land–atmosphere feedbacks using different metrics and both observation and simulation data (Seneviratne et al., 2006; Koster et al., 2009; Teuling et al., 2009). They found that water and energy limit shifts may be further strengthened by the interaction between the land and atmosphere, particularly with positive feedback mechanisms, known as the “dry-gets-dryer” effect (Seneviratne et al., 2010; Gentine et al., 2019). Over extended temporal periods, this phenomenon may lead to the persistence of arid and high-temperature conditions (Zhang et al., 2020). Consequently, it is necessary to quantify the CSM characteristics and the influencing environmental factors of water and energy limit shifts.

Traditionally, under the framework based on the ratio of LE to the total of LE and H (Haghighi et al., 2018; Fu et al., 2022a), sparse eddy covariance observations (Feldman et al., 2019; Fu et al., 2022b) pose challenges in adequately capturing comprehensive regional- or continental-scale CSM and its variations (Dong et al., 2023; Hsu and Dirmeyer, 2023a). In recent years, the feasibility of conducting large-scale analyses has been enhanced by the growing accessibility of multi-source remotely sensed datasets (Liu et al., 2012). Globally, some studies used the ratio of LE to net radiation (Seneviratne et al., 2010; Schwingshackl et al., 2017), surface temperature diurnal amplitude (Feldman et al., 2019; Fu et al., 2024), and LE (Hsu and Dirmeyer, 2023b; Duan et al., 2023). In addition, the advancement of global remote sensing technology has facilitated the generation of reliable GPP products (Yuan et al., 2014; Li and Xiao, 2019; Zhang et al., 2019; Bi et al., 2022; He et al., 2022; F. Li et al., 2023) and ET products (Yao et al., 2013; Yao et al., 2014; Zhang et al., 2019; Cheng et al., 2021; He et al., 2022; C. Li et al., 2022; F. Li et al., 2023) for CSM detection. Denissen et al. (2020) proposed a new tipping-point metric, the difference between the correlation between SM and ET and the correlation between VPD and ET, to directly determine CSM at continental scales. Fu et al. (2022b) first demonstrated that the covariance between GPP and VPD indirectly quantifies CSM. The point at which the covariance between GPP and VPD transitions from positive to negative during a period of soil drying is denoted as CSM. However, a source of considerable uncertainty when considering only a single data source and estimation approach exists at a large spatial scale. There are significant differences among satellite-based ET and GPP datasets, and CSM varies with different methods, leading to uncertainty as to whether the CSM of carbon flux is the same as that of water flux.

Chinese land surfaces frequently experience water and energy limit shifts (Xiao, 2014; Zhu et al., 2023). Diagnosing large-scale CSM helps to understand water- and energy-limited regimes determined by distinct flora and soil types (Homaee et al., 2002; Hsu and Dirmeyer, 2023b). The association between water, energy, and flux helps to define water and energy limit shifts. As such, this study uses two innovative metrics and eight satellite-based products to diagnose CSM and water and energy limit shifts across China. The goal of this study is to (1) assess the consistency of different methods in calculating CSM at flux sites; (2) examine CSM variations across land cover types, soil textures, and water resource subregions; and (3) investigate dominant factors from hydrological, meteorological, and ecological variables that influence water and energy limit shifts.

2 Material and methods

2.1 Data

Eddy covariance flux datasets were compared with eight satellite-based ET and GPP datasets in Sect. 3.1. Then, CSM derived from the relationship between SM and evaporative fraction (EF) was used to evaluate the performance of CSM derived from the covariance and correlation difference methods in Sect. 3.2. Layer-wise SM and satellite-based ET and GPP were used for the large-scale detection of CSM. Land cover types, soil textures, and water resource subregions were used to examine CSM variations in Sect. 3.3. SM, ET, GPP, and meteorological data were used to investigate dominant factors influencing water and energy limit shifts in Sect. 3.4. All energy, vegetation, and water variables were resampled or combined to 0.1°-8 d resolution. The period, limited by the temporal availability of data sources, covered 2001–2018.

2.1.1 Evapotranspiration and gross primary production

Figure 1 illustrates the locations of 21 flux sites, and Table 1 shows the detailed information on flux sites. Eddy-covariance-derived measurements were applied to evaluate the performance of satellite-based ET and GPP. Given the fact that Huazhaizi, Dashalong, Luodi, Arou, Guantao, Huailai, Miyun, and Daxing did not have GPP data, the REddyProc website (https://www.bgc-jena.mpg.de/5622399/REddyProc/, last access: 16 July 2023) was used to calculate GPP. REddyProc imported half-hourly net ecosystem exchange, LE, H, and meteorological measurements to partition net ecosystem exchange into GPP and ecosystem respiration.

Figure 1 (a) Flux site locations, land cover types (2020), and water resource subregions of China. Distributions of (b) clay, (c) silt, and (d) sand content (1995). ZGE: Zhungaer basin, PR: Pearl River basin, YTR: Yangtze River basin, SWR: Southwestern River basin, TR: Tarim basin, SR: Songhua River basin, CT: Changthang region, NM: Inner Mongolian Plateau region, LR: Liaohe River basin, YR: Yellow River basin, HR: Huaihe River basin, HC: Hexi Corridor region, HAR: Haihe River basin, SER: Southeastern River basin, QB: Qaidam basin.

Table 1 Flux site information used in this study.

TPDC: National Tibetan Plateau Data Center.

Table 2 Spatial gridded datasets used in this study.

Table 2 contains a list of all spatial datasets used in this study. Eight satellite-based ET and GPP products are included. The Penman–Monteith–Leuning (PML) product integrates the stomatal conductance theory to relate ET and GPP processes (Zhang et al., 2019; He et al., 2022) and applies daily meteorological data, land surface temperature from ERA5, enhanced Whittaker-filtered MODIS leaf area index (LAI), albedo, and emissivity. The interdependency and mutual restrictions between GPP and ET considerably improve the simulation accuracy. The Global LAnd Surface Satellite (GLASS) ET integrates the MOD16, a revised remote-sensing-based Penman–Monteith product; the Priestley–Taylor Jet Propulsion Laboratory, a modified satellite-based Priestley–Taylor product; and the Semi-Empirical Algorithm of the University of Maryland using the Bayesian model averaging approach (Yao et al., 2013, 2014). The GLASS GPP algorithm incorporates the effects of atmospheric carbon dioxide content, radiation components, and VPD based on the eddy covariance–light use efficiency model introduced by Yuan et al. (2007). It is founded on two underlying assumptions: the fraction of absorbed photosynthetically active radiation has a linear relationship with the normalized difference vegetation index, and constant light use efficiency is governed by either air temperature or soil moisture, depending on which component imposes the greatest limitation.

In addition, Collocation-Analyzed Multi-source Ensembled Land Evapotranspiration (CAMELE) provides long-term (1981–2020) ET, employing ERA5, FLUXCOM, PML, GLDAS, and GLEAM (C. Li et al., 2022), at 0.1° 8 d and 0.25° daily resolutions. Surface Energy Balance Algorithm for Land (SEBAL) ET focuses on 1 km daily resolution during 2001–2018. This product integrates GMAO's meteorological data and NASA's MOD43A1 daily surface albedo, the MOD11A1 daily surface temperature, and the MOD13 vegetation index (Cheng et al., 2021). Two-leaf light use efficiency (TL-LUE) GPP offers comprehensive worldwide GPP, shaded GPP, and sunlit GPP covering the period from 1992 to 2020. This model applies recent data inputs such as the GLOBMAP LAI, CRUJRA meteorological data, and ESA-CCI land cover information (Bi et al., 2022). Global Orbiting Carbon Observatory-2 SIF (GOSIF) GPP spans from 2000 to 2020. A total of eight SIF–GPP relationships, including both universal and biome-specific formulations, are used to estimate GPP from SIF on a per-pixel basis and are examined with and without intercept terms to account for the uncertainty in converting SIF into GPP estimates (Li and Xiao, 2019).

2.1.2 Layer-wise soil moisture and meteorological data

Given the recent availability of state-of-the-art gridded SM in China as released by Q. Li et al. (2022), CSM can now be investigated in the context of the SM state. Gridded SM reaches 100 cm soil depth with 10 cm intervals at 1 km daily resolution during 2000–2020. It is trained by predictors of ERA5-Land time series, leaf area index (LAI), land cover type, topography, and in situ observed soil attributes at 1789 stations throughout China using the robust random forest machine learning technique. Based on the findings of Q. Li et al. (2022), the product demonstrates notable benefits over both ERA5-Land and SMAP-L4 datasets, especially in terms of a superior quality level compared to the SoMo.ml dataset at soil depths of 10, 20, 80, and 100 cm. Thus, this study utilized the SM of these layers.

Yang et al. (2010) and He et al. (2020) put forth a comprehensive dataset for Chinese regional surface meteorological forcing. This dataset encompasses air temperature, air pressure, specific humidity, wind speed, downward shortwave radiation, downward longwave radiation, and precipitation. It is presented in NetCDF format with a spatiotemporal resolution of 0.1° and 3 h during 1979–2018. The primary input includes Princeton University's global land surface model data, GLDAS, GEWEX-SRB radiation, TRMM precipitation, and China Meteorological Administration observations. Data quality control techniques include the elimination of physically implausible values and statistical interpolation using ANU-Spline. This dataset demonstrates precision levels that lie between those of site-based observation and satellite-based estimation and, therefore, exceed the accuracy of current international reanalysis datasets. In this study, VPD was computed based on specific humidity and air temperature; VPD, air temperature, precipitation, and downward shortwave radiation were employed in water–energy limit shift attribution.

2.1.3 Land cover types, soil textures, and water resource subregions

Land cover types, soil textures, and water resource subregions influence CSM. Land cover types (2020) are created by human visual interpretation relying on Landsat satellite remote sensing images. The categorization scheme includes cropland, forests, grassland, water, ice, urban, and barren. Soil textures, expressed as sand, silt, and clay content within each grid cell, are compiled from the 1:1 000 000 soil type map and the Second National Soil Survey. Water resource subregions are divided according to the China Geological Survey, including the Zhungaer basin, Pearl River basin, Yangtze River basin, Southwestern River basin, Tarim basin, Songhua River basin, Changthang region, Inner Mongolian Plateau region, Liaohe River basin, Yellow River basin, Huaihe River basin, Hexi Corridor region, Haihe River basin, Southeastern River basin, and Qaidam basin. Water resource subregions are based on the principles of groundwater systems and water cycles and are focused on the inherent features of groundwater resources within distinct natural units.

2.2 Determination of CSM

CSM derived by the SM and EF method was used to assess the CSM using ET and GPP at the flux site. There must be both positive and negative metrics from the covariance and correlation difference methods. For each grid cell and the entire period per year, negative metrics are displayed when SM is less than CSM, and positive metrics are shown when SM is greater than CSM. The data are taken into account only when the temperature surpasses 10 °C (Denissen et al., 2020) to avoid the influence of ice and snow, and the covariance between VPD and GPP must exhibit a minimum of seven covariance values within 9 d moving windows, with a minimum of 15 data points (Fu et al., 2022b). Hence, we concentrated on the warm season, June–September, which includes 16 data points each year, with eight covariance values within 9 d moving windows. CSM was considered in each grid cell using satellite-based ET and GPP over the period 2001–2018.

2.2.1 Soil moisture–evaporative fraction method

Investigating the relationship between SM and EF in the dry period can isolate the transition from energy limitation to water limitation (Feldman et al., 2019). CSM captures the interconnectedness between SM and EF. If SM is greater than or less than CSM, the relationship between SM and EF appears as a flat line or a positive slope line. A linear-plus-plateau model characterizes the relationship precisely (Seneviratne et al., 2010; Schwingshackl et al., 2017):

$$\begin{array}{}\text{(1)}& \text{EF}=\left\{\begin{array}{ccc}{\text{EF}}_{max}+S\left(\text{SM}-\text{CSM}\right),& \text{if}& \text{SM}<\text{CSM}\\ {\text{EF}}_{max},& \text{if}& \text{SM}\ge \text{CSM}\end{array}\right.,\end{array}$$

where EF is the evaporative fraction, defined as LE $/$ (LE+H); EF_max represents the maximum EF in the energy-limited stage; and S is the gradient in the water-limited stage. Here, specific estimated CSM was simultaneously estimated by the Monte Carlo method. For a set of optimal parameters, a Nash–Sutcliffe efficiency (Nash and Sutcliffe, 1970) above 0.5 was considered to be satisfactory (Herman et al., 2018). Thus, only eight sites, including Xilingele in 2004, Damshung in 2004, CN-Sw2 in 2011, CN-Du2 in 2007, CN-Cng in 2010, Miyun in 2009, Huailai in 2015, and Qianyanzhou in 2010, were chosen for CSM detection. In addition, the Bayesian information criterion (BIC) (Schwarz, 1978) was used to select the best fit among three segmented regression candidates (the flat line, the positive slope line, and the linear-plus-plateau). If the flat-line regression or the positive-slope regression outperformed the linear-plus-plateau regression, CSM was considered to be unidentified.

2.2.2 Covariance method

The covariance method presents a novel method for assessing ecosystem water stress in direct correlation with GPP, as illustrated by Fu et al. (2022b). It serves to quantify CSM over large areas. Positive covariances between VPD and GPP indicate that energy limits GPP. Negative covariances indicate that water limitation has a larger impact on GPP. VPD is determined by the disparity between the saturation vapor pressure (e_s) and the actual vapor pressure (e_a). Bolton (1980) posits that the calculation of e_a involves specific humidity (SH) and surface pressure (Pr):

$$\begin{array}{}\text{(2)}& e_{\mathrm{a}}={\displaystyle \frac{\text{SH}\times \text{Pr}}{\text{SH}\times \mathrm{0.378}+\mathrm{0.622}}}.\end{array}$$

2.2.3 Correlation difference method

Another novel correlation difference metric, proposed by Denissen et al. (2020), evaluates water-limited conditions versus energy-limited conditions using the detrended anomaly of VPD, ET, and SM:

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }\text{corr}=\text{corr}(\text{ET},\text{VPD})-\text{corr}(\text{ET},\text{SM}).\end{array}$$

MATLAB's corr tool calculates this metric, which uses Kendall's rank correlation (corr) rather than assuming linear correlations between variables (van Doorn et al., 2018). If Δcorr > 0 then the grid cell is under an energy-limited regime, and vegetation anomalies (i.e., ET) correlate more strongly with energy anomalies (i.e., VPD) than with water anomalies (i.e., SM). In contrast, if Δcorr < 0 then the grid cell is under a water-limited regime. When Δcorr ≈ 0, SM is labeled as CSM, indicating that water- and energy-limited regimes are transitioning.

2.3 Evaluation criteria

The correlation coefficient was applied to evaluate the performance of satellite-based ET from CAMELE, GLASS, PML, and SEBAL and GPP from GOSIF, GLASS, PML, and TL compared to the eddy-covariance-observed in situ ET and GPP. A point-to-pixel evaluation was carried out to evaluate the overestimation or underestimation of ET and GPP for each land cover type from all 21 flux sites. We summed 8 d ET and GPP in grassland, evergreen broadleaf forests, evergreen needleleaf forests, mixed forests, cropland, wetland, and barren land.

The alignment of CSM obtained by different methods was determined using the chi-square test (McHugh, 2013; Hsu and Dirmeyer, 2023a). CROSSTAB in MATLAB was used to perform the chi-square test. SM values were divided into two groups, below and above CSM. In this case, categorical data were tagged as a binary variable of 0 for drier than CSM and 1 for wetter than CSM. If there were significant differences with a 95 % confidence level, CSM was different.

2.4 Partial least square regression

Partial least square regression has been widely acknowledged as a viable approach for mitigating collinearity issues among independent variables (Karthikeyan et al., 2020) and has been used extensively in quantifying interannual impacts of climate and plant growth variations on water and energy dynamics. Here, performances of the partial least square regression model were assessed by 5-fold cross-validation using the mean absolute percentage error. The dominant factor (precipitation, temperature, incoming shortwave radiation, VPD, ET, GPP, or SM) of Δcorr was identified by the largest variable importance in projection scores.

Figure 2 Correlation coefficients between daily eddy covariance observations and satellite-based (a) ET and (b) GPP products and average (c) ET and (d) GPP across land cover types. The bars show 95 % confidence intervals. GL: grassland, EBF: evergreen broadleaf forests, ENF: evergreen needleleaf forests, MF: mixed forests, CL: cropland, WL: wetland.

Figure 3 Variations in soil moisture and evaporative fraction at (a) Xilingele in 2004, (c) Damshung in 2004, (e) CN-Sw2 in 2011, (g) CN-Du2 in 2007, (i) CN-Cng in 2010, (k) Miyun in 2009, (m) Huailai in 2015, (o) Qianyanzhou in 2010. Variations in covariance (referred to as Cov) between vapor pressure deficit (VPD) and gross primary production (GPP) and correlation difference metric (referred to as Δcorr) at (b) Xilingele in 2004, (d) Damshung in 2004, (f) CN-Sw2 in 2011, (h) CN-Du2 in 2007, (j) CN-Cng in 2010, (l) Miyun in 2009, (n) Huailai in 2015, (p) Qianyanzhou in 2010.

3 Results

3.1 Consistency of ET and GPP

Figure 2a and b show good agreement between daily satellite-based products and site observations in most land cover types. Across all sites, correlation coefficients obtained from CAMELE, GLASS, PML, and SEBAL ET were 0.74, 0.65, 0.78, and 0.59, respectively; correlation coefficients obtained from GLASS, TL, GOSIF, and PML GPP were 0.75, 0.71, 0.77, and 0.74, respectively. For ET, the highest correlation coefficient occurred between GLASS and eddy covariance observations in mixed forests (0.96), while the lowest value was between SEBAL and site observations in barren land (0.47). For GPP, the highest correlation coefficient was found between TL and site measurements in mixed forests (0.97), while the lowest value was between GLASS and site-based data in barren land (0.32). In general, no single product consistently outperformed others over all land cover types. As shown in Fig. 2c and d, ET had the highest value in evergreen needleleaf forests and was the lowest in barren land, while GPP peaked in evergreen broadleaf forests and was the lowest in wetland. In these land cover types, ET and GPP derived from satellite-based products were also substantially different and varied quite a bit between different products. Especially in evergreen broadleaf forests, ET derived from GLASS (3.37 mm) and CAMELE (3.05 mm) and GPP derived from GOSIF (7.55 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) and TL (7.59 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) were higher than site observations of 1.74 mm and 5.26 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, respectively. If satellite-based ET and GPP were around ± 10 % of site-observed values, they were termed as satisfactory; otherwise, they were either overestimated or underestimated. CAMELE, GLASS, PML, and SEBAL ET and GLASS, TL, GOSIF, and PML GPP met satisfactory values in one, one, three, and one and two, one, three, and two land cover types, respectively. PML ET provided the most satisfactory estimates in evergreen broadleaf forests, cropland, and barren land, with an average bias of 1.05 %, 1.13 %, and 1.34 %, respectively; GOSIF GPP provided the most satisfactory estimates in grassland, evergreen needleleaf forests, and mixed forests, with an average bias of 4.31 %, 9.14 %, and 4.29 %, respectively. Although discrepancies existed among multi-source remotely sensed products across flux sites, they offered an opportunity to quantify characteristics of large-scale CSM and to examine uncertainties from single-source data.

3.2 Consistency of CSM

Variations in SM and EF were depicted in Fig. 3 for eight sites. Fitted lines represented controlling mechanisms in various evaporative regimes. Overall, the linear-plus-plateau regression with the lowest BIC outperformed the flat line and the positive-slope line in the study period. Specifically, the CN-Du2 and Qianyanzhou sites showed a great slope at low SM values, with BIC values of −80.29 and −98.64, respectively. We also found that grassland CSM in different regions varied greatly. For example, grassland at Xilingele had the lowest CSM of 0.079 m³ m⁻³, with SM ranging from 0.06 to 0.20 m³ m⁻³; CSM at Damshung in southwestern China was 0.175 m³ m⁻³, with SM ranging from 0.14 to 0.26 m³ m⁻³; CSM at CN-Cng in northeastern China was 0.457 m³ m⁻³, with high SM ranging from 0.30 to 0.70 m³ m⁻³. Ranges of SM determined CSM value. Moreover, vertical lines of different colors represented CSM derived from Δcorr using the correlation difference method and covariance between VPD and GPP. To explore the performance of both methods at different sites and whether they can be used on a large scale, the data applied to both methods were averaged for 8 d, consistently with gridded data. For the CN-Du2 and Qianyanzhou sites, only positive or negative VPD–GPP covariance and Δcorr were found. For the Damshung, CN-Cng, and Huailai sites, we found more than one SM value where VPD–GPP covariance or Δcorr was zero. Along with surface soil wetting, there was a change in VPD–GPP covariance and Δcorr from positive to negative at these sites, inconsistently with the transition from water to energy limitation, indicating that CSM was not identifiable. Differently from the above, VPD–GPP covariance had the optimal CSM value that agreed best with the EF–SM-derived CSM at the Xilingele and CN-Sw2 sites. Through another technique, Δcorr was better than VPD–GPP covariance at the Miyun site. At these sites, VPD–GPP covariance and Δcorr changed from negative (water limit) to positive (energy limit). Therefore, VPD–GPP covariance and Δcorr had the potential to obtain large-scale CSM.

Figure 4 Spatial pattern of wet binary bit number at 10 cm depth using covariance between vapor pressure deficit (VPD) and gross primary production (GPP) from (a) GOSIF, (b) GLASS, (c) PML, and (d) TL. Spatial pattern of wet binary bit number at 10 cm depth using correlation difference metric with Kendall's rank correlation between the detrended anomaly of soil moisture (SM) and evapotranspiration (ET) from (e) CAMELE, (f) GLASS, (g) PML, and (h) SEBAL and the correlation between the detrended anomaly of VPD and those ET products. (i) Disparity of soil moisture regimes among all methods and (j) percentage of area with p < 0.05. ZGE: Zhungaer basin, PR: Pearl River basin, YTR: Yangtze River basin, SWR: Southwestern River basin, TR: Tarim basin, SR: Songhua River basin, CT: Changthang region, NM: Inner Mongolian Plateau region, LR: Liaohe River basin, YR: Yellow River basin, HR: Huaihe River basin, HC: Hexi Corridor region, HAR: Haihe River basin, SER: Southeastern River basin, QB: Qaidam basin.

3.3 Spatial pattern of CSM

The number of wet binary bits was used to quantify the agreement among eight ET- and GPP-based models at 10 cm soil depth. If CSM was identified, SM wetter than CSM was represented as 1, and 0 was used for others. If CSM was not identified within a year, digits of the mode were treated as 0. If CSM was not detected for all 18 years, it was displayed as empty. The intercomparison provided helpful insights to examine consistency and discrepancies between multi-source ET and GPP products in depicting the spatial distribution of CSM. Figure 4 shows strong disparities in northern and central China, especially in the Inner Mongolian Plateau region, Songhua River basin, Yangtze River basin, and Yellow River basin. In these regions, the chi-square test showed significant differences among GPP-based models due to their large number of wet binary bits. In addition, TL GPP displayed no CSM value in northwestern China. Note that the SM wetter than CSM showed agreement in eastern and southern basins, such as the Huaihe River basin, Liaohe River basin, Southeastern River basin, and Pearl River basin, indicating that ET- and GPP-based models were consistent in these basins.

Figure 5 Spatial pattern of critical soil moisture (CSM) at 10 cm depth using covariance between vapor pressure deficit (VPD) and gross primary production (GPP) from (a) GOSIF, (b) GLASS, (c) PML, and (d) TL and CSM using correlation difference metric with Kendall's rank correlation between the detrended anomaly soil moisture (SM) and evapotranspiration (ET) from (e) CAMELE, (f) GLASS, (g) PML, and (h) SEBAL and the correlation between the detrended anomaly VPD and those ET products. (i–w) The basin-average values. ZGE: Zhungaer basin, PR: Pearl River basin, YTR: Yangtze River Basin, SWR: Southwestern River Basin, TR: Tarim Basin, SR: Songhua River basin, CT: Changthang region, NM: Inner Mongolian Plateau region, LR: Liaohe River basin, YR: Yellow River basin, HR: Huaihe River basin, HC: Hexi Corridor region, HAR: Haihe River basin, SER: Southeastern River basin, QB: Qaidam basin.

Table 3 Site-based CSM from EF–SM and grid-based CSM using satellite-based ET and GPP and 10 cm depth SM.

Figure 5 shows the spatial distribution of CSM obtained from the covariance between VPD and GOSIF, GLASS, PML, and TL GPP; the correlation difference metric with Kendall's rank correlation between the detrended anomaly of CAMELE, GLASS, PML, and SEBAL ET; and 10 cm soil depth SM and the correlation between the detrended anomaly of ET and VPD. Geographically, these spanned large swaths of land through water-scarce desert regions and lush, rainy forests. Overall, spatial patterns of CSM obtained through four ET products were consistent with those from four GPP products, showing a decreasing variation from south to north. Specifically for water resource subregions, CSM in the semi-humid Huaihe River basin, Haihe River basin, and Yellow River basin was about 0.3 m³ m⁻³ and increased to approximately 0.4 m³ m⁻³ in the Southeastern River basin and Pearl River basin. In addition, Table 3 shows the comparison of site-based CSM from EF–SM-based and grid-based CSM. It was found that gridded CSM values at the CN-Cng, Miyun, and Huailai sites were generally consistent with site-based values. Gridded data had spatial continuity, while site observations showed significant differences in CSM, even between adjacent sites (e.g., CN-Du2 value of 0.113 m³ m⁻³ and Miyun value of 0.274 m³ m⁻³), resulting in inconsistent CSM between satellite- and site-based values.

Figure 6 Soil moisture (SM) at 10, 20, 80, and 100 cm soil depths and critical soil moisture (CSM) derived from CAMELE, GLASS, PML, and SEBAL ET at corresponding soil depths for (a) cropland; (b) forests; (c) grassland; (d) water; (e) ice; (f) urban; and (g) barren soils with a majority of (h) clay, (i) silt, and sand and with contents (j) less than 60, (k) between 60 % and 70 %, (l) between 70 % and 80 %, (m) between 80 % and 90 %, and (n) higher than 90 %.

Furthermore, large-scale CSM depended on roots pulling water out of the unsaturated soil matrix (Feldman et al., 2019) and varied across vegetation types and soil textures at four soil layers (Fig. 6). With shorter root systems and less vegetation, barren areas showed low CSM. Forest regions displayed a relatively high CSM (e.g., 0.18 m³ m⁻³ using PML ET and 10 cm depth SM). As for soil textures, sand covering the large area was further divided into contents of less than 60 %, 60 %–70 %, 70 %–80 %, 80 %–90 %, and higher than 90 %. Soil with a majority of clay had a wetter CSM than others (e.g., 0.38 m³ m⁻³ using PML ET and 10 cm depth SM); this was to be expected given the fact that clay had a larger negative matric potential compared to coarse soil textures dominated by sand and silt. In summary, fine soils and luxuriant vegetation had wetter CSM. Additionally, a layer-wise CSM analysis was conducted to highlight variations in SM properties for different soil layers. It was evident that there were variations in the CSM behavior across layers with higher SM and CSM at 20 cm soil depth. We also found that CSM for grassland and clay was higher than average SM at all four layers, which led to the identification of a large range of SM within water-limited regimes. However, for cropland and forests, differences existed in CSM among four ET-based methods, with higher CSM from GLASS and SEBAL than from others.

Figure 7 Spatial pattern of (a) Δcorr derived from PML ET and 10 cm soil depth soil moisture and (b) the mean absolute percentage error based on partial least square regression for Δcorr estimations. ZGE: Zhungaer basin, PR: Pearl River basin, YTR: Yangtze River basin, SWR: Southwestern River basin, TR: Tarim basin, SR: Songhua River basin, CT: Changthang region, NM: Inner Mongolian Plateau region, LR: Liaohe River basin, YR: Yellow River basin, HR: Huaihe River basin, HC: Hexi Corridor region, HAR: Haihe River basin, SER: Southeastern River basin, QB: Qaidam basin.

3.4 Attribution of water and energy limit shifts

We assessed the spatial pattern of multi-annual average Δcorr at 10 cm soil depth over the period 2001–2018. PML ET was used for Δcorr given the fact that it had the best performance (Sect. 3.1). As shown in Fig. 7a, water-limited regimes were most common in dry and semi-arid areas. Western and northern regions were generally water-limited, while southern regions were energy-limited. The cross-validation using partial least square regression shows that the variance in Δcorr was explained by precipitation, temperature, incoming shortwave radiation, VPD, ET, GPP, and SM, ranging from 73.34 % in the Yangtze River basin to 99.95 % in the Haihe River basin (Fig. 7b).

Figure 8 Spatial pattern of significance (p < 0.05) of (a) Δcorr, (b) precipitation (P), (c) incoming shortwave radiation (SRa), (d) GOSIF gross primary production (GPP), (e) PML evapotranspiration (ET), (f) soil moisture (SM), (g) temperature (Ta), and (h) vapor pressure deficit (VPD) during the period of Δcorr detection using the Mann–Kendall test (Mann, 1945; Kendall, 1948). “De.” means “decreasing” and “In.” means “increasing”. (i) Attribution of Δcorr variations. Colors indicate the variable that best predicts Δcorr dynamics. ZGE: Zhungaer basin, PR: Pearl River basin, YTR: Yangtze River basin, SWR: Southwestern River basin, TR: Tarim basin, SR: Songhua River basin, CT: Changthang region, NM: Inner Mongolian Plateau region, LR: Liaohe River basin, YR: Yellow River basin, HR: Huaihe River basin, HC: Hexi Corridor region, HAR: Haihe River basin, SER: Southeastern River basin, QB: Qaidam basin.

Variations in dominant factors underlined the relevance of hydrological, meteorological, and ecological variables in inducing interannual changes in Δcorr. As shown in Fig. 8, blue pixels represented the significant decrease in Δcorr, indicating increased water stress and correlation between ET and SM. Several typical regions had relatively large areas of significant decreases in Δcorr, such as the Changthang region (2.62 %) and the Tarim basin (3.49 %). ET was the most important predictor across 42 % of the Changthang region and 24 % of the Tarim basin, which confirmed that increasing ET pushed increased water stress in these regions. In addition, a significant increase in VPD in the west might cause drought, especially in the Tarim basin, where VPD was the most important predictor across 19 % of the area. For the Haihe River basin, decreasing SM contributed to increased water limitation. On the contrary, 4.29 % of the Hexi Corridor region showed significant increases in Δcorr; increasing SM contributed to decreased water stress. A total of 16.65 % of the Songhua River basin showed significant increases in Δcorr; decreased water limitation was associated with increasing GPP (greening) in these regions. Moreover, ET and VPD played the most important role in 30 % and 24 % of the Pearl River basin, respectively; the significant decrease in VPD mitigated drought in these regions.

4 Discussion

Analysis of the spatial patterns of CSM using multi-source satellite-based water and carbon fluxes (Fig. 2) derived from different methods (Fig. 3) further enables us to effectively reflect variations from energy to water limitation in spatiotemporally continuous grid cells. To address the question of how soil textures and plant features define constraints imposed by water supply and energy availability, there has been a growing focus on CSM from site to continental scales. For instance, northern California exhibits CSM of 0.15 m³ m⁻³ in semi-arid grassland at the site scale (Baldocchi et al., 2004); CSM using satellite-based surface temperature diurnal amplitude in the semi-arid grassland of Africa has been reported to be 0.12 m³ m⁻³ at the continental scale (Feldman et al., 2019). For specific plants, CSM is around 0.238 m³ m⁻³ using PML ET in the Inner Mongolian Plateau region (Fig. 5), where grass is abundant. This is in line with the grassland CSM of 0.214 m³ m⁻³ from the covariance approach across 195 global sites from the Integrated Carbon Observation System, AmeriFlux, and FLUXNET2015 (Fu et al., 2022a). Another study based on the correlation difference method using SM from the European Space Agency Climate Change Initiative program and ET from FLUXCOM reported large-scale CSM of around 0.21 m³ m⁻³ throughout Europe across all grid cells (Denissen et al., 2020). Researchers also found that plants exhibit a great vertical water uptake range to alleviate the impact of water stress (Gallego-Elvira et al., 2016), with water uptake extending to below 50 cm (Case et al., 2020) or 1–2 m (Tumber-Davila et al., 2022). Low CSM may be attributed to shorter rooting systems in water-limited environments (Konings and Gentine, 2017), while locations with high humidity, such as tropical West Africa and the southern part of the Congo Basin (Feldman et al., 2022), exhibit high CSM. Deep-rooted forests can better regulate their response to drought with high CSM among soil layers, which means that root systems of plants play a key part in determining water- and energy-limited regimes and may help us understand regional- or continental-scale water- and energy-limited regimes that arise from different vegetation and soil conditions. To comprehend the underlying factors driving CSM, it is necessary to do a more comprehensive analysis of climate and ecosystem conditions: CSM detection shows that grassland had a large range of SM within water-limited regimes (Fig. 6), where CSM was higher than average SM, probably because of shallow root systems affected by moisture; therefore, facing decreased Δcorr (Fig. 8), the grassland located in the northwestern arid region was more vulnerable. Further, water-limited regions exhibit great sensitivity in hydrologic cycles to variations in vegetation functioning, climate variability, and catchment physical conditions. Consequently, water-limited vegetation exhibits a higher degree of sensitivity to surface disturbances compared to locations with higher levels of precipitation. In this scenario, the effect of ET is more pronounced, resulting in a decline in energy limitation, such as in the Tarim basin. However, this study focusing on a specific time of year may not be enough to explain the critical value that may be shown throughout the rest of the year. Since CSM values in some grids are not detected by the eight products considered here, further research is needed for the CSM that may appear throughout the rest of the year in different regions. In addition, to compare the performance of multi-source remotely sensed water and carbon fluxes, all data are unified into the 8 d resolution. Therefore, a more refined timescale, such as a 1 d scale study, is also needed.

Multiple factors contribute to inherent constraints in identifying different regimes associated with the utilization of multi-source satellite-based ET and GPP. For example, ET and GPP exhibit great uncertainties (Liu et al., 2021) in areas with barren land, as indicated in Sect. 3.1. In eastern and southern regions, where satellite-based methods were more reliable, the eight considered satellite-based SM regimes were in good agreement (Fig. 4). Since the CAMELE ET was combined with PML ET, they showed consistency in cropland and forests, with a lower CSM than GLASS and SEBAL (Sect. 3.3). By considering variations in energy and water limitations in terrestrial ecosystems (Sect. 3.4), there is potential to improve the water and carbon flux estimation in turn. In addition, SM from ground samplings and gridded sources (Koster et al., 2009) contributes to the uncertainty in characterizing CSM, as discussed in Sect. 3.3. For gridded SM, surface climate shows a significant effect on the upper-soil-layer SM modeling, while the background aridity leads to low variability in the deeper-layer SM (Q. Li et al., 2022). Besides, external forcings seem to be responsible for a shift towards enhanced land–atmosphere coupling (Zhang et al., 2020). It should be noted that the South–North Water Transfer Project and the Pinglu Canal Project in China will result in significant modifications to SM characteristics, which are fundamental components of the concept known as CSM. Water management measures may reduce water stress in grasslands affected by climate change and make southern coastal clay areas more resistant to possible disturbances. Overall, our research can inform large-scale water conservancy projects for better water resource allocation from the perspective of the critical effect of SM.

5 Conclusion

Our main accomplishment is observing and identifying water and energy limit shifts using multi-source satellite-based water and carbon fluxes over China. These shifts show which areas are more likely to be affected by climate change. To do this, we first examined the consistency of ET and GPP derived from the site- and satellite-based grid observations and the consistency of CSM derived from the EF–SM, covariance, and correlation difference methods. CSM detected by the covariance between VPD and GPP and CSM using the correlation difference metric using VPD, ET, and SM matched well with CSM using the EF–SM method at the site scale, suggesting that these methods could detect large-scale CSM. According to satellite-based CSM from four ET products, four GPP products, and the latest SM dataset, surface water- and energy-limited regimes varied among land cover types, soil textures, and water resource subregions; soil textures of clay and land cover types of grassland had a large range of SM within water-limited regimes. Based on the spatial pattern of CSM, we further identified the dominant factor of Δcorr and discovered that VPD was the most important predictor across 24 % of the Pearl River basin and 19 % of the Tarim basin. However, unlike the declining VPD in the Pearl River basin, the increasing VPD aggravated the water stress in the Tarim basin, especially for the more fragile grassland in these areas. As environmental change and extreme disturbances affect CSM, future research directions will aim to consider the impact of hydraulic projects such as inter-basin water transfers on CSM, the impact of extreme disturbances such as tropical cyclones and wildfires on CSM, and possible changes in CSM.

This study used multi-source satellite-based water and carbon fluxes and different methods to detect CSM, and more efforts were put into the evaluation and validation of CSM. Datasets covering 18 years, as used for CSM, were quite typical of the long-term climatology of continental wetness. Since CSM, an emerging property, is generated by multiple processes occurring on the land surface, in the atmosphere, and at the interface between them, uncertainties in ET and GPP from the algorithm, uncertainties in SM from ground sampling, and enhanced land–atmosphere coupling due to external forcing all contribute to CSM uncertainties. We emphasize that SM behavior below and above CSM determines ET and GPP and that water-limited regimes of the SM range depend on CSM. Water and carbon fluxes are vulnerable to the sensitivity of Δcorr to hydrological, meteorological, and ecological predictors. Accordingly, the water and carbon algorithm should consider water–energy limit shifts to improve the simulation accuracy. Thus, applying our new understanding of Δcorr and CSM under changing land–atmosphere conditions will provide a more complete perspective of the evolution of regional terrestrial ecosystems over extended periods.