Inter-annual variability of the global terrestrial water cycle

Variability of the terrestrial water cycle, i.e. precipitation (P), evapotranspiration (E), runoff (Q) and water storage change (ΔS) is the key to understanding hydro-climate extremes. However, a comprehensive global assessment for the partitioning of variability in P between E, Q and ΔS is still not available. In this study, we use the recently released global monthly hydrologic reanalysis product known as the Climate Data Record (CDR) to conduct an initial investigation of the inter-annual variability of the global terrestrial water cycle. We first examine global patterns in partitioning the long-term mean $\stackrel{\mathrm{‾}}{P}$ between the various sinks $\stackrel{\mathrm{‾}}{E}$, $\stackrel{\mathrm{‾}}{Q}$ and $\stackrel{\mathrm{‾}}{\mathrm{\Delta }S}$ and confirm the well-known patterns with $\stackrel{\mathrm{‾}}{P}$ partitioned between $\stackrel{\mathrm{‾}}{E}$ and $\stackrel{\mathrm{‾}}{Q}$ according to the aridity index. In a new analysis based on the concept of variability source and sinks we then examine how variability in the precipitation ${\mathit{\sigma }}_P^{\mathrm{2}}$ (the source) is partitioned between the three variability sinks ${\mathit{\sigma }}_E^{\mathrm{2}}$, ${\mathit{\sigma }}_Q^{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}$ along with the three relevant covariance terms, and how that partitioning varies with the aridity index. We find that the partitioning of inter-annual variability does not simply follow the mean state partitioning. Instead we find that ${\mathit{\sigma }}_P^{\mathrm{2}}$ is mostly partitioned between ${\mathit{\sigma }}_Q^{\mathrm{2}}$, ${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}$ and the associated covariances with limited partitioning to ${\mathit{\sigma }}_E^{\mathrm{2}}$. We also find that the magnitude of the covariance components can be large and often negative, indicating that variability in the sinks (e.g. ${\mathit{\sigma }}_Q^{\mathrm{2}}$, ${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}$) can, and regularly does, exceed variability in the source (${\mathit{\sigma }}_P^{\mathrm{2}}$). Further investigations under extreme conditions revealed that in extremely dry environments the variance partitioning is closely related to the water storage capacity. With limited storage capacity the partitioning of ${\mathit{\sigma }}_P^{\mathrm{2}}$ is mostly to ${\mathit{\sigma }}_E^{\mathrm{2}}$, but as the storage capacity increases the partitioning of ${\mathit{\sigma }}_P^{\mathrm{2}}$ is increasingly shared between ${\mathit{\sigma }}_E^{\mathrm{2}}$, ${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}$ and the covariance between those variables. In other environments (i.e. extremely wet and semi-arid–semi-humid) the variance partitioning proved to be extremely complex and a synthesis has not been developed. We anticipate that a major scientific effort will be needed to develop a synthesis of hydrologic variability.

1 Introduction

In describing the terrestrial branch of the water cycle, the precipitation (P) is partitioned into evapotranspiration (E), runoff (Q) and change in water storage (ΔS). With averages taken over many years, $\stackrel{\mathrm{‾}}{\mathrm{\Delta }S}$ is usually assumed to be zero and it has long been recognized that the partitioning of the long-term mean annual precipitation ($\stackrel{\mathrm{‾}}{P}$) between $\stackrel{\mathrm{‾}}{E}$ and $\stackrel{\mathrm{‾}}{Q}$ was jointly determined by the availability of both water ($\stackrel{\mathrm{‾}}{P}$) and energy (represented by the net radiation expressed as an equivalent depth of water and denoted $\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}$). Using data from a large number of watersheds, Budyko (1974) developed an empirical relation relating the evapotranspiration ratio ($\stackrel{\mathrm{‾}}{E}/\stackrel{\mathrm{‾}}{P}$) to the aridity index ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}$). The resultant empirical relation and other Budyko-type forms (e.g. Fu, 1981; Choudhury, 1999; Yang et al., 2008; Roderick and Farquhar, 2011; Sposito, 2017) that partition P between E and Q have proven to be extremely useful in both understanding and characterizing the long-term mean annual hydrological conditions in a given region.

However, the long-term mean annual hydrologic fluxes rarely occur in any given year. Instead, society must (routinely) deal with variability around the long-term mean. The classic hydro-climate extremes are droughts and floods but the key point here is that hydrologic variability is expressed on a full spectrum of time and space scales. To accommodate that perspective, we need to extend our thinking beyond the long-term mean to ask how the variability of P is partitioned into the variability of E, Q and ΔS (e.g. Orth and Destouni, 2018).

Early research on hydrologic variability focussed on extending the Budyko curve. In particular, Koster and Suarez (1999) used the Budyko curve to investigate inter-annual variability in the water cycle. In their framework, the evapotranspiration standard deviation ratio (defined as the ratio of standard deviation for E to P, σ_E∕σ_P) was (also) estimated using the aridity index ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}$). The classic Koster and Suarez framework has been widely applied and extended in subsequent investigations of the variability in both E and Q, using catchment observations, reanalysis data and model outputs (e.g. McMahon et al., 2011; Wang and Alimohammadi, 2012; Sankarasubramanian and Vogel, 2002; Zeng and Cai, 2015). However, typical applications of the Koster and Suarez framework have previously been at regional scales and there is still no comprehensive global assessment for partitioning the variability of P into the variability of E, Q and ΔS. One reason for the lack of a global comprehensive assessment is the absence of gridded global hydrologic data. Interestingly, the atmospheric science community have long used a combination of observations and model outputs to construct gridded global-scale atmospheric re-analyses and such products have become central to atmospheric research. Those atmospheric products also contain estimates of some of the key water cycle variables (e.g. P, E), such as in the widely used interim ECMWF Re-Analysis (ERA-Interim; Dee et al. 2011). Though efforts have been taken to develop land-based products from atmospheric reanalyses, e.g. ERA-Land (Balsamo et al., 2015) and MERRA-Land (Reichle et al., 2011) databases, the central aim of atmospheric re-analysis is to estimate atmospheric variables. That atmospheric-centric aim, understandably, ignores many of the nuances of soil water infiltration, vegetation water uptake, runoff generation and many other processes of central importance in hydrology.

