Influence of irrigation on root zone storage capacity estimation

Vegetation plays a crucial role in regulating the water cycle through transpiration, which is the water flux from the subsurface to the atmosphere via roots. The amount and timing of transpiration is controlled by the interplay of seasonal energy and water supply. The latter strongly depends on the size of the root zone storage capacity (S_r), which represents the maximum accessible volume of water that vegetation can use for transpiration. S_r is primarily influenced by hydroclimatic conditions, as vegetation optimizes its root system in such a way that it guarantees water uptake and overcomes dry periods. S_r estimates are commonly derived from root zone water deficits that result from the phase shift between the seasonal signals of root zone water inflow (i.e., precipitation) and outflow (i.e., evaporation). In irrigated croplands, irrigation water serves as an additional input into the root zone. However, this aspect has been ignored in many studies, and the extent to which irrigation influences S_r estimates has never been comprehensively quantified. In this study, our objective is to quantify the influence of irrigation on S_r and identify the regional differences therein. To this end, we integrated two irrigation methods, based on the respective irrigation water use and irrigated area fractions, into the S_r estimation. We evaluated the effects compared with S_r estimates that do not consider irrigation for a sample of 4856 catchments globally with varying degrees of irrigation activity. Our results show that S_r consistently decreased when considering irrigation, with a larger effect in catchments with a larger irrigated area. For catchments with an irrigated area fraction exceeding 10 %, the median decrease in S_r was 19 and 23 mm for the two methods, corresponding to decreases of 12 % and 15 %, respectively. S_r decreased the most for catchments in tropical climates. However, the relative decrease was the largest in catchments in temperate climates. Our results demonstrate, for the first time, that irrigation has a considerable influence on S_r estimates over irrigated croplands. This effect is as strong as the effects of snowmelt that have previously been documented in catchments that have a considerable amount of precipitation falling as snow.

1 Introduction

Vegetation strongly influences the water cycle, as it controls the partitioning of precipitation into discharge and evaporation by mediating soil evaporation, interception evaporation, and transpiration (Milly, 1994). Transpiration is defined as the water transport from the subsurface back to the atmosphere via the roots of vegetation, and it is, on average, the largest terrestrial water flux globally (Schlesinger and Jasechko, 2014). The amount and timing of vegetation transpiration at catchment scales is largely controlled by the interplay between seasonal energy and water availability signals (Gentine et al., 2012). At the individual plant scale, plants also regulate transpiration by root biomass adjustments, anatomical alterations, and physiological acclimation (e.g., Brunner et al., 2015), depending on the vegetation species (Zhang et al., 2020). However, at the ecosystem scale, which represents the collective of individual plants, the subsurface water removal by transpiration is regulated by the liquid water input and by the available subsurface water buffer. This water buffer, the root zone storage capacity (S_r), is defined as the maximum volume per unit square of subsurface moisture that is accessible to the roots of vegetation for uptake (Gao et al., 2014). S_r is an essential property of hydrological systems – and parameter in land surface models and hydrological models – regulating terrestrial water, carbon, and energy balances at all scales, from the plot scale to the global scale (Seneviratne et al., 2010; Wang and Dickinson, 2012; Dralle et al., 2020a; Singh et al., 2022). Increasing evidence suggests that the extent of root systems, and consequently the magnitude of S_r, is primarily controlled by climate conditions (Kleidon and Heimann, 1998; Gao et al., 2014; De Boer-Euser et al., 2016; Kuppel et al., 2017). More specifically, the results of many studies suggest that the extent of root systems is a manifestation of vegetation (i.e., the collective of all individual plants within a specified spatial domain) having efficiently adapted to past hydroclimatic conditions. In other words, individual plants within an ecosystem have survived in competition with other plants because they found a more efficient (or optimal) balance between aboveground and belowground resource allocation (Kleidon and Heimann, 1998; Collins and Bras, 2007; Guswa, 2008; Sivandran and Bras, 2013; Fan et al., 2017; Singh et al., 2020). Direct observations of S_r at scales larger than the plot scale do not exist; therefore, several indirect methods have been developed to estimate S_r from other observable ecosystem properties considering optimality principles (Kleidon, 2004; Gao et al., 2014; Speich et al., 2018; Dralle et al., 2020a).

One of these methods is the memory method (a term coined by Van Oorschot et al., 2021), which is also referred to as the water balance method (Nijzink et al., 2016; Hrachowitz et al., 2021) or mass curve technique (Gao et al., 2014; Zhao et al., 2016). This method allows one to estimate S_r based on root zone water deficits arising from the phase shift between the seasonal signals of precipitation and evaporation – here defined as the total of transpiration, soil evaporation, and interception evaporation, following the terminology proposed by Savenije (2004) and Miralles et al. (2020). This approach is based on evidence that root systems of present-day vegetation are a legacy that reflects the memory of past water deficits during dry spells. Vegetation has efficiently adapted the extent of its root system to past water deficits with a specific memory (i.e., the dry-spell return period) to guarantee continuous access to water to satisfy canopy water demand, although no more than that (Savenije and Hrachowitz, 2017). Numerous studies have successfully demonstrated the potential of the memory method to provide estimates of climate-controlled S_r for river catchments based on discharge data (Gao et al., 2014; De Boer-Euser et al., 2016; Van Oorschot et al., 2021) as well as on larger scales based on remotely sensed estimates of evaporation (Wang-Erlandsson et al., 2016; Singh et al., 2020; McCormick et al., 2021; Stocker et al., 2023). In addition, the method proved valuable to track the temporal evolution of S_r due to changing hydroclimatic conditions (Bouaziz et al., 2022) and human interventions, such as forest management (Nijzink et al., 2016; Hrachowitz et al., 2021).

