Substantial root-zone water storage capacity observed by GRACE and GRACE/FO

Root-zone water storage capacity (S_r) – the maximum water volume available for vegetation uptake – bolsters ecosystem resilience to droughts and heatwaves, influences land–atmosphere exchange, and controls runoff and groundwater recharge. In land models, S_r serves as a critical parameter to simulate water availability for vegetation and its impact on processes like transpiration and soil moisture dynamics. However, S_r is difficult to measure, especially at large spatial scales, hindering an accurate understanding of many biophysical processes, such as photosynthesis, evapotranspiration, tree mortality, and wildfire risk. Here, we present a global estimate of S_r using measurements of total water storage (TWS) anomalies from the Gravity Recovery and Climate Experiment (GRACE) and GRACE Follow-On satellite missions. We find that the median S_r value for global vegetated regions is at least 220±40 mm, which is over 50 % larger than the latest estimate derived from tracking storage change via water fluxes and 380 % larger than that calculated using a typical soil and rooting-depth parameterization. These findings reveal that plant-available water stores exceed the storage capacity of 2 m deep soil in nearly half of Earth's vegetated surface, representing a notably larger extent than previous estimates. Applying our S_r estimates in a global hydrological model improves evapotranspiration simulations compared to other S_r estimates across much of the globe, particularly during droughts, highlighting the robustness of our approach. Our study highlights the importance of continued refinement and validation of S_r estimates and provides a new observational approach for further exploring the impacts of S_r on water resource management and ecosystem sustainability.

1 Introduction

During periods of insufficient precipitation, vegetation relies on water stored underground to survive (Miguez-Macho and Fan, 2021). The larger the root-zone water storage capacity (S_r), the more water the root zone can store during wet periods for use in droughts (Teuling et al., 2006). S_r, therefore, plays an important role in regulating ecosystem resilience to droughts and heatwaves and in affecting wildfire outbreaks and mortality risk (Callahan et al., 2022; Chen et al., 2013; Goulden and Bales, 2019; Hahm et al., 2019; Humphrey et al., 2018; Stocker et al., 2023). It is also an essential parameter for modeling plant carbon uptake, transpiration, soil evaporation, streamflow, and groundwater (Maxwell and Condon, 2016; Zhao et al., 2022; Peterson et al., 2021). Despite its critical role in modulating the carbon and water cycles, global patterns of S_r remain poorly characterized.

S_r is typically calculated as the integration of plant rooting depth and soil-texture-dependent water-holding capacity (Seneviratne et al., 2010; Vereecken et al., 2022; Speich et al., 2018; Federer et al., 2003). However, this approach (hereafter referred to as the rooting-depth-based estimation) suffers from uncertainties associated with plant rooting depth and substrate hydraulic properties, particularly at depth, both of which undermine the accuracy of the calculated S_r (Vereecken et al., 2022; Novick et al., 2022). Moreover, this approach assumes a static root zone confined to the near-surface unsaturated soil layer. However, recent studies have shown that this assumption is not always accurate. In many ecosystems, plant roots can penetrate beyond the shallow soil layer into weathered bedrock, accessing rock moisture and tapping into groundwater, especially during prolonged dry periods (Li et al., 2015; Hahm et al., 2020; McCormick et al., 2021; Rempe and Dietrich, 2018; Maxwell and Condon, 2016; Fan et al., 2017; Baldocchi et al., 2021). Thus, the rooting-depth-based estimation may significantly underestimate S_r.

More recently, Earth observations of precipitation (P) and evapotranspiration (ET) have been used to estimate S_r. Several studies (Stocker et al., 2023; Wang-Erlandsson et al., 2016; Gao et al., 2014; McCormick et al., 2021) have proxied S_r using the maximum cumulative difference in ET and P during dry periods (when ET > P), which reflects the largest water volume that an ecosystem has withdrawn from its root zone. This method (hereafter referred to as the water-deficit-based estimation) is based on mass balance and thus eliminates the need for assumptions about plant access to rock moisture and groundwater, rooting depth, and soil and bedrock hydraulics. However, obtaining accurate P and ET data is challenging at scale (Sun et al., 2018; Miralles et al., 2016), and errors in these data can accumulate and deteriorate S_r calculations. Here, to avoid this shortcoming, we estimated root-zone storage dynamics directly from total water storage (TWS) anomalies measured by the Gravity Recovery and Climate Experiment (GRACE) and GRACE Follow-On (GRACE-FO) satellite missions (hereafter GRACE/FO). With these observations, we characterized global patterns of S_r and found that both the rooting-depth-based estimate and the water-deficit-based estimate have significantly underestimated S_r.

2 Materials and methods

2.1 GRACE/FO TWS