Hydrologists have only recently accepted the challenge of developing their own re-analysis-type products with perhaps the first serious hydrologic re-analysis being published as recently as a few years ago (Rodell et al., 2015). More recently, the Princeton University group has extended this early work by making available a gridded global terrestrial hydrologic re-analysis product known as the Climate Data Record (CDR) (Zhang et al., 2018). Briefly, the CDR was constructed by synthesizing multiple in situ observations, satellite remote sensing products and land surface model outputs to provide gridded estimates of global land precipitation P, evapotranspiration E, runoff Q and total water storage change ΔS ($\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$, monthly, 1984–2010). In developing the CDR, the authors adopted local water budget closure as the fundamental hydrologic principle. That approach presented one important difficulty. Global observations of ΔS start with the GRACE satellite mission from 2002. Hence before 2002 there is no direct observational constraint on ΔS and the authors made the further assumption that the mean annual ΔS over the full 1984–2010 period was zero at every grid box. That is incorrect in some regions (e.g. Scanlon et al., 2018) and represents an observational problem that cannot be overcome. However, our interest is in the year-to-year variability and for that application, the assumption of no change in the mean annual ΔS over the full 1984–2010 period is unlikely to lead to major problems since we are not looking for subtle changes over time. With that caveat in mind, the aim of this study is to use this new 27-year gridded hydrologic re-analysis product to conduct an initial investigation of the inter-annual variability of the terrestrial branch of the global water cycle.

The paper is structured as follows. We begin in Sect. 2 by describing the various climate and hydrologic databases used in this study and also include a further assessment of the suitability of the CDR database for this initial variability study. In Sect. 3, we examine relationships between the mean and variability in the four water cycle variables (P, E, Q and ΔS). In Sect. 4, we first relate the variabilities to the classical aridity index and then use those results to evaluate the theory of Koster and Suarez (1999). Subsequently we examine how the variance of P is partitioned into the variances (and relevant covariances) of E, Q and ΔS and undertake an initial survey that investigates some of the factors controlling the variance partitioning. We conclude the paper with a discussion summarizing what we have learnt about water cycle variability over land by using the CDR database.

2 Methods and data

2.1 Methods

The water balance is defined by

$$\begin{array}{}\text{(1)}& P\left(t\right)=E\left(t\right)+Q\left(t\right)+\mathrm{\Delta }S\left(t\right),\end{array}$$

with P the precipitation, E the evapotranspiration, Q the runoff and ΔS the total water storage change in time step t (annual in this study). By the usual variance law, we have

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\mathit{\sigma }}_P^{\mathrm{2}}& ={\mathit{\sigma }}_E^{\mathrm{2}}+{\mathit{\sigma }}_Q^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}+\mathrm{2}\mathrm{cov}\left(E,Q\right)+\mathrm{2}\mathrm{cov}\left(E,\mathrm{\Delta }S\right)\\ & +\mathrm{2}\mathrm{cov}(Q,\mathrm{\Delta }S),\end{array}\end{array}$$

which includes all relevant variances (denoted σ²) and covariances (denoted cov). Equation (2) can be thought of as the hydrologic variance balance equation.

2.2 Hydrologic and climatic data

We use the CDR database (Zhang et al., 2018), which is a recently released global land hydrologic re-analysis. This product includes global precipitation P, evapotranspiration E, runoff Q and water storage change ΔS ($\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$, monthly, 1984–2010). In this study we focus on the inter-annual variability and the monthly water cycle variables (P, E, Q and ΔS) are aggregated to annual totals. The CDR does not report additional radiation variables and we use the NASA/GEWEX Surface Radiation Budget (SRB) Release-3.0 (monthly, 1984–2007, $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$) database (Stackhouse et al., 2011) to calculate E_o (defined as the net radiation expressed as an equivalent depth of liquid water, Budyko, 1974). We then calculate the aridity index ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}$) using P from the CDR and E_o from the SRB databases (see Fig. S1a in the Supplement).

In general, we anticipate two important factors, i.e. the water storage capacity and the presence of ice/snow at the surface, which are most likely to have influence on the partitioning of hydrologic variability. For the storage, the active range of the monthly water storage variation was used to approximate the water storage capacity (S_max). In more detail, the water storage S(t) at each time step t (monthly here) was first calculated from the accumulation of ΔS(t), i.e. $S\left(t\right)=S(t-\mathrm{1})+\mathrm{\Delta }S\left(t\right)$, where we assumed zero storage at the beginning of the study period (i.e. S(0)=0). With the resulting time series available, S_max was estimated as the difference between the maximum and minimum S(t) during the study period at each grid box (see Fig. S1b in the Supplement). The estimated S_max shows a large range from 0 to 1000 mm with the majority of values from 50 to 600 mm (Fig. S1b), which generally agrees with global rooting depth estimates assuming that water occupies from 10 % to 30 % of the soil volume at field capacity (Jackson et al., 1996; Wang-Erlandsson et al., 2016; Yang et al., 2016). To characterize snow/ice cover, and to distinguish extremely hot and cold regions, we also make use of a gridded global land air temperature dataset from the Climatic Research Unit (CRU TS4.01 database, monthly, 1901–2016, $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$) (Harris et al., 2014) (see Fig. S1c in the Supplement).

2.3 Spatial mask to define study extent

The CDR database provides an estimate of the uncertainty (±1σ) for each of the hydrologic variables (P, E, Q, ΔS) in each month. We use those uncertainty estimates to identify and remove regions with high relative uncertainty in the CDR data. The relative uncertainty is calculated as the ratio of root mean square of the uncertainty (±1σ) to the mean annual P, E and Q at each grid box following the procedure used by Milly and Dunne (2002). Note that the long-term mean ΔS is zero by construction in the CDR database, and for that reason we did not use ΔS to calculate the relative uncertainty. Grid boxes with a relative uncertainty (in P, E and Q) of more than 10 % are deemed to have high relative uncertainty (Milly and Dunne, 2002) and were excluded from the study extent. The excluded grid boxes were mostly in the Himalayan region, the Sahara Desert and in Greenland. The final spatial mask is shown in Fig. S2 and this has been applied throughout this study.

