Climate-controlled root zone parameters show potential to improve water flux simulations by land surface models

Abstract

The root zone storage capacity ($S_{\mathrm{r}}$) is the maximum volume of water in the subsurface that can potentially be accessed by vegetation for transpiration. It influences the seasonality of transpiration as well as fast and slow runoff processes. Many studies have shown that $S_{\mathrm{r}}$ is heterogeneous as controlled by local climate conditions, which affect vegetation strategies in sizing their root system able to support plant growth and to prevent water shortages. Root zone parameterization in most land surface models does not account for this climate control on root development and is based on lookup tables that prescribe the same root zone parameters worldwide for each vegetation class. These lookup tables are obtained from measurements of rooting structure that are scarce and hardly representative of the ecosystem scale. The objective of this research is to quantify and evaluate the effects of a climate-controlled representation of $S_{\mathrm{r}}$ on the water fluxes modeled by the Hydrology Tiled ECMWF Scheme for Surface Exchanges over Land (HTESSEL) land surface model. Climate-controlled $S_{\mathrm{r}}$ is estimated here with the “memory method” (MM) in which $S_{\mathrm{r}}$ is derived from the vegetation's memory of past root zone water storage deficits. $S_{\text{r,MM}}$ is estimated for 15 river catchments over Australia across three contrasting climate regions: tropical, temperate and Mediterranean. Suitable representations of $S_{\text{r,MM}}$ are implemented in an improved version of HTESSEL (Moisture Depth – MD) by accordingly modifying the soil depths to obtain a model$S_{\text{r,MD}}$ that matches $S_{\text{r,MM}}$ in the 15 catchments. In the control version of HTESSEL (CTR), $S_{\text{r,CTR}}$ is larger than$S_{\text{r,MM}}$ in 14 out of 15 catchments. Furthermore, the variability among the individual catchments of $S_{\text{r,MM}}$ (117–722 $\mathrm{mm}$) is considerably larger than of $S_{\text{r,CTR}}$ (491–725 $\mathrm{mm}$). The climate-controlled representation of $S_{\mathrm{r}}$ in the MD version results in a significant and consistent improvement of the modeled monthly seasonal climatology (1975–2010) and interannual anomalies of river discharge compared with observations. However, the effects on biases in long-term annual mean river discharge are small and mixed. The modeled monthly seasonal climatology of the catchment discharge improved in MD compared to CTR: the correlation with observations increased significantly from 0.84 to 0.90 in tropical catchments, from 0.74 to 0.86 in temperate catchments and from 0.86 to 0.96 in Mediterranean catchments. Correspondingly, the correlations of the interannual discharge anomalies improve significantly in MD from 0.74 to 0.78 in tropical catchments, from 0.80 to 0.85 in temperate catchments and from 0.71 to 0.79 in Mediterranean catchments. The results indicate that the use of climate-controlled $S_{\text{r,MM}}$ can significantly improve the timing of modeled discharge and, by extension, also evaporation fluxes in land surface models. On the other hand, the method has not been shown to significantly reduce long-term climatological model biases over the catchments considered for this study.

1 Introduction

Vegetation controls the partitioning of precipitation into evaporation and runoff by transporting water through its roots to the atmosphere and is thereby key in the representation of land surface–atmosphere interactions . The moisture flow from the land surface to the atmosphere through vegetation root water uptake is defined as transpiration and is globally the largest water flux from terrestrial ecosystems . The contribution of transpiration to total land evaporation is regulated by the interplay between the atmospheric water demand and the soil moisture within the reach of vegetation's roots. The root zone is defined as the part of the subsurface where vegetation has developed roots and can be characterized by parameters such as root depth and root density. The importance of the root zone in land surface and climate modeling is widely acknowledged, and multiple studies emphasize the climate sensitivity to changes in the vegetation's root zone . However, the parameterization of the root zone in state-of-the-art land surface models (LSMs) is a possible cause for the large uncertainties in water flux representations in these models , which is particularly true for land evaporation simulations .

The hydrologically relevant magnitude of the vegetation's root zone can be described by the root zone water storage capacity, $S_{\mathrm{r}}$, that represents the maximum subsurface moisture volume that can be accessed by the vegetation's roots. The size of $S_{\mathrm{r}}$ controls the variability and timing of water fluxes and specifically the ability of vegetation to maintain transpiration during the dry season when there is little to no recharge . It is important to note that $S_{\mathrm{r}}$ is not necessarily proportional to the depth of roots. While root depth only describes the vertical root profile, $S_{\mathrm{r}}$ also accounts for lateral root extent and root density. For example, an ecosystem covered by deep-rooting vegetation with roots with low density likely has a smaller $S_{\mathrm{r}}$ than one covered by vegetation with shallow, high-density roots .

However, most global LSMs do not have the explicit objective to estimate $S_{\mathrm{r}}$ and rather aim for a description of root zone parameters (e.g., root depth, root density and root distribution) for different vegetation classes combined with soil type information and a model-dependent fixed soil depth. The generally shallow ($<\mathrm{2}$ $\mathrm{m}$) fixed soil depth limits the size of $S_{\mathrm{r}}$ and, as a consequence, also the moisture extraction by roots from deep soil layers . LSMs use lookup tables that prescribe the same root zone parameters worldwide for each combination of vegetation and soil class as obtained from a very limited number of point-scale observations of rooting structure . The spatial distribution of the root zone parameterization in LSMs is obtained by combining these lookup table values with maps of vegetation cover and soil texture. The limitations of this approach are as follows: the root observations (1) are uncertain due to the fact that they mostly vertically extrapolate root measurements while excavating only the first meter or less ; (2) do not adequately represent global distributions of root structures because observations are extremely scarce – e.g., the dataset includes 475 root profiles in 209 geographical locations; (3) are observations of individual plants that do not represent spatial variations in ecosystem composition at larger scales than the plot scale; and (4) are snapshots in time and therefore do not represent their evolution over time due to continuous adaptation of ecosystems to changing environmental conditions.

An alternative to the lookup tables based on point-scale root observations for describing the vegetation's root zone is a climate-controlled approach. The only LSM to our knowledge in which climate-controlled root zone parameters are used is the JSBACH3.2 model in which rooting depths are based on the optimization model of net primary production from . Yet, there is generally strong evidence that climate is the dominant control of root development in many environments, as vegetation tends to optimize its aboveground and belowground carbon investment in order to optimally function by avoiding water shortages and maintaining transpiration and productivity . For example, it is likely that trees in a dry climate will develop a larger $S_{\mathrm{r}}$ than trees in a wet climate because trees in a dry climate need to invest more in growing roots to sustain their water demand .

A widely applied climate-controlled approach in catchment hydrological studies to describe $S_{\mathrm{r}}$ is the “memory method”. In this method $S_{\mathrm{r}}$ is derived from water storage deficit calculations in the root zone at the catchment scale, assuming vegetation is able to keep memory of past deficit conditions to size roots in such a way to guarantee continuous access to water (hereinafter $S_{\text{r,MM}}$) . Recent studies have demonstrated that this method provides plausible catchment-scale estimates of $S_{\mathrm{r}}$ e.g., that result in improvements in modeling catchment discharge compared to soil-derived $S_{\mathrm{r}}$ estimates . However, climate-controlled root zone parameters have not yet been widely incorporated in LSMs.

Table 1

Average hydrological characteristics of the catchments in the three climate regions for the time period 1973–2010 with long-term mean annual discharge $\stackrel{\mathrm{‾}}{Q}$, long-term mean annual precipitation $\stackrel{\mathrm{‾}}{P}$, long-term mean annual potential evaporation ${\stackrel{\mathrm{‾}}{E}}_{\mathrm{p}}$, aridity index $I_{\mathrm{A}}={\stackrel{\mathrm{‾}}{E}}_{\mathrm{p}}/\stackrel{\mathrm{‾}}{P}$, and the seasonality index of precipitation $I_{\mathrm{S}}=\frac{\mathrm{1}}{\stackrel{\mathrm{‾}}{P_{\mathrm{a}}}}{\sum }_{m=\mathrm{1}}^{m=\mathrm{12}}|\stackrel{\mathrm{‾}}{P_m}-\frac{\stackrel{\mathrm{‾}}{P_{\mathrm{a}}}}{\mathrm{12}}|$, where $\stackrel{\mathrm{‾}}{P_{\mathrm{a}}}$ is the annual mean precipitation and $\stackrel{\mathrm{‾}}{P_m}$ the monthly mean precipitation in month $m$ ; $\mathit{\varphi }$ is the time lag between long-term mean maximum monthly precipitation ($P$) and potential evaporation ($E_{\mathrm{p}}$). Values for all individual catchments are provided in Table S2.

Climate region	$\stackrel{\mathrm{‾}}{Q}$	$\stackrel{\mathrm{‾}}{P}$	$\stackrel{\mathrm{‾}}{E_{\mathrm{p}}}$	$I_{\mathrm{A}}$	$I_{\mathrm{S}}$	$\mathit{\varphi }$
	($\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$)	($\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$)	($\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$)	(–)	(–)	($\mathrm{months}$)
Tropical (seven catchments)	302	1101	1869	2	0.9	2.3
Temperate (five catchments)	57	651	1488	2.5	0.2	0.6
Mediterranean (three catchments)	53	879	1276	1.7	0.3	5.7