We use monthly measurements of the TWS anomaly from GRACE for the years 2002–2017 and from GRACE-FO for the years 2018–2022. These measurements were obtained from the Jet Propulsion Laboratory (JPL) RL06 solutions (Watkins et al., 2015; Wiese et al., 2016), which provide monthly average anomalies of the gravity field over an equal-area 3° × 3° mass concentration block (mascon). We opted for the JPL mascon solutions because each JPL mascon is relatively uncorrelated with neighboring mascons and thus offers more localized spatial variations than other mascon solutions and the spherical harmonic solutions (Watkins et al., 2015; Wiese et al., 2016). We did not fill the 11-month gap (July 2017 to May 2018) between GRACE and GRACE-FO. However, we linearly interpolated other missing months from the nearest previous and subsequent non-missing values (Rodell et al., 2018; Zhao et al., 2021). Because we aimed to estimate root-zone storage capacity S_r, we only included mascon locations with over 50 % fractional vegetation cover based on the land cover product (MCD12Q1) version 6.1 from the Moderate Resolution Imaging Spectroradiometer (MODIS) (Sulla-Menashe and Friedl, 2018).

2.2 S_r from TWS drawdown and uncertainty estimate

Ecosystem use of land water storage for ET is reflected in TWS drawdowns, which are consecutive declines in TWS despite seasonal or intermittent recharge and after accounting for long-term trends due to anthropogenic groundwater use. An example is illustrated in Fig. 1 at a mascon location in southern Idaho, where the largest TWS drawdowns are annotated. From the water balance, a TWS drawdown over a time period Δt is equal to

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }\mathrm{TWS}=\sum P-\sum \mathrm{ET}-\sum R,\end{array}$$

where ∑P, ∑ET, and ∑R are the total precipitation, total evapotranspiration, and net runoff out of the mascon over Δt, respectively. Based on Eq. (1), when precipitation exceeds runoff ($\sum P-\sum R\mathit{>}\mathrm{0}$), any TWS drawdown (or negative ΔTWS) must be influenced by a change in storage due to ET. To determine if precipitation exceeds runoff during GRACE/FO-observed TWS drawdowns, we compared R estimates from a multi-forcing observation-based global runoff reanalysis (Ghiggi et al., 2021) to P estimates from the Global Precipitation Climatology Project (Gebremichael et al., 2003). We found that, in nearly all analyzed mascon locations, the average P−R is positive during, at least, the five largest TWS drawdowns (Fig. A1 in the Appendix), confirming that these TWS drawdowns reflect root-zone water storage transpired by ecosystems and not loss of water in the mascon due to runoff.

Figure 1 TWS time series showing the three largest drawdowns at a mascon location in southern Idaho.

We estimated root-zone water storage capacity S_r to be the largest TWS drawdown during the record period of GRACE/FO (denoted as $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$). To avoid overestimating S_r, we removed the impact of groundwater pumping, snow, and surface water on TWS drawdowns. Anthropogenic groundwater use often manifests as a negative long-term trend in the TWS time series (Rodell et al., 2018; Rodell et al., 2009; Feng et al., 2013). For example, regions showing significant decreasing trends in TWS largely coincide with well-known groundwater irrigation areas identified in AQUASTAT data (Fig. A2). To avoid conflating this drawdown with S_r, we first calculated the TWS trend by simultaneously fitting an annual and a semi-annual signal, a linear trend, and a constant to the GRACE/FO time series (Fig. A2). Then, we assumed that any negative trend was attributable to groundwater pumping and removed the negative trend from the original GRACE/FO time series before calculating the TWS drawdowns. In high-latitude and mountainous regions, the maximum TWS anomaly during drawdowns may include snow. To avoid attributing snow storage to root-zone water storage, we first determined the largest drawdown from the full GRACE/FO time series and then calculated S_r using the maximum and minimum TWS anomaly with a monthly mean air temperature above 5 °C. We obtained air temperature data from the fifth-generation European Centre for Medium-Range Weather Forecasts atmospheric reanalysis of the global climate (ERA5) (Hersbach et al., 2020). Following Wang et al. (2023a), we used total runoff from Ghiggi et al. (2021), which includes both surface runoff and subsurface runoff, as a proxy for surface water storage change (i.e., ΔSW =R) and removed it from TWS drawdowns to isolate Δ SW contributions (water stored in rivers, lakes, and reservoirs) to the GRACE/FO signal. This approach assumes that (1) R contributes directly to an increase in surface water levels within the drainage network, and (2) it takes approximately 1 month for R to exit the drainage system, aligning with the monthly time step of GRACE/FO data. Note that total runoff from Ghiggi et al. (2021) stopped in 2019, and we used monthly climatology values between 2002 and 2019 to extend the data to 2022 and to align with the GRACE/FO record length. Other contributions to TWS drawdowns, such as changes in water intercepted by leaf and branch surfaces and internal plant water storage, are too small to be detected by GRACE/FO (Rodell et al., 2005) and are unlikely to significantly affect our estimates. Our method also implicitly includes moisture stored in the topsoil for soil evaporation (Stoy et al., 2019). However, the contribution of soil evaporation to ET decreases quickly as TWS draws down (Stocker et al., 2023), and we expect that the magnitude of the largest drawdown will be determined by the root-zone depletion magnitude reflected at the end of the drawdown.