2.4 Further evaluation of CDR data for variability analysis

In the original work, the CDR database was validated by comparison with independent observations including (i) mean seasonal cycle of Q from 26 large basins (see Fig. 8 in Zhang et al., 2018), (ii) mean seasonal cycle of ΔS from 12 large basins (Fig. 10 in Zhang et al., 2018), (iii) monthly runoff from 165 medium size basins and a further 862 small basins (Fig. 14 in Zhang et al., 2018), and (iv) summer E from 47 flux towers (Fig. 16 in Zhang et al., 2018). Those evaluations did not directly address variability in various water cycle elements. With our focus on the variability we decided to conduct further validations of the CDR database beyond those described in the original work. In particular, we focussed on further independent assessments of E and we use monthly (as opposed to summer) observations of E from FLUXNET to evaluate the variability in E. We also compare the evapotranspiration E in the CDR with two other gridded global E products that were not used to develop the CDR including the LandFluxEval database ($\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$, monthly, 1989–2005) (Mueller et al., 2013) and the Max Planck Institute database (MPI, $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$, monthly, 1982–2011) (Jung et al., 2010). The runoff Q in the CDR is further compared with the gridded European Q product E-RUN ($\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$, monthly, 1951–2015) (Gudmundsson and Seneviratne, 2016a).

For the comparison to FLUXNET observations (Baldocchi et al., 2001; Agarwal et al., 2010) we identified 32 flux tower sites (site locations are shown in Fig. S3 and details are shown in Table S1) with at least 3 years of continuous (monthly) measurements using the FluxnetLSM R package (v1.0) (Ukkola et al., 2017). The monthly totals and annual climatology of P and E from CDR generally follow FLUXNET observations, with high correlations and reasonable root mean square error (Figs. S4–S5, Table S1). Comparison of the point-based FLUXNET (∼100 m–1 km scale) with the grid-based CDR (∼50 km scale) is problematic since the CDR represents an area that is at least 2500 times larger than the area represented by the individual FLUXNET towers and we anticipate that the CDR record would be “smoothed” relative to the FLUXNET record. With that in mind, we chose to compare the ratio of the standard deviation of E to P between the CDR and FLUXNET databases and this normalized comparison of the hydrologic variability proved encouraging (Fig. S6).

The comparison of E between the CDR and the LandFluxEval and MPI databases also proved encouraging. We found that the monthly mean E from the CDR database is slightly underestimated compared with the LandFluxEVAL database (Fig. S7a) but agrees closely with the MPI database (Fig. S8a). In terms of variability, the standard deviations of monthly E from the CDR are in very close agreement with the LandFluxEVAL database (Fig. S7c), but there is a bias and scaling offset for the comparison with the MPI database, particularly for the grid cells with low standard deviation of E (Fig. S8c). The comparison of runoff Q between the E-RUN and CDR databases shows that the two databases have very similar spatial patterns of both the long-term mean ($\stackrel{\mathrm{‾}}{Q}$) and standard deviation (σ_Q) of the monthly Q (Fig. S10). The grid-by-grid comparison results are also encouraging, showing slight bias of both the long-term mean and standard deviation of monthly Q in the CDR database compared with the E-RUN database (Fig. S11).

We concluded that while the CDR database was unlikely to be perfect, it was nevertheless suitable for an initial exploratory survey of inter-annual variability in the terrestrial branch of the global water cycle.

3 Mean and variability of water cycle components

3.1 Mean annual P, E and Q and the Budyko curve

The global pattern of mean annual P, E and Q using the CDR data (1984–2007) is shown in Fig. 1. The mean annual P ($\stackrel{\mathrm{‾}}{P}$) is prominent in tropical regions, southern China, and eastern and western North America (Fig. 1a). The magnitude of mean annual E ($\stackrel{\mathrm{‾}}{E}$) more or less follows the pattern of $\stackrel{\mathrm{‾}}{P}$ in the tropics (Fig. 1b) while the mean annual Q ($\stackrel{\mathrm{‾}}{Q}$) is particularly prominent in the Amazon, South and Southeast Asia, tropical parts of western Africa, and in some other coastal regions at higher latitudes (Fig. 1c).

Figure 1 Mean annual (1984–2010) (a) P, (b) E and (c) Q. Note that the mean annual ΔS in the CDR database is zero by construction and is not shown.

Figure 2 Relationship of mean annual (a) evapotranspiration ($\stackrel{\mathrm{‾}}{E}/\stackrel{\mathrm{‾}}{P}$) and (b) runoff ($\stackrel{\mathrm{‾}}{Q}/\stackrel{\mathrm{‾}}{P}$) ratios to the aridity index ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}$) from the CDR and SRB databases. For comparison, the Budyko (1974) curve is shown on panel (a). The curve on panel (b) is calculated assuming a steady state ($\stackrel{\mathrm{‾}}{Q}/\stackrel{\mathrm{‾}}{P}=\mathrm{1}-\stackrel{\mathrm{‾}}{E}/\stackrel{\mathrm{‾}}{P}$).

We relate the grid-box level ratio of $\stackrel{\mathrm{‾}}{E}$ to $\stackrel{\mathrm{‾}}{P}$ in the CDR database to the classical Budyko (1974) curve using the aridity index ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}$) (Fig. 2a). As noted previously, in the CDR database, $\stackrel{\mathrm{‾}}{\mathrm{\Delta }S}$ is forced to be zero and this enforced steady state (i.e. $\stackrel{\mathrm{‾}}{P}=\stackrel{\mathrm{‾}}{E}+\stackrel{\mathrm{‾}}{Q}$) allowed us to also predict the ratio of $\stackrel{\mathrm{‾}}{Q}$ to $\stackrel{\mathrm{‾}}{P}$ using the same Budyko curve (Fig. 2b). The Budyko curves follow the overall pattern in the CDR data, which agrees with previous studies showing that the aridity index can be used to predict water availability (e.g. Gudmundsson et al., 2016). However, there is substantial scatter due to, for example, regional variations related to seasonality, water storage change and the landscape characteristics (Milly, 1994a, b; Padrón et al., 2017). With that caveat in mind, the overall patterns are as expected, with $\stackrel{\mathrm{‾}}{E}$ following $\stackrel{\mathrm{‾}}{P}$ in dry environments ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}>\mathrm{1.0}$) while $\stackrel{\mathrm{‾}}{E}$ follows $\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}$ in wet environments ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}\le \mathrm{1.0}$) (Fig. 2).