It is important to note that the memory method is based on liquid water input to the root zone. As such, solid-phase precipitation and storage as transient, seasonal, or perennial snowpacks introduces time lags between the moment of precipitation and the release of liquid water (i.e., meltwater) into the subsurface. These time lags can lead to considerable temporal shifts in liquid water supply, thereby affecting the development of seasonal water deficits and the associated magnitudes of S_r. Various models with different levels of complexity have previously been integrated into the memory method to account for the time lags due to snow accumulation and melt dynamics (de Boer-Euser et al., 2019; Dralle et al., 2021; Stocker et al., 2023). Dralle et al. (2021) recently showed that explicitly accounting for snow accumulation and associated time lags in meltwater release in the memory method does generally lead to lower values of S_r in regions where a significant fraction of precipitation occurs in the form of snow.

Irrigation similarly affects the timing of water input to the soil. Besides its effect on timing, irrigation during the growing season leads to the input of additional water besides precipitation that otherwise would not be accessible for roots and, thus, unavailable for vegetation uptake. Irrigation thereby also affects the magnitude of water input and actively shapes the root development of crops. Irrigation leads to shallower roots and higher root densities in the upper soil compared with nonirrigated vegetation, as it reduces the need for resource allocation for root growth, instead allowing increased resource allocation for aboveground growth (Klepper, 1991; Engels et al., 1994; Bakker et al., 2009; Maan et al., 2023). The strength of this signal is variable and depends on the irrigation method applied (Lv et al., 2010; Jha et al., 2017; Wang et al., 2020). Currently, approximately 20 % of global croplands are irrigated (FAO, 2022); with the increasing demand for crop production, irrigation requirements are expected to increase in the future (Alexandratos and Bruinsma, 2012). In spite of some exceptions (e.g., Roodari et al., 2021), irrigation is rarely systematically represented in hydrological and biogeophysical models (McDermid et al., 2023), mostly due to a lack of sufficient data (e.g., Meier et al., 2018). This also holds for the memory method, as most studies using the memory method for S_r estimation have not accounted for irrigation, which has likely led to an overestimation of S_r in irrigated areas (Gao et al., 2014; De Boer-Euser et al., 2016; Stocker et al., 2023). To our knowledge, only Wang-Erlandsson et al. (2016) explicitly accounted for irrigation when estimating S_r, by adding irrigation estimates simulated by the Lund–Potsdam–Jena managed Land (LPJmL) dynamic global vegetation model to the precipitation input (Jägermeyr et al., 2015). However, the extent to which irrigation influences the magnitudes of S_r estimates as well as the regions in which it is most relevant to take irrigation into account at the global scale remain unknown.

Our objective here is to quantify the influence of irrigation on the root zone storage capacity estimated with the memory method and to identify the regional differences therein. We do so by using a sample of 4856 catchments globally with varying degrees of irrigation activity. Specifically, we test the hypothesis that irrigation considerably reduces the root zone storage capacity (S_r) and, therefore, needs to be accounted for in the estimation of S_r. To this end, we introduce two methods that represent irrigation based on catchment water balances to the memory method using irrigation data from two different sources. The first method explicitly uses estimates of irrigation water use from Zhang et al. (2022) in the S_r calculation with the memory method. The second method is a simpler parameterization based on the irrigated area fraction (Siebert et al., 2015b).

2 Methods

2.1 Data

For this study, we used station-based discharge (Q) data from the following sources: the Global Streamflow Indices and Metadata Archive (GSIM; Do et al., 2018a; Gudmundsson et al., 2018a), the Australian edition of the Catchment Attributes and Meteorology for Large-sample Studies (CAMELS-AUS) dataset (Fowler et al., 2021), the LArge-SaMple DAta for Hydrology and Environmental Sciences for Central Europe (LamaH-CE; Klingler et al., 2021), and the Italian Hydrological Portal (Lendvai, 2020). We used annual mean discharge ($\stackrel{\mathrm{‾}}{Q}$) for the catchment-specific available time period. For the 1981–2010 period, we obtained the catchment average daily precipitation (P) and daily mean temperature (T_a) from the Global Soil Wetness Project Phase 3 (GSWP3; Dirmeyer et al., 2006) and daily potential evaporation (E_p) from the Global Land Evaporation Amsterdam Model version 3.5a (GLEAMv3.5a), which is based on the Priestley–Taylor approach (Martens et al., 2017; Miralles et al., 2011). We selected 4856 catchments based on the following criteria: (1) at least 10 years of Q data during the 1981–2010 period; (2) catchment area < 10 000 km² to limit the heterogeneity within catchments; (3) annual mean discharge ($\stackrel{\mathrm{‾}}{Q}$) smaller than annual mean precipitation ($\stackrel{\mathrm{‾}}{P}$) for the specific catchment.

For each catchment, we obtained irrigation estimates from two different data sources. Firstly, we used the average irrigated area fraction I_a (–), which is the areal fraction of land equipped with infrastructure for irrigation. I_a was obtained from the “AEI_HYDE_FINAL_IR” dataset developed by Siebert et al. (2015b), which is representative of the irrigation extent in the year 2005 (Fig. 1a). This dataset was based on subnational irrigation statistics and the History Database of the Global Environment (HYDE), version 3.1, land use data (Klein Goldewijk et al., 2011; Siebert et al., 2015b). Secondly, we used estimates of annual mean irrigation water use representative of the 2011–2018 period ($\stackrel{\mathrm{‾}}{I_{\mathrm{w}}}$ (mm yr⁻¹)) from Zhang et al. (2022), who developed an algorithm to estimate irrigation from multiple satellite-based products and the Priestley–Taylor Jet Propulsion Laboratory (PT-JPL) model (Fig. 1b).

Figure 1 Global irrigation characteristics. (a) Irrigated area fraction (I_a, –) representative of 2005 based on subnational irrigation statistics and the HYDE 3.1 land use data (Siebert et al., 2015b). (b) Annual mean irrigation water use ($\stackrel{\mathrm{‾}}{I_{\mathrm{w}}}$, mm yr⁻¹) for the 2011–2018 period based on multiple satellite-based products and the PT-JPL model (Zhang et al., 2022). White areas indicate I_a or $\stackrel{\mathrm{‾}}{I_{\mathrm{w}}}$ equal to zero.