The objective of this study is to quantify and evaluate the effects of a climate-controlled representation of $S_{\mathrm{r}}$ on the water fluxes modeled by the Hydrology Tiled ECMWF Scheme for Surface Exchanges over Land (HTESSEL) land surface model. Specifically, we will test the hypothesis that implementing $S_{\text{r,MM}}$ in HTESSEL can improve the modeled magnitude and timing of catchment discharge and evaporation fluxes. By applying the memory method for estimating ecosystem-scale $S_{\mathrm{r}}$ for use in LSMs, the first three limitations of using sparse root observations mentioned above can be overcome, but it should be acknowledged that, although the memory method in principle allows for adaptively updating $S_{\mathrm{r}}$, in this work we use a fixed value in time. In this study, $S_{\text{r,MM}}$ values representative for the 1973–2010 time period are estimated for 15 Australian catchments across different climate regions (Sect.  and Appendix A). The $S_{\text{r,MM}}$ estimates are then used to constrain the $S_{\mathrm{r}}$ in HTESSEL (Sect. ). Section  evaluates the effects on discharge and evaporation in HTESSEL by performing offline simulations with and without the improved representation of $S_{\mathrm{r}}$. Finally, in Sects.  and the potential for a wider application of climate-controlled root zone parameters is discussed.

2 Methods

2.1 Study area

Australia is characterized by large spatial differences in precipitation (Fig. ), vegetation coverage and temperatures, varying from hot and dry deserts in the interior to tropical forests with a monsoon season in the north. We have selected 15 Australian river catchments with station observations of river discharge at the outlet of the catchment to estimate $S_{\mathrm{r}}$ by applying the memory method (Fig. ; Table S1 in the Supplement) . The catchments are selected based on available discharge data (at least 30 years of station observations), size (at least one-third of the land surface model grid cell area of approximately 5500 ${\mathrm{km}}^{\mathrm{2}}$ in order to spatially extrapolate catchment characteristics to grid cells) and differences in climate (spatial spread of the catchments across Australia for the analysis of different climate zones). The catchments are classified in three climate regions based on their hydrological characteristics (Table ; Fig. ; Table S2). The tropical catchments are characterized by pronounced seasonality of rainfall with a seasonality index of precipitation ($I_{\mathrm{S}}$) of 0.7 or higher, while temperate and Mediterranean catchments have year-round rainfall ($I_{\mathrm{S}}<\mathrm{0.7}$). The Mediterranean catchments are characterized by a time lag $\mathit{\varphi }$ between long-term mean maximum monthly potential evaporation $E_{\mathrm{p}}$ and precipitation $P$ of 5 or 6 months, while in tropical and temperate catchments mean maximum monthly $E_{\mathrm{p}}$ and $P$ occur within 3 months.

Figure 1

Location of the 15 study catchments within Australia. The green, red and orange markers indicate the climate region, and the blue shades indicate long-term mean annual precipitation . A list of the catchments and their characteristics is provided in Table S1.

[Figure omitted. See PDF]

Figure 2

Monthly seasonal climatology of precipitation ($P$) and potential evaporation ($E_{\mathrm{p}}$) for the (a) tropical, (b) temperate and (c) Mediterranean catchments, with the solid lines representing $P$ and the dashed lines $E_{\mathrm{p}}$, for the time series 1973–2010. The different shades indicate the 15 individual study catchments.

[Figure omitted. See PDF]

2.2 Data

For this study we use daily discharge data from station observations in the catchments for the time period 1973–2010 . For the same time period we use daily precipitation and daily mean temperature data from the GSWP-3 dataset on a regular 0.5${}^{\circ }$grid . Daily $E_{\mathrm{p}}$ is calculated by applying the Hargreaves and Samani formulation based on temperature and radiation . The FLUXCOM RS $+$ METEO dataset is used as a reference dataset to benchmark modeled actual evaporation. FLUXCOM provides a gridded product of interpolated monthly evaporation as a fusion of FLUXNET eddy covariance towers, satellite observations and meteorological data (GSWP-3) for the time period 1975–2010 . This dataset has shown plausible estimates of mean annual and seasonal evaporation and is generally considered a suitable tool for global land model evaluations . However, we found considerable differences between the long-term annual mean evaporation ${\stackrel{\mathrm{‾}}{E}}_{\text{FLUXCOM}}$ and $\stackrel{\mathrm{‾}}{E}$ derived from the catchment water balance (${\stackrel{\mathrm{‾}}{E}}_{\text{WB}}$) based on observed $Q$ and GSWP-3 $P$ (${\stackrel{\mathrm{‾}}{E}}_{\text{WB}}=\stackrel{\mathrm{‾}}{P}-\stackrel{\mathrm{‾}}{Q}$) (Fig. ). Figure  clearly illustrates that the ${\stackrel{\mathrm{‾}}{E}}_{\text{FLUXCOM}}$ is consistently lower than ${\stackrel{\mathrm{‾}}{E}}_{\text{WB}}$, with an average difference of 150 $\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$, which is equivalent to about 20 $\mathrm{\%}$ of the long-term water balances. ${\stackrel{\mathrm{‾}}{E}}_{\text{WB}}$ is likely to be more reliable than ${\stackrel{\mathrm{‾}}{E}}_{\text{FLUXCOM}}$ because ${\stackrel{\mathrm{‾}}{E}}_{\text{WB}}$ provides an integrated catchment-scale estimate as it is derived from observations of $Q$ assuming that the catchments are large enough to neglect deep groundwater drainage to or from other catchments . In addition, ${\stackrel{\mathrm{‾}}{E}}_{\text{FLUXCOM}}$ is based on point-scale estimates of FLUXNET stations that do not coincide with and are mostly located far from the study catchments . The discrepancy between the FLUXCOM and catchment water balance is addressed by scaling the monthly FLUXCOM evaporation:

1$${\displaystyle }{\displaystyle }{E}_{\text{FLUXCOM-WB}}=E_{\text{FLUXCOM}}{\displaystyle \frac{{\stackrel{\mathrm{‾}}{E}}_{\text{WB}}}{{\stackrel{\mathrm{‾}}{E}}_{\text{FLUXCOM}}}},$$ with $E_{\text{FLUXCOM-WB}}$ as the monthly reference evaporation representative for the catchment scale, $E_{\text{FLUXCOM}}$ from in the catchment corresponding grid cells and $\frac{{\stackrel{\mathrm{‾}}{E}}_{\text{WB}}}{{\stackrel{\mathrm{‾}}{E}}_{\text{FLUXCOM}}}$ the catchment-specific scaling factor.

Figure 3

Long-term mean annual evaporation ($\stackrel{\mathrm{‾}}{E}$) as estimated from long-term water balance data (${\stackrel{\mathrm{‾}}{E}}_{\text{WB}}$) compared to the FLUXCOM dataset (${\stackrel{\mathrm{‾}}{E}}_{\text{FLUXCOM}}$) for the 1975–2010 period.

[Figure omitted. See PDF]

We use gridded data on vegetation type and coverage derived from the GLCC1.2 and soil texture data from the FAO/UNESCO Digital Soil Map of the World . Characteristics of the different soil textures are based on the Van Genuchten soil parameters . These data are needed as input to the HTESSEL model and for the estimation of $S_{\mathrm{r}}$.

2.3 Memory method for estimating root zone storage capacity

$S_{\text{r,MM}}$ is estimated based on catchment hydrometeorological data according to the methodology described in the studies of , and . $S_{\text{r,MM}}$ is based on an extreme value analysis of the annual maximum water storage deficits in the vegetation's root zone ($S_{\mathrm{d}}$). $S_{\mathrm{d}}$ maximizes during dry periods, and therefore $S_{\mathrm{r}}$ represents an upper limit of root zone storage assuming that vegetation has sufficient access to water to overcome these dry periods. The cumulative water storage deficit $S_{\mathrm{d}}$ ($\mathrm{mm}$) in the root zone is based on daily time series of effective precipitation $P_{\mathrm{e}}$ ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$) and transpiration $E_{\mathrm{t}}$ ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$) for the time period 1973–2010 and is described by