3.2 Inter-annual variability in P, E, Q and ΔS

We use the variance balance equation (Eq. 2) to partition the inter-annual ${\mathit{\sigma }}_P^{\mathrm{2}}$ into separate components due to ${\mathit{\sigma }}_E^{\mathrm{2}}$, ${\mathit{\sigma }}_Q^{\mathrm{2}}$, ${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}$ along with the three covariance components [2cov(EQ), 2cov(E,ΔS), 2cov(Q,ΔS)] (Fig. 3). The spatial pattern of ${\mathit{\sigma }}_P^{\mathrm{2}}$ (Fig. 3a) is very similar to that of $\stackrel{\mathrm{‾}}{P}$ (Fig. 1a), which implies that the ${\mathit{\sigma }}_P^{\mathrm{2}}$ is positively correlated with $\stackrel{\mathrm{‾}}{P}$. In contrast the partitioning of ${\mathit{\sigma }}_P^{\mathrm{2}}$ to the various components is very different from the partitioning of $\stackrel{\mathrm{‾}}{P}$ (cf. Figs. 1 and 3). First we note that while the overall spatial pattern of ${\mathit{\sigma }}_E^{\mathrm{2}}$ more or less follows ${\mathit{\sigma }}_P^{\mathrm{2}}$, the overall magnitude of ${\mathit{\sigma }}_E^{\mathrm{2}}$ is much smaller than ${\mathit{\sigma }}_P^{\mathrm{2}}$ and ${\mathit{\sigma }}_Q^{\mathrm{2}}$ in most regions, and in fact ${\mathit{\sigma }}_E^{\mathrm{2}}$ is also generally smaller than ${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}$. The prominence of ${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}$ (compared to ${\mathit{\sigma }}_E^{\mathrm{2}}$) surprised us. The three covariance components [cov(E,Q), cov(E,ΔS), cov(Q,ΔS)] are also important in some regions. In more detail, the cov(E,Q) term is prominent in regions where ${\mathit{\sigma }}_Q^{\mathrm{2}}$ is large and is mostly negative in those regions (Fig. 3e), indicating that years with lower E are associated with higher Q and vice versa. There are also a few regions with prominent positive values for cov(EQ) (e.g. the seasonal hydroclimates of northern Australia) indicating that in those regions, years with a higher E are associated with higher Q. The cov(E,ΔS) term (Fig. 3f) has a similar spatial pattern to the cov(E,Q) term (Fig. 3e) but with a smaller overall magnitude. Finally, the cov(Q,ΔS) term shows a more complex spatial pattern, with both prominent positive and negative values (Fig. 3g) in regions where ${\mathit{\sigma }}_Q^{\mathrm{2}}$ (Fig. 3c) and ${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}$ (Fig. 3d) are both large.

Figure 3 Water cycle variances (${\mathit{\sigma }}_P^{\mathrm{2}}$, ${\mathit{\sigma }}_E^{\mathrm{2}}$, ${\mathit{\sigma }}_Q^{\mathrm{2}}$, ${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}$) and covariances [cov(EQ), cov(E,ΔS), cov(Q,ΔS)]. Note that we have multiplied the covariances by two (see Eq. 2).

These results show that the spatial patterns in variability are not simply a reflection of patterns in the long-term mean state. On the contrary, we find that of the three primary variance terms, the overall magnitude of (inter-annual) ${\mathit{\sigma }}_E^{\mathrm{2}}$ is the smallest, implying the least (inter-annual) variability in E. This is very different from the conclusions based on spatial patterns in the mean P, E and Q (see Sect. 3.1). Further, while ${\mathit{\sigma }}_Q^{\mathrm{2}}$ more or less follows ${\mathit{\sigma }}_P^{\mathrm{2}}$ as expected, we were surprised by the magnitude of ${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}$ which, in general, substantially exceeds the magnitude of ${\mathit{\sigma }}_E^{\mathrm{2}}$. Further, the magnitude of the covariance terms can be important, especially in regions with high ${\mathit{\sigma }}_Q^{\mathrm{2}}$. However, unlike the variances, the covariance can be both positive and negative and this introduces additional complexity. For example, with a negative covariance it is possible for the variance in Q (${\mathit{\sigma }}_Q^{\mathrm{2}}$) to exceed the variance in P (${\mathit{\sigma }}_P^{\mathrm{2}}$). To examine that in more detail we calculated the equivalent frequency distribution for each of the plots in Fig. 3. The results (Fig. S9) further emphasize that in general, ${\mathit{\sigma }}_E^{\mathrm{2}}$ is the smallest of the variances (Fig. S9b). We also note that the frequency distributions for the covariances (Fig. S9e, f, g) are not symmetrical. In summary, it is clear that spatial patterns in the inter-annual variability of the water cycle (Fig. 3) do not simply follow the spatial patterns for the inter-annual mean (Fig. 1).

3.3 Relation between variability and the mean state for P, E, Q

Differences in the spatial patterns of the mean (Fig. 1) and inter-annual variability (Fig. 3) in the global water cycle led us to further investigate the relation between the mean and the variability for each separate component. Here we relate the standard deviation (σ_P, σ_E, σ_Q) instead of the variance to the mean of each water balance flux (Fig. 4) since the standard deviation has the same physical units as the mean, making the results more comparable. As inferred previously, we find σ_P to be positively correlated with $\stackrel{\mathrm{‾}}{P}$ but with substantial scatter (Fig. 4a). The same result more or less holds for the relation between σ_Q and $\stackrel{\mathrm{‾}}{Q}$ (Fig. 4c). In contrast the relation between σ_E and $\stackrel{\mathrm{‾}}{E}$ is very different (Fig. 4b). In particular, σ_E is a small fraction of $\stackrel{\mathrm{‾}}{E}$ and this complements the earlier finding (Fig. 4b) that the inter-annual variability for E is generally smaller than for the other physical variables (P, Q and ΔS). (The same result was also found using both LandFluxEVAL and MPI databases; see Fig. S12 in the Supplement.) Importantly, unlike P and Q, E is constrained by both water and energy availability (Budyko, 1974) and the limited inter-annual variability in E presumably reflects limited inter-annual variability in the available (radiant) energy (E_o). This is something that could be investigated in a future study.