To identify the effects of irrigation for different regions, we used the Köppen–Geiger climate classes as a climate indicator. For each catchment, we selected the predominant Köppen–Geiger climate class based on a global map at a 1 km resolution representing the 1980–2016 period (Beck et al., 2018a). The gridded data products for P, E_p, I_a, and $\stackrel{\mathrm{‾}}{I_{\mathrm{w}}}$ were converted to catchment estimates using area-weighted averages of the grid cells with more than 50 % of their area located inside the catchment. Before area weighting, the gridded products were resampled to a spatial resolution of 0.05° using nearest-neighbor interpolation. This way, all gridded products were treated similarly and problems with small catchments with no matching grid cells were avoided.

Figure 2 (a) Schematic bucket model representation of the memory method including an irrigation model showing storages (mm) – interception storage (S_i), snow storage (S_sn), storage deficit (S_d), surplus storage (S_s), and root zone storage capacity (S_r) – and fluxes (mm d⁻¹) – total precipitation (P), liquid precipitation (P_l), precipitation falling as snow (P_sn), interception evaporation (E_i), snowmelt (P_m), effective precipitation (P_e), transpiration (E_t), precipitation surplus (P_s), and irrigation (I). (b) An example time series of S_s, S_d, and I based on Eqs. (1–6), where Δt_d is the length of the deficit period (days) and S_s(t_s1) is the surplus storage at the end of the surplus period. Note that this time series represents only 2 years to illustrate the method, while all catchments have at least 10 years of data.

2.2 Memory method with irrigation methods

Figure 2a shows a conceptualization of the memory method based on four storage components (mm): interception storage (S_i), snow storage (S_sn), “surplus” storage (S_s), and storage deficit (S_d). S_d is initially conceptualized as an infinite deficit storage volume and its temporal evolution can be described by

$$\begin{array}{}\text{(1)}& S_{\mathrm{d}}\left(t\right)=\underset{t_{\mathrm{0}}}{\overset{\mathit{\tau }}{\int }}\left(P_{\mathrm{e}}-E_{\mathrm{t}}+I-P_{\mathrm{s}}\right)\mathrm{d}t,\end{array}$$

where P_e represents effective precipitation (mm d⁻¹), E_t is transpiration (mm d⁻¹), I is irrigation (mm d⁻¹), and P_s is surplus precipitation (mm d⁻¹) (Fig. 2a). In Eq. (1), t₀ corresponds to the first day of the first hydrological year and τ corresponds to the daily time steps ending on the last day of the last hydrological year. Our hydrological year starts the first day of the month after the wettest month, which is defined as the month with the largest positive difference between monthly mean P and E_p on average. At t₀, the starting point of the analysis, S_d=0. In Eq. (1), P_e (mm d⁻¹) is calculated from the water balance of the interception storage S_i (Fig. 2a), and E_t (mm d⁻¹) is described as a fraction of daily potential evaporation E_p (mm d⁻¹) based on the catchment water balance. We used a simple snow model based on the degree-day method (e.g., Bergstrom, 1975; Gao et al., 2017) to account for the delay in liquid water input to the soil by describing liquid precipitation (P_l, mm d⁻¹), precipitation falling as snow (P_sn, mm d⁻¹), and snowmelt (P_m, mm d⁻¹). The equations for the interception storage, snow storage, and transpiration calculation are described in Appendix A.

Surplus precipitation P_s (mm d⁻¹) in Eq. (1) is described by Eq. (2), in which we used the following notation for the sum of the fluxes between two time steps: $F_t=\underset{t-\mathrm{1}}{\overset{t}{\int }}F\mathrm{d}t$, where F is either P_e, E_t, I, or P_s. Thus P_s,t is described by

$$\begin{array}{}\text{(2)}& P_{\mathrm{s},t}=\max \left(\mathrm{0},S_{\mathrm{d}}+P_{\mathrm{e},t}-E_{\mathrm{t},t}+I_t\right),\end{array}$$

with S_d and I_t approaching zero during periods of abundant precipitation. Hence, it then holds that $P_{\mathrm{s},t}\approx P_{\mathrm{e},t}-E_{\mathrm{t},t}$.

For the computation of applied irrigation (I), we split the time series into surplus and deficit periods (Fig. 2b). For each hydrological year, we defined one deficit period, which is the longest deficit period with the largest S_d in the hydrological year. Surplus periods were defined as the periods in between the deficit periods. For each surplus period, the surplus precipitation P_s (Eq. 2) accumulates in the surplus store S_s:

$$\begin{array}{}\text{(3)}& S_{\mathrm{s}}\left(t\right)=\max \left(\mathrm{0},\underset{t_{\mathrm{s}\mathrm{0}}}{\overset{t_{\mathrm{s}\mathrm{1}}}{\int }}\left(P_{\mathrm{s}}-I\right)\mathrm{d}t\right),\end{array}$$

where t_s0 is the first day of the surplus period and t_s1 is the last day of the surplus period (Fig. 2b). S_s does not have a maximum storage capacity, but it is reset to zero each year, after each deficit period. This storage conceptualizes any water buffers that can be used for irrigation in the consecutive deficit period and may encompass ditches, lakes, and aquifers. This method assumes that irrigated water only originates from water inside the catchment boundaries and that it is sustainably extracted so that the long-term water balance is closed. During a deficit period, the fraction of S_s that is used for irrigation is defined by irrigation factor f (–), which determines how much of the surplus water stored during the surplus period is used for irrigation during the consecutive deficit period. f represents both the water evaporated or discharged during the irrigation process before recharging the soil and the spatial extent of the irrigation. It is assumed that daily irrigation I is equally distributed over the deficit period (Fig. 2b), so that I (mm d⁻¹) is defined as follows:

$$\begin{array}{}\text{(4)}& I\left(t\right)={\displaystyle \frac{f{S}_{\mathrm{s}}\left(t_{\mathrm{s}\mathrm{1}}\right)}{\mathrm{\Delta }{t}_{\mathrm{d}}}},\end{array}$$

where Δt_d is the length of the deficit period (t_d1−t_d0) in days (Fig. 2b). Here, based on the two irrigation data sources used (Sect. 2.1), we have developed two methods to estimate f in Eq. (4):