2$${\displaystyle }{\displaystyle }{S}_{\mathrm{d}}\left(t\right)=\text{max}\left(\mathrm{0},-\underset{t_{\mathrm{0}}}{\overset{\mathit{\tau }}{\int }}(P_{\mathrm{e}}-E_{\mathrm{t}})\mathrm{d}t\right),$$ with an integration from $t_{\mathrm{0}}$ that corresponds to the first day in the hydrological year 1973 to $\mathit{\tau }$ that corresponds to the daily time steps ending on the last day of the hydrological year 2010. $P_{\mathrm{e}}$ ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$) is derived from the water balance of the interception storage $S_{\mathrm{i}}$: 3$${\displaystyle }{\displaystyle \frac{\mathrm{d}{S}_{\mathrm{i}}}{\mathrm{d}t}}=P-E_{\mathrm{i}}-P_{\mathrm{e}},$$ with $P$ representing the precipitation ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$) and $E_{\mathrm{i}}$ the interception evaporation ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$). Equation () can be solved by Eqs. ()–(). Herein, for the sum of fluxes between two time steps the following notation is used: $F_t={\int }_{t-\mathrm{1}}^{t}F\mathrm{d}t$, where $F$ is either $P$, $E_{\mathrm{i}}$, $P_{\mathrm{e}}$ or $E_{\mathrm{p}}$ (potential evaporation, $\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$). The numerical solution was then thus obtained as follows using daily time steps. $$\begin{array}{}\text{4}& {\displaystyle }& {\displaystyle }{P}_{\mathrm{e},t}=\left\{\begin{array}{ll}\mathrm{0}& \text{if }{P}_t+S_{\mathrm{i},t-\mathrm{1}}\le S_{\text{i,max}}\\ P_t+S_{\mathrm{i},t-\mathrm{1}}-S_{\text{i,max}}& \text{if }{P}_t+S_{\mathrm{i},t-\mathrm{1}}>S_{\text{i,max}}\end{array}\right.\text{5}& {\displaystyle }& {\displaystyle }{S}_{\mathrm{i},t}^{*}=S_{\mathrm{i},t-\mathrm{1}}+P_t-P_{\mathrm{e},t}\text{6}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{i},t}=\left\{\begin{array}{ll}{E}_{\mathrm{p},t}& \text{if }{E}_{\mathrm{p},t}<S_{\mathrm{i},t}^{*}\\ S_{\mathrm{i},t}^{*}& \text{if }{E}_{\mathrm{p},t}\ge S_{\mathrm{i},t}^{*}\end{array}\right.\end{array}$$ Here, $S_{\text{i,max}}$ is the maximum interception storage ($\mathrm{mm}$) that depends on the land cover and is estimated between 2–8 $\mathrm{mm}$ for a tropical forest and between 0–3 $\mathrm{mm}$ for a temperate forest . However, found that the sensitivity of $S_{\mathrm{r}}$ to the value of $S_{\text{i,max}}$ is small, and therefore a value of 2.5 $\mathrm{mm}$ is used here in all catchments for simplicity.

Daily $E_{\mathrm{t}}$ ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$) in Eq. () was calculated by 7$${\displaystyle }{\displaystyle }{E}_{\mathrm{t}}=c{E}_{\mathrm{p}},$$ where $c$ ($-$) is a coefficient that represents the ratio between transpiration and potential evaporation $c=\stackrel{\mathrm{‾}}{E_{\mathrm{t}}}/\stackrel{\mathrm{‾}}{E_{\mathrm{p}}}$. $\stackrel{\mathrm{‾}}{E_{\mathrm{t}}}$ ($\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$) is the long-term mean transpiration derived from the water balance ($\stackrel{\mathrm{‾}}{E_{\mathrm{t}}}=\stackrel{\mathrm{‾}}{P_{\mathrm{e}}}-\stackrel{\mathrm{‾}}{Q}$) and $\stackrel{\mathrm{‾}}{E_{\mathrm{p}}}$ ($\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$) the long-term mean potential evaporation. $E_{\mathrm{t}}$ considered here includes both transpiration and soil evaporation, but as the latter is much smaller, we use the term transpiration for simplicity. The subtle interactions between atmospheric water demand and vegetation-available water supply can lead to interannual variability in $c$. The above-described approach that provided constant estimates of $c$ is therefore extended by an iterative procedure to estimate annually varying values of the coefficient $c$ as described in Appendix .

Catchment $S_{\text{r,MM}}$ ($\mathrm{mm}$) is estimated based on the assumption that a catchment's ecosystem designs its rooting system while keeping memory of water stress events with certain return periods. Previous studies provide evidence that these return periods are likely to be longer for high vegetation (e.g., forest) than for low vegetation (e.g., grass). Based on the results of , and drought return periods (RP) for high and low vegetation are set to 40 and 2 years, respectively. The $S_{\text{r,MM}}$ corresponding to these drought return periods is calculated by applying the Gumbel extreme value distribution to annual maximum storage deficits. Theoretically we could treat $S_{\mathrm{r}}$ separately for high and low vegetation in HTESSEL. However, this would require changing the root distributions (see Sect. 2.4), which we decided not to do as we did not want to change multiple parameters at the same time. Therefore, for the implementation of $S_{\text{r,MM}}$ in HTESSEL, catchment $S_{\text{r,MM}}$ is estimated as a weighted sum of the high and low vegetation $S_{\mathrm{r}}$ based on the coverage fraction of high ($C_{\mathrm{H}}$) and low ($C_{\mathrm{L}}$) vegetation in the corresponding grid cell of that specific catchment, described by 8$${\displaystyle }{\displaystyle }{S}_{\text{r,MM}}=C_{\mathrm{L}}{S}_{\text{r,L,2yr}}+C_{\mathrm{H}}{S}_{\text{r,H,40yr}}.$$

2.4 HTESSEL model description

In this study we use the Hydrology Tiled ECMWF Scheme for Surface Exchanges over Land (HTESSEL) land surface model . This section presents the model parameterization of vegetated areas in the HTESSEL control model version (hereinafter CTR) based on the IFS documentation of cycle CY43R1 and the model code itself . The core structure of this model is described by , and major changes in the hydrology parameterization were made by with the implementation of a global soil texture map instead of a single soil type and a runoff scheme accounting for sub-grid variability, which resulted in improvements in global water budget simulations .

Figure 4

Root zone parameterization in the HTESSEL CTR version with the directly changed parameters in the HTESSEL MD version highlighted in red. (a) 3D overview of a single grid cell. (b) Schematic image of the four-layer subsurface. (c) Scheme of equations for the calculation of soil moisture, discharge and evaporation. The symbols in this figure are as follows, with $i$ representing high ($\mathrm{H}$) and low ($\mathrm{L}$) vegetation and $k$ layers 1–4: $C$ ($-$) vegetation coverage, $T$ dominant vegetation type, $z$ ($\mathrm{m}$) layer depth, $P$ ($\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) precipitation, $P_{\mathrm{t}}$ ($\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) precipitation through-fall, $M$ ($\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) snowmelt, $Q$ ($\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) total discharge, $Q_{\mathrm{s}}$ ($\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) surface runoff, $Q_{\text{sb}}$ ($\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) subsurface runoff, $I_{\text{max}}$ ($\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) maximum infiltration rate, $b$ ($-$) variable representing sub-grid orography, $E$ ($\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) total evaporation, $E_{\mathrm{t}}$ ($\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) transpiration, $E_{\mathrm{s}}$ ($\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) soil evaporation, $E_{\mathrm{c}}$ ($\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) canopy evaporation, $R$ ($\mathrm{\%}$) root distribution, $\mathit{\theta }$ (${\mathrm{m}}^{\mathrm{3}}{\mathrm{m}}^{-\mathrm{3}}$) unfrozen soil moisture, ${\mathit{\theta }}_{\text{pwp}}$ (${\mathrm{m}}^{\mathrm{3}}{\mathrm{m}}^{-\mathrm{3}}$) soil moisture at permanent wilting point, ${\mathit{\theta }}_{\text{cap}}$ (${\mathrm{m}}^{\mathrm{3}}{\mathrm{m}}^{-\mathrm{3}}$) soil moisture at field capacity, ${\mathit{\theta }}_{\text{sat}}$ (${\mathrm{m}}^{\mathrm{3}}{\mathrm{m}}^{-\mathrm{3}}$) soil moisture at saturation, $S_{\mathrm{r}}$ ($\mathrm{m}$) the root zone storage capacity, ${\mathit{\theta }}_{\text{roots}}$ (${\mathrm{m}}^{\mathrm{3}}{\mathrm{m}}^{-\mathrm{3}}$) the root extraction efficiency, $r_{\mathrm{c}}$ ($\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$) canopy resistance, $r_{\mathrm{a}}$ ($\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$) atmospheric resistance, $R_{\mathrm{s}}$ ($\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$) downward shortwave radiation, $D_{\mathrm{a}}$ ($\mathrm{hPa}$) atmospheric water vapor deficit, $q$ specific humidity ($\mathrm{kg}{\mathrm{kg}}^{-\mathrm{1}}$), $r_{\text{s,min}}$ ($\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$) minimum canopy resistance, LAI ($-$) leaf area index, $S_{\mathit{\theta }}$ (${\mathrm{m}}^{\mathrm{3}}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$) root extraction rate, $\mathit{\gamma }$ ($\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) hydraulic conductivity, $\mathit{\lambda }$ (${\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) hydraulic diffusivity and ${\mathit{\rho }}_{\mathrm{w}}$ ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$) density of water.

[Figure omitted. See PDF]

Figure a represents a simplified 3D view of a single grid cell. The HTESSEL model describes eight different surface fractions within a grid cell , but we only considered the vegetation-covered fractions (high and low vegetation) because of the presence of roots. Considering exclusively vegetated areas, the grid cell surface is subdivided into high and low vegetation-covered area ($C_{\mathrm{H}}$ and $C_{\mathrm{L}}$) with a dominant type of vegetation ($T_{\mathrm{H}}$ and $T_{\mathrm{L}}$) based on the GLCC1.2 vegetation database. This database distinguishes 18 different vegetation types (e.g., evergreen broadleaf, tall grass, crops), each described with vegetation-specific parameters based on experiments and literature (e.g., minimum canopy resistance, root distribution). The subsurface has a single soil texture based on and is subdivided into four model layers with a total depth $z$ of 2.89 $\mathrm{m}$ that is kept uniformly constant in the global domain.

