Global root zone storage capacity from satellite-based evaporation

. This study presents an “Earth observation-based” method for estimating root zone storage capacity – a critical, yet uncertain parameter in hydrological and land surface modelling. By assuming that vegetation optimises its root zone storage capacity to bridge critical dry periods, we were able to use state-of-the-art satellite-based evaporation data computed with independent energy balance equations to derive gridded root zone storage capacity at global scale. This approach does not require soil or vegetation information, is model independent, and is in principle scale independent. In contrast to a traditional look-up table approach, our method captures the variability in root zone storage capacity within land cover types, including in rainforests where direct measurements of root depths otherwise are scarce. Implementing the estimated root zone storage capacity in the global hydrological model STEAM (Simple Terrestrial Evaporation to Atmosphere Model) improved evaporation simulation overall, and in particular during the least evaporating months in sub-humid to humid regions with moderate to high seasonality. Our results suggest that several forest types are able to create a large storage to buffer for severe droughts (with a very long return period), in contrast to, for example, savannahs and woody savannahs (medium length return period), as well as grasslands, shrublands, and croplands (very short return period). The presented method to estimate root zone storage capacity eliminates the need for poor resolution soil and rooting depth data that form a limitation for achieving progress in the global land surface modelling community. horizontal heterogeneities, and uncertainties in soil proﬁle data including hard pans.

However, root zone storage capacity is very difficult to measure and observe in the field, especially at the larger scales that are relevant for many modelling needs. Rooting profiles measurements are also scarce, and difficult to generalise since vegetation rooting systems naturally adapt to pre-L. Wang-Erlandsson et al.: Global root zone storage capacity vailing climates and soil heterogeneities (e.g. Gentine et al., 2012;Sivandran and Bras, 2013). Even when rooting profiles are available, difficulties arise in translating them to root zone storage capacity, due to variations in root densities, hydrological activity, horizontal spatial heterogeneities, and uncertainties in soil profile data including hard pans.