• The first of these techniques is the “irrigation water use” method (IWU), in which f_d,IWU (–) is defined for each deficit period d for each catchment by

  $$\begin{array}{}\text{(5)}& \begin{array}{rl}& f_{\mathrm{d},\mathrm{IWU}}=\max \left(\mathrm{1},{\displaystyle \frac{\stackrel{\mathrm{‾}}{I_{\mathrm{w}}}\mathrm{d}t}{S_{\mathrm{s}}\left(t_{\mathrm{s},\mathrm{1}}\right)}}\right),\mathrm{so}\mathrm{that}\\ & I\left(t\right)={\displaystyle \frac{\stackrel{\mathrm{‾}}{I_{\mathrm{w}}}\mathrm{d}t}{\mathrm{\Delta }{t}_{\mathrm{d}}}}\mathrm{if}\mathrm{sufficient}\mathrm{water}\mathrm{is}\mathrm{available}\mathrm{in}{S}_{\mathrm{s}}.\end{array}\end{array}$$

  Here, $\stackrel{\mathrm{‾}}{I_{\mathrm{w}}}$ (mm yr⁻¹) is the catchment annual mean irrigation water use, dt=1 yr, and S_s(t_s,1) (mm) is the surplus storage at the end of the surplus storage accumulation period, i.e., the amount of water stored in S_s at the start of the deficit period. In this method, f_d,IWU is different for each deficit period d, as S_s also varies. For each catchment, f_IWU is defined as the average f_d,IWU.

• The second of these techniques it the “irrigated area fraction” method (IAF), in which f_IAF (–) is temporally non-varying and is defined for each catchment by

  $$\begin{array}{}\text{(6)}& f_{\mathrm{IAF}}=\mathit{\beta }{I}_{\mathrm{a}},\end{array}$$

  where I_a (–) is the catchment irrigated area fraction and β (–) is a correction factor that is constant in space and time for all catchments. β was chosen as a constant to create a relatively simple approach that does not directly rely on irrigation water use data, which is beneficial for application in time periods (both historical and future) without irrigation water use data as well as in regions where no reliable irrigation water use data are available. We estimated β by minimizing the difference between f_IAF and f_IWU in terms of the root-mean-square error (RMSE). We generated 1000 linearly spaced values for β between 0 and 2.5 and computed f_IAF for all the catchments. For all of these cases, the RMSE of catchment f_IAF and f_IWU was computed (Fig. 3). The RMSE was lowest for β=0.9 (RMSE=0.042), which is applied for all catchments in Eq. (6).

Figure 3 The root-mean-square error (RMSE) between the catchment irrigation factors f_IWU (Eq. 5) and f_IAF (Eq. 6) for 4856 catchments for 1000 linearly spaced values of β between 0 and 2.5. β_opt represents the value for β where the RMSE minimizes. (b) Scatterplot of f_IWU (Eq. 5) and f_IAF (Eq. 6) for β=β_opt=0.9, with lighter colors indicating a higher point density.

To evaluate the effect of these methods on estimated S_r, we tested a third case, referred to as “no irrigation” (NI), in which f_NI=0. A priori, we cannot and do not consider either of the two methods, i.e., the IWU or the IAF method, to be more representative than the other. While the IWU method uses the irrigation data more directly than the IAF method, the latter directly takes the interannual variability in surplus water into account.

Figure 4 Catchment S_r for the no irrigation (NI) case, with dots representing catchment outlets. Similar panels for the IWU and IAF cases are presented in Fig. S3.

2.2.1 Root zone storage capacity calculation

Here, the catchment-scale root zone storage capacity S_r was derived from the catchment-scale storage deficit S_d time series for the three different irrigation cases NI, IWU, and IAF (Table 1). For each catchment, the annual maximum storage deficits (S_d,M) were defined for each hydrological year as follows:

$$\begin{array}{}\text{(7)}& S_{\mathrm{d},\mathrm{M}}=\max \left(S_{\mathrm{d}}\left(t\right)\right)-\min \left(S_{\mathrm{d}}\left(t\right)\right),\end{array}$$

with the min(S_d) occurring earlier in the hydrological year than the max(S_d). Previous studies (e.g., Gao et al., 2014; Wang-Erlandsson et al., 2016) applied a Gumbel distribution to the S_d,M values to estimate S_r for different return periods T. Wang-Erlandsson et al. (2016) found that the best evaporation simulations for croplands, and thus irrigated land, with a global hydrological model were achieved with an S_r based on a return period of 2 years, as croplands adapt to survive droughts with relatively short return periods. Here, we directly used the observed S_d,M values with occurrences closest to T=2 years instead of a fitted extreme value distribution, because fitting an extreme value distribution is ambiguous for return periods of interest (here 2 years) much smaller than the time series length (here >10 years). For all catchments, for each irrigation case separately, the S_r was estimated as the mean of the three observed S_d,M values with occurrences closest to T=2 years, as represented by the cross-markers closest to the vertical dashed line at T=2 years in Fig. 5b, d, f, and h.

Table 1 Details of the irrigation cases considered in this study.

2.2.2 Evaluation

To visualize the effects of irrigation on S_d and S_r, we selected four example catchments with different irrigation magnitudes (i.e., I_a and $\stackrel{\mathrm{‾}}{I_{\mathrm{w}}}$) located on four different continents and in different climate zones. For quantification of the effects of irrigation on S_r, we computed absolute (Δ) and relative (Δ_r) differences between the S_r estimates for the NI, IWU, and IAF cases (Table 1). Catchments were stratified based on (1) four different ranges of irrigated area I_a, namely, I_a≤0.01, 0.01$<I_{\mathrm{a}}\le \mathrm{0.05}$, $\mathrm{0.05}<I_{\mathrm{a}}\le \mathrm{0.1}$, and I_a>0.1 (Fig. S1 in the Supplement); (2) regions – South America, North America, Europe, and Asia; and (3) climate zones based on the Köppen–Geiger climate classification – subdivided into tropical (Af, Am, and Aw), arid (BWh, BWk, BSh, and BSk), temperate (Cfa, Cfb, and Cfc), Mediterranean (Csa and Csb), and continental (Dfa, Dfb, Dfc, and Dfd) climates (Beck et al., 2018a) (Fig. S2). Uncertainty in the differences in S_r were represented by the interquartile range (IQR).