Figure 4 Relation between inter-annual mean and standard deviation for (a) P, (b) E and (c) Q from the CDR database. Note that the mean annual ΔS is zero by construction and is not shown.

4 Relating the variability of water cycle components to aridity

In the previous section, we investigated spatial patterns of the mean and the variability in the global water cycle. In this section, we extend that by investigating the partitioning of ${\mathit{\sigma }}_P^{\mathrm{2}}$ to the three primary physical terms (${\mathit{\sigma }}_E^{\mathrm{2}}$, ${\mathit{\sigma }}_Q^{\mathrm{2}}$, ${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}$) along with the three relevant covariances. For that, we begin by comparing the Koster and Suarez (1999) theory against the CDR data and then investigate how the partitioning of the variance is related to the aridity index $\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}$ (see Fig. S1a in the Supplement). Following that, we investigate variance partitioning in relation to both our estimate of the storage capacity S_max (see Fig. S1b in the Supplement) as well as the mean annual air temperature $\stackrel{\mathrm{‾}}{T_{\mathrm{a}}}$ (see Fig. S1c in the Supplement) that we use as a surrogate for snow/ice cover. We finish this section by examining the partitioning of variance at three selected study sites that represent extremely dry/wet, high/low water storage capacity and the hot/cold spectrums.

4.1 Comparison with the Koster and Suarez (1999) theory

We first evaluate the classical empirical curve of Koster and Suarez (1999) by relating ratios σ_E∕σ_P and σ_Q∕σ_P to the aridity index (Fig. 5). The ratio σ_E∕σ_P in the CDR database is generally overestimated by the empirical Koster and Suarez curve, especially in dry environments (e.g. $\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}>\mathrm{3}$) (Fig. 5a). The inference here is that the Koster and Suarez theory predicts σ_E∕σ_P to approach unity in dry environments while the equivalent value in the CDR data is occasionally unity but is generally smaller. With σ_E∕σ_P generally overestimated by the Koster and Suarez theory we expect, and find, that σ_Q∕σ_P is generally underestimated by the same theory (Fig. 5b). The same overestimation was found based on the other two independent databases for E (LandFluxEVAL and MPI) (Fig. S13). This overestimation is discussed further in Sect. 5.

Figure 5 Relationship of inter-annual standard deviation of (a) evapotranspiration (σ_E∕σ_P) and (b) runoff (σ_Q∕σ_P) ratios to aridity ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}$). The curves represent the semi-empirical relations from Koster and Suarez (1999).

4.2 Relating inter-annual variability to aridity

Here we examine how the fraction of the total variance in precipitation accounted for by the three primary variance terms along with the three covariance terms varies with the aridity index ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}$) (Fig. 6). (Also see Fig. S14 for the spatial maps.) The ratio ${\mathit{\sigma }}_E^{\mathrm{2}}/{\mathit{\sigma }}_P^{\mathrm{2}}$ is close to zero in extremely wet regions and has an upper limit noted previously (Fig. 5a) that approaches unity in extremely dry regions (Fig. 6a). The ratio ${\mathit{\sigma }}_Q^{\mathrm{2}}/{\mathit{\sigma }}_P^{\mathrm{2}}$ is close to zero in extremely dry regions but approaches unity in extremely wet regions but with substantial scatter (Fig. 6b). The ratio ${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}/{\mathit{\sigma }}_P^{\mathrm{2}}$ is close to zero in both extremely dry/wet regions (Fig. 6c) and shows the largest range at an intermediate aridity index ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}\sim \mathrm{1.0}$).

Figure 6 Relation between water cycle variances and covariances (see Fig. 3b–g) as a fraction of the variance of P (${\mathit{\sigma }}_P^{\mathrm{2}}$) and the aridity index ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}$) coloured by density. Note that we have multiplied the covariance components by two (see Eq. 2).

The covariance ratios are all small in extremely dry (e.g. $\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}\ge \mathrm{6.0}$) environments and generally show the largest range in semi-arid and semi-humid environments. The peak magnitudes for the three covariance components consistently occur when $\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}$ is close to 1.0, which is the threshold often used to separate wet and dry environments.

4.3 Further investigations on the factors controlling partitioning of the variance

Results in the previous section demonstrated that spatial variation in the partitioning of ${\mathit{\sigma }}_P^{\mathrm{2}}$ into ${\mathit{\sigma }}_E^{\mathrm{2}}$, ${\mathit{\sigma }}_Q^{\mathrm{2}}$, ${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}$ and the three covariance components is complex (Fig. 6). To help further understand inter-annual variability of the terrestrial water cycle, we conduct further investigations in this section using two factors likely to have a major influence on the variance partitioning of ${\mathit{\sigma }}_P^{\mathrm{2}}$. The first is the storage capacity S_max (see Fig. S1b in the Supplement). The second is the mean annual air temperature $\stackrel{\mathrm{‾}}{T_{\mathrm{a}}}$ (see Fig. S1c in the Supplement), which is used here as a surrogate for snow/ice presence.

4.3.1 Relating inter-annual variability to storage capacity