We calculated the random error of $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ by adding the errors of the two GRACE/FO measurements and the uncertainty of groundwater pumping and surface water signals in quadrature. To calculate the GRACE/FO measurement error, we used the formal error product provided by the JPL mascon solutions (Watkins et al., 2015; Wiese et al., 2016). For the uncertainty of groundwater pumping and surface water signals, we assumed a ±50 % error based on the magnitude of our calculated signals, following Zhao et al. (2021). This assumption implies that the uncertainty range is equal to the signals themselves, leading to what is likely to be a conservative error estimate.

2.3 Comparison to other S_r estimates

We compared our $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ estimate to two other S_r datasets. These datasets represent the typical approach covering rooting depth and soil-texture-dependent water-holding capacity (referred to as $S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$) and the water deficit accumulation approach (referred to as $S_{\mathrm{r}}^{\mathrm{accum}}$). We chose the $S_{\mathrm{r}}^{\mathrm{accum}}$ estimate from Stocker et al. (2023) because it used the latest Earth-observation-constrained estimates of precipitation and evapotranspiration. We used their S_CWDX80 product, which was estimated based on cumulative water deficit extremes occurring with a return period of 80 years. We calculated $S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$ using existing datasets on rooting depths and soil texture. The $S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$ approach requires knowing effective rooting depths (Federer et al., 2003; Speich et al., 2018; Stocker et al., 2023; Bachofen et al., 2024). We obtained effective rooting depths from Yang et al. (2016), who derived them using an analytical model that balances the marginal carbon costs and benefits of deeper roots. While such model-based datasets are valuable for providing comprehensive coverage and insights into complex processes, they do not incorporate direct observational data for validation or correction. Soil water-holding capacity, defined as the difference between field capacity and permanent wilting point, is calculated based on soil texture information from the Harmonized World Soil Database version 1.2 (Wieder et al., 2014) and pedo-transfer functions based on Balland et al. (2008). The Harmonized World Soil Database provides information for depths of 0–0.3 and 0.3–1 m. For depths greater than 1 m, we assume texture values from the 0.3–1 m depth following Stocker et al. (2023). For consistency, we spatially averaged both $S_{\mathrm{r}}^{\mathrm{accum}}$ and $S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$ estimates to match the GRACE/FO spatial scale (3° × 3°).

2.4 Evaluation using the USGS monthly hydrologic model

Validating large-scale S_r remains inherently difficult because direct measurement of S_r is challenging. Previous studies have primarily employed two indirect methods: comparison to measured rooting depths and hydrological modeling. Stocker et al. (2023) converted their deficit-based S_r estimates (∼ 5 km resolution) into rooting depths using soil texture and water-holding-capacity parameters and then compared them to field rooting-depth measurements aggregated at biome levels to mitigate the scale mismatch. However, this approach is not suitable for our study as GRACE/FO-derived S_r (∼ 300 km) encompasses multiple biome types within the effective resolution of GRACE/FO data, making biome-level aggregation less meaningful. Additionally, the rooting-depth method overlooks groundwater and rock moisture contributions to S_r, which Stocker et al. (2023) found to be significant in over half of their root measurement sites. This omission will likely become more critical at the spatial scale of GRACE/FO, which averages larger areas and includes more diverse biome types. These factors make the rooting-depth comparison unsuitable for evaluating GRACE/FO-derived S_r. Wang-Erlandsson et al. (2016) used deficit-based S_r estimates in a simple hydrological model and assessed improvements in simulating hydrologic time series. While this approach better aligns with the scale of GRACE/FO, it is constrained by the limited availability of high-quality global hydrologic data. This can lead to a circular use of the same data for both S_r estimation and model evaluation, as seen when Wang-Erlandsson et al. (2016) used satellite-based ET data for both S_r estimation and model evaluation, reducing the independence of the validation process. We also note that existing gridded ET products generally have assumed ecosystem responses to water stress in their algorithms and thus are highly uncertain (Miralles et al., 2016).