3 Results

3.1 Irrigation influence on root zone storage capacity

Globally, the S_r estimates without accounting for irrigation ranged from 0 to 800 mm, with larger values in semiarid regions with high rainfall seasonality, such as northeastern Brazil (median S_r≈250 mm), or monsoon regions, such as northeastern Indian (median S_r≈450 mm), than in regions with temperate climates with year-round rainfall, such as western Europe (median S_r≈70 mm), or continental, colder, climates, such as Canada (median S_r≈40 mm) (Fig. 4).

Figure 5 (a, c, e, g) Time series of storage deficits S_d (mm) (Eq. 1) for four illustrative catchments with increasing irrigation from top to bottom for the three irrigation cases, NI, IWU, and IAF (Table 1); for each catchment, the associated annual mean irrigation water use ($\stackrel{\mathrm{‾}}{I_{\mathrm{w}}}$), irrigated area fraction (I_a), and root zone storage capacity (S_r) values are shown. (b, d, f, h) Return-level plot of annual maximum storage deficits (S_d,M) (Eq. 7) for the three irrigation cases (NI, IWU, and IAF), with the vertical dashed line corresponding to a return period (T) of 2 years (Sect. 2.2.1). The locations of the catchments are shown in Fig. 6. The catchment identity, continent, and Köppen–Geiger climate zone are as follows: (a, b) br_0002356, South America, temperate (Cfb); (c, d) ca_0000689, North America, continental (Dfb); (e, f) es_0000742, Europe, Mediterranean (Csa); (g, h) in_0000252, Asia, tropical (Aw).

The storage deficits S_d (Eq. 1) generally decreased when accounting for irrigation effects according to the IWU and IAF cases compared with the case without irrigation (NI). These overall effects of the method are illustrated by four selected example catchments in Fig. 5. More pronounced effects of irrigation on S_d are visible for the example catchments in Europe (Fig. 5e and f) and Asia (Fig. 5g and h), with larger $\stackrel{\mathrm{‾}}{I_{\mathrm{w}}}$ and I_a, than in the example catchments in South America (Fig. 5a and b) and North America (Fig. 5c and d). As S_d decreased, the annual maximum storage deficits S_d,M, as determined by Eq. (7), decreased as well. Consequently, the estimated S_r decreased for the IWU and IAF cases compared with NI, with more pronounced effects in the example catchments with larger $\stackrel{\mathrm{‾}}{I_{\mathrm{w}}}$ and I_a (Fig. 5). Globally, S_r consistently decreased for IWU and IAF (Fig. 6), although the magnitudes varied to a considerable extent. Nevertheless, relatively clear regional patterns of the effects of irrigation on S_r emerged. The most pronounced effects cluster in catchments in regions that are characterized by widespread and intense crop cultivation, and thus high irrigation water use, such as northern Spain, France, and parts of India (Fig. 1).

Figure 6 Relative difference in S_r (Δ_rS_r, –) for (a) IWU compared with NI ($(\mathrm{NI}-\mathrm{IWU})/\mathrm{NI}$) and (b) IAF compared with NI ($(\mathrm{NI}-\mathrm{IAF})/\mathrm{NI}$). Red markers indicate the selected catchments from Fig. 5. See Table 1 for details on the irrigation cases.

3.2 Regional differences in irrigation influence on root zone storage capacity

Figure 7 shows that the effects of irrigation on S_r increased with increasing irrigated area fraction I_a for both the IWU and IAF cases. We found the largest effects in catchments with I_a>0.1, such as the example catchment in Asia (Fig. 5g). For these catchments, the median ΔS_r was 19 mm (IQR of 10–31 mm) for IWU and 23 mm (IQR of 11–42 mm) for IAF (Fig. 7), corresponding to decreases of 12 % and 15 %, respectively (Table 2). These effects were considerably larger than the effects of irrigation in catchments with $\mathrm{0.05}<I_{\mathrm{a}}\le \mathrm{0.1}$, which reached a median Δ_rS_r of 6 %, corresponding to ΔS_r≈9 mm (Fig. 7, Table 2). Although the median effects of irrigation on S_r for catchments with I_a≤0.05 were relatively small, the effects can be considerable for specific individual catchments, as shown by the outliers in Fig. 7.

Figure 7 Box plots of the absolute S_r difference (ΔSr, m) between the irrigation cases (IWU and IAF) and the no irrigation case (NI) (Table 1). Catchments are stratified in four groups based on the irrigated area fraction I_a (Fig. S2), where n is the number of catchments in each group. The black line represents the median, the box represents the interquartile range (IQR), and the whiskers represent the 5th and 95th percentiles. White markers denote the points presented in Fig. 5. Median and IQR values for relative S_r differences (Δ_rSr, %) are presented in Table 2.

The strongest irrigation influence on S_r for catchments with I_a>0.05 was found in Asia, followed by South America, North America, and Europe, for both IWU and IAF (Fig. 8a). For the catchments in Asia, we found median ΔS_r values of 21 mm (IQR of 13–41 mm) and 27 mm (IQR of 12–56 mm) for IWU and IAF, respectively . However, the relative differences in S_r (Δ_rS_r=9 %–10 %) were smaller in Asia than in other regions, reaching up to 14 % in South America, because the initial S_r without accounting for irrigation was considerably larger in Asia than in other regions (Fig. 4, Table 2). Figure 8b shows that S_r decreased the most in tropical catchments, with a median ΔS_r of 19 mm for IWU and 24 mm for IAF. These findings are in line with the results presented in Fig. 8a, as most of the tropical catchments that we evaluated were located in Asia (Fig. S3). For catchments in the arid, Mediterranean, temperate, and continental climate zones, the median ΔS_r was smaller and varied between 5 and 15 mm. However, catchments in temperate climates exhibited the largest relative influence of irrigation on S_r: median Δ_rS_r=14 % for IWU and 15 % for IAF (Table 2).