We first relate the partitioning of ${\mathit{\sigma }}_P^{\mathrm{2}}$ to water storage capacity (S_max) by repeating Fig. 6 but instead we use a logarithmic scale for the x axis and we distinguish S_max via the background colour (Fig. 7). To eliminate the possible overlap of grid cells in the colouring process, all the grid cells over land are further separated using different latitude ranges (as shown in the four columns of Fig. 7), i.e. 90–60^∘ N, 60–30^∘ N, 30–0^∘ N and 0–90^∘ S. We find that S_max is relatively high in wet environments ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}\le \mathrm{1.0}$, Fig. 7a) but shows no obvious relation to the partitioning of ${\mathit{\sigma }}_P^{\mathrm{2}}$. However, in dry environments ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}>\mathrm{1.0}$) the ratio ${\mathit{\sigma }}_E^{\mathrm{2}}/{\mathit{\sigma }}_P^{\mathrm{2}}$ apparently decreases with the increase in S_max (Fig. 7a–d). That relation is particularly obvious in extremely dry environments ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}\ge \mathrm{6.0}$) at equatorial latitudes where there is an upper limit of ${\mathit{\sigma }}_E^{\mathrm{2}}/{\mathit{\sigma }}_P^{\mathrm{2}}$ close to 1.0 when S_max is small (blue grid cells in Fig. 7c). The interpretation for those extremely dry environments is that when S_max is small, ${\mathit{\sigma }}_P^{\mathrm{2}}$ is almost completely partitioned into ${\mathit{\sigma }}_E^{\mathrm{2}}$ (Fig. 7b, c) with the other variance and covariance components close to zero. While for those same extremely dry environments, as S_max increases, the partitioning of ${\mathit{\sigma }}_P^{\mathrm{2}}$ is shared between ${\mathit{\sigma }}_E^{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}$ and their covariance (Fig. 7c, k, s), while ${\mathit{\sigma }}_Q^{\mathrm{2}}$ and its covariance components remain close to zero (Fig. 7g, o, w). However, at polar latitudes in the Northern Hemisphere (panels in the first and second columns of Fig. 7) there are variations that could not be easily associated with variations in S_max, which led us to further investigate the role of snow/ice on the variance partitioning in the following section.

Figure 7 Relation between water cycle variances and covariances (see Fig. 3b–g) as a fraction of the variance for P (${\mathit{\sigma }}_P^{\mathrm{2}}$) and the aridity index ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}$) for grid cells over different latitude ranges (i.e. 90–60^∘ N, 60–30^∘ N, 30–0^∘ N and 0–90^∘ S). The colours relate to the water storage capacity S_max. Note that we have multiplied the covariances by two (see Eq. 2). The vertical grey dashed lines represent thresholds used to separate extremely dry ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}\ge \mathrm{6.0}$) and wet ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}\le \mathrm{0.5}$) environments. Note the use of a logarithmic x axis and scale bar for S_max.

4.3.2 Relating inter-annual variability to mean air temperature

To understand the potential role of snow/ice in modifying the variance partitioning, we repeat the previous analysis (Fig. 7) but here we use the mean annual air temperature ($\stackrel{\mathrm{‾}}{T_{\mathrm{a}}}$) to colour the grid cells to (crudely) indicate the presence of snow/ice (Fig. 8). The results are complex and not easy to simply understand. The most important difference revealed by this analysis is in the hydrologic partitioning between cold (first column) and hot (third column) conditions in wet environments ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}\le \mathrm{0.5}$). In particular, when $\stackrel{\mathrm{‾}}{T_{\mathrm{a}}}$ is high, ${\mathit{\sigma }}_P^{\mathrm{2}}$ is almost completely partitioned into ${\mathit{\sigma }}_Q^{\mathrm{2}}$ in wet environments (e.g. $\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}\le \mathrm{0.5}$, Fig. 8g). In contrast, when $\stackrel{\mathrm{‾}}{T_{\mathrm{a}}}$ is low in a wet environment ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}\le \mathrm{0.5}$ in first column of Fig. 8), there are substantial variations in the hydrologic partitioning. That result reinforces the complexity of variance partitioning in the presence of snow/ice.

Figure 8 Relation between water cycle variances and covariances (see Fig. 3b–g) as a fraction of the variance for P (${\mathit{\sigma }}_P^{\mathrm{2}}$) and the aridity index ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}$) for grid cells over different latitude ranges (i.e. 90–60^∘ N, 60–30^∘ N, 30–0^∘ N and 0–90^∘ S). The colours relate to the mean air temperature ($\stackrel{\mathrm{‾}}{T_{\mathrm{a}}}$). Note that we have multiplied the covariances by two (see Eq. 2). The vertical grey dashed lines represent thresholds used to separate extremely dry ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}\ge \mathrm{6.0}$) and wet ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}\le \mathrm{0.5}$) environments.

4.4 Case studies

The previous results (Sect. 4.3) have demonstrated that the partitioning of ${\mathit{\sigma }}_P^{\mathrm{2}}$ is influenced by the water storage capacity (S_max) in extremely dry environments ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}\ge \mathrm{6.0}$) and that the presence of snow/ice is important (as indicated by mean air temperature, $\stackrel{\mathrm{‾}}{T_{\mathrm{a}}}$) in extremely wet environments ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}\le \mathrm{0.5}$). In this section, we examine, in greater detail, several sites to gain deeper understanding of the partitioning of ${\mathit{\sigma }}_P^{\mathrm{2}}$. For that purpose, we selected three sites based on extreme values for the three explanatory parameters, i.e. $\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}$ (Fig. S1a), S_max (Fig. S1b) and $\stackrel{\mathrm{‾}}{T_{\mathrm{a}}}$ (Fig. S1c). The criteria to select three climate sites are as follows – Site 1: dry ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}\ge \mathrm{6.0}$) and small S_max ($S_{max}\approx \mathrm{0}$), Site 2: dry ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}\ge \mathrm{6.0}$) and relatively large S_max ($S_{max}\gg \mathrm{0}$) and Site 3: wet ($\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}\le \mathrm{0.5}$) and hot ($\stackrel{\mathrm{‾}}{T_{\mathrm{a}}}>\mathrm{25}$ ^∘C). For each of the three classes, we use a representative grid cell (Fig. 9) to show the original time series (Fig. 10) and the partitioning of the variability (Fig. 11).

Figure 9 Locations of three representative grid cells used as case study sites.

Figure 10 Inter-annual time series (P, E, Q and ΔS) and the associated variance-covariance matrix (E, Q and ΔS) for case study Sites 1–3. Left column shows time series for (a) Site 1, (c) Site 2 and (e) Site 3, with right column, i.e. panels (b), (d) and (f), the associated variance–covariance matrix for three sites. Note that the covariance values in the tables should be multiplied by two to agree with the variance–covariance balance in Eq. (2).

Figure 11 Location of three case study sites in the water cycle variability space. The grey background dots are from Fig. 6.