Figure b presents the connection of the subsurface with the surface through roots and transpiration fluxes ($E_{\mathrm{t}}$) in more detail. $S_{\mathrm{r}}$ is not explicitly described in the model parameterization, and therefore it is formulated based on our own understanding of its relation to the HTESSEL vegetation and root zone parameterizations (Eq. ). Vegetation has roots in all four model soil layers (except for the vegetation types desert and tundra that can only access the upper layer and the upper three layers, respectively; ). There is a variable root distribution across the layers that is different for each vegetation type. The vegetation-specific root distribution ($R_k$) describes the root fraction with respect to the total amount of roots in each model soil layer. At a single time step, the capability of roots to extract soil moisture (${\mathit{\theta }}_{k,\text{roots}}$, represented by the brown boxes in b) is a function of $R_k$ and the layer of unfrozen soil moisture content (${\mathit{\theta }}_k$). Thus, the more roots we have in a soil layer, the more moisture can be extracted at each time step. In the long term, however, the vegetation is able to extract all the plant-available soil moisture in the layers where roots are present. Therefore, $S_{\text{r,CTR}}$, represented in blue in Fig. b, is described by 9$${\displaystyle }{\displaystyle }{S}_{\text{r,CTR}}=z({\mathit{\theta }}_{\text{cap}}-{\mathit{\theta }}_{\text{pwp}}),$$ with $z$ representing the hydrologically active depth, which corresponds to the combined depth of all soil layers with roots ($z=\mathrm{2.89}$ $\mathrm{m}$ is a default value in HTESSEL for all vegetation types except for desert and tundra), and ${\mathit{\theta }}_{\text{cap}}-{\mathit{\theta }}_{\text{pwp}}$ the plant-available moisture, which is constant over the four soil layers. The plant-available moisture is bounded by the soil-texture-specific moisture contents at field capacity (${\mathit{\theta }}_{\text{cap}}$), above which soil moisture drains by gravity, and at the wilting point (${\mathit{\theta }}_{\text{pwp}}$), below which soil moisture is not accessible to roots. It should be noted that we aimed for a physical definition of $S_{\text{r,CTR}}$ but that the effective water used by vegetation may be different. We come back to this point more elaborately in the Discussion section (Sect. ).