Figure 8 Box plots of absolute S_r difference (ΔS_r, mm) between the irrigation cases (IWU and IAF) and the no irrigation case (NI) (Table 1). In panel (a) catchments are stratified regionally, similar to the maps in Fig. 6, whereas catchments are stratified based on climate zone (Sect. 2.2.2, Fig. S3) in panel (b); for both panels (a) and (b), only catchments with an irrigated area fraction I_a>0.05 are shown. The total number of catchments (n) in each group is given, with the numbers in parentheses representing n for $\mathrm{0.05}<I_{\mathrm{a}}\le \mathrm{0.1}$ and I_a>0.1, respectively. The black line represents the median, the box represents the interquartile range (IQR), and the whiskers represent the 5th and 95th percentiles. Median and IQR values for relative S_r differences (Δ_rSr, %) are presented in Table 2.

3.3 Comparison of the IWU and IAF methods

Figure 6 shows similar spatial patterns of Δ_rS_r for IWU and IAF, but the magnitudes differ. For most groups of catchments, IAF had a more pronounced effect on S_r than IWU (Table 2; Fig. S4). The different results for IWU and IAF can be explained by the different methodologies (Table 1). The IWU method directly used annual mean irrigation water use ($\stackrel{\mathrm{‾}}{I_{\mathrm{w}}}$) from Zhang et al. (2022) as an estimate for I (if sufficient water was available in the surplus store S_s). On the other hand, in the IAF method I was defined as a fraction of S_s based on the irrigated area fraction (I_a) and the constant β. Therefore, the estimated I in IAF directly reflected the interannual variability in surplus water. Another cause of the different results for IWU and IAF lies in the estimation of β in IAF, which was based on minimization of the differences between f_i,IWU and f_i,IAF (Sect. 2.2; Fig. 3). In spite of this optimization, differences between f_i,IWU and f_i,IAF remained, which partially explain the differences in ΔS_r between the two methods.

Table 2 Median and interquartile range (IQR) of the relative S_r difference (Δ_rSr, %) between the irrigation cases (IWU and IAF) and the no irrigation case (NI) with the catchments stratified based on irrigated area fraction (I_a; Fig. 7) in the top four rows, based on region (only catchments with I_a>0.05; Fig. 8a) in the middle four rows, and based on climate zones (only catchments with I_a>0.05; Fig. 8b) in the bottom five rows. The IQR is given as the 25th percentile–75th percentile.

4 Discussion

4.1 Synthesis of results

Our results show that the effect of irrigation on S_r is discernible in all regions, but the magnitude of the effect depends on the amount of irrigation applied (Figs. 5–8). For many parts of the world, the integration of irrigation in the S_r estimation did not have a large influence (Fig. 6). However, S_r considerably decreased for catchments with irrigated area fractions I_a>0.05, and ignoring irrigation in these regions would lead to biased estimates of S_r and, as a consequence, inadequate modeling of vegetation transpiration (Fig. 7). The reduction in S_r in catchments with irrigation was expected based on the fact that the memory method is founded on the theory that vegetation will invest less in roots if sufficient water is available (Guswa, 2008). Here, the observed changes in S_r are attributed to changes in the roots of vegetation, as they are directly related to the size of S_r. Additionally, adaptations at the plant scale associated with irrigation, such as adjustments in stomatal aperture (Chaves et al., 2016) and root hydraulic conductance (Lo Gullo et al., 1998), are also implicitly related to changes in S_r. The influence of irrigation on S_r estimates, as presented in Fig. 6, resembles the spatial pattern found in global assessments of irrigation water withdrawal (Huang et al., 2018) and the extent of irrigation activities (McDermid et al., 2023). This was expected because we used similar underlying irrigation data in the irrigation methods developed here. Given the ongoing irrigation expansion as presented by McDermid et al. (2023), it is expected that larger irrigation water volumes will lead to further reductions in S_r at catchment scales in the near-future compared with the reductions reported in this study. At the same time, irrigation efficiency is also improving (McDermid et al., 2023), but this effect on S_r is less straightforward. Improved irrigation efficiency (i.e., reduced soil evaporation) reduces the irrigation water volumes needed, which, at the catchment scale, leads to increased long-term mean discharge, and thus reduced long-term mean evaporation. This would result in reduced S_r in the memory method compared with a situation with lower irrigation efficiency. However, it has been shown that increased efficiency does not necessarily lead to reduced irrigation water use, as the saved water due to increasing irrigation efficiency is often applied elsewhere (Grafton et al., 2018; Lankford et al., 2020).

Previous studies using the memory method have not considered irrigation in S_r estimates (e.g., De Boer-Euser et al., 2016; Stocker et al., 2023), nor have they evaluated its effects (Wang-Erlandsson et al., 2016). To put our results into perspective, we looked at the effects of snow accumulation and melt on S_r estimates from Dralle et al. (2021) for the continental United States, as this process alters the S_d time series in a similar way to irrigation by temporally shifting liquid water input into the system. Dralle et al. (2021) estimated that integrating snow accumulation and melt in the memory method led to an average reduction in S_r of 6 mm (2 %) for areas with >10 % winter snow coverage and of 28 mm (17 %) for areas with >80 % winter snow coverage (Dralle et al., 2020b). These magnitudes are broadly consistent with our findings for irrigation (Table 2). Our results indicate that the effects of snow and irrigation on S_r are comparable. In our study, 27 catchments have both considerable snowfall (snow days >5 % of the total days) and irrigation (I_a>0.05). For these catchments, the snow model (Appendix A) led to an S_r reduction of 7 mm (7 %) on average for the NI case compared with a setup without the snow model. With irrigation, S_r further decreased by 6 mm (7 %) and 11 mm (12 %) for IAF in these catchments.