We show the P, E, Q and ΔS time series along with the relevant variances and covariances in Fig. 10. Starting with the two dry sites, at the site with low storage capacity (Site 1), the time series shows that E closely follows P, leaving annual Q and ΔS close to zero (Fig. 10a). The variance of P (${\mathit{\sigma }}_P^{\mathrm{2}}=\mathrm{206.9}$ mm²) is small and almost completely partitioned into the variance of E (${\mathit{\sigma }}_E^{\mathrm{2}}=\mathrm{196.9}$ mm²), leaving very limited variance for Q, ΔS and all three covariance components (Fig. 10b). At the dry site with larger storage capacity (Site 2), E, Q and ΔS do not simply follow P (Fig. 10c). As a consequence, the variance of P (${\mathit{\sigma }}_P^{\mathrm{2}}=\mathrm{2798.0}$ mm²) is shared between E (${\mathit{\sigma }}_E^{\mathrm{2}}=\mathrm{1150.2}$ mm²), ΔS (${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}=\mathrm{800.5}$ mm²) and their covariance component ($\mathrm{2}\mathrm{cov}(E,\mathrm{\Delta }S)=\mathrm{538.4}$ mm², Fig. 10d). Switching now to the remaining wet and hot site (Site 3), we note that Q closely follows P, with ΔS close to zero and E showing little inter-annual variation (Fig. 10e). The variance of P (${\mathit{\sigma }}_P^{\mathrm{2}}=\mathrm{57}\mathrm{374.4}$ mm²) is relatively large and almost completely partitioned into the variance of Q (${\mathit{\sigma }}_Q^{\mathrm{2}}=\mathrm{57}\mathrm{296.4}$ mm²), leaving very limited variance for E and ΔS and the three covariance components (Fig. 10f). We also examined numerous other sites with similar extreme conditions as the three case study sites and found the same basic patterns as reported above.

To put the data from the three case study sites into a broader variability context we position the site data onto a backdrop of Fig. 6. As noted previously, at Site 1, the ratio ${\mathit{\sigma }}_E^{\mathrm{2}}/{\mathit{\sigma }}_P^{\mathrm{2}}$ is very close to unity (Fig. 11a), and under this extreme condition, we have the following approximation:

$$\begin{array}{}\text{(3)}& {\mathit{\sigma }}_P^{\mathrm{2}}\approx {\mathit{\sigma }}_E^{\mathrm{2}}(\mathrm{Site}\mathrm{1},\mathrm{dry}\mathrm{and}{S}_{max}\approx \mathrm{0})\end{array}$$

In contrast, for Site 2 with the same aridity index but higher S_max, we have

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\mathit{\sigma }}_P^{\mathrm{2}}\approx {\mathit{\sigma }}_E^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}+\mathrm{2}\mathrm{cov}\left(E,\mathrm{\Delta }S\right)\\ & (\mathrm{Site}\mathrm{2},\mathrm{dry}\mathrm{and}{S}_{max}\gg \mathrm{0}).\end{array}\end{array}$$

Finally, at Site 3, we have

$$\begin{array}{}\text{(5)}& {\mathit{\sigma }}_P^{\mathrm{2}}\approx {\mathit{\sigma }}_Q^{\mathrm{2}}(\mathrm{Site}\mathrm{3},\mathrm{wet}\mathrm{and}\mathrm{hot}).\end{array}$$

4.5 Synthesis

The above simple examples demonstrate that aridity $\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}$, storage capacity S_max and, to a lesser extent, air temperature $\stackrel{\mathrm{‾}}{T_{\mathrm{a}}}$, all play some role in the partitioning of ${\mathit{\sigma }}_P^{\mathrm{2}}$ to the various components. Our synthesis of the results for the partitioning of ${\mathit{\sigma }}_P^{\mathrm{2}}$ is summarized in Fig. 12. In dry environments with low storage capacity ($S_{max}\approx \mathrm{0}$) we have minimal runoff and expect that ${\mathit{\sigma }}_P^{\mathrm{2}}$ is more or less completely partitioned into ${\mathit{\sigma }}_E^{\mathrm{2}}$ (Fig. 12a). In those environments, (inter-annual) variations in storage ${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}$ play a limited role in setting the overall variability. However, in dry environments with larger storage capacity ($S_{max}\gg \mathrm{0}$), ${\mathit{\sigma }}_E^{\mathrm{2}}$ is only a small fraction of ${\mathit{\sigma }}_P^{\mathrm{2}}$ (Fig. 12a), leaving most of the overall variance in ${\mathit{\sigma }}_P^{\mathrm{2}}$ to be partitioned to ${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}$ and the covariance between E and ΔS (Fig. 12c and e). This emphasizes the hydrological importance of water storage capacity in buffering variations of the water cycle under dry conditions.

Figure 12 Synthesis of factors controlling variance partitioning. The arrows denote trends with increasing S_max. The grey background dots are from Fig. 6.

Under extremely wet conditions, the largest difference in variance partitioning is not due to differences in storage capacity but is instead related to differences in mean air temperature. In wet and hot environments, we have maximum runoff and find that ${\mathit{\sigma }}_P^{\mathrm{2}}$ is more or less completely partitioned into ${\mathit{\sigma }}_Q^{\mathrm{2}}$ (Fig. 12b) while the partitioning to ${\mathit{\sigma }}_E^{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}$ is small. However, in wet and cold environments, the variance partitioning shows great complexity, with ${\mathit{\sigma }}_P^{\mathrm{2}}$ being partitioned into all possible components. We suggest that this emphasizes the hydrological importance of thermal processes (melting/freezing) under extremely cold conditions.

However, the most complex patterns to interpret are those for semi-arid to semi-humid environments (i.e. $\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}$ ∼ 1.0). Despite a multitude of attempts over an extended period we were unable to develop a simple useful synthesis to summarize the partitioning of variability in those environments. We found that the three covariance terms all play important roles and we also found that simple environmental gradients (e.g. dry/wet, high/low storage capacity, hot/cold) could not easily explain the observed patterns. We anticipate that vegetation-related processes (e.g. phenology, rooting depth, gas exchange characteristics, disturbance) may prove to be important in explaining hydrologic variability in these biologically productive regions that support most of the human population. This result implies that a major scientific effort will be needed to develop a synthesis of the controlling factors for variability of the water cycle in these environments.