Figure c presents the equation scheme of HTESSEL for calculating soil moisture as well as discharge and evaporation fluxes, with $i$ representing high ($\mathrm{H}$) or low vegetation cover ($\mathrm{L}$) and $k$ the four soil layers. The relative soil moisture content $\mathit{\theta }$ controls the calculations of discharge and evaporation fluxes. The surface runoff ($Q_{\mathrm{s}}$) is defined by the precipitation through-fall ($P_{\mathrm{t}}$), snowmelt ($M$) and maximum infiltration rate ($I_{\text{max}}$) (Eq. 10). $I_{\text{max}}$ is a function of $P_{\mathrm{t}}$, $M$, a spatially variable parameter ($b$) that is defined by the standard deviation in sub-grid orography, and the vertically integrated (top 0.5 $\mathrm{m}$) soil moisture ($\mathit{\theta }$) and saturation soil moisture (${\mathit{\theta }}_{\text{sat}}$) (Eq. 11) . The subsurface runoff ($Q_{\text{sb}}$) consists of two components: free drainage from layer 4, which is a function of hydraulic conductivity in this layer (${\mathit{\gamma }}_{\mathrm{4}}$) and water density (${\mathit{\rho }}_{\mathrm{w}})$ (Eq. 12), and the excess absolute soil moisture when ${\mathit{\theta }}_k>{\mathit{\theta }}_{\text{sat}}$ (Eq. 13). Total discharge ($Q$) is the sum of $Q_{\mathrm{s}}$ and $Q_{\text{sb}}$ (Eq. 14), and as typical in-stream travel times through the catchments are about 1 $\mathrm{d}$ at most, we did not consider routing to be important at the monthly timescale for which we analyze the results. The average root extraction efficiency in all layers (${\mathit{\theta }}_{\text{roots}}$) is described by Eqs. (15) and (16) as the weighted sum of the vegetation-specific $R_k$ and ${\mathit{\theta }}_k$. The canopy resistance ($r_{\mathrm{c}}$) (Eq. 17) describes the resistance of vegetation to transpiration and is a function of vegetation-specific values for minimum canopy resistance ($r_{\text{s,min}}$) and LAI, a function of shortwave radiation ($f_{\mathrm{1}}\left(R_{\mathrm{s}}\right)$), a function of atmospheric water vapor deficit ($f_{\mathrm{2}}\left(D_{\mathrm{a}}\right)$) and a function of the root extraction efficiency ($f_{\mathrm{3}}({\mathit{\theta }}_{\text{roots},i},{\mathit{\theta }}_{\text{pwp}},{\mathit{\theta }}_{\text{cap}}$). The canopy resistance defines $E_{\text{t,i}}$ together with specific humidity ($q$) and an atmospheric resistance term ($r_{\mathrm{a}}$) (Eq. 18). Total $E_{\mathrm{t}}$ is a weighted sum of the separate transpiration products based on the sub-grid coverage $C_{\mathrm{L}}$ and $C_{\mathrm{H}}$ (Eq. 19), and total evaporation ($E$) is a sum of transpiration ($E_{\mathrm{t}}$), soil evaporation ($E_{\mathrm{s}}$) and canopy evaporation ($E_{\mathrm{c}}$) fluxes (Eq. 20). The detailed formulations of the latter two fluxes are not relevant in this study and are therefore not included in this model description. $E_{\text{t,i}}$ is attributed to the different soil layers in the calculation of the root extraction ($S_{\mathit{\theta }}$) based on the layer depth ($z_k$) and ${\mathit{\theta }}_{k,\text{roots}}$ (Eq. 21). The change in soil moisture over time ($\partial \mathit{\theta }/\partial t$) is calculated by applying the Darcy–Richards equation with $\mathit{\gamma }$ and $\mathit{\lambda }$ representing hydraulic conductivity and diffusivity (Eq. 22). This equation is solved with a top soil boundary condition of $P-E-Q_{\mathrm{s}}$ and a bottom soil boundary condition of free drainage.

2.5 Implementation of memory method root zone storage capacity estimates in HTESSEL

Here we develop an approach to implement the climate-controlled $S_{\text{r,MM}}$ (results in Sect. ) in HTESSEL, while maintaining the modeling framework of the CTR model described in Sect. . We found that $S_{\text{r,CTR}}$ is exclusively defined by the soil type and the hydrologically active model soil depth ($z$) (Eq. ). In our modified version of HTESSEL, hereafter referred to as the Moisture Depth (MD) model, the soil depth for moisture calculations is changed to satisfy the following equation:

23$${\displaystyle }{\displaystyle }{S}_{\text{r,MM}}=z_{\text{MD}}({\mathit{\theta }}_{\text{cap}}-{\mathit{\theta }}_{\text{pwp}}),$$ with $z_{\text{MD}}$ as the total soil depth in the MD model modified to satisfy $S_{\text{r,MD}}=S_{\text{r,MM}}$. This depth change is achieved by changing model layer 4, except in the case that this would cause the model depth of layer 4 to approach zero ($z_{\mathrm{4}}\approx \mathrm{0}$). In this case a minimum threshold (0.2 $\mathrm{m}$) is set for $z_{\mathrm{4}}$, and the depth of layer 3 is further changed to obtain $S_{\text{r,MD}}=S_{\text{r,MM}}$ as required in Eq. (). This is necessary because $z_{\mathrm{4}}\approx \mathrm{0}$ in the moisture calculation would cause inconsistencies in the thermal diffusion calculations as the layer soil temperature is a function of the layer soil moisture. The layer depths for thermal diffusion calculations are not modified in the MD model, and we found that the soil layer temperatures are insensitive to depth changes in MD. The directly changed parameters in MD are highlighted in red in Fig. . Also, the root distribution is not modified in MD because we aimed for a physical representation of $S_{\mathrm{r}}$ (Eq. ) and we did not want to change multiple model parameters at the same time. Furthermore, we would like to reiterate that the soil depth in the model should be interpreted neither as actual soil depth nor rooting depth, but merely as a way to represent the plant-accessible water volume.

2.6 Model simulations

Simulations are performed in a stand-alone version of HTESSEL as it was implemented in the framework of version 3 of the EC-EARTH Earth system model (http://www.ec-earth.org, last access: August 2020) for both the CTR (Sect. ) and MD (Sect. ) model versions. The model is forced with 3-hourly GSWP-3 atmospheric boundary conditions for the historical time series 1970–2010, with the first 5 years used for spin-up. The spatial resolution of the HTESSEL model is a reduced Gaussian grid (N128), with the grid cells over Australia being approximately 5500 ${\mathrm{km}}^{\mathrm{2}}$.

Table 2

Long-term annual mean modeled discharge ($\stackrel{\mathrm{‾}}{Q}$) and evaporation ($\stackrel{\mathrm{‾}}{E}$) in the HTESSEL CTR and MD versions for the tropical, temperate and Mediterranean climate regions (catchment averages) as well as reference $\stackrel{\mathrm{‾}}{Q}$ (station observations) and $\stackrel{\mathrm{‾}}{E}$ (${\stackrel{\mathrm{‾}}{E}}_{\text{WB}}$, Sect. ). The $p$ biases of the modeled climate region average $\stackrel{\mathrm{‾}}{Q}$ and $\stackrel{\mathrm{‾}}{E}$ are presented between brackets. Similar values for the individual catchments are shown in Tables S3 and S4.

	$\stackrel{\mathrm{‾}}{Q}$ ($\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$)			$\stackrel{\mathrm{‾}}{E}$ ($\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$)
Climate region	Observations	HTESSEL CTR	HTESSEL MD	WB	HTESSEL CTR	HTESSEL MD
Tropical	291	331 ($+\mathrm{13.7}$ $\mathrm{\%}$)	345 ($+\mathrm{18.6}$ $\mathrm{\%}$)	834	790 ($-$5.3 $\mathrm{\%}$)	776 ($-$7.1 $\mathrm{\%}$)
Temperate	56	59 ($+\mathrm{4.9}$ $\mathrm{\%}$)	70 ($+\mathrm{24.4}$ $\mathrm{\%}$)	626	624 ($-$0.4 $\mathrm{\%}$)	611 ($-$2.4 $\mathrm{\%}$)
Mediterranean	49	171 ($+\mathrm{249.9}$ $\mathrm{\%}$)	177 ($+\mathrm{263.8}$ $\mathrm{\%}$)	836	717 ($-$14.2 $\mathrm{\%}$)	709 ($-$15.2 $\mathrm{\%}$)

2.7 Model evaluation

Most study catchments are smaller than single HTESSEL grid cells (Table S1). For catchments completely falling within a single HTESSEL grid cell, this cell is selected for analysis. In the case that a catchment falls within more than one grid cell, the average of the model output in the separate grid cells is used for analysis. The model performances of CTR and MD are compared based on modeled monthly discharge and evaporation fluxes for 1975–2010: long-term annual means, monthly seasonal climatology and interannual anomalies of monthly fluxes (monthly fluxes minus monthly climatology) are evaluated. Modeled $Q$ is compared to station observations and modeled $E$ to the FLUXCOM-WB evaporation (Sect.  and Eq. ). For long-term annual means, the percent bias between the reference and modeled fluxes is calculated (evaporation $p\text{ bias}=({\stackrel{\mathrm{‾}}{E}}_{\text{mod}}-{\stackrel{\mathrm{‾}}{E}}_{\text{ref}})/{\stackrel{\mathrm{‾}}{E}}_{\text{ref}}$). For the monthly seasonal climatology and interannual anomalies, the model performance is quantified by using the Pearson correlation coefficient ($r$) and a variability performance metric ($v=(\mathrm{1}-\mathit{\alpha }{)}^{\mathrm{2}}$) that depends on the ratio of modeled and reference standard deviation ($\mathit{\alpha }={\mathit{\sigma }}_{\text{mod}}/{\mathit{\sigma }}_{\text{ref}}$). These performance metrics are calculated for the individual catchments and then averaged to evaluate model performance over tropical, temperate and Mediterranean climate regions.

To test the significance of the improvement in model performance of MD compared to CTR, a Monte Carlo bootstrap method (1000 repetitions) is employed. The 1000 samples are taken by randomly resampling with replacement among CTR and MD values at each time step. The null hypothesis of getting as high or higher performance parameters simply by chance is tested at the 5 $\mathrm{\%}$ and 10 $\mathrm{\%}$ significance levels for the individual catchments as well as for the performance averages over the tropical, temperate and Mediterranean climate regions. $P$ values of the model improvements are provided in Tables S5 and S6.

3 Results

3.1 Root zone storage capacity estimates

Figure  shows that there is no relation between $S_{\text{r,MM}}$ and $S_{\text{r,CTR}}$. The range of $S_{\text{r,MM}}$ (125–722 $\mathrm{mm}$) in the study catchments is much larger than the range of $S_{\text{r,CTR}}$ (491–725 $\mathrm{mm}$), indicating that HTESSEL may not adequately represent the spatial heterogeneity of $S_{\mathrm{r}}$ (Table S2). The range of $S_{\text{r,MM}}$ in the catchments is consistent with , who found similar ranges of $S_{\mathrm{r}}$ (approximately 100–600 $\mathrm{mm}$) over Australia by using gridded products of $S_{\mathrm{r}}$ based on rooting depths from observations and optimized inverse modeling, and they also found similar ranges of global $S_{\text{r,MM}}$ estimated based on satellite evaporation products. $S_{\text{r,MM}}$ estimates are on average smaller in the five temperate (194 $\mathrm{mm}$) catchments than in the three Mediterranean (321 $\mathrm{mm}$) and the seven tropical (437 $\mathrm{mm}$) catchments. In the tropical and Mediterranean regions vegetation needs to bridge extensive dry seasons as rainfall seasonality is high (Fig. , Table ), resulting in larger $S_{\text{r,MM}}$ than in temperate regions with year-round precipitation. In the Mediterranean, the average time lag between $P$ and $E_{\mathrm{p}}$ of 5.7 months results in large root zone storage deficits in the hot and dry summers and therefore larger $S_{\text{r,MM}}$ than in the temperate catchments.

Figure 5

Catchment $S_{\mathrm{r}}$ as estimated from the memory method ($S_{\text{r,MM}}$) compared to the HTESSEL CTR parameterization ($S_{\text{r,CTR}}$) in the catchment corresponding grid cells.

[Figure omitted. See PDF]

3.2 Long-term mean annual climatology

The HTESSEL CTR version overestimates observed $\stackrel{\mathrm{‾}}{Q}$ in 9 out of 15 catchments with on average 40 $\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$ (tropical), 3 $\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$ (temperate) and 122 $\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$ (Mediterranean) (Tables , S3 and S4). This overestimation of observed $Q$ goes together with an average underestimation of ${\stackrel{\mathrm{‾}}{E}}_{\text{WB}}$ by CTR. As $S_{\text{r,MM}}$ is generally smaller than $S_{\text{r,CTR}}$ (Fig. ), the MD version results in reduced $\stackrel{\mathrm{‾}}{E}$ and increased $\stackrel{\mathrm{‾}}{Q}$ compared to CTR, but the changes are quite small (Table ). The MD increase in modeled $\stackrel{\mathrm{‾}}{Q}$ compared to CTR results on average in larger $p$ biases in tropical ($+\mathrm{16.9}$ $\mathrm{\%}$ vs. $+\mathrm{13.7}$ $\mathrm{\%}$), temperate ($+\mathrm{24.4}$ $\mathrm{\%}$ vs. $+\mathrm{4.9}$ $\mathrm{\%}$) and Mediterranean ($+\mathrm{263.8}$ $\mathrm{\%}$ vs. $+\mathrm{249.9}$ $\mathrm{\%}$) catchments, but the results are largely variable among the individual catchments (Table S4).

Figure 6

Monthly seasonal climatology of observed discharge ($Q$) (a–c) and FLUXCOM-WB evaporation ($E_{\text{FLUXCOM-WB}}$) (d–f) as well as modeled values in the HTESSEL CTR and MD versions, averaged for the tropical (a, d), temperate (b, e) and Mediterranean (c, f) catchments for the time series 1975–2010. Labels (b1) and (c1) represent the same data as in (b2) and (c2), but with a different $y$ axis. Similar illustrations for the individual catchments are shown in Figs. S1 ($Q$) and S2 ($E$).

[Figure omitted. See PDF]

3.3 Monthly seasonal climatology

Although $\stackrel{\mathrm{‾}}{Q}$ does not considerably change in MD compared to CTR (Sect. ), MD reproduces the seasonal variations in $Q$ considerably better than CTR (Fig. a–c and Table ). In the tropical and Mediterranean catchments, MD increases $Q$ in the wet months, while it decreases $Q$ in the dry months compared to CTR and hence improves the seasonal timing of observed $Q$ (Fig. a, c and Table ). In the temperate catchments, MD increases $Q$ in the wet months (July–September) compared to CTR in accordance with observations, although in the other months the changes in MD compared to CTR are mixed (Fig. b). In terms of the correlation between modeled and observed monthly seasonal climatology, $Q$ improved in MD compared to CTR in 12 out of 15 catchments, with 7 catchments passing the 5 $\mathrm{\%}$ significance level for improvement (Table S5). For the climate region averages, the correlation significantly improved in MD from 0.84 to 0.90 (tropical), from 0.74 to 0.86 (temperate) and from 0.86 to 0.96 (Mediterranean) compared to CTR (Table ). On average, MD resulted in larger variations in monthly $Q$ than CTR (Fig. a–c). The variability term $v=(\mathrm{1}-{\mathit{\sigma }}_{\text{mod}}/{\mathit{\sigma }}_{\text{obs}}{)}^{\mathrm{2}}$ improved from 0.17 to 0.06 (tropical) and from 0.17 to 0.10 (temperate) in MD compared to CTR, but in the Mediterranean catchments the models strongly overestimate the observed variations in $Q$ (Fig. c), with the variability term increasing from 2.80 in CTR to 8.73 in MD (Tables  and S5).

Table 3

Model performance parameters of monthly seasonal discharge ($Q$) and evaporation ($E$) climatologies (1975–2010), with $r$ representing Pearson correlation and $v=(\mathrm{1}-\mathit{\alpha }{)}^{\mathrm{2}}$ variability (where $\mathit{\alpha }={\mathit{\sigma }}_{\text{mod}}/{\mathit{\sigma }}_{\text{obs}}$), in tropical, temperate and Mediterranean climate regions for the HTESSEL CTR and MD versions (catchment averages). Modeled $Q$ is compared to station observations and modeled $E$ to FLUXCOM-WB (Eq. ). For $r$, a value of 1 represents a perfect model; for $v$ a value of 0 represents a perfect model. The significance test of the MD improvements compared to CTR is represented by ${}^{**}$ (passing the 5 $\mathrm{\%}$ level) and ${}^{*}$ (passing the 10 $\mathrm{\%}$ level). Values of $r$ and $\mathit{\alpha }$ for the individual catchments and $p$ values of improvement are shown in Tables S5 ($Q$) and S6 ($E$).

		Discharge		Evaporation
Climate region	HTESSEL version	$r$ (–)	$v$ (–)	$r$ (–)	$v$ (–)
Tropical	CTR	0.84	0.17	0.98	0.07
	MD	0.90${}^{**}$	0.05${}^{**}$	0.98	0.07
Temperate	CTR	0.74	0.17	0.99	0.04
	MD	0.86${}^{**}$	0.10${}^{**}$	0.98	0.05
Mediterranean	CTR	0.86	2.80	0.81	0.08
	MD	0.96${}^{*}$	8.73	0.80	0.07

In contrast to the improvement in the monthly seasonal climatology of $Q$ in MD, the monthly seasonal cycle of $E$ appears not to be significantly affected, as shown in Fig. d–f and Table .

Figure 7

Interannual monthly anomalies of observed discharge ($Q$) (a, c, e) and FLUXCOM-WB evaporation ($E$) (b, d, f) fluxes as well as modeled values in the HTESSEL CTR and MD versions in an individual representative tropical (catchment Mi) (a, b), temperate (catchment Na) (c, d) and Mediterranean (catchment K) (e, f) catchment based on the time series for 1975–2010. Similar illustrations for the individual catchments are shown in Figs. S1 ($Q$) and S2 ($E$).

[Figure omitted. See PDF]

3.4 Interannual monthly anomalies

Figure a and c show that MD is better in capturing the variations in interannual $Q$ anomalies than CTR in the presented tropical and temperate catchments, while in the Mediterranean catchment both models strongly overestimate the interannual $Q$ anomalies compared to observations (Fig. e). In 14 out of 15 catchments, the variability in the interannual $Q$ anomalies increases in MD compared to CTR (Fig. S1 in the Supplement and Table S5). This results in an average improvement in the interannual anomaly variability ($v$) from 0.12 to 0.11 (tropical) and from 0.09 to 0.06 (temperate) in MD compared to CTR (Table ). However, in the Mediterranean catchments, the increased variability in the $Q$ anomalies leads to a strong overestimation of $Q$ anomalies with respect to observations (Figs. e and S1m–o), with $v$ increasing from 0.99 in CTR to 4.26 in MD. Figure a, c and e also show that the timing of the $Q$ anomalies improves in MD compared to CTR; in particular, the improved timing of the falling limbs is clearly visible in Fig. a and e. The interannual $Q$ anomaly correlation (corresponding to the timing) improves in 14 out of 15 catchments, with 9 catchments passing the 5 $\mathrm{\%}$ significance level for improvement (Table S5). On average, the correlation ($r$) increases from 0.74 to 0.78 (tropical), from 0.80 to 0.85 (temperate) and from 0.71 to 0.79 (Mediterranean) in MD compared to CTR. In contrast to the improvement in the interannual $Q$ anomalies in MD, the interannual $E$ anomalies do not considerably change compared to CTR (Fig. b, d and f; Tables  and S6).

4 Discussion

4.1 Synthesis of results

$S_{\text{r,MM}}$ is lower than $S_{\text{r,CTR}}$ in 14 out of 15 catchments (Fig. ). This is seemingly in contrast to literature suggesting that the root depth in land surface models is too low and that the absence of deep roots is a cause for uncertainties in simulated evaporation . However, $S_{\mathrm{r}}$ represents a conceptual water volume that is accessible to roots without defining where this volume is in reality. Therefore, it is not necessarily proportional to root depth as a small $S_{\mathrm{r}}$ does not preclude the presence of deep roots, as illustrated in Fig. 4 in .

The modeling results show that the difference in long-term mean $\stackrel{\mathrm{‾}}{Q}$ and $\stackrel{\mathrm{‾}}{E}$ fluxes between CTR and MD are small (Table ), whereas the differences between monthly (climatological and interannual) variations are clearly visible (Figs.  and ). This corresponds to other studies on catchment hydrology that suggest that the root zone storage mainly affects the fast hydrological response of a catchment . Furthermore, previous studies found larger improvements in modeled discharge using $S_{\text{r,MM}}$ in humid regions with large rainfall seasonality . This is not found in our study, as we obtain slightly smaller improvements in the discharge correlation for the tropical catchments than for the temperate and Mediterranean ones. This is at least partly related to the smaller difference between $S_{\text{r,MM}}$ and $S_{\text{r,CTR}}$ in the tropical catchments than in temperate and Mediterranean ones (Fig. ). The Mediterranean catchments have large climatological biases and overly large discharge variability in the seasonal cycle and interannual anomalies in CTR, and MD further degrades the performance with respect to bias and variability (Tables –). On the other hand, the correlation of seasonal climatology and interannual anomalies consistently improves in all climate regions with the implementation of $S_{\text{r,MM}}$. Therefore, it is suggested that aspects of the hydrology parameterization other than $S_{\mathrm{r}}$ (e.g., the lack of a groundwater layer) could be primarily leading to the large climatological biases and overly large discharge variability in the seasonal cycle and interannual anomalies in the Mediterranean. On the other hand, uncertainties in the GSWP-3 forcing could also partly cause the large biases in the Mediterranean. In this climate region, it is found that GSWP-3 $\stackrel{\mathrm{‾}}{P}$ (0.5${}^{\circ }$grid) is considerably larger than $\stackrel{\mathrm{‾}}{P}$ from the SILO dataset, which provides $P$ on a 0.05${}^{\circ }$grid directly derived from ground-based observational data .

Table 4

Model performance parameters of interannual monthly discharge ($Q$) and evaporation ($E$) anomalies (1975–2010), with $r$ representing Pearson correlation and $v=(\mathrm{1}-\mathit{\alpha }{)}^{\mathrm{2}}$ variability (where $\mathit{\alpha }={\mathit{\sigma }}_{\text{mod}}/{\mathit{\sigma }}_{\text{obs}}$), in tropical, temperate and Mediterranean climate regions for the HTESSEL CTR and MD versions (catchment averages). Modeled $Q$ is compared to station observations and modeled $E$ to FLUXCOM-WB (Eq. ). For $r$, a value of 1 represents a perfect model; for $v$ a value of 0 represents a perfect model. The significance test of the MD improvements compared to CTR is represented by ${}^{**}$ (passing the 5 $\mathrm{\%}$ level) and ${}^{*}$ (passing the 10 $\mathrm{\%}$ level). Values of $r$ and $\mathit{\alpha }$ for the individual catchments and $p$ values of improvement are shown in Tables S5 ($Q$) and S6 ($E$).

		Discharge		Evaporation
Climate region	HTESSEL version	$r$ (–)	$v$ (–)	$r$ (–)	$v$ (–)
Tropical	CTR	0.74	0.12	0.79	1.39
	MD	0.78${}^{**}$	0.11	0.80${}^{**}$	1.52
Temperate	CTR	0.80	0.09	0.81	1.12
	MD	0.85${}^{**}$	0.06${}^{*}$	0.82${}^{**}$	1.46
Mediterranean	CTR	0.71	0.99	0.78	1.17
	MD	0.79${}^{**}$	4.26	0.78	1.31

Although we found significant differences in modeled $Q$ between CTR and MD, the discrepancy in $E$ was very limited in all climate regions (Tables , and S6; Fig. S2). As stated before, the reliability of the FLUXCOM $E$ is questionable in our study catchments (Fig. ). Although the model performance with respect to $E$ fluxes is uncertain, the lack of evaporation sensitivity to $S_{\mathrm{r}}$ was unexpected and requires more in-depth evaluation of the results in view of the HTESSEL model parameterization.

Figure 8

Modeled transpiration and soil moisture content in the HTESSEL CTR and MD versions in an individual representative tropical (catchment Mi) (a), temperate (catchment Na) (b) and Mediterranean (catchment K) (c) catchment. From top to bottom: transpiration, relative difference between CTR and MD transpiration ($\frac{E_{\text{t,CTR}}-E_{\text{t,MD}}}{E_{\text{t,CTR}}}$), soil moisture layer 1, soil moisture layer 2, soil moisture layer 3, soil moisture layer 4. Additionally, the vegetation coverage ($C_{\mathrm{L}}$ and $C_{\mathrm{H}}$) and the relative rooting distribution ($R_k$) for the dominant high and low vegetation types are presented.

[Figure omitted. See PDF]

In order to further explain the evaporation (in)sensitivity, we analyzed the modeled soil moisture and specifically looked at a wet period (mid-1990) and a dry period (beginning of 1991) in a temperate catchment, as shown in Fig. b. During the wet period, soil moisture in the upper three layers is above or close to ${\mathit{\theta }}_{\text{cap}}$ for both MD and CTR, while in the fourth layer MD has larger soil moisture than CTR. In this case evaporation is not moisture-limited and controlled by the top three layers because of the larger root distribution in these layers (Eqs. 14 and 15). Therefore, the modeled transpiration is not sensitive to the increase in layer 4 soil moisture in MD compared to CTR. During the transition from wet to dry periods, the upper three layers dry out first as there is a reduction in precipitation input. As these layers are relatively dry, evaporation is controlled by the fourth layer in which $\mathit{\theta }$ is reduced to values close to ${\mathit{\theta }}_{\text{pwp}}$ in MD, while it remains relatively wet in CTR. It is this difference in ${\mathit{\theta }}_{\mathrm{4}}$ that causes the sensitivity of transpiration in MD during the wet-to-dry transition. However, most of the time the modeled soil moisture is in the wet and insensitive regime, and therefore the overall effect of MD on modeled evaporation tends to be small in the catchments considered in this study. To further analyze the evaporation sensitivity to $S_{\mathrm{r}}$ changes, it would be useful to evaluate to what extent it is model-dependent and compare HTESSEL behavior with other LSMs in a multi-model context (e.g., ). On the other hand, we also expect the evaporation sensitivity to $S_{\mathrm{r}}$ to be related to the methodology applied, which will be further discussed in Sect. .

4.2 Methodological uncertainty

Although the catchments were selected carefully, their location and sizes do not completely match the HTESSEL grid cells. Thus, assuming a one-to-one relation between precipitation, evaporation, river discharge and root zone storage capacities at the catchment and the grid cell is a potential source of error. However, this configures as random error and is therefore likely to cancel out in multiple catchment settings as is done in this study. Another source of uncertainty is the parameterization of the memory method for estimating catchment $S_{\mathrm{r}}$. This method requires estimations of maximum interception storage, seasonal and interannual transpiration signals, and return periods, which lead to differences in $S_{\text{r,MM}}$ when other values are chosen. A sensitivity analysis of $S_{\text{r,MM}}$ with a high $S_{\text{r,MM}}$ ($S_{\text{i,max}}=\mathrm{1.5}$ $\mathrm{mm}$, ${\text{RP}}_{\text{low}}=\mathrm{3}$ years, ${\text{RP}}_{\text{high}}=\mathrm{60}$ years, $f=\mathrm{0.15}$; see Appendix A) and a low $S_{\text{r,MM}}$ ($S_{\text{i,max}}=\mathrm{3.5}$ $\mathrm{mm}$, ${\text{RP}}_{\text{low}}=\mathrm{1.5}$ years, ${\text{RP}}_{\text{high}}=\mathrm{20}$ years, $f=\mathrm{0.35}$) on average deviated 45 $\mathrm{mm}$ from the average $S_{\text{r,MM}}$ estimates used in this study ($S_{\text{i,max}}=\mathrm{2.5}$ $\mathrm{mm}$, ${\text{RP}}_{\text{low}}=\mathrm{2}$ years, ${\text{RP}}_{\text{high}}=\mathrm{40}$ years, $f=\mathrm{0.25}$). This deviation is small considering the average $S_{\text{r,MM}}$ of 319 $\mathrm{mm}$. In addition, irrigation, as a possible external water source in catchments with crops (Table S1), and deep groundwater, as a water source for deep-rooting vegetation, are not accounted for in the approach. However, we think that the estimation of transpiration is the main uncertainty in the approach. The assumption that the seasonal variations in $E_{\mathrm{t}}$ and $E_{\mathrm{p}}$ are in phase may not hold in Mediterranean regions where $E_{\mathrm{p}}$ and $P$, and thereby the water available for transpiration, tend to be out of phase. Applying the seasonal pattern of transpiration modeled by CTR to the memory method in Mediterranean catchments results in smaller $S_{\text{r,MM}}$ estimates in these catchments (average: 292 $\mathrm{mm}$) than with the initial approach whereby the seasonality of $E_{\mathrm{t}}$ was based on $E_{\mathrm{p}}$ (average: 321 $\mathrm{mm}$). The relatively low deviation for both the parameter uncertainty and the uncertainty in the timing of $E_{\mathrm{t}}$ leads us to conclude that these assumptions have a small impact on the general finding that $S_{\text{r,MM}}$ is lower than $S_{\text{r,CTR}}$ and that HTESSEL does not represent the spatial heterogeneity of $S_{\mathrm{r}}$.

Station observations of river discharge are used in both the $S_{\text{r,MM}}$ estimation and the model evaluation. However, because the memory method is only based on observations of long-term annual mean discharge ($\stackrel{\mathrm{‾}}{Q}$) and the model evaluation is mainly based on the monthly seasonal and interannual variations in $Q$, we consider model evaluation based on these data appropriate.

4.3 Root zone storage capacity implementation

The HTESSEL CTR version does not explicitly formulate $S_{\mathrm{r}}$, and therefore we formulate $S_{\text{r,CTR}}$ based on the root zone parameterization as presented in Sect.  in order to modify the model parameters in a way to make the model consistent with the $S_{\text{r,MM}}$ estimates. This formulation represents the theoretical $S_{\text{r,CTR}}$, but it may not fully correspond to the soil moisture in the four layers that is actually used by the modeled vegetation. The effective $S_{\mathrm{r}}$ ($S_{\text{r,CTR,eff}}$) can be derived a posteriori from the modeled soil moisture storage deficits and an extreme value analyses as done in the memory method (Sect. ). $S_{\text{r,CTR,eff}}$ is smaller than $S_{\text{r,CTR}}$ based on depths (Fig. S3c), which is likely related to the relatively small root percentage in layer 4 compared to the other layers for most vegetation types . On the other hand, the $S_{\text{r,MM}}$ we implemented in MD by changing soil depths is close the $S_{\text{r,MD,eff}}$ based on modeled soil moisture deficits in MD (Fig. S3d).

In MD the depths for soil moisture calculations are changed, directly resulting in changes in absolute soil moisture and thereby in indirect changes in discharge and transpiration. This modification is relatively simple and flexible, and there is no limitation on the possible range of soil depths for moisture calculations, so it could therefore similarly be implemented in other land surface models. However, it should be noted that this strategy chosen for changing the HTESSEL $S_{\mathrm{r}}$ is not the only one possible. As follows from Eq. (), the plant-available soil moisture (${\mathit{\theta }}_{\text{cap}}-{\mathit{\theta }}_{\text{pwp}}$) also defines the $S_{\mathrm{r}}$. However, modifications in the model's ${\mathit{\theta }}_{\text{cap}}$ or ${\mathit{\theta }}_{\text{pwp}}$ are not desired as these parameters are soil-texture-specific properties. Moreover, modifications in the formulations of the root-available moisture for each time step (${\mathit{\theta }}_{\text{roots}}$) appear not to be conceptually meaningful.

There are several alternative hypotheses that may potentially explain the limited sensitivity of modeled $E$ to the modified $S_{\mathrm{r}}$. First, the resistance of vegetation to transpiration is a function of the moisture supply (soil moisture) and the moisture demand (atmospheric condition) (Eqs. 14–16). The atmospheric conditions, which define moisture demand and thereby constrain transpiration, are similar in both CTR and MD because the models are run in an offline version. Therefore, the soil moisture–atmosphere feedback is not represented and the moisture demand side dominates the moisture supply side in the evaporation calculations. This issue could be overcome by using coupled climate simulations. Second, although $S_{\mathrm{r}}$ is changed in MD compared to CTR, the parameterization of the vegetation water stress is kept constant. found that different formulations of root water uptake considerably influence modeled water budgets, and therefore it is likely that changes in evaporation in MD compared to CTR are constrained by the vegetation water stress formulations (Eqs. 14–16). Third, the insensitivity of evaporation to the changes in model soil depth is probably also related to the fact that the resistance of vegetation to transpiration is a function of the relative soil moisture ($\mathit{\theta }$), which is not directly affected by changing the soil depth. On the other hand, soil depth changes directly affect the modeled $Q$, as modeled surface ($Q_{\mathrm{s}}$) and subsurface runoff ($Q_{\text{sb}}$) directly depend on the absolute moisture storage capacity of the soil (see Eqs. 10 and 12), with $Q_{\mathrm{s}}$ being a function of the absolute moisture in the top 50 $\mathrm{cm}$ of soil and $Q_{\text{sb}}$ a function of the total excess soil moisture when the layer's moisture content exceeds saturation moisture content. Fourth, monthly fluxes of $Q$ are often a full order of magnitude smaller than $E$. Hence, small changes in the partitioning simply add up to larger relative changes for $Q$.

5 Conclusions

This study is an attempt to overcome major limitations in the representation of the vegetation's root zone in land surface models. Specifically, we looked at the HTESSEL land surface model and found that the root zone storage capacity $S_{\mathrm{r}}$ is only a function of soil texture and soil depth, the latter being kept constant over the modeled global domain (in HTESSEL $z=\mathrm{2.89}$ $\mathrm{m}$), while from the state-of-the-art literature e.g., it is indicated that $S_{\mathrm{r}}$ is, to a large extent, climate-controlled. We found that the HTESSEL control version (CTR) indeed does not adequately represent the spatial heterogeneity of $S_{\mathrm{r}}$, with the range of $S_{\text{r,CTR}}$ (491–725 $\mathrm{mm}$) much narrower than the range obtained for the climate-controlled estimate $S_{\text{r,MM}}$ (125–722 $\mathrm{mm}$) in 15 Australian catchments with contrasting climate characteristics considered in this study. Furthermore, $S_{\text{r,CTR}}$ was found to be considerably larger than the climate-controlled estimate $S_{\text{r,MM}}$ in 14 out of 15 catchments. It is noted that these findings could be different for other LSMs when they have shallower soil depths.

We developed a new version of HTESSEL by suitably modifying the soil depths (MD) to obtain modeled $S_{\text{r,MD}}$ that matches $S_{\text{r,MM}}$ over the 15 catchments considered over Australia, while maintaining the overall HTESSEL model setup (Fig. ). This strategy to modify the model's $S_{\mathrm{r}}$ is relatively simple and could similarly be implemented in other land surface models. Moreover, the applied methodology could allow for a time-varying $S_{\mathrm{r}}$ in LSMs, and hence all four limitations of using sparse root observations mentioned in Sect.  could be overcome.

The comparison of the offline simulations with original (CTR) and modified (MD) versions of HTESSEL shows that the difference of the biases in the modeled long-term mean climatology of discharge and evaporation fluxes is generally small. On the other hand, the seasonal timing of the discharge flux is significantly improved in MD, indicating the beneficial effect of the climate-controlled representation of $S_{\mathrm{r}}$. Consistently, MD improves the correlation with observations for the monthly seasonal climatology of discharge fluxes in 12 out of 15 catchments (with 7 catchments passing the 5 $\mathrm{\%}$ significance level) and for the interannual monthly discharge anomalies in 14 out of 15 catchments (with 9 catchments passing the 5 $\mathrm{\%}$ significance level) (Table S5). Considering the climate region averages, the correlations of monthly seasonal climatology significantly improve in MD compared to CTR from 0.843 to 0.902 (tropical), from 0.741 to 0.855 (temperate) and from 0.860 to 0.951 (Mediterranean). The averaged correlations of the interannual monthly anomalies significantly improve in MD compared to CTR from 0.741 to 0.778 (tropical), from 0.795 to 0.847 (temperate) and from 0.705 to 0.785 (Mediterranean). Surprisingly, the modeled evaporation is shown to be relatively insensitive to changes in $S_{\mathrm{r}}$. In HTESSEL evaporation only depends on the relative moisture content in each soil layer, which in the model is not directly affected by the depth of the soil. Investigation of this insensitivity showed that it is only sensitive during dry periods when evaporation is dominated by transpiration from the fourth layer (Fig. ). On the other hand, surface runoff and subsurface runoff in HTESSEL depend on the total moisture content of the soil at any given time. Other than the relative moisture content this depends on the absolute moisture storage capacity of the soil that will vary together with the change in soil depth. Moreover, small changes in absolute fluxes translated to larger relative changes for runoff compared to evaporation (Fig. ).

As a final conclusion, we believe that a global application of climate-controlled root zone parameters has the potential to improve the timing of modeled water fluxes by land surface models, but from the results of this study a significant reduction of annual mean climatological biases cannot be expected. More work will be needed in the future to improve long-term mean simulations of discharge and evaporation fluxes by exploiting station-based and latest-generation satellite observations. Toward this aim the use of coordinated multi-model frameworks for the intercomparison of state-of-the-art LSMs could be fundamental.

Appendix A Iterative procedure for transpiration estimation

Daily transpiration is estimated by Eq. () with $c$, a coefficient that represents the ratio between transpiration and potential evaporation (Sect. ). With $c=\stackrel{\mathrm{‾}}{E_{\mathrm{t}}}/\stackrel{\mathrm{‾}}{E_{\mathrm{p}}}$ as a constant value, we do not account for interannual variability in transpiration caused by the interplay between atmospheric water demand and vegetation-available water supply. Therefore, we add an iterative procedure to estimate annually varying values for $c$, which is described here.

Steps 1 to 6 describe the procedure used to estimate $c$, with step 1 providing the initial estimates and steps 2 to 6 executed iteratively. $i$ represents the iterations (0–9) and $a$ the hydrological years (1973–2010). $P_{\mathrm{e}}$, $E_{\mathrm{t}}$, $E_{\mathrm{p}}$ and $S_{\mathrm{d}}$ are daily values. After 10 iterations ($i=$ 9) the resulting annual transpiration estimates stabilized and the corresponding storage deficits were used for the Gumbel $S_{\mathrm{r}}$ analysis as described in Sect. .

• Create initial estimates ($i=\mathrm{0}$) of $E_{\mathrm{t}}$ and $S_{\mathrm{d}}$ with a constant $c_{\mathrm{0},a}={\stackrel{\mathrm{‾}}{E}}_{\mathrm{t}}/{\stackrel{\mathrm{‾}}{E}}_{\mathrm{p}}$ for $a=$ 1973–2010. $$\begin{array}{}\text{A1}& {\displaystyle }& {\displaystyle }{E}_{\text{t,0}}=c_{\mathrm{0},a}{E}_{\mathrm{p}}\text{A2}& {\displaystyle }& {\displaystyle }{S}_{\mathrm{d},\mathrm{0}}=\text{max}\left(\mathrm{0},-\underset{\mathrm{1973}}{\overset{\mathrm{2010}}{\int }}(P_{\mathrm{e}}-E_{\mathrm{t},\mathrm{0}})\mathrm{d}t\right)\end{array}$$

• Calculate the annual change in storage in the root zone ($S$) with $t_{\mathrm{0}}$ and $t_{\mathrm{1}}$ representing the start and end of a hydrological year. A3$${\displaystyle }{\displaystyle }\mathrm{\Delta }{S}_{i,a}=S_{\mathrm{d},i}\left(t_{\mathrm{0}}\right)-S_{\mathrm{d},i}\left(t_{\mathrm{1}}\right)$$

• Calculate annual transpiration following the water balance. A4$${\displaystyle }{\displaystyle }{\stackrel{\mathrm{‾}}{E}}_{\mathrm{t},i,a}={\stackrel{\mathrm{‾}}{P}}_{\mathrm{e},i,a}-{\stackrel{\mathrm{‾}}{Q}}_{i,a}-{\displaystyle \frac{\mathrm{\Delta }{S}_{i,a}}{t_{\mathrm{1}}-t_{\mathrm{0}}}}$$

• Calculate $c_a$ for each hydrological year based on the annual $E_{\mathrm{t}}$ estimate from step 3 and calculate daily $E_{\mathrm{t}}$. $$\begin{array}{}\text{A5}& {\displaystyle }& {\displaystyle }{c}_{i,a}={\displaystyle \frac{{\stackrel{\mathrm{‾}}{E}}_{\mathrm{t},i,a}}{{\stackrel{\mathrm{‾}}{E}}_{\mathrm{p},a}}}\text{A6}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{t},i}=c_{i,a}{E}_{\mathrm{p}}\end{array}$$

• Calculate storage deficits based on daily $E_{\mathrm{t}}$ from step 4. A7$${\displaystyle }{\displaystyle }{S}_{\mathrm{d},i}=\text{max}\left(\mathrm{0},-\underset{\mathrm{1973}}{\overset{\mathrm{2010}}{\int }}(P_{\mathrm{e}}-E_{\mathrm{t},i})\mathrm{d}t\right)$$

• The input storage deficit of iteration $i+\mathrm{1}$ in step 2 is the average of iteration $i$ and $i-\mathrm{1}$. A8$${\displaystyle }{\displaystyle }{S}_{\mathrm{d},i+\mathrm{1}}={\displaystyle \frac{S_{\mathrm{d},i}+S_{\mathrm{d},i-\mathrm{1}}}{\mathrm{2}}}$$

The following three constraints are set to the iterations.

• The long-term water balance closes (${\stackrel{\mathrm{‾}}{P}}_{\mathrm{e}}-\stackrel{\mathrm{‾}}{Q}-{\stackrel{\mathrm{‾}}{E}}_{\mathrm{t}}\approx \mathrm{0}$).

• Annual transpiration is always larger than zero and smaller than the annual potential evaporation.

• Variations in $c$ are limited by $c_{\mathrm{0},a}-f{c}_{\mathrm{0},a}<c_{i,a}<c_{\mathrm{0},a}+f{c}_{\mathrm{0},a}$, with $f$ as a coefficient set to 0.25.

Figure  illustrates the iterative approach for storage deficit calculations. Daily $P$, $E_{\mathrm{p}}$ and $E_{\mathrm{t}}$ based on Eq. () are presented in Fig. a. Figure b shows annual variations of $P_{\mathrm{e}}$ and $E_{\mathrm{t}}$. During the years 1980–1984, $P_{\mathrm{e}}$ is clearly less than average, and the $E_{\text{t,0}}$ estimate is likely too high in these years because vegetation has less water available for transpiration. The final iteration $E_{\text{t,9}}$ provides a more realistic interannual pattern of transpiration. Initial and final iteration storage deficits are presented in Fig. c.

Figure A1

Storage deficit iteration approach in a temperate catchment for the time period 1977–1987. (a) Daily water fluxes with $P$ representing precipitation, $E_{\mathrm{p}}$ potential evaporation and $E_{\mathrm{t}}$ the initial transpiration calculation based on Eq. (). (b) Annual water fluxes with $P_{\mathrm{e}}$ representing effective precipitation, $E_{\text{t,0}}$ the initial transpiration estimate and $E_{\text{t,9}}$ the final iteration transpiration estimate. Mean $P_{\mathrm{e}}$ is based on the full time period (1973–2010). (c) Daily storage deficit with $S_{\text{d,0}}$ representing the initial calculation and $S_{\text{d,9}}$ the final iteration.

[Figure omitted. See PDF]

Code and data availability

Catchment discharge observations were taken from the Australian Bureau of Meteorology and can be downloaded from http://www.bom.gov.au/water/hrs/ (last access: September 2019) . FLUXCOM evaporation data were taken from the FLUXCOM initiative and can be downloaded from http://www.fluxcom.org/EF-Download/ (last access: December 2019) . Top-of-atmosphere radiation data were taken from Mines ParisTech and can be downloaded from http://www.soda-pro.com/web-services/radiation/extraterrestrial-irradiance-and-toa/ (last access: September 2019) . The offline HTESSEL model was provided by EC-EARTH, together with the GSWP-3 forcing data as well as vegetation and soil data. The adapted modules, model output and analysis codes are available upon request. The Python scripts used for $S_{\mathrm{r}}$ calculation and statistical significance of the results can be downloaded from https://github.com/fvanoorschot/Python-scripts-van-Oorschot-2021/ (last access: June 2021) .

The supplement related to this article is available online at: https://doi.org/10.5194/esd-12-725-2021-supplement.

Author contributions

The study was conceived by RJvdE and amended with input from all authors. FvO carried out the study, analyzed the results and wrote the paper with input and feedback from RJvdE, MH and AA. Specific knowledge and support for the $S_{\mathrm{r}}$ calculations were provided by MH, and specific knowledge and support for the EC-EARTH model were provided by AA.

Competing interests

The authors declare that they have no conflict of interest.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acknowledgements

Acknowledgement is given for the use of ECMWF's computing and archive facilities for this research, which were provided by the KNMI and by ECMWF in the framework of the special project SPITALES. This work was supported by the European Union's Horizon 2020 research and innovation program under grant agreement no. 101004156 (CONFESS project). Ruud J. van der Ent acknowledges funding from the Netherlands Organization for Scientific Research (NWO) under project number 016.Veni.181.015.

Financial support

This research has been supported by the Horizon 2020 project (CONFESS project no. 101004156) and the NWO (project no. 016.Veni.181.015).

Review statement

This paper was edited by Axel Kleidon and reviewed by Andrew Guswa, Stefan Hagemann, and one anonymous referee.