To address these challenges, we evaluated the relative accuracy of $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$, $S_{\mathrm{r}}^{\mathrm{accum}}$, and $S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$ by separately parameterizing a hydrological model with each estimate; these are referred to as HydroModel($S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$), HydroModel($S_{\mathrm{r}}^{\mathrm{accum}}$), and HydroModel($S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$). We then assessed their accuracy in simulating ET using an independent dataset: version 4.1 of the Global Land Evaporation Amsterdam Model (GLEAM) ET (Miralles et al., 2025a). This dataset was not involved in the calculation of $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$, $S_{\mathrm{r}}^{\mathrm{accum}}$, or $S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$, ensuring independence and avoiding the circular validation that affected previous studies (e.g., Wang-Erlandsson et al., 2016). Furthermore, the GLEAM ET product provides several key improvements over other gridded ET products. For example, it combines hybrid learning from eddy covariance and sap flow to better capture vegetation responses to drought (Koppa et al., 2022) and explicitly accounts for plant access to groundwater (Hulsman et al., 2023). The atmospheric-forcing data and model parameters were identical across simulations, with S_r being the only variable parameter. Therefore, differences in model performance reflect the relative accuracy of the three S_r estimates. A monthly hydrologic model developed by the United States Geological Survey (USGS) (McCabe and Markstrom, 2007) was used due to its simplicity and transparency with regard to physical processes. Specifically, the model relies on a straightforward specification of S_r as a “water bucket” depth rather than an indirect specification through prescribed rooting depth, soil texture, and pedo-transfer functions across the profile. This allows us to parameterize the model directly with $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$, $S_{\mathrm{r}}^{\mathrm{accum}}$, and $S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$. The USGS model was run at each GRACE mascon location with air temperature forcing from ERA5 and precipitation forcing from the GPCP. We used climate forcing from 1993 to 2001 to spin up the model and performed water cycle simulations for the study period from 2002 to 2022. No calibrations were carried out.

To mitigate the impact of possible biases embedded in GLEAM ET and in the forcing data and those caused by model uncertainty (as the USGS model is uncalibrated), we used standardized ET anomalies (i.e., z scores) as the target of validation and focused on assessing whether S_r improves the temporal dynamics of ET simulations (i.e., seasonal and interannual variations) rather than the absolute values of ET. The Nash–Sutcliffe model efficiency (NSE) coefficient was used to assess the predictive skill of each USGS hydrologic model, defined as follows:

$$\begin{array}{}\text{(2)}& \mathrm{NSE}=\mathrm{1}-{\displaystyle \frac{{\sum }_{t=\mathrm{1}}^T(X_{\mathrm{o}}^t-X_{\mathrm{m}}^t{)}^{\mathrm{2}}}{{\sum }_{t=\mathrm{1}}^T(X_{\mathrm{o}}^t-\stackrel{\mathrm{‾}}{X_{\mathrm{o}}}{)}^{\mathrm{2}}}},\end{array}$$

where X represents the standardized ET anomaly; $\stackrel{\mathrm{‾}}{X_{\mathrm{o}}}$ is the mean of the observed X; and $X_{\mathrm{o}}^t$ and $X_{\mathrm{m}}^t$ are the observed and modeled X at time t, respectively (Nash and Sutcliffe, 1970). An NSE value closer to 1 indicates a better model performance in simulating X, while an NSE value less than 0 indicates that the mean observed value is a better predictor than the simulated value, suggesting an unsatisfactory model performance (Nash and Sutcliffe, 1970). If HydroModel($S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$), HydroModel($S_{\mathrm{r}}^{\mathrm{accum}}$), and HydroModel($S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$) all yield negative NSE values, the efficacy of using the USGS hydrologic model to evaluate the relative accuracy of the three S_r estimates is compromised.

3 Results

3.1 S_r from GRACE/FO ($S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$)

We find a substantial root-zone water storage capacity worldwide. Across the global vegetated domain, $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ (or the largest TWS drawdown) spans from 22 to 2131 mm (Fig. 2a). The distribution of $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ is positively skewed, with a median value of 221 mm (129–389 mm interquartile range; note that, hereafter, values in parentheses always refer to the interquartile range). Larger $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ values are associated with densely vegetated regions like tropical rainforests, the southeastern US, the Pacific Northwest, and the southern part of China. By contrast, smaller $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ values are found in sparsely vegetated regions like central Asia, much of Australia, and some Arctic regions (Fig. 2a). Figure 2b shows the duration of the maximum TWS drawdown, with a global median of 2.8 years (1.6–5.2 years). We find no correlation between the duration and the magnitude of the largest TWS drawdown across different regions (Fig. 2a–b). The impact of random error sources on our $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ estimate remains moderate, with a global median relative error of 18 % (13 %–26 %) (Fig. 2c).