5 Discussion and conclusions

Importantly, hydrologists have long been interested in hydrologic variability, but without readily available databases it has been difficult to quantify water cycle variability. For example, we are not aware of maps showing global spatial patterns in variance for any terms of the water balance (except for P). In this study, we describe an initial investigation of the inter-annual variability of the terrestrial branch in the global water cycle that uses the recently released global monthly Climate Data Record (CDR) database for P, E, Q and ΔS. The CDR is one of the first dedicated hydrologic reanalysis databases and includes data for a 27-year period. Accordingly, we could only examine hydrologic variability over this relatively short period. Further, we expect future improvements and modifications as the hydrologic community seeks to further develop and refine these new reanalysis databases. With those caveats in mind, we started this analysis by first investigating the partitioning of P in the water cycle in terms of long-term mean and then extended that to the inter-annual variability using a theoretical variance balance equation (Eq. 2). Despite the initial nature of this investigation we have been able to establish some useful general principles.

The mean annual P is mostly partitioned into mean annual E and Q, as is well known, and the results using the CDR were generally consistent with the earlier Budyko framework (Fig. 2). Having established that, the first general finding is that the spatial pattern in the partitioning of inter-annual variability in the water cycle is not simply a reflection of the spatial pattern in the partitioning of the long-term mean. In particular, with the variance calculations, the annual anomalies are squared and hence the storage anomalies do not cancel out like they do when calculating the mean. With that in mind, we were surprised that the inter-annual variability of water storage change (${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}$) is typically larger than the inter-annual variability of evapotranspiration (${\mathit{\sigma }}_E^{\mathrm{2}}$) (cf. Fig. 3b and d). The consequence is that ${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}$ is more important than ${\mathit{\sigma }}_E^{\mathrm{2}}$ for understanding inter-annual variability of the global water cycle. A second important generalization is that unlike the variance components which are all positive, the three covariance components in the theory (Eq. 2) can be both positive and negative. We report results here showing both large positive and negative values for the three covariance terms (Fig. 3e, f, g). This was especially prevalent in biologically productive regions ($\mathrm{0.5}<\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}<\mathrm{1.5}$, Fig. 3e, g). When examining the mean state, we are accustomed to think that P sets a limit to E, Q and ΔS, as per the mass balance (Eq. 1). But the same thinking does not extend to the variance balance since the covariance terms on the right-hand side of Eq. (2) can be both large and negative, leading to circumstances where the variability in the sinks (${\mathit{\sigma }}_E^{\mathrm{2}}$, ${\mathit{\sigma }}_Q^{\mathrm{2}}$, ${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}$) could actually exceed variability in the source (${\mathit{\sigma }}_P^{\mathrm{2}}$). These general principles of variance partitioning in the water cycle above may vary at different timescales (e.g. monthly, daily), and we expect more details of the variability partitioning across various temporal scales to be investigated in future studies.

Our initial attempt to develop deeper understanding of variance partitioning was based on a series of case studies located in extreme environments (wet/dry vs. hot/cold vs. high/low water storage capacity). The results offered some further insights about hydrologic variability. For example, under extremely dry (water-limited) environments, with limited storage capacity (S_max), we found that E follows P and ${\mathit{\sigma }}_E^{\mathrm{2}}$ follows ${\mathit{\sigma }}_P^{\mathrm{2}}$, with ${\mathit{\sigma }}_Q^{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}$ both approaching zero. However, as S_max increases, the partitioning of ${\mathit{\sigma }}_P^{\mathrm{2}}$ progressively shifts to a balance between ${\mathit{\sigma }}_E^{\mathrm{2}}$, ${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}$ and cov(E,ΔS) (Figs. 10–12). This result explains the overestimation of σ_E∕σ_P by the empirical theory of Koster and Suarez (1999) which implicitly assumed no inter-annual change in storage. The Koster and Suarez empirical theory is perhaps better described as an upper limit that is based on minimal storage capacity, and that any increase in storage capacity would promote the partitioning of ${\mathit{\sigma }}_P^{\mathrm{2}}$ to ${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}$, particularly under dry conditions (Figs. 10–12).

In extremely wet/hot environments (i.e. no snow/ice presence) we found ${\mathit{\sigma }}_P^{\mathrm{2}}$ to be mostly partitioned to ${\mathit{\sigma }}_Q^{\mathrm{2}}$ (with both ${\mathit{\sigma }}_E^{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathrm{\Delta }S}^{\mathrm{2}}$ approaching zero, Fig. 10). In contrast, in extremely wet/cold environments, the partitioning of ${\mathit{\sigma }}_P^{\mathrm{2}}$ was highly (spatially) variable, presumably because of spatial variability in the all-important thermal processes (freeze/melt).

The most complex results were found in mesic biologically productive environments ($\mathrm{0.5}<\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}<\mathrm{1.5}$), where all three covariance terms (Eq. 2) were found to be relatively large and therefore they all played critical roles in the overall partitioning of variability (Fig. 6). As noted above, in many of these regions, the (absolute) magnitudes of the covariances were actually larger than the variances of the water balance components E, Q and ΔS (e.g. Fig. 3). That result demonstrates that deeper understanding of the process-level interactions that are embedded within each of the three covariance terms (e.g. the role of seasonal vegetation variation) will be needed to develop process-based understanding of variability in the water cycle in these biologically productive regions ($\mathrm{0.5}<\stackrel{\mathrm{‾}}{E_{\mathrm{o}}}/\stackrel{\mathrm{‾}}{P}<\mathrm{1.5}$).

The syntheses of the long-term mean water cycle originated in 1970s (Budyko, 1974), and it took several decades for those general principles to become widely adopted in the hydrologic community. The hydrologic data needed to understand hydrologic variability are only now becoming available. With those data we can begin to develop a process-based understanding of hydrologic variability that can be used for a variety of purposes; e.g. deeper understanding of hydro-climatic behaviour, hydrologic risk analysis, climate change assessments and hydrologic sensitivity studies are just a few applications that spring to mind. The initial results presented here show that a major intellectual effort will be needed to develop a general understanding of hydrologic variability.