Both the results of IWU and IAF showed considerable effects of irrigation on S_r (Figs. 6–8), and both are suitable for use in the memory method, keeping in mind the individual uncertainties related to data and methodological assumptions. We think that the IWU method is more suitable for regional applications for periods with available I_w data (Zhang et al., 2022) than IAF, as I_w was derived from water balances, which strongly depend on the evaluated period. However, for spatial and temporal extrapolation, the direct use of the I_w data in the IWU method is more uncertain than the simpler IAF method, as the irrigated area fraction (I_a) used in IAF is expected to be temporally less variable than the water used for irrigation (I_w). Therefore, we think that the simpler parameterized IAF method is more suitable for use in the memory method for global applications and varying time periods. Moreover, IAF has the potential to be integrated dynamically in hydrological or land surface models employed for global Earth system model studies and future predictions.

4.2 Methodological limitations

By using several data sources, we obtained a large sample of 4856 catchments on different continents, characterized by a wide spectrum of climates and, in particular, regions with various levels of irrigation activity. However, the global coverage is not entirely balanced because Africa and large parts of Asia were undersampled. A further limitation may arise from the assumption in the memory method with irrigation methods proposed here that catchments are hydrologically closed systems. However, inter-catchment lateral flows, such as groundwater and irrigation water, can significantly alter catchment water balances (e.g., Bouaziz et al., 2018; Fan, 2019; Condon et al., 2020). Moreover, the extraction of fossil groundwater for irrigation (Siebert et al., 2010; Grogan et al., 2017; de Graaf et al., 2019) can violate the assumption of closing water balances for the irrigation methods in the memory method developed here. Our methodology, based on a sustainable water use assumption, provides a lower boundary of S_r reduction in irrigated catchments. It is expected that irrigation exceeding sustainable use would lead to larger S_r reductions than reported here, as more water is available to crops than derived from the water balance in this case. Furthermore, the methodology assumes single succession of excess and deficit periods within a year (Fig. 2b), which is not necessarily representative in regions with double-cropping systems or bimodal monsoons (Biradar and Xiao, 2011). Another limitation was the availability and quality of irrigation data (Sect. 2.1, Fig. 1). The annual mean $\stackrel{\mathrm{‾}}{I_{\mathrm{w}}}$ used in IWU was based on the 2011–2018 period, while the catchment time series varied between 1981 and 2010. Similarly, the I_a that we used represented the 2005 irrigated area fraction (Siebert et al., 2015b). The temporal mismatch between catchment hydrological time series and irrigation data may have led to an overestimation of I for the catchment-specific period, as irrigated area and irrigation techniques and efficiency have developed over the evaluated period (McDermid et al., 2023). Although this inconsistency in the temporal data influences the catchment-specific outcomes, we believe that it did not have a major influence on the quantification of the general patterns of the effects of irrigation on S_r, which was the aim of this study.

An additional source of uncertainty in the application of the memory method, as used in this study, relates to the derivation of S_r from the S_d time series (Fig. 5). Given that an ecosystem has developed its S_r such that it functions optimally and can overcome dry periods (e.g., Guswa, 2008), S_r for a specific time period would correspond to the maximum S_d value observed during that same time period (S_r=max(S_d)). However, it is important to note that the memory method represents a simplified approximation of real ecosystem behavior and has inherent limitations. The most important limitation is that our application of the memory method did not account for the feedback between S_d and $E_{\mathrm{t}}/E_{\mathrm{p}}$, which likely led to an overestimation of S_d (Van Oorschot et al., 2021). In this study, we primarily focused on crops that do not exhibit a multiyear root adaptation for survival, as there is no remaining S_r after each year's harvest. However, the catchments used here are not entirely covered by crops in any of the cases; therefore, we used a return period of 2 years for the S_r estimation (following Wang-Erlandsson et al., 2016).

5 Conclusions

Using a large sample of catchments globally, the presented results support the hypothesis that irrigation considerably reduces the root zone storage capacity (S_r) estimated with the memory method. We found a median S_r reduction of 12 % (IQR of 7 %–21 %) for the IWU method and of 15 % (IQR of 6 %–33 %) for the IAF method for catchments with an irrigated area fraction I_a>10 %. In general, these effects were less pronounced in catchments with a smaller irrigated area, although the S_r for individual catchments could also be considerably influenced by irrigation. S_r decreased the most for catchments in tropical climates, with a median decrease of 19–24 mm (for I_a>5 %). The reductions in S_r found in this study are of the same order of magnitude as the snow effects on S_r estimated by Dralle et al. (2021). Of paramount relevance for regional-to-global hydrological and climate modeling studies, this study demonstrates the relevance of irrigation for adequately estimating S_r. The irrigation water use can be expected to further increase over the next decades; therefore, the related effects on S_r should be represented in the Earth system models that are used for the next climate projections. The methodological approach developed in this study could be profitably used in this respect.

Appendix A: Memory method equations

These equations follow Van Oorschot et al. (2021), who based their methods on Gao et al. (2014), De Boer-Euser et al. (2016), Nijzink et al. (2016), and Wang-Erlandsson et al. (2016). Based on temperature, total precipitation (P, mm d⁻¹) was split into liquid precipitation (P_l, mm d⁻¹) and precipitation falling as snow (P_sn, mm d⁻¹) (Fig. 2). As temperature varies with altitude, we divided each catchment into elevation zones of 250 m. For each elevation zone (z), the daily temperature (T_z) was calculated as follows:

$$\begin{array}{}\text{(A1)}& T_{\mathrm{z}}\left(t\right)=T_{\mathrm{a}}\left(t\right)+\mathit{\lambda }\mathrm{\Delta }H,\end{array}$$

where T_a (°C) is the catchment average temperature (GSWP3), λ is the lapse rate of 0.0064 °C m⁻¹, and ΔH (m) is the elevation difference between the elevation zone and the mean elevation. For each elevation zone (z), daily P_l and P_sn were defined as follows:

$$\begin{array}{}\text{(A2)}& {\displaystyle }& {\displaystyle }{P}_{\mathrm{l},\mathrm{z}}\left(t\right)=\left\{\begin{array}{ll}P\left(t\right)& \mathrm{if}{T}_{\mathrm{z}}\left(t\right)>T_{\mathrm{t}}\\ \mathrm{0}& \mathrm{if}{T}_{\mathrm{z}}\left(t\right)<T_{\mathrm{t}},\end{array}\right.\text{(A3)}& {\displaystyle }& {\displaystyle }{P}_{\mathrm{sn},\mathrm{z}}\left(t\right)=\left\{\begin{array}{ll}P\left(t\right)& \mathrm{if}{T}_{\mathrm{z}}\left(t\right)<T_{\mathrm{t}}\\ \mathrm{0}& \mathrm{if}{T}_{\mathrm{z}}\left(t\right)>T_{\mathrm{t}}.\end{array}\right.\end{array}$$

The water balance of the snow storage (S_sn; Fig. 2) for each elevation zone (z) was described by

$$\begin{array}{}\text{(A4)}& {\displaystyle \frac{\mathrm{d}{S}_{\mathrm{sn},\mathrm{z}}}{\mathrm{d}t}}=P_{\mathrm{sn},\mathrm{z}}-P_{\mathrm{m},\mathrm{z}}.\end{array}$$

Equation (A4) can be solved by Eqs. (A5) and (A6), in which we used $F_t=\underset{t-\mathrm{1}}{\overset{t}{\int }}F\mathrm{d}t$, where F is either P_sn or P_m, for the sum of fluxes between two time steps. Thus, the numerical solution using daily time steps can be described as follows:

$$\begin{array}{}\text{(A5)}& {\displaystyle }& {\displaystyle }{S}_{\mathrm{sn},\mathrm{z},t}=S_{\mathrm{sn},\mathrm{z},t-\mathrm{1}}+P_{\mathrm{sn},\mathrm{z},t}-P_{\mathrm{m},\mathrm{z},t},\text{(A6)}& {\displaystyle }& {\displaystyle }{P}_{\mathrm{m},\mathrm{z},t}=\left\{\begin{array}{ll}\max \left(M\left(T_{\mathrm{z},t}-T_{\mathrm{t}}\right),S_{\mathrm{sn},\mathrm{z},t}\right)& \mathrm{if}{T}_{\mathrm{z},\mathrm{t}}<T_{\mathrm{t}}\\ \mathrm{0}& \mathrm{if}{T}_{\mathrm{z},\mathrm{t}}>T_{\mathrm{t}},\end{array}\right.\end{array}$$

where T_t is the threshold temperature for snowfall of 0 °C and M is the snowmelt factor of 2 mm d⁻¹ °C⁻¹. Total catchment P_l, P_sn, and P_m were calculated as an area-weighted sum of the values for the different elevation zones.

The calculation of effective precipitation P_e (mm d⁻¹) and transpiration (E_t) in Eq. (1) is similar to that in Van Oorschot et al. (2021). The water balance of the interception store (S_i; Fig. 2a) is described by

$$\begin{array}{}\text{(A7)}& {\displaystyle \frac{\mathrm{d}{S}_{\mathrm{i}}}{\mathrm{d}t}}=P_{\mathrm{l}}-P_{\mathrm{e}}-E_{\mathrm{i}},\end{array}$$

where P_l is the liquid precipitation (mm d⁻¹) and E_i is the interception evaporation (mm d⁻¹). Equation (A7) can be solved by Eqs. (A8)–(A10), in which we used $F_t=\underset{t-\mathrm{1}}{\overset{t}{\int }}F\mathrm{d}t$, where F is either P_l, E_i, P_e, or E_p (potential evaporation, mm d⁻¹), for the sum of fluxes between two time steps. Thus, the numerical solution using daily time steps can be described as follows:

$$\begin{array}{}\text{(A8)}& {\displaystyle }& {\displaystyle }{P}_{\mathrm{e},t}=\left\{\begin{array}{ll}\mathrm{0}& \mathrm{if}{P}_{\mathrm{l},t}+S_{\mathrm{i},t-\mathrm{1}}\le S_{\mathrm{i},\max }\\ P_{\mathrm{l},t}+S_{\mathrm{i},t-\mathrm{1}}-S_{\mathrm{i},\max }& \mathrm{if}{P}_{\mathrm{l},t}+S_{\mathrm{i},t-\mathrm{1}}>S_{\mathrm{i},\max },\end{array}\right.\text{(A9)}& {\displaystyle }& {\displaystyle }{S}_{\mathrm{i},t}^{*}=S_{\mathrm{i},t-\mathrm{1}}+P_{\mathrm{l},t}-P_{\mathrm{e},t},\text{(A10)}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{i},t}=\left\{\begin{array}{ll}{E}_{\mathrm{p},t}& \mathrm{if}{E}_{\mathrm{p},t}<S_{\mathrm{i},t}^{*}\\ S_{\mathrm{i},t}^{*}& \mathrm{if}{E}_{\mathrm{p},t}\ge S_{\mathrm{i},t}^{*},\end{array}\right.\end{array}$$

where E_p is potential evaporation (mm d⁻¹) and S_i,max is the maximum interception storage (mm). The size of S_i,max has a minor influence on estimates of S_r as shown by, e.g., Hrachowitz et al. (2021) and Bouaziz et al. (2020), and was therefore set to a constant value of 2.5 mm. Daily transpiration (E_t) in Eq. (1) was calculated as a fraction of daily E_p by

$$\begin{array}{}\text{(A11)}& E_{\mathrm{t}}=\left(E_{\mathrm{p}}-E_{\mathrm{i}}\right){\displaystyle \frac{\stackrel{\mathrm{‾}}{E_{\mathrm{t}}}}{\stackrel{\mathrm{‾}}{E_{\mathrm{p}}}-\stackrel{\mathrm{‾}}{E_{\mathrm{i}}}}},\end{array}$$

where $\stackrel{\mathrm{‾}}{E_{\mathrm{t}}}$ is the long-term mean E_t derived from the water balance ($\stackrel{\mathrm{‾}}{E_{\mathrm{t}}}=\stackrel{\mathrm{‾}}{P_{\mathrm{e}}}-\stackrel{\mathrm{‾}}{Q}$) and $\stackrel{\mathrm{‾}}{E_{\mathrm{p}}}$ is the long-term mean E_p.