Figure 2 S_r estimated from GRACE/FO total water storage (TWS) anomaly. (a) Global patterns of $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ for Earth's vegetated regions. (b) The duration of the maximum TWS drawdown. (c) Global patterns in the random error of $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$. Insets in (a)–(c) show the histograms of corresponding mapping variables across our study area. White spaces on land represent mascon locations with less than 50 % vegetation cover.

To characterize the utilization of root-zone water storage capacity, we compared the second- and third-largest TWS drawdowns to $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$. We find that, on average, the second-largest TWS drawdown consumes 83 % (71 %–92 %) of the $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ estimate (Fig. 3a), while the third-largest uses 68 % (54 %–82 %) (Fig. 3b). The average duration of the second- and third-largest TWS drawdowns decreases from 1.6 years (1.1–3.2 years) to 1.2 years (0.5–1.7 years) (Fig. 3c–d). In about 40 % of our analyzed mascons, the longest TWS drawdown period does not coincide with the largest drawdown magnitude. These findings underscore the nuanced dynamics of water storage use within the root zone, suggesting variability in both magnitude and duration across different regions.

Figure 3 Utilization of root-zone water storage capacity. Panels (a) and (b) are consumption percentages of $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ during the second- and third-largest TWS drawdowns. Panels (c) and (d) are the durations of the second- and third-largest TWS drawdowns. Insets in (a)–(d) show the histograms of corresponding mapped variables.

3.2 Comparison with other S_r estimates

Our $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ estimate is larger than those of $S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$ and $S_{\mathrm{r}}^{\mathrm{accum}}$ over much of the globe. Figure 4a–b show the difference of $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ compared to $S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$ and $S_{\mathrm{r}}^{\mathrm{accum}}$, respectively. Across the global vegetated domain, $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ surpasses $S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$ in over 90 % of mascon locations, with a median value that is 175 mm (or 380 %) higher than that of $S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$. The $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ exceeds $S_{\mathrm{r}}^{\mathrm{accum}}$ in over 70 % of the study area, with a median value that is 77 mm (or 53 %) higher than that of $S_{\mathrm{r}}^{\mathrm{accum}}$, despite exhibiting lower values in many regions of Africa, India, Mexico, and northeastern Brazil (Fig. 4b). Notably, these differences are greater than the random error of $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$, emphasizing that the underestimations by $S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$ and $S_{\mathrm{r}}^{\mathrm{accum}}$ are significant.

Figure 4 $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ values are notably larger than those of other datasets over much of the globe. (a) The difference between $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ and $S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$. (b) The difference between $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ and $S_{\mathrm{r}}^{\mathrm{accum}}$.

3.3 Implementation in the USGS hydrologic model

To assess whether $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ is an improvement over $S_{\mathrm{r}}^{\mathrm{accum}}$ and $S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$, we used each of them to separately parameterize the USGS hydrologic model. We first evaluated the accuracy of HydroModel($S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$), HydroModel($S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$), and HydroModel($S_{\mathrm{r}}^{\mathrm{accum}}$) in replicating the full time series of standardized GLEAM ET anomalies. For over 95 % of the global vegetated domain, at least one model achieved a positive NSE value. In these regions, the average NSE for HydroModel($S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$) is 0.73 (0.65–0.89), that for HydroModel($S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$) is 0.69 (0.63–0.86), and that for HydroModel($S_{\mathrm{r}}^{\mathrm{accum}}$) is 0.72 (0.64–0.87). These relatively high NSE values indicate that the USGS model is effective in simulating ET. While the global average NSE values for the three models are similar, HydroModel($S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$) demonstrates slightly superior performance, outperforming HydroModel($S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$) in 66 % of the vegetated regions and HydroModel($S_{\mathrm{r}}^{\mathrm{accum}}$) in 59 % of these regions (Fig. 5).

Figure 5 $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ improves overall model performance in simulating standardized ET anomalies over much of the globe. (a) The NSE difference between HydroModel($S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$) and HydroModel($S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$) for full time series. (b) The NSE difference between HydroModel($S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$) and HydroModel($S_{\mathrm{r}}^{\mathrm{accum}}$) for full time series. Gray areas indicate regions where all models fail to achieve a positive NSE value.

We hypothesized that a more accurate S_r would have a greater impact on improving ET simulations during drought periods when ET is more dependent on deep subsurface water storage. To test this, we calculated NSE values specifically for drought periods, defined as when the 3-month standardized precipitation index was less than −1.2, indicative of severe drought conditions (McKee et al., 1993). Across 87 % of the global vegetated domain, at least one model achieved a positive NSE value. In these regions, the average NSE for HydroModel($S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$) is 0.65 (0.52–0.86), that for HydroModel($S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$) is 0.52 (0.40–0.81), and that for HydroModel($S_{\mathrm{r}}^{\mathrm{accum}}$) is 0.61 (0.48–0.84). These lower NSE values compared to the full ET time series reflect the challenges faced by the USGS model in simulating ET during droughts, consistently with previous findings (e.g., Zhao et al., 2022). However, HydroModel($S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$) showed notable improvement over the other two models, particularly in high-latitude regions (Fig. 6). These results suggest that, while $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ provides only marginal improvements for the full time series of standardized GLEAM ET anomalies, its superiority over $S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$ and $S_{\mathrm{r}}^{\mathrm{accum}}$ becomes more pronounced during drought conditions.