Background
Broadly, six types of approaches to estimate the root zone storage capacity have been suggested or are in use in hydrological and land surface models: the field observation based approach, the look-up table approach, the optimisation approach, the inverse modelling approach, the calibration approach, and the mass balance based approach. These approaches are described below and compared in Table S1. Some of these approaches estimate rooting depth or root profiles, and can be translated to root zone storage capacity through combination with soil plant-available water (Sect. 3.2., Eq. 10), even though it is a simplification.
The field observation based approach provide estimates of rooting depths based on rooting depth measurements (Doorenbos and Pruitt, 1977;Dunne and Willmott, 1996;Jackson et al., 1996;Schenk and Jackson, 2002;Zeng, 2001) and has the advantage of being constructed from actual observations of vertical rooting distribution Jackson et al., 1996). To scale up rooting depth to the global scale, Schenk and Jackson (2002) used the mean biome rooting depth and Schenk and Jackson (2009) employed an empirical regression model based on reported root profiles from literature. However, this method suffers from data scarcity and location bias, and risks unlikely vegetation and soil combinations due to data uncertainty (Feddes et al., 2001). Moreover, it requires assumptions on water uptake from a certain fraction of the entire observed root profile. Observations show that many woody and herbaceous vegetation species are able to access very deep layers in a variety of soil conditions Stone and Kalisz, 1991), up to 18 m in Amazonian tropical forest (Nepstad et al., 1994), 53 m in the desert of the south-western USA (Phillips, 1963), and 68 m (possibly 140 m) in the central Kalahari dry savannah (Jennings, 1974). However, isolated roots that go very deep do not necessarily mean that vegetation across the landscape can exploit the full soil to that depth.
The look-up table approach is used in hydrological and land surface modelling to parameterise root zone storage capacity based on literature values of mean biome rooting depth and soil texture data (e.g. Müller Schmied et al., 2014;Wang-Erlandsson et al., 2014). This approach facilitates land cover change experiments and is grounded in literature, but assumes root zone storage capacity to be a function of merely land cover and soil type, with little consideration for climatic adjustments. This is a major oversight, as plants within the same vegetation type can exhibit a large span of root zone storage capacities in different climates and landscapes by adaptation to environmental conditions (Collins and Bras, 2007;Feldman, 1984;Gentine et al., 2012;Nepstad et al., 1994). Moreover, an incompatibility issue may arise if the literature based rooting depths employs a land cover classification different from that of the land surface model (Zeng, 2001).
The optimisation approach predicts vertical rooting depth based on soil, climate, and vegetation data, and assumptions about the soil hydraulic properties and root distribution behaviour. Often, optimal root profiles are derived based on maximised net primary production (Kleidon and Heimann, 1998a), carbon or transpiration gain (e.g. Collins and Bras, 2007;Schwinning and Ehleringer, 2001;van Wijk and Bouten, 2001), and sometimes also while being as shallow as possible (e.g. Laio et al., 2006;Schenk, 2008). The optimisation techniques used differ widely, including genetic algorithms (Schwinning and Ehleringer, 2001;van Wijk and Bouten, 2001), physical ecohydrological modelling (Collins and Bras, 2007;Hildebrandt and Eltahir, 2007), simple analytical modelling (Laio et al., 2006), and stochastic modelling (Schenk, 2008). This approach is powerful for improving the understanding of root profile development and can be useful for land surface models with explicit root distribution description (Smithwick et al., 2014). Nevertheless, further model development is needed to handle all types of environments (e.g. additional routines to handle groundwater uptake, acidic soil horizons, or low soil temperature) (Schenk, 2008).
The inverse modelling approach estimate rooting depth using a model to iteratively simulate a variable available from satellite data (e.g. net or gross primary production, absorbed photosynthetically active radiation, or total terrestrial evaporation) with different rooting depth parameterisations (Ichii et al., 2007(Ichii et al., , 2009Kleidon, 2004). This approach not only has a large spatial coverage while being indirectly observation-based, but also is dependent on soil information as well as the land surface model performance. Recently, this approach has also been applied at the local scale to approximate the root zone storage capacity by minimising differences between evaporation modelled from water balance and evaporation from remote sensing (Campos et al., 2016).
The calibration approach is widely used in hydrology, whereby a hydrological model is calibrated on the root zone storage capacity, using hydrological records on precipitation, runoff and evaporation, sometimes in combination with expert knowledge (e.g. Feddes et al., 1993;Fenicia et al., 2008;Jhorar et al., 2004;Winsemius et al., 2009;Gharari et al., 2014). However, the parameters derived are tied to the model used for calibration and are not necessarily comparable to measurable variables in nature, since they tend to compensate for uncertainties in model structure and data. In addition, since discharge is often the only observed variable (or one of only a few), the calibration approach is only suitable for applications at the catchment scale. For global hydrological models, parameters can be calibrated separately for a selection of gauged river basins and transferred to neighbour-L. Wang-Erlandsson et al.: Global root zone storage capacity 1461 ing ungauged catchments (Döll et al., 2003;Güntner, 2008;Hunger and Döll, 2008;Nijssen et al., 2001;Widén-Nilsson et al., 2007). This procedure, known as regionalisation, has (to our knowledge) only been performed for other parameter values than the root zone storage capacity, although the principle does not change with the parameters tuned. Nevertheless, challenges remain with discharge data uncertainty and parameter equifinality (Beven, 2006).
Recently, Gao et al. (2014) used a mass balance approachmore specifically, the mass curve technique -to estimate the root zone storage capacity at the catchment scale in the USA and in Thailand. The underlying assumption is based on the tested hypothesis that plants will not root deeper than necessary (Milly and Dunne, 1994;Milly, 1994;Schenk, 2008). The water demand during the dry season equaled a constant transpiration rate, which was obtained through a water balance approach together with a normalised difference vegetation index (NDVI). Their results suggested that ecosystems develop their root zone storage capacity to deal with droughts with specific return periods, beyond which the costs of carbon allocation to roots are too high from the perspective of the plants. This resonates well with past economic analyses of plant behaviour and traits (e.g. Givnish , 2014). Yet another mass balance approach was applied by de Boer-Euser et al. (2016) to catchments in New Zealand, using an interception and a root zone storage reservoir to record soil moisture storage deficit from variations in precipitation and transpiration. They derived mean annual transpiration from annual water balances, and seasonality of transpiration was added through estimate of potential transpiration and assumption about vegetation dormancy. The largest storage deficit of individual years were then used to derive catchment representative root zone storage capacity from the Gumbel extreme value distribution assuming dry spell return periods of 10 years. These two applications of the mass balance approach have the advantage of being both model independent and indirectly observation based. In addition, no land cover or soil information is needed, making the method parsimonious and flexible. Irrigation was, however, not considered and their assumption of ecosystem adaptation does not apply very well to seasonal crops (de Boer-Euser et al., 2016).
In a similar cumulative mass balance approach, van Dijk et al. (2014) combined a satellite evapotranspiration product with monthly precipitation data to estimate a mean seasonal storage range (MSSR) at 250 m resolution across Australia, as one of the inputs into national-scale mapping of groundwater-dependent ecosystems (http://www. bom.gov.au/water/groundwater/gde/). MSSR expresses the estimated mean seasonal range in the amount of water stored in all water stores combined (surface, soil and groundwater). A large range was considered likely to indicate a large use of water from storage during low rainfall periods from, for example, root water uptake from deeper soil or groundwater storages. Separate mapping of areas subject to irrigation or flood inundation was used to identify areas likely to rely on groundwater. The main conceptual drawback of this method is that the longer-term average seasonal pattern is likely to underestimate rooting depth in general, and even more so in regions without a strong seasonality in rainfall. The method also proved sensitive to any bias in evaporation and rainfall estimates and, in some conditions, simplifying assumptions about runoff and drainage rates (van Dijk et al., 2014).

Research aims
This study constitutes a first attempt to estimate global root zone storage capacity from satellite-based evaporation and precipitation data using a mass balance approach, which is possible thanks to recent development, testing, and validation of remote sensing evaporation products (e.g. Anderson et al., 2011;Guerschman et al., 2009;Hofste, 2014;Mu et al., 2011). Similar to the other mass-balancebased approaches, we assume that all hydrologically active roots are being used during the driest time and is not deeper than necessary. While we make use of the same mass balance principle as applied by Gao et al. (2014) and de Boer-Euser et al. (2016), our algorithm is based on indirect measurements of every unique pixels. Methodologically, in contrast to these two studies, the analyses here are carried out on global gridded data rather than by catchment and use total evaporation instead of interception and transpiration estimates.
Our aims are to (1) present a method for estimating root zone storage capacity using remote sensing evaporation and precipitation data at global scale that include the influence of irrigation; (2) evaluate how the new method influences evaporation simulation in a global hydrological model, in comparison to a classical look-up table approach; and (3) investigate the drought return periods different land cover types adjust to. This study, thus, provides an Earth observation-based and model-independent estimate of global root zone storage capacity that can be useful in models without the need for root distribution and soil information.

Estimating root zone storage capacity
The root zone storage capacity S R is estimated from soil moisture deficit D constructed from time series of water outflow F out and inflow F in from the root zone storage system. The algorithm is explained in this section and conceptually illustrated in Fig. 1.
First, we define the inflows and outflows from the system. The drying F out of the system is the total daily evaporation E: (1) Note that the total evaporation E is defined as the sum of transpiration, interception evaporation, soil moisture evaporation, and open-water evaporation. The shaded areas represent the accumulated differences A that are positive when outflow F out > inflow F in , and negative when F out < F in . Moisture deficit D is increased by positive A and decreased by negative A. Note that D never becomes negative.
The wetting F in of the system is the total daily precipitation P and the effective irrigation water F irr (i.e. additional evaporation from surface, wet soil, and ponding water at the tail end of irrigation borders that originates from irrigation): We need the term F irr in order to prevent S R from becoming overestimated in irrigated regions. This is because irrigation is captured in satellite-based evaporation data, but obviously not in precipitation data. Without correction, the irrigation evaporation in the satellite evaporation data would erroneously contribute to accumulation of soil moisture deficit in our computations. Beside irrigation, additional evaporation from natural non-soil water storages (e.g. floodplains, wetlands, and groundwater) may contribute to overestimation of soil storage dynamics (see also Sect. 4.5). In regions (see Appendix A) where the annual accumulated evaporation exceeds annual accumulated precipitation, also the long-term average of the difference of E − (P + F irr ) is added to F in in order to compensate for lateral inflow or estimation errors in evaporation or precipitation. Second, the difference between inflow and outflow is calculated at the daily scale. The accumulated difference A is represented by the shaded areas in Fig. 1 and can be defined as where t n is either the start of the accounting period or a point in time when F out = F in . Third, we calculate the moisture deficit D, being the shortage of water from rainfall: (4) The accumulation of D will occur in our algorithm only during periods where F out > F in , and reductions of D will occur when F out < F in . However, D never becomes negative by definition, since it can be considered a running estimate of the root zone storage reservoir size (see Fig. 1 at t 2 ). Not allowing negative D also means that any excess precipitation is assumed to be runoff or deep drainage. In this way, for every hydrological year, one maximum accumulated moisture deficit can be determined, representing the largest annual drought. A long time series of these maximum annual values creates the opportunity to study the return period of the maximum moisture deficits. Extreme values analysis, such as with Gumbel's method (Gumbel, 1935), then yield estimates of extreme moisture deficits with different probabilities of exceedance (see Sect. 2.3).
Finally, the root zone storage capacity (S R ) is defined as the maximum of the obtained D values: S R estimate based on an evaporation and precipitation time series would (in the absence of additional water supply) theoretically constitute a minimum root zone storage capacity (see Fig. S1 in the Supplement). If the root water uptake by plants does not abstract water until wilting point, the root zone storage may not utilise its full capacity. Note also that the S R computed is not to be confused with time-variable moisture availability. The time-variable water availability can be inferred from hydrological models using S R as the water holding capacity. During dry periods, the magnitude of surface runoff and deep drainage is usually small, and therefore is assumed to not affect root zone storage capacity calculations.

Implementation in a hydrological model
The newly derived root zone storage capacity is used in the global hydrological model STEAM (Simple Terrestrial Evaporation to Atmosphere Model) (Wang-Erlandsson et al., 2014) to evaluate its influence on evaporation simulation. STEAM is a process-based model that partitions evaporation into five fluxes (i.e. vegetation interception, floor interception, transpiration, soil moisture evaporation, and open-water evaporation). Potential evaporation is computed using the Penman-Monteith equation (Monteith, 1965), surface stomatal resistance is based on the Jarvis-Stewart equation (Stewart, 1988), and phenology is expressed as a function of minimum temperature, soil moisture content, and daylight (Jolly et al., 2005). The model operates at 1.5 • and 3 h resolution.
In STEAM, root zone storage capacity is originally calculated as the product of soil plant-available water (depending on soil texture) and rooting depth (depending on land cover type), using volumetric soil moisture as input to the stress function (here, the formulation of van Genuchten, 1980): where θ is the actual volumetric soil moisture content (dimensionless), θ wp is the volumetric soil moisture content at wilting point, and θ fc at field capacity. (This soil moisture stress function departs from the original formulation in STEAM (Matsumoto et al., 2008;Wang-Erlandsson et al., 2014), which is described in Sect. 2 in the Supplement.) However, the root zone storage capacity S R is simply location bound (depending on climatic variables alone) and no longer considered a land-cover-and soil-based parameter. Thus, to use S R directly, we do not account for soil moisture below wilting point and assume where h is the rooting depth (m). The reformulated stress function of soil moisture becomes where S is the actual root zone storage (m). This reformulation is possible since the stress function retains its shape. Thus, S R can in similar ways be implemented in other hydrological models.
To measure improvement, the root mean square error (ε RMS ) for simulated evaporation is calculated using the original look-up-table-based root zone storage capacity S R,STEAM and the newly derived root zone storage capacity S R,new (i.e. S R,CRU-SM or S R,CHIRPS-CSM ), respectively. The root mean square error improvement (ε RMS,imp ) is positive if the E simulated using S R is closer to a benchmark evaporation data set than the E simulated using S R,STEAM . The equation below shows the ε RMS,imp of S R,new : The remote-sensing-based ensemble evaporation product E SM (and E CSM , see Fig. S7) was used as benchmark E benchmark . This use may seem circular when E benchmark is used to derive S R,new , but is in fact valid due to differences in algorithms, precipitation input data, model types, and time span covered. First, the algorithms for estimating S R,new , and for estimating E in STEAM are very different. While S R,new is derived based on the E overshoot over P , STEAM is a process-based model where evaporation originates from five different compartments, each constrained by potential evaporation and related stress functions. This means that it is impossible to reproduce E benchmark simply by inserting S R,new to STEAM. Second, the precipitation products (CRU and CHIRPS respectively) used to derive S R,new differ from the precipitation forcing (ERA-I) used in STEAM. Third, E benchmark and STEAM are truly independent to each other as well. Whereas STEAM is process and water balance based, the ensemble E product is based on a combination of two (E SM ) or three (E CSM ) energy balance methods. Last, S R,new is based on a single year value of E benchmark (i.e. the year of maximum storage deficit), whereas the analyses of improvements are based on the entire available time series of 10-11 years. The only difference of the new STEAM simulations is the inclusion of updated information on root zone storage so that during longer periods of drought, more realistic estimations of continued evaporation processes can be expected. Thus, if S R,new dimensioned on 1 year of E benchmark nevertheless improves E simulation in STEAM with regard to 10-11 years of E benchmark (i.e. the overall RMS decreases when S R,new is used in STEAM) is a strong indication that the storage capacity correction was implemented for the right reason.
To investigate where the performance increases are most significant, improvements in mean annual, mean maximum monthly and mean minimum monthly E is calculated separately. ε RMS,imp by climate are done for bins of precipitation seasonality index and aridity index (defined in Appendix B) containing more than 200 grid cells. ε RMS,imp by land cover types are analysed for grid cells where single land cover occupancy exceeds 90 % in a 1.5 • grid cell. ε RMS analyses are carried out on area-weighted evaporation values to avoid bias caused by differences in grid cell areas. Results are shown in Sect. 4.4.

Frequency analysis
We calculate S R for 10 to 11 years (2003-2012 and 2003-2013, respectively; see Sect. 3.1) depending on data availability. However, different ecosystems may adapt their root system depths to different return periods of drought, which may or may not correspond to the available data time series length. Thus, we also determine the S R,L yrs for different return periods of drought L (see Sect. 4.4) based on Gumbel's distribution (Gumbel, 1935). The resulting S R,L yrs is a function of the mean and standard deviation of the extremes in the data series: where y n is the reduced mean as a function of the number of available years n (y 10 = 0.4952 and y 11 = 0.4996), σ n is the reduced standard deviation as a function of n (σ 10 = 0.9496 and σ 11 = 0.9676), σ S R is the standard deviation of S R , while y L is the reduced variate of the Gumbel distribution:

Evaporation and precipitation input for estimating S R
We present two S R data sets, one covering the latitudes 50 • N-50 • S (S R,CHIRPS-CSM ), and one with global coverage 80 • N-56 • S (S R,CRU-SM ). See Table 1 for an overview of the data input for each S R data set.  (Senay et al., 2013), and the MODIS evapotranspiration (MOD16) at 0.05 • (Mu et al., 2011). These three different evaporation models are all based on MODIS satellite data, but they use different parts of the electro-magnetic spectrum. CMRSET combines a vegetation index, which estimates vegetation photosynthetic activity, and shortwave infrared spectral data to estimate vegetation water content and presence of standing water. SSEBop relies on the thermal infrared data for determination of the latent heat flux, and MOD16 relies on the visible and near-infrared data to account for leaf area index variability. Hence, their input data, model structure, and output data are not necessarily similar, which makes them attractive for deriving an ensemble evaporation product. S R,CHIRPS-CSM is based on data covering the years 2003-2012 as CMRSET was not available for 2013.
For the global coverage S R,CRU-SM map, we used the 0.5 • Climatic Research Unit Time Series version 3.22 (CRU TS3.22) precipitation data (P CRU ) (Harris et al., 2014) together with the ensemble mean (E SM ) of only SSEBop and MOD16, since we found CMRSET to overestimate evaporation at high latitudes, possibly due to the effect of snow cover on estimates. In addition, the irrigation effect was analysed for S R,CRU-SM by including evaporation originating from ir-rigation water simulated at 0.5 • and at the daily scale by the dynamic global vegetation model LPJmL (Jägermeyr et al., 2015). S R,CRU-SM is computed based on evaporation data covering the years 2003-2013. Irrigation data cover the years 2003-2009 (monthly mean irrigation evaporation were used for years after 2009).
We present S R,CHIRPS-CSM , because P CHIRPS is the lead precipitation product and we can make use of three evaporation data sets. However, P CHIRPS is unfortunately not available at the global scale, and CMRSET is not reliable in high latitudes. Thus, we added the global-scale S R,CRU-SM to this study. This allows for application in global-scale models as well as investigations at the global scale (e.g. climate and land cover based analyses).
The input precipitation and evaporation data are shown in Figs. 2 and S2. This study required global coverage data at a grid cell resolution for both evaporation and precipitation. Importantly, these products must not be produced using assumptions on root zone storage capacity, to prevent circularity (since we are estimating root zone storage capacity). In other words, there should be no water balance type of computation process involved in the determination of S R . We used satellite-based evaporation products because they are the only options available that fulfill these criteria, (i.e. reanalyses and land surface model evaporation contain soil depth information, whereas FLUXNET data are too sparse for acquiring consistently good quality global coverage). The monthly satellite-based evaporation data used in the manuscript were those available at the time of this research. Conversely, precipitation data do not need to be satellite based, but can also be ground based. Intercomparisons of precipitation products show that both CRU and CHIRPS are good quality precipitation products. In particular, CHIRPS performance stands out in a comprehensive intercomparison of 13 different precipitation products in the Nile basin (Hessels, 2015). Nevertheless, data uncertainties still persist. The mean annual accumulated evaporation of E CSM and E SM is sometimes higher than the mean annual accumulated precipitation P CHIRPS and P CRU , which is discussed in Appendix A. The use of three evaporation data sets decreases uncertainties related to individual evaporation products, because there is simply not one single preferred model. To compare the effect of different input data, we also present results of S R based on the separate evaporation and precipitation data (Figs. S4 and S5).
In addition, ECMWF re-analysis interim (ERA-I) (Dee et al., 2011) daily 0.5 • evaporation and precipitation data were used to temporally downscale the monthly evaporation and precipitation data. In the temporal downscaling, we first established the ratios between daily values to the mean monthly ERA-I, and second, used the relationship to estimate daily values from monthly E SM or E CSM values. This allows for daily products of evaporation and precipitation, which was necessary in order to incorporate also short drought periods.

Other data used in analyses
The following data sets were compared with our S R estimates: -the estimated 1 • rooting depth for 95 % of the roots from Schenk and Jackson (2009); -the 1 • rooting depth estimated by the optimised inverse modelling from Kleidon (2004), (where the minimum rooting depth producing the long-term maximum net primary production is selected as the best estimate); -the 1 • rooting depth estimated by the assimilated inverse modelling from Kleidon (2004), (where the rooting depth that minimises the difference between the modelled and the satellite-derived absorbed photosynthetically active radiation is selected as the best estimate); -the root zone storage capacity look-up table-based parameterisation used in a global hydrological model, i.e. the STEAM (Wang-Erlandsson et al., 2014).

L. Wang-Erlandsson et al.: Global root zone storage capacity
In order to enable comparison between rooting depth h and root zone storage capacity S R , we assumed that the root zone reaches its wilting point and converted between h and S R using soil properties: where θ paw is the maximum plant available soil moisture, θ fc is the volumetric soil moisture content at field capacity, and θ wp is the volumetric soil moisture content at wilting point. Soil texture data at 30 is taken from the Harmonised World Soil Database (HWSD) (FAO/IIASA/ISRIC/ISSCAS/JRC, 2012), and field capacity and wilting point information is based on the US Department of Agriculture (USDA) soil classification (Saxton and Rawls, 2006). To analyse if and how the inferred S R may improve simulations in a hydrological model, we applied S R,CRU-SM to the evaporation simulation model STEAM. Input ERA-I data to STEAM were at 3 h and 1.5 • resolution and include: precipitation, snowfall, snowmelt, temperature at 2 m height, dew point temperature at 2 m height, wind speed vector fields (zonal and meridional components) at 10 m height, incoming shortwave radiation, net long-wave radiation, and evaporation (only used to scale potential evaporation from daily to 3 h). To analyse the improvements in simulated evaporation by using S R,CRU-SM as input to STEAM (see Sect. 4.3), we used an aridity index based on precipitation and reference evaporation from CRU TS3.22 (Harris et al., 2014).
For land cover-based analyses, we used the 0.05 • landcover-type climate modelling grid (CMG) MCD12C1 created from Terra and Aqua MODIS data (Friedl et al., 2010) for the year 2008, based on the land cover classification according to the International Geosphere-Biosphere Programme (IGBP) (shown in Fig. S3). Land cover fractions are preserved in upscaling to 0.5 • . Only 0.5 • grid cells containing at least 95 % of a single land cover type are used in the land cover-based analyses (see Sect. 4.2.2) and grid cells containing at least 95 % water are removed from all S R analyses.
Data with finer resolution than 0.5 • have been upscaled to 0.5 • by simple averaging (i.e. assuming that the value of a 0.5 • grid cell correspond to the mean of the overlapping finer grid cell values). Data with coarser resolution than 0.5 • were downscaled by oversampling (i.e. transferring grid cell values assuming that the finer 0.5 • grid cell values correspond to those overlapped by the coarser degree grid cell values). Figure 3 shows the S R,CHIRPS-CSM (clipped, based on E CSM and P CHIRPS ) and S R,CRU-SM (global, based on E SM and P CRU ) estimates adjusted for irrigation (provided in the Supplements as ASCII files). Independent of the input data used, large root zone storage capacities are observed in the semiarid Sahel, South American and African savannah, central USA, India, parts of Southeast Asia, and northern Australia. The lowest root zone storage capacities are observed in the most arid and barren areas, and in the humid and densely vegetated tropics. The largest differences between S R,CHIRPS-CSM and S R,CRU-SM are observed over the Amazon, along the Andes, in Central Asia, and in the Sahara. Along mountain ridges (for example along the Andes and Himalaya), the S R estimates are generally large, possibly due to data uncertainty in these transition regions or evaporation in foothills sustained by lateral water fluxes from the mountains in addition to precipitation. The positive values of S R,CHIRPS-CSM in the Sahara desert are caused by overestimation of evaporation in the CMRSET evaporation product (see also Figs. S4 and S5). Notably, Fig. 3 shows contrasting root zone storage capacity over the South American and African tropical forests, although they belong to the same ecological class (i.e. evergreen broadleaf forests). This variability is purely due to temporal fluctuations between precipitation and evaporation and is independent of soil properties.  Fig. 4g) (Wang-Erlandsson et al., 2014). When the different data sets are compared to S R,CRU-SM (Fig. 4b, d, f and h), we see both agreements and significant differences. All data sets appear to more or less agree on the approximate range of root zone storage capacity in large parts of the Northern Hemisphere. Around the Equator, all data sets indicate root zone storage capacity to be lower or similar to that of S R,CRU-SM in the tropical forests of the Amazon and the Indonesian islands. In the Congo region and in Central America, S R,Kleidon,O and S R,Kleidon,A are larger than both S R,CRU-SM and the other. In the south temperate zone, S R,CRU-SM appears to correspond to or be lower than the other data sets. Figure 4 also reveals patterns specific to the different data sets that can be explained by the underlying method used for estimating rooting depth or root zone storage. For example, both S R,Schenk and S R,STEAM contain spuriously large values in the deserts (such as the Sahara and the Gobi) where vegetation is non-existent or extremely sparse. The methods based on satellite data (S R,CRU-SM , S R,Kleidon,O and S R,Kleidon,A ) appear to reflect reality in deserts more accurately. The S R,Kleidon,O presents the largest root zone storage capacities (most pronounced over Africa, India, parts of South America), since this data set represents an idealised and optimised case. On the contrary, the smallest root zone storage capacities are presented in the Amazon rainforest by S R,Schenk . These smaller values could be due to the lack of observations, since S R,Schenk is derived from rooting depth field measurements. But any difference between rooting depth and root zone storage capacity could also be due to discrepancies between actual rooting depth and hydrologically active rooting depth (see also Sect. 3.2). In contrast to the other data sets, S R,STEAM is relatively homogenous and does not contain any large values (basically all < 400 mm) (Fig. 4g). This is natural, since the other data sets are based on more heterogeneous observations, whereas S R,STEAM is based on a homogenous look-up table. Nevertheless, different input data were also used in the different studies. Thus, it is difficult to attribute the variations in root zone storage capacity estimates to differences in methods or differences in input data. Additional comparisons in scatter plots and root mean square error are shown in Fig. S6 and Table S2. Figure 5 shows the root zone storage capacity distribution for different land cover types and S R data sets, (S R,CHIRPS-CSM is not shown since it does not have global coverage). Except for deciduous broadleaf forests, the S R,CRU-SM of forests (Fig. 5a-e) are closer to S R,Kleidon,O and S R,Kleidon,A than to S R,Schenk . Interestingly, the range of S R is large in the ever-green forest types for the "adaptive" estimates (S R,CRU-SM , S R,Kleidon,O , and S R,Kleidon,A ), but small for the literature based methods (S R,Schenk and S R,STEAM ). In open shrublands and grasslands ( Fig. 5f and i) root zone storage capacities are similar across all estimates, except for the higher S R,STEAM . In savannahs, croplands, and natural/vegetation mosaic areas (Fig. 5h, j, and k), S R,Kleidon,O , and S R,Kleidon,A appear to have higher values than others. In woody savannahs (Fig. 5g), S R,Kleidon,O has a notably large range as well as high mean root zone storage capacity. In barren land (Fig. 5l), S R,Schenk and S R,STEAM are counter-intuitively high.

Implementation in a hydrological model
We implemented S R,CRU-SM , S R,CHIRPS-CSM , and S R,STEAM in the hydrological model STEAM (see Sect. 2.2 for methods) in order to analyse how the new root zone storage capacities might improve model performance. This section shows the performance analyses using S R,CRU-SM as input, since it has global coverage. A comparison in E simulation performance between using S R,CHIRPS-CSM and S R,CRU-SM as input to STEAM is shown in Fig. S7, and discussed in the Supplement. Figure 6 compares the STEAM-simulated evaporation when using, on the one hand, S R,CRU-SM and, on the other hand, the look-up table-based S R,STEAM . In general, S R,CRU-SM estimated higher evaporation rates in the tropics and lower evaporation in the subtropics and temperate zone. In particular, the differences are pronounced during the warm and dry seasons. For example, the evaporation reductions with S R,CRU-SM is widespread in the Northern Hemisphere during July. During the dry seasons (e.g. January in the Sahel, July in Congo south of the Equator), the evaporation increase is the most significant. Moreover, the changes in evaporation also depend on land cover type. In South America, evaporation increases in the seasonal tropical forests of the Amazon, whereas evaporation decreases in the savannahs and shrublands in the south. These results suggest that S R,CRU-SM has the greatest potential to influence model simulations for the hot and dry seasons, in regions where the root zone storage varies strongly. Figure 7 shows the ε RMS improvements of simulated mean annual, mean maximum monthly, and mean minimum monthly E sorted by seasonality and aridity, using S R,CRU-SM as input and E SM as a benchmark. The analysis reveals that our S R,CRU-SM estimate has the greatest potential to improve model simulations for minimum monthly evaporation. In particular, the improvements become significant with increased seasonality of rainfall, and in sub-humid to humid regions, resonating the findings of de Boer-Euser et al. (2016).

The effect of different drought return periods
Vegetation may adapt to a different time period than the 10-11 years of data that were available for this study. Thus, Figure 5. Comparison of root zone storage capacity estimates by land cover type using Tukey box plots. The central markers of the boxes mark the median, and the box edges mark the 25th and 75th percentile. The whiskers extend to 1.5 times the interquartile range.
we normalised S R using the Gumbel distribution in order to assess the effect of different drought return periods (see Sect. 2.3). Normalised S R are provided in the Supplements as ASCII files. Figure 8 shows the mean latitudinal S R,CHIRPS-CSM,L yrs and S R,CRU-SM,L yrs for different drought return periods L based on the Gumbel distribution. As may be expected, both S R,CHIRPS-CSM and S R,CRU-SM based on the 10-11 years where data were available correspond most closely to the S R,L yrs for L = 10 years (S R,10 yrs ). S R,L yrs always increases with L, but more strongly for small L and less so for large L following the Gumbel distribution. The largest spans are seen in the northern latitudes and around the Equator. Figure 9 shows a comparison of how Gumbel normalised S R,CRU-SM,L yrs affect the evaporation simulation ε RMS improvements by land cover type. Interestingly, a drought return period of 2 years (S R,CRU-SM, 2 yrs ) offers the best evaporation simulation performance in deciduous broadleaf forests, open shrublands, grasslands, croplands, and barren lands, whereas S R,CRU-SM,10 yrs or S R,CRU-SM,20 yrs are best in evergreen needleleaf forests, woody savannahs, and sa-vannahs, and S R,CRU-SM,60 yrs is best in evergreen broadleaf forests, deciduous needleleaf forests, and mixed forests.
A short drought return period of 2 years improves evaporation simulation the most in short vegetation types probably because these land cover types adapt to average years rather than to extreme drought years. In extreme years, they survive by going dormant. Evergreen broadleaf forests, on the other hand, adapt to 40-60 years of drought return period since they deal with droughts by accessing deeper soil moisture storages and thus invest in root growth (Brunner et al., 2015). The performance increases in deciduous needleleaf forests by using 60 years of drought return period could be explained by their need to cater to dry periods during their most active summer months. Shedding the leaves during the wet season (semi-arid tropics) or the growing season (summer in temperate climates) is not attractive because it prevents reproduction. Interestingly, deciduous broadleaf forests appear to adapt to a 2-year drought return period -i.e. radically different to deciduous needleleaf forests. This is possibly due to their younger age (Poulter, 2012;Hicke et al., 2007) and considerably shorter longevity (Larson, 2001;Loehle, 1988). Longevity could be explained by strong defence mechanisms against fungi and insects, lack of physical environmental damage, as well as low occurrence of environmental stress such as drought (Larson, 2001). Thus, it seems logical that the older and longer living deciduous needleleaf forests have developed their root zone storage capacities to stand against more extreme droughts. Analysing the performance of each land cover type not only reveals interesting patterns (such  as the contrast between deciduous needleleaf and broadleaf forests), but also leads to small sample sizes (particularly for evergreen needleleaf forests and the deciduous forest types) that should be considered when interpreting the results.
Based on the best performing drought return periods for each land cover type, we created a Gumbel normalised root zone storage capacity map (Fig. S9), which is shown and analysed in Sect. 3 in the Supplement. In addition, we also analyse how S R of different land cover types can be associated with climatic indicators in Appendix B.

Limitations
Although research indicates that most ecosystem rooting depth are limited by water rather than other resources (Schenk, 2008), other factors may still cause S R to be larger than what is considered here. A minimum rooting depth of 0.3-0.4 m are for example considered in Schenk and Jackson (2009). Although we are comparing others' rooting depth estimates to S R,CRU-SM , they are not directly comparable. Our approach deals with the accessible water volume in the root zone, which is not always related to root zone depth since the root density can vary over the depth. Our S R estimates  Wang-Erlandsson et al., 2014) in the global hydrological model STEAM, where the satellite-based E SM was used as the benchmark for improvements (see methods described in Sect. 2.2). The improvements for root zone storage capacities with different return periods L 2-60 years (i.e. S R,CRU-SM,2 yrs , S R,CRU-SM,5 yrs , S R,CRU-SM,10 yrs , S R,CRU-SM,20 yrs , S R,CRU-SM,40 yrs , and S R,CRU-SM,60 yrs ) are shown for the different land cover types that have > 90 % grid cell coverage. The number of represented grid cells are provided in the parenthesis following each land cover type label along the x axis.
implicitly capture the root density that is active in water uptake.
The S R,CHIRPS-CSM and S R,CRU-SM have been derived using evaporation and precipitation data from recent years (i.e. the 2000s), and should be used with caution if applied to past or future model simulations. Land cover change during the years 2003-2013 have not been taken into account. This has a potential impact on the computation of additional evaporation from irrigated areas with fast changing acreage.
Wetlands and groundwater-dependent ecosystems produce additional evaporation that cannot be ascribed to local rainfall (van Dijk et al., 2014). Bastiaanssen et al. (2014) recently demonstrated for the Nile basin that in some areas, natural withdrawals exceed man-made withdrawals to the irrigation sector. Since satellite evaporation data captures all types of evaporation, and we only corrected for irrigation, natural additional evaporation sources are implicitly included in S R,CHIRPS-CSM and S R,CRU-SM . Thus, our S R estimates may not strictly represent the root zone storage capacities in regions where water uptake from groundwater is significant (see Fig. A1).
The choice of remotely sensed evaporation products influenced the resulting S R more than the choice of precipitation product in this study (see Fig. S4). In particular, the largest standard deviations in the ensemble evaporation products are located in central South America, the Sahel, India, and northern Australia (see Fig. 2e and f). To reduce uncertainty, the presented method is preferably applied using ensemble products based on reliable evaporation and precipitation data sets identified in comparison and evaluation studies (e.g. Hofste, 2014;Yilmaz et al., 2014;Trambauer et al., 2014;Bitew and Gebremichael, 2011;Herold et al., 2015;Hessels, 2015;Moazami et al., 2013;Trenberth et al., 1991).
Finally, while the S R estimates are model independent, the analyses of the best performing drought return periods of different land cover types will depend on the hydrological model used, given the large variations of evaporation estimates (and in particular transpiration/evaporation ratios) among land surface models (e.g. Wang and Dickinson, 2012). Thus, although the contrasting return periods for woody land cover types and annual short vegetation types reported here are supported by current knowledge about ecohydrological response to droughts, the calculated values are subject to assumptions. Uncertainties are probably largest for heterogeneous land cover types (such as savannahs) because they tend to be challenging to parameterise and simulate. Therefore, implementation of S R in other hydrological or land surface models would require model-specific analyses of optimal return periods.

Summary and conclusion
This study presents a method to estimate root zone storage capacity in principle from remotely sensed evaporation and observation-based precipitation data, by assuming that plants do not invest more in their roots than necessary to bridge a dry period. Two global root zone storage estimates (S R,CRU-SM and S R,CHIRPS-CSM ) are presented based on different precipitation and evaporation data sets, but show in general similar patterns globally. S R,CRU-SM and S R,CHIRPS-CSM both improve mean annual E simulation in STEAM (see Fig. S7), and there is not a preferred product.
Different ecosystems have evolved to survive droughts of different return periods with different strategies. Our analyses showed that whereas long drought return periods increased performance for many forest types, short drought return periods increased performance for many short vegetation types. The best E simulation results were achieved when normalising the S R estimate using a very short drought return period (∼ 2 years) for deciduous broadleaf forests, grasslands, shrublands, croplands, and barren or sparsely vegetated lands, a medium length drought return period (∼ 10-20 years) for evergreen needleleaf forests, woody savannahs, and savannahs, and a very long drought return period (∼ 60 years) for evergreen broadleaf, deciduous needleleaf, and mixed forests. This is probably because grasslands survive extreme droughts by going dormant, whereas forests invest in root growth (Brunner et al., 2015). Thus, the root zone storage capacities of short vegetation types seem to adapt to average years, whereas those of forests adapt to extreme years. Differences among forest types are thought to be related to forest age and drought coping strategy. Normalisation to longer drought return periods should not be done for short-lived annual plants such as two-third of the world's croplands (Cox et al., 2006), nor beyond the age of the ecosystem of concern, because vegetation can not be assumed to adapt beyond their age.
The S R estimates presented here are both globally gridded and observation based. They have an advantage over the field-study-based and statistically derived S R,Schenk (Schenk and Jackson, 2009) by being directly based on gridded data and by covering regions where observational studies are limited (e.g. the evergreen broadleaf forests). In comparison to the inverse modelling approaches of Kleidon (2004), the method presented in this study is independent of model simulations and therefore closer to direct observations.
The new S R estimates can be used in hydrological and land surface modelling to improve simulation results, particularly in the dry season and in seasonal tropical forests where variations of root zone storage capacity are large. Using the new S R as input to the hydrological model STEAM improved the evaporation simulation considerably in subhumid to humid regions with high seasonality. In particular, the most significant improvements occurred in the months with the least evaporation. Normalisation of S R to different drought return periods for different land cover types could further improve evaporation simulation in STEAM, suggesting that Gumbel normalisation is a viable method to optimise the S R estimates prior to implementation in global hydrological or land surface models.
The presented method is simple to apply and in principle scale independent. For researchers working at regional or local scales, root zone storage capacities can easily be derived using available evaporation and precipitation data. Moreover, when information on irrigation and groundwater use is available, they can be used to adjust S R , as was done for example by van Dijk et al. (2014). Satellite-based evaporation data sets are also quickly being developed and improved. New global-scale evaporation products such as ALEXI (Atmosphere-Land Exchange Inverse) (Anderson et al., 2011) and ETMonitor  are underway based on 375 and 1000 m pixels. More sophisticated twolayer surface energy balance models also have the capacity to distinguish transpiration from other forms of evaporation. This implies that local root zone storage capacity can be computed, based on transpiration fluxes, which is preferred from a bio-physical point of view (although it would require estimate of interception evaporation to calculate effective precipitation). As new evaporation data sets become available, the S R estimates can easily be updated. In addition, this method can be used to diagnose and compare different evaporation products, in particular for identifying variations in seasonality. With longer time series of land cover and climate data, this method can possibly also be used to infer the effect of climate change on root zone storage capacity as a function of the adaptability of vegetation to altered conditions.

L. Wang-Erlandsson et al.: Global root zone storage capacity Appendix A: Evaporation exceedance over precipitation
The mean annual accumulated evaporation of E CSM and E SM is sometimes higher than the mean annual accumulated precipitation P CHIRPS and P CRU (see Fig. A1). In these areas, overestimation of S R may be expected because it is unlikely that the 10 or 11 year accumulation of E is more than rainfall, except for hydrological situations with lateral inflow through inundation, irrigation, or groundwater inflow. The evaporation data set E CSM exhibits larger and more widely spread exceedance over P CHIRPS in comparison to the E SM − P CRU combination. Most notably, the exceedance is high and potentially spurious in arid and semi-arid zones (e.g. the Sahara, western USA, and Central Asia), which suggests that the evaporation from deserts is not accurate. Regions where both S R,CRU-SM and S R,CHIRPS-CSM show high accumulated evaporation exceedance are along the Andes, patches in western USA, East Africa, Ivory Coast, Central Asia, northwest China and spots in Australia. These are essentially irrigated areas, lakes, reservoirs, wetlands, and coastal deltas. Possibly, overestimation of S R can also be caused for example by vegetation tapping into groundwater. Uncertainty in evaporation and precipitation products also propagates to errors in S R . The uncertainty of evaporation is location specific, (grid cells with a large standard deviation between the individual E products are shown in Fig. 2e and f).
Interestingly, the high evaporation exceedance appears to be much more pronounced during drier years. In Fig. A2, we sort every grid cell by the annual precipitation amount, from dry to wet, and plot the mean latitudinal E exceedance for the regions where the long-term accumulated E − P is positive. The figure clearly shows that E exceedance decreases with increase in rainfall, indicating that increased water demand during dry years is satisfied by withdrawing moisture from the soil matrix that is bounded with more potential (higher pF), or from underlying groundwater through deeply rooting vegetation.

B1 Methods and data
We analyse how S R,CRU-SM of different land cover types can be associated with climatic indicators. Stepwise multiple regression method based on the Akaike information criterion (AIC) is used to analyse how these climatic indicators may explain variations in S R within a land cover type. The climatic indicators used are precipitation seasonality (I s ), aridity (I a ), and interstorm duration (I isd ) (as these were found to be important by Gao et al., 2014): and where P m is the mean precipitation of the month, P a is the mean annual precipitation, and E p is the potential evaporation. We defined I isd as the mean continuous number of days per year without precipitation. Interaction effects between the variables are taken into account. The climate variables interstorm duration, aridity, and precipitation seasonality are developed based on monthly 0.5 • reference evaporation from CRU TS3.22 (Harris et al., 2014) and monthly 0.5 • precipitation for 1982-2009 from the Global Precipitation Climatology Centre (GPCC) (Schneider et al., 2011). Here, GPCC data (instead of CRU) are used in order to prevent false correlation with the CRU-based S R,CRU-SM . Figure B1. Coefficient of determination R 2 of the multiple linear regression model of S R,CRU-SM (based on P CRU and E SM ) based on the climate variables interstorm duration I isd , precipitation seasonality I s , and aridity I a . The green bars are forests or wooded land, the yellow bars represent croplands, and the teal bars represent short vegetation types.

B2 Results and discussion
We use multiple linear regression to correlate S R,CRU-SM values to climatic indicators, with the aim to investigate how well climate indicators can predict root zone storage capacities in different land cover types. It appears that climate indicators predict root zone storage capacities much better in evergreen forests than in short vegetation types. Figure B1 shows high R 2 in mostly evergreen forests; moderate R 2 in other forest types and croplands; and low R 2 in savannahs, shrublands, and grasslands. This is probably because of their different drought survival strategies. While evergreen forests bridge droughts with water uptake from storage in their root zone, deciduous forests shed their leaves, and short vegetation types such as grasslands go dormant and decrease their transpiration to a minimum. The multiple linear regression model for S R in croplands is moderately explained by climate indicators, potentially due to human management. All climate variables were selected by AIC in the multiple linear regression model (Table B1). Table B1. Predictor variables selected by Akaike information criterion (AIC) for the different land cover types. The predictor variables are interstorm duration (I isd ), precipitation seasonality (I s ), and aridity index (I a ).

Land cover type
Predictor variables 02: evergreen needleleaf forest S R = I isd + I s + I a + I isd : I a + I s : I a 03: evergreen broadleaf forest S R = I isd + I s + I a + I isd : I s + I s : I a 04: deciduous needleleaf forest S R = I isd + I s + I a + I s : I a + I isd : I a + I isd : I s + I isd : I s : I a 05: deciduous broadleaf forest S R = I isd + I s + I a + I s : I a + I isd : I s 06: mixed forests S R = I isd + I s + I a + I isd : I a + I isd : I s + I s : I a + I isd : I s : I a 08: open shrublands S R = I isd + I s + I a + I isd : I a + I isd : I s + I s : I a + I isd : I s : I a 09: woody savannas S R = I isd + I s + I a + I isd : I a + I s : I a 10: savannas S R = I isd + I s + I a + I isd : I a + I s : I a + I isd : I s + I isd : I s : I a 11: grasslands S R = I isd + I s + I a + I isd : I s + I s : I a 13: croplands S R = I isd + I s + I a + I isd : I a + I isd : I s + I s : I a + I isd : I s : I a 15: cropland/natural veg. mosaic S R = I isd + I s + I a + I isd : I s 17: barren or sparsely vegetated S R = I isd + I s + I a + I s : I a + I isd : I a + I isd : I s L. Wang-Erlandsson et al.: Global root zone storage capacity The Supplement related to this article is available online at doi:10.5194/hess-20-1459-2016-supplement.