Figure 6 $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ improves model performance notably in simulating standardized ET anomalies during drought periods across much of the globe. Panels (a) and (b) are similar to Fig. 5a and b, respectively, except they are for drought time periods. Gray areas indicate regions where all models fail to achieve a positive NSE value.

4 Discussion

4.1 Limitations of and uncertainty in $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$

Our $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ estimate provides a conservative lower bound for S_r because the largest TWS drawdown during the GRACE/FO record period may not cover a period during which ET from storage exhausts the entire root-zone water storage capacity, particularly in areas experiencing water accumulation in the root zone due to increases in precipitation. This likely explains why our $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ estimate is lower than that of $S_{\mathrm{r}}^{\mathrm{accum}}$ in northern and eastern Africa, where strong increasing TWS trends were observed (Figs. 3b and A2) in response to increasing precipitation trends (e.g., Rodell et al., 2018). Additionally, our approach to account for groundwater pumping and surface water may overestimate these signals' actual magnitudes and, thus, likely contributes to underestimating S_r. Specifically, we assumed all negative TWS trends to be caused by groundwater withdrawals and removed them from $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$. However, intense groundwater withdrawals are concentrated in specific regions such as northwestern India, California's Central Valley, and the North China Plain (Rodell et al., 2009; Feng et al., 2013; Liu et al., 2022). Consequently, we may have removed TWS depletion trends caused by natural variability, as seen in the drought-stricken southeast of Brazil (Rodell et al., 2018). This likely explains why $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ is lower than $S_{\mathrm{r}}^{\mathrm{accum}}$ there (Fig. 3b). Furthermore, we used total runoff (which includes surface runoff, snowmelt, and groundwater discharge) as a proxy to remove surface water storage change from the TWS drawdown. We used total runoff – as opposed to surface runoff alone (Wang et al., 2023a) – due to observational data availability, though doing so may lead to an overestimation of surface water storage change and, therefore, an underestimation of S_r.

4.2 $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ vs. $S_{\mathrm{r}}^{\mathrm{accum}}$: impact of multi-year drawdowns and drought return period assumption

Despite being conservative, $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ reveals a substantially larger volume of root-zone water storage capacity than $S_{\mathrm{r}}^{\mathrm{accum}}$. One reason for this discrepancy may be the lack of interannual storage variability considered in the $S_{\mathrm{r}}^{\mathrm{accum}}$ calculation (Stocker et al., 2023). Although Stocker et al. (2023) used a cumulative water deficit approach to infer root-zone water storage drawdown, akin to our TWS drawdown approach, they found that the annual totals of P exceeded those of ET at almost all locations. The fact that their method resets the calculation whenever accumulated P − ET is positive suggests that their method was generally unable to account for carryover storage and multi-year drawdowns of root-zone storage. Our use of GRACE/FO TWS, which allows for multi-year drawdowns, is supported by recent observations (Goulden and Bales, 2019; McCormick et al., 2021; Pérez-Ruiz et al., 2022; Peterson et al., 2021; Scott and Biederman, 2019) and modeling efforts (Miguez-Macho and Fan, 2021; Livneh and Hoerling, 2016) suggesting widespread carryover storage effects. Our calculations of $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ found that the largest TWS drawdown period lasted a median of 2.8 years, with an interquartile range between 1.6 and 5.2 years (Fig. 2c). Even the second- and third-largest TWS drawdowns had a median duration of more than 1 year globally (Fig. 3c–d). These findings align with the results reported in the previously referenced studies on carryover storage effects.

The discrepancy between $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ and $S_{\mathrm{r}}^{\mathrm{accum}}$ is further influenced by differences in how drought return periods are defined. $S_{\mathrm{r}}^{\mathrm{accum}}$ values reported in Stocker et al. (2023) are statistically scaled to represent an 80-year drought return period rather than being directly derived from observed drought events. This extrapolation assumes a fixed probability distribution of drought occurrence, which may not fully capture real-world hydrological variability. In contrast, $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ is based on observed multi-year TWS drawdowns from the period of 2002–2022, directly reflecting droughts that occurred over the past 2 decades. A more comparable approach would require using the unscaled $S_{\mathrm{r}}^{\mathrm{accum}}$, which was derived using Earth observations of precipitation and ET from the period of 2003–2018 without statistical adjustments. However, the unscaled $S_{\mathrm{r}}^{\mathrm{accum}}$ is not publicly available. Given that the 80-year return period scaling inflates $S_{\mathrm{r}}^{\mathrm{accum}}$ values and that $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ is still substantially larger, it follows that the unscaled $S_{\mathrm{r}}^{\mathrm{accum}}$ would be even lower, further supporting the conclusion that the $S_{\mathrm{r}}^{\mathrm{accum}}$ approach underestimates root-zone storage capacity compared to $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$.

4.3 Groundwater and rock moisture in $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ and differences compared to $S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$

The $S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$ estimate falls notably below that of both $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ and $S_{\mathrm{r}}^{\mathrm{accum}}$. This discrepancy may be attributed to the $S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$ approach ignoring plant access to moisture stored beneath the soil, such as in weathered and fractured bedrock and groundwater. These deep moisture sources are known to affect ET significantly, thus contributing to S_r (e.g., Fan et al., 2017; Rempe and Dietrich, 2018; McCormick et al., 2021). Unlike $S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$, the definitions of $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ and $S_{\mathrm{r}}^{\mathrm{accum}}$ incorporate natural variability in these deep moisture reserves, broadening the traditional root-zone concept beyond the unsaturated soil layer. This expanded definition acknowledges the dynamic nature of the root zone, with plants accessing deep groundwater and rock moisture during prolonged droughts and periods of high transpiration demand (Gao et al., 2024). Indeed, root-accessible water does not require roots to physically occupy the entire storage domain. Processes like capillary rise can move deep water upward to the traditional root zone for vegetation transpiration, especially during dry seasons and droughts.

The importance of including groundwater and rock moisture in S_r is well-supported by recent evidence. Studies using in situ groundwater (Fan et al., 2017; Thompson et al., 2011; Baldocchi et al., 2021; Li et al., 2015), remote sensing observations (Koirala et al., 2017; Rohde et al., 2024), and modeling efforts (Miguez-Macho and Fan, 2021; Hain et al., 2015) have demonstrated that plants can access these deep moisture sources and have highlighted their critical role in sustaining ET, especially during extreme droughts. In many ecosystems, water stress can stimulate root growth into deep subsurface layers through the capillary rise effect, with roots extending to the capillary fringe and the water table, as observed in both field and laboratory studies (Naumburg et al., 2005; Orellana et al., 2012; Fan et al., 2017; Kuzyakov and Razavi, 2019). Although individual shallow-rooted plants (e.g., grassland sites) may not directly tap into groundwater or rock moisture, the large spatial scale of GRACE/FO likely captures water uptake across diverse vegetation types. Even in areas primarily covered by shallow-rooted vegetation, deeper-rooted plants within the same GRACE/FO mascon may redistribute water upward through hydraulic redistribution, making it available for shallow-rooted plants to use (e.g., Espeleta et al., 2004; Orellana et al., 2012). In fact, satellite observations have confirmed widespread plant–groundwater interactions at large spatial scales (Koirala et al., 2017), even in dryland regions dominated by grasslands (Rohde et al., 2024; Wang et al., 2023b). Recognizing and incorporating groundwater and bedrock moisture into root-zone storage capacity can enhance our understanding of land–atmosphere interactions (Maxwell and Condon, 2016; Schlemmer et al., 2018; Dong et al., 2022), improve runoff simulations (Hahm et al., 2019), and provide a more accurate representation of vegetation resilience to droughts and heatwaves (Jiménez-Rodríguez et al., 2022; Esteban et al., 2021).

The $S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$ approach, while useful for simplifying root-zone complexity, overlooks critical aspects of root density, its vertical and lateral distribution, and the ability of plants to access deep water stores – factors that have significant implications for understanding ecosystem water uptake and improving land models. For instance, this approach reduces the complexity of rooting systems into a single effective rooting-depth parameter (Federer et al., 2003; Speich et al., 2018), which tends to be shallower than both the maximum rooting depth (Federer et al., 2003) and the depth that contains the upper 95 % of the root biomass (Yang et al., 2016). These deeper layers, however, often play a disproportionately important role in ecosystem water uptake (Fan et al., 2017; Jackson et al., 1999; Bachofen et al., 2024). Additionally, when dividing $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ with the same water-holding capacity (WHC) used in $S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$ to calculate effective rooting depth, this depth exceeds 2 m in nearly 50 % of global vegetated areas, which is in contrast to Yang et al.'s (2016) estimate of 10 % and Stocker et al.'s (2023) estimate of 37 %. These results indicate that the potential for plants to tap into deep water stores is more prevalent than previously understood. For land models that do not explicitly incorporate S_r as a variable, this suggests that models with a soil depth of less than 2 m (e.g., the Noah model within the Global Land Data Assimilation System (GLDAS)) may be unable to accurately simulate these deeper water drawdowns. Consequently, this limitation could impact studies of groundwater that rely on GLDAS to separate soil moisture from TWS (e.g., Rodell et al., 2009).

4.4 Strengths and limitations of $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ validation

Although direct observations of S_r at large spatial scales are limited, our validation effort for $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ shows two notable strengths. First, we used an independent dataset for the validation of USGS models parameterized by different S_r estimates, unlike a previous study (Wang-Erlandsson et al., 2016), which relied on a dataset already used in their S_r calculation. Second, the GLEAM ET dataset used here for validation addresses key limitations of other gridded ET products by using a data-driven embedding of plant–water relationships (rather than explicitly assuming these a priori as most ET products do) and by explicitly accounting for groundwater contributions to ET (Miralles et al., 2025a).

Despite these strengths, our validation effort is not without limitations. First, the mechanistic linkage between S_r and commonly used hydrological indicators (e.g., ET and streamflow) is complex. Identifying decisive indicators that are highly sensitive to S_r is an ongoing research challenge. In this context, our findings provide an initial step towards understanding this relationship, demonstrating that a more accurate S_r improves simulations of drought-time ET anomalies more effectively than all-time variations (Figs. 5 and 6). However, resolving such a complex relationship is further complicated by model structural errors or uncertainties in other model parameters, which can obscure the true impact of accurate S_r parameterization on ecohydrological processes. For example, in our study, streamflow simulated by the USGS model is mainly driven by precipitation and shows limited sensitivity to S_r (results not shown). This aligns with the findings of another simple hydrologic model used by Wang-Erlandsson et al. (2016), as discussed in their open peer review file, where streamflow measurements were also not used for model evaluation. Second, we used standardized ET anomalies (z scores) as the validation target, focusing on temporal dynamics such as seasonal and interannual variations rather than absolute ET values. While this approach effectively mitigates the impact of data biases and ensures consistency, it narrows the scope of the validation.

4.5 Implications for high-resolution land surface models

Despite the coarse resolution of GRACE/FO observations, $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ and our proposed approach remain valuable for improving the operational configuration of higher-resolution land models. First, $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ can be used to evaluate and refine default S_r parameterizations within models once aggregated to the coarse scale of GRACE/FO data in conjunction with other diagnostic analyses. For instance, if a model underestimates ET during droughts in a region where its S_r value is significantly lower than that of $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$, the default S_r value may be increased based on $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ even if the model's resolution is much higher than that of $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$. Second, in the future, our methodology can be extended to downscaled GRACE/FO products, leveraging techniques such as data assimilation systems or artificial intelligence to improve the spatial resolution of $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ (Li et al., 2019; Gou and Soja, 2024).

5 Conclusions

We used GRACE/FO TWS observations to estimate root-zone water storage capacity (S_r), an essential yet challenging-to-observe variable. The overall improved performance of HydroModel($S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$) in simulating ET, particularly during droughts, implies that $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$ more accurately reflects the real-world root-zone water storage capacity compared to $S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$ and $S_{\mathrm{r}}^{\mathrm{accum}}$. Overall, our results suggest that S_r is, on average, at least 50 % larger than the water-deficit-based estimate and a staggering 380 % larger compared to the rooting-depth-based estimate. The underestimations by $S_{\mathrm{r}}^{\mathrm{accum}}$ and $S_{\mathrm{r}}^{\mathrm{RD}\times \mathrm{WHC}}$ exceed the random error of $S_{\mathrm{r}}^{\mathrm{GRACE}/\mathrm{FO}}$, underscoring the need for continued refinement and validation of S_r. Underestimating S_r may lead to overestimating ecosystem sensitivity to water stress, potentially biasing predictions of future carbon cycles (Ukkola et al., 2021; Giardina et al., 2023). Given the strong coupling between the carbon and water cycles, underestimating S_r may also lead to underestimating ecosystem water consumption and overestimating human-available water resources, particularly during droughts and heatwaves, with important implications for water resource planning (Zhao et al., 2022; Mastrotheodoros et al., 2020).

Appendix A

Figure A1 The average P−R during the largest (a), the second-largest (b), the third-largest (c), the fourth-largest (d), and the fifth-largest (e) TWS drawdowns.

Figure A2 (a) Trends in TWS obtained from GRACE/FO observations from 2002 to 2022. (b) Percentage of area equipped for irrigation that is actually irrigated. (c) Map quality marks assigned to each country for area equipped for irrigation in (b). Panels (b)–(c) are from the Global Map of Irrigation Areas version 5.0 by AQUASTAT.