Assessing Global Water Storage Variability from GRACE: Trends, Seasonal Cycle, Subseasonal Anomalies and Extremes

Surveys in Geophysics, Feb 2016

Throughout the past decade, the Gravity Recovery and Climate Experiment (GRACE) has given an unprecedented view on global variations in terrestrial water storage. While an increasing number of case studies have provided a rich overview on regional analyses, a global assessment on the dominant features of GRACE variability is still lacking. To address this, we survey key features of temporal variability in the GRACE record by decomposing gridded time series of monthly equivalent water height into linear trends, inter-annual, seasonal, and subseasonal (intra-annual) components. We provide an overview of the relative importance and spatial distribution of these components globally. A correlation analysis with precipitation and temperature reveals that both the inter-annual and subseasonal anomalies are tightly related to fluctuations in the atmospheric forcing. As a novelty, we show that for large regions of the world high-frequency anomalies in the monthly GRACE signal, which have been partly interpreted as noise, can be statistically reconstructed from daily precipitation once an adequate averaging filter is applied. This filter integrates the temporally decaying contribution of precipitation to the storage changes in any given month, including earlier precipitation. Finally, we also survey extreme dry anomalies in the GRACE record and relate them to documented drought events. This global assessment sets regional studies in a broader context and reveals phenomena that had not been documented so far.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://link.springer.com/content/pdf/10.1007%2Fs10712-016-9367-1.pdf

Assessing Global Water Storage Variability from GRACE: Trends, Seasonal Cycle, Subseasonal Anomalies and Extremes

Surv Geophys Assessing Global Water Storage Variability from GRACE: Trends, Seasonal Cycle, Subseasonal Anomalies and Extremes Vincent Humphrey 0 Sonia I. Seneviratne 0 Lukas Gudmundsson 0 0 Institute for Atmospheric and Climate Science, ETH Zurich , Universitaetstrasse 16, 8092 Zurich , Switzerland Throughout the past decade, the Gravity Recovery and Climate Experiment (GRACE) has given an unprecedented view on global variations in terrestrial water storage. While an increasing number of case studies have provided a rich overview on regional analyses, a global assessment on the dominant features of GRACE variability is still lacking. To address this, we survey key features of temporal variability in the GRACE record by decomposing gridded time series of monthly equivalent water height into linear trends, inter-annual, seasonal, and subseasonal (intra-annual) components. We provide an overview of the relative importance and spatial distribution of these components globally. A correlation analysis with precipitation and temperature reveals that both the inter-annual and subseasonal anomalies are tightly related to fluctuations in the atmospheric forcing. As a novelty, we show that for large regions of the world high-frequency anomalies in the monthly GRACE signal, which have been partly interpreted as noise, can be statistically reconstructed from daily precipitation once an adequate averaging filter is applied. This filter integrates the temporally decaying contribution of precipitation to the storage changes in any given month, including earlier precipitation. Finally, we also survey extreme dry anomalies in the GRACE record and relate them to documented drought events. This global assessment sets regional studies in a broader context and reveals phenomena that had not been documented so far. GRACE decomposition; Water storage Precipitation Temperature Drought Signal 1 Introduction Land water resources are essential for human society and are affected by climate variability and human water use (Jime´nez Cisneros et al. 2014). It is thus important to monitor changes in land water storage, as well as the underlying processes leading to their variations in space and time. The Gravity Recovery and Climate Experiment (GRACE), launched in 2002, constitutes an essential tool for such analyses, as was demonstrated in a wealth of studies (Tapley et al. 2004a; Wahr et al. 2004; Rodell et al. 2004; Andersen et al. 2005; Velicogna and Wahr 2006; Gu¨ ntner et al. 2007a; Ramillien et al. 2008; Zaitchik et al. 2008; Rodell et al. 2009; Chen et al. 2010a; Houborg et al. 2012; Sasgen et al. 2012; Gardner et al. 2013; Wouters et al. 2014; Reager et al. 2014; Famiglietti 2014; Chen et al. 2015; Wahr 2015). After more than a decade of observations, the GRACE mission has resulted in an unprecedented view on global water storage variability, with a great diversity in terms of temporal scales, ranging from long-term trends to short-lived deviations from the seasonal cycle. These different scales of temporal variability often constitute a common denominator between GRACE studies, either implicitly—as when the discussion focuses on specific aspects like the seasonal cycle, trends or extremes—or explicitly—as when water storage time series are decomposed into subseries. Since the earliest GRACE studies, it has been, for instance, very common to refer to the phasing and amplitude of the seasonal cycle when comparing GRACE terrestrial water storage with other datasets such as model simulations (e.g., Tapley et al. 2004b; Wahr et al. 2004). As the temporal coverage of the GRACE record extended, more comprehensive studies also identified secular trends and inter-annual anomalies by separating the GRACE signal into long-term trends, periodical components and residual noise (Ramillien et al. 2005; Schmidt et al. 2008b; Steffen et al. 2009). However, there is still no global overview on the relative magnitude and distribution of these features of temporal variability. In addition, while certain of these features (e.g., seasonal cycles and trends) are relatively well described, others (e.g., high-frequency residuals and extremes) have typically attracted much less attention so far and remain more difficult to interpret. From a global perspective, terrestrial water storage anomalies derived from GRACE are dominated by a seasonal signal in most parts of the world. Consequently, the earliest studies comparing GRACE with hydrological models have primarily focused on the seasonal component. Most often, the seasonal cycle in GRACE was shown to compare relatively well with model simulations, both with respect to the signal’s amplitude and phase (Wahr et al. 2004; Swenson and Milly 2006; Syed et al. 2008; Schmidt et al. 2008b; Do¨ ll et al. 2014a). Reviews (Ramillien et al. 2008; Gu¨ ntner 2008; Schmidt et al. 2008a) showed that seasonal disagreement between GRACE and model data was usually attributed to deficiencies in the modelling of water storage compartments and to errors in the precipitation forcing, but also to signal leakage and inaccuracies of the GRACE data. Multiple studies have shown that long-term variability in the GRACE record over land can be related to long-term trends in groundwater (Rodell et al. 2009; Voss et al. 2013; Do¨ ll et al. 2014b; Chen et al. 2015; Richey et al. 2015a, b) and surface water (Swenson and Wahr 2009; Singh et al. 2012), teleconnections (Phillips et al. 2012) and mass variations in the cryosphere (Sasgen et al. 2012; Velicogna and Wahr 2013). The hydrological signal extracted from GRACE can also be contaminated by glacial isostatic adjustment (Wu et al. 2010) and crustal deformations caused by major earthquakes (Han et al. 2006, 2011, 2013). While the seasonal cycle, long-term anomalies and secular trends are arguably well documented, fewer studies have focused on subseasonal variability and extreme events at a global scale. So far only case studies have shown that major drought and flood events can be observed in the GRACE record (e.g., Andersen et al. 2005; Seitz et al. 2008; Frappart et al. 2012; Long et al. 2013; Abelen et al. 2015). Only recently, the potential of GRACE for monitoring drought conditions (Houborg et al. 2012; Thomas et al. 2014) and predicting flood potential (Reager et al. 2014) was investigated globally. However, large challenges remain since month-to-month variability in GRACE is highly contaminated with outliers, measurement errors and uncertainties arising from data processing (Bonin et al. 2012). The overarching goal of this study is to provide a global and comprehensive survey of the dominant features of temporal variability in terrestrial water storage observed from GRACE. Our approach is to decompose the total signal at each grid point into (1) linear trends, (2) inter-annual variability, (3) seasonal cycle and (4) subseasonal variability. We first assess the contribution of each component to the total signal at the global scale (Sect. 4.1). In Sect. 4.2, the magnitude and significance of the linear trends are discussed in the context of previous regional studies. Subsequently, the decomposed subseries of terrestrial water storage are compared with decomposed precipitation and temperature fields. Starting with the inter-annual anomalies, regions of high correlation between GRACE and these atmospheric drivers are identified (Sect. 4.3). Section 4.4 provides global maps of the maximum and minimum seasonal water storage and identifies phase shifts with respect to the seasonal cycles of both precipitation and temperature (Sect. 4.4). In Sect. 4.5, we focus on the subseasonal residuals and show that a careful averaging of the daily atmospheric data to the monthly resolution reveals excellent correlations with the high-frequency component of the GRACE signal. Finally, we use the decomposition approach to identify and analyse drought events in the GRACE record (Sect. 4.6). 2 Data 2.1 GRACE Data Monthly grids of terrestrial water storage anomalies used in this study are based on the spherical harmonic coefficients (Release 05) provided by the Center for Space Research (CSR), the Jet Propulsion Laboratory (JPL) and the GeoForschungsZentrum Potsdam (GFZ) for the period April 2002–August 2015. For more information on the GRACE mission, the gravity recovery process and the derivation of water storage anomalies from the spherical harmonic coefficients, we refer the reader to the reviews from Wouters et al. (2014) or Wahr (2015) and the references therein. The gridded product used in this study is the GRACE Tellus dataset (available at http://grace.jpl.nasa.gov). This dataset provides mass grids in units of equivalent water height for the three different sets of harmonic coefficients mentioned above, at a temporal resolution of approximately 1 month and with a grid resolution of 1 . It is worth noting that although 1 (or even finer) grids are commonly used in global analyses of terrestrial water storage anomalies, this does not reflect the actual spatial resolution of the GRACE measurements. Due to the truncation of spherical harmonics, the effective spatial resolution is by construction limited to a few hundreds of kilometres (Landerer and Swenson 2012). Additionally, postprocessing filters that are used to reduce spatially correlated errors further degrade the spatial resolution of the GRACE signal (Swenson and Wahr 2006; Duan et al. 2009; Longuevergne et al. 2010; Frappart et al. 2011b; Wouters et al. 2014). This causes spatial autocorrelation in the gridded dataset, as can be seen in Fig. 1, which also provides a general overview of the regions where hydrological variability, as detected by GRACE, has the largest magnitude. At the time of writing, the GRACE Tellus product is obtained through the following processing: the degree one harmonic coefficients (Earth’s geocenter) are estimated from Swenson et al. (2008), the coefficients of degree-order 2–0 (related to Earth oblateness) are replaced with more reliable solutions from Satellite Laser Ranging (Cheng et al. 2011) and correction for glacial isostatic adjustment is applied following Geruo et al. (2013). A known issue is that GRACE maps are heavily contaminated with correlated noise; hence, several spatial filtering techniques have been proposed that aim at restoring the geophysical signal (Kusche 2007; Ramillien et al. 2008; Werth et al. 2009; Frappart and Ramillien 2012). In the Tellus product, the destriping filter of Swenson and Wahr (2006) is applied to correct for North–South oriented stripes in GRACE maps and a 300 km Gaussian filter is additionally applied to the data to reduce residual noise. Finally, it is worth mentioning that GRACE time series are not evenly spaced in time. GRACE ‘‘months’’ most often do not correspond to calendar months due to instrument issues and solutions for several months can be missing, in particular after 2011. Instead, GRACE months represent approximately 1 month long periods with varying numbers of days. Sources of errors in GRACE include measurement errors, aliasing errors originating from the inaccurate correction of atmospheric and oceanic mass redistribution, and spatial leakage (Seo et al. 2006). Spatial leakage is caused both by the truncation of spherical harmonics and the postprocessing filters applied to the data (Chen et al. 2007a; Landerer and Swenson 2012). Since there are no other large-scale observations of terrestrial water storage that could be used as ground truth, estimating errors and confidence intervals for GRACE data is a major challenge (Gu¨ ntner 2008). One possibility to reduce uncertainty in the GRACE data is to use the ensemble mean of the solutions obtained by different processing centres (Werth et al. 2009; Sakumura et al. 2014). In this study, we use the mean of the three solutions from CSR, JPL and GFZ provided in the GRACE Tellus dataset. In order to correct for the amplitude attenuation caused by the postprocessing filters applied to the GRACE data, the Tellus dataset also provides the scaling factors proposed April 2002 - August 2015 Equivalent Water Height Standard Deviation [mm] 0 50 100 150 200 250 300 by Landerer and Swenson (2012). These scaling factors are derived by first applying the complete GRACE processing to modelled estimates of terrestrial water storage and subsequently comparing the agreement between the original and processed model data. A disadvantage of these scaling factors is that they can depend on the hydrological model used as a reference, especially in semi-arid and arid regions as well as over irrigated areas (Long et al. 2015). Long et al. (2015) also mention that scaling factors found in some regions should be interpreted carefully. For these reasons, scaling factors were not applied to the GRACE data in this study. 2.2 Filtered Grids of Atmospheric Reanalysis The atmospheric reanalysis ERA-Interim, from the European Centre for Medium-Range Weather Forecasts (ECMWF), is used to derive daily fields of mean 2 m air temperature and precipitation totals (Dee et al. 2011; available at http://apps.ecmwf.int/datasets/data/ interim-full-daily/). This dataset is obtained at a 0.25 resolution and averaged to the 1 resolution of the GRACE Tellus dataset. However, the effective spatial resolution of the hydrological signal observed in GRACE is still coarser than 1 , due to the resolution of the GRACE measurements (see Sect. 2.1). For the Tellus product, this effective spatial resolution is approximately 300 km (3 at the equator). In order to make the atmospheric data comparable with GRACE, we apply a 300 km Gaussian filter to the atmospheric grids. Without this filter, the atmospheric fields would show much more detailed patterns than the GRACE data. It is important to note that when GRACE solutions are compared with modelled estimates of terrestrial water storage, a common practice is to apply the whole GRACE processing to the model data, including an expansion of the modelled mass distribution into spherical harmonics and the subsequent postprocessing (e.g., Wahr et al. 2004; Schmidt et al. 2006; Swenson and Milly 2006; Syed et al. 2008). However, this latter approach cannot be applied to global fields of temperature and precipitation, which is why we only apply a Gaussian filter. We also note that the correlations between GRACE and filtered atmospheric fields are expected to increase as a consequence of this filtering. This effect has already been documented in a similar setting by Abelen and Seitz (2013) when comparing GRACE results with both modelled and remotely sensed soil moisture. 3 Methods 3.1 Signal Decomposition 3.1.1 Background and Previous Approaches Decomposition of the GRACE hydrological signal is common practice in the recent literature, and different methods have been used to address different objectives. One possibility is to aim at isolating the contribution of specific water storage compartments such as groundwater, soil moisture or snow mass to the total GRACE signal. This leads to highly underdetermined inversion problems of blind signal separation and gives rise to non-unique solutions as the contributing geophysical signals are most often not statistically independent. To account for this, inversion methods have been proposed that can use higher-order statistical information derived from model data to decompose the total signal (Ramillien et al. 2004, 2005; Frappart et al. 2006; Schmeer et al. 2012). Assimilation of GRACE data into a land surface scheme could also be seen as another approach relating GRACE variability to water storage compartments that are already partitioned in a model structure (Zaitchik et al. 2008; Eicker et al. 2014). In groundwater studies, a common strategy is to directly subtract model estimates of snow storage, soil moisture and surface water from the total GRACE signal and use the remainder as an estimate of groundwater changes (Rodell and Famiglietti 2002; Rodell et al. 2007, 2009; Chen et al. 2015). Another decomposition approach is based on extracting the dominant spatio-temporal patterns of long-term trends and periodic GRACE signals by means of dimensionality reduction methods. This has been done, for instance, with principal component analysis (Schrama et al. 2007; Rangelova et al. 2007; Schmidt et al. 2008b), independent component analysis (Forootan and Kusche 2012; Frappart et al. 2011b) or multichannel singular spectrum analysis (Rangelova et al. 2010). A last option is based on extracting temporal components (i.e. at each grid cell) using time series decomposition techniques. This approach has been used to assess the properties and the relative importance of the resulting features of temporal variability (Barletta et al. 2012; Frappart et al. 2013). Occasionally, the employed decomposition also assumes that the data follows a predefined pattern, as, for instance, when the seasonal cycle is represented by fitted harmonic functions (Steffen et al. 2009). In this paper, we aim at a temporal decomposition of the time series, making as few assumptions as possible and accounting for the irregular spacing of the GRACE ‘‘months’’. This additive decomposition is summarized in Eq. 1, where the original signal (Xtot) is represented as the sum of a long-term component (Xlong), a seasonal cycle (Xseas) and the remaining subseasonal residuals (Xsub). These high-frequency residuals are expected to be a combination of both a real signal representing subseasonal water storage variability and the noise that is present in the GRACE data. The long-term component (Xlong) is further divided into linear trends (Xlin) and the anomalies with respect to this linear trend, being here referred to as inter-annual variability (Xinter). Xtot ¼ |Xffl {loznfflg} þXseas þ Xsub XlinþXinter ð1Þ 3.1.2 Seasonal Trend Decomposition Using Loess (STL) The Seasonal Trend Decomposition using Loess procedure (STL) introduced by Cleveland et al. (1990) is a robust decomposition method that is used to extract the mean seasonal cycle and to separate the remaining deseasonalized signal into a low- and a high-frequency component, where the low-frequency component should contain only periodicities larger than 12 months. This procedure was already used with GRACE data by Baur (2012) and Hassan and Jin (2014) as a method to derive the long-term component, in Bergmann et al. (2012) to robustly deseasonalize GRACE time series, and in Frappart et al. (2013) to compare terrestrial water storage with monthly rainfall time series in the Amazon basin. It has also been successfully applied, for instance, in a hydro-climatological setting (Gudmundsson et al. 2011) or to extract temperature trends (Dufresne et al. 2013). The STL procedure is based on locally weighted smoothing of the deseasonalized time series in which the smoothing parameters are analytically optimized to minimize spectral leakage between the high- and the low-frequency components. We introduce here an adaptation of the original algorithm allowing us to apply this method to unevenly spaced time series, accounting for the irregular temporal spacing of the GRACE data. The STL procedure consists of passes of different smoothing filters and includes the calculation of robustness California - 121.5°W/38.5°N GRACE (Xtot) Seasonal (Xseas) Linear trend (Xlin) taeenw t][mm 500 t ilavuqE ieghh-1-5000 IInntteerr aannnnuuaall (+XSinutbesre)asonal (Xinter+Xsub) 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 weights in order to account for the possible influence of outliers in the time series. A detailed description of the modified algorithm is presented in Appendix 1. The STL procedure decomposes the time series into the three components: Xseas, Xsub and Xlong (Eq. 1). The latter component Xlong is the long-term (or low-frequency) component of the time series and is further decomposed into the components Xlin and Xinter (Eq. 1). The linear trend Xlin is first estimated from the long-term component (Xlong) using the Theil–Sen estimator (Sen 1968), and Xinter is computed as the deviation from this linear trend (Xinter = Xlong - Xlin). Compared to classical linear regression, using the Theil–Sen slope provides an estimate of the trend that is more robust and less sensitive to large anomalies occurring near the beginning or the end of the time series. This procedure is applied to each grid cell of both the monthly GRACE data and the daily atmospheric forcing so that we obtain decomposed time series for each of these datasets. In Fig. 2, we illustrate how the presented approach decomposes the GRACE signal into the different subcomponents for the case of a specific grid cell located in California. 3.2 Monthly Averaging of the Daily Decomposed Forcing Time Series The decomposed daily atmospheric forcing data need to be averaged to monthly values in order to enable a comparison with the GRACE time series. The common approach for this is to use the monthly arithmetic mean (e.g., Frappart et al. 2013; Forootan et al. 2014a; Ahmed et al. 2014). As a reference method, we use the arithmetic mean of the days exactly covered by each GRACE monthly solution. We thus obtain monthly series for each component of the atmospheric daily series. In addition, we present hereafter a more sophisticated averaging method that accounts for storage processes that specifically influence the high-frequency component (Xsub). 3.2.1 Limitations of the Arithmetic Mean for the Comparison of High-Frequency Anomalies When comparing averaged time series of water storage with precipitation, some systematic errors are introduced simply because of the arbitrarily chosen averaging intervals (e.g., monthly intervals in the present case). As water storage is a state and precipitation a flux variable, temporal averages can at times cause a mismatch of the two monthly time series, especially in the case of high-frequency anomalies. A typical example is when a very large precipitation event occurs just at the end of a given month: this extreme event will have a large effect on the precipitation average of the given month but its influence on water storage will be most relevant for the subsequent months. Such artefacts are often falsely attributed to observational errors. In order to address this issue, we propose an alternative to the arithmetic mean that takes the effect of earlier precipitation into account. 3.2.2 Comparing Flux and State Variables at Different Temporal Resolutions Hereafter, precipitation anomalies correspond to a time-dependent flux variable, denoted f ðtÞ where t ¼ ft1; . . .; ti; . . .; tng is an evenly spaced time vector of length n, with units of days. Similarly, daily water storage anomalies correspond to a time-dependent state variable denoted sðtÞ. In our case, the state variable sðtÞ is not observed at the daily time scale; however, the GRACE product provides average values of sðtÞ for arbitrary time intervals which approximatenly correspond tooa month. We define this new averaged time series as s ðt Þ, where t ¼ t1; . . .; tj ; . . .; tm is an unevenly spaced time vector of length m corresponding to the GRACE ‘‘months’’. The relation between sðtÞ and s ðt Þ can be represented by the arithmetic mean (see Fig. 3 for a schematic illustration of the presented relations): s tj 1 ¼ nj X ti2½aj;bj f tj 1 ¼ nj ti2½aj;bj ... ... sðtiÞ where aj and bj correspond to the edges of the jth time interval (e.g., of the jth GRACE month) and nj is the number of days falling within this interval (nj ¼ bj aj). Our main concern is now to determine the relation between f ðtÞ and f ðt Þ. As mentioned above, a common approach is to compute the mean of the daily values over the given time intervals. Analogously to Eq. 2, this simply corresponds to: X Here, we suggest the use of a weighted mean of f ðtÞ as an alternative approach: ð2Þ ð3Þ : Daily TWS Eq. 2 : Monthly TWS : Daily forcing Eq. 3 or 4 : Monthly forcing from the 1st to the nth day. For quasi-monthly time steps, t ¼ t1; . . .; tj ; . . .; tm going from the 1st to the mth month, aj and bj denoting the first and last days of a given month X n i¼1 f tj ¼ where the normalized weights W^ tj ; ti , which will be defined in the next section, depend both on tj and ti and have the property that: X n i¼1 For illustrative purposes, this summation is shown in Fig. 4. The two examples correspond to the case of a flux anomaly f ðtiÞ occurring either before (Fig. 4a) or during (Fig. 4b) the given time interval aj; bj . The last step is to ensure that the property set by Eq. 5 is fulfilled by normalizing the weights (Eq. 8): When this is done with a fixed tj and for all values of ti, we obtain the averaging filter illustrated in Fig. 5 for different values of s—the free parameter controlling the rate of the exponential decay. From this figure, we see that weights are assigned to flux anomalies including to those occurring before the time interval tj . Additionally, we provide a more 3.2.3 Weights Based on Integrated Exponential Decay Functions A simple way to represent the effect of a short-term precipitation anomaly (e.g., a daily precipitation event) on the subsequent state of water storage is the exponential decay function. This is equivalent to assuming linear storage components (bucket models), which is common practice in conceptual hydrological modelling (Beven 2012). Here, we assume that the influence of a flux anomaly (e.g., a precipitation event) on the state variable (e.g., water storage) will decrease exponentially with time. Formally, we define wðt; tiÞ as the influence of a given flux anomaly f ðtiÞ observed at time ti on the subsequent values of the state variable sðtÞ at time t [ ti. wðt; tiÞ ¼ 0; e 1sðt tiÞ; if t\ti if t ti where s is a free parameter controlling the rate of exponential decay and is expressed in units of time (e.g., in days). The influence of the given flux anomaly f ðtiÞ on the earlier values of the state variable (i.e. when t\ti) is of course zero. However, wðt; tiÞ only represents the influence of f ðtiÞ on the subsequent daily values of sðtÞ, but we are in fact interested in the influence of f ðtiÞ on the values of s ðt Þ—the monthly values. For a given tj , summing wðt; tiÞ over the corresponding time interval t 2 aj; bj yields: t2½aj;bj ð4Þ ð5Þ ð6Þ ð7Þ ð8Þ Fig. 4 Illustration of Eqs. 6 (red) and 7 (blue). The red curve depicts the exponentially decaying influence of a given daily flux anomaly (precipitation) occurring at time ti on the state variable (water storage) at subsequent time steps. The summation of this influence over the interval covered by a given GRACE month corresponds to the relative weight (blue area) assigned to the flux at time ti. See the text for a description of the different symbols (a) (b) GRACE month GRACE month practical formulation of this weighting function obtained after integrating and normalizing Eq. 6 analytically (Appendix 2): 3.2.4 Shape and Properties of the Weighting Function The parameter s controls the rate of exponential decay and will hereafter be referred to as the decay time scale of the weighting function. Inverting Eq. 6 for s shows that s corresponds to the number of time steps (e.g., days) after which the influence of a given flux anomaly f ðtiÞ will have reduced to 1/e & 37 % of its initial influence at time ti. An interesting property is that when s tends to small values (see Fig. 5 for s = 1), W^ tj ; ti converges very quickly to a weighting function that is almost equivalent to the arithmetic mean performed over the interval aj; bj (i.e. Eq. 3). Small values of s Fig. 5 Illustration of the shape of the weighting function (Eq. 9) for different values of the decay time scale s. The y-axis corresponds to the normalized weight (W^ ) that is applied to the daily flux time series when it is averaged to the approximately monthly GRACE resolution 330 300 270 240 210 180 150 120 Days preceding the end of a GRACE month ( ) 90 60 GRACE month 30 0 correspond to small decay time scales, indicating that a single flux anomaly will not have a prolonged effect on the state variable. A hydrological interpretation of this feature suggests a short mean residence time of the water store. Inversely, large values of s imply longer residence times and, therefore, more weight is given to anomalies occurring before the time interval of interest. In such a case, it is interesting to note that anomalies occurring near the end of the given time interval are assigned small weights. Hence, the difference between the presented weighting scheme and an arithmetic mean becomes more important for larger s. Since the value of s at each grid cell is unknown in our application, it needs to be estimated from the data. Here, we optimize the agreement between the monthly averaged subseasonal forcing (i.e. Xsub of f ðtÞ) and the subseasonal monthly GRACE (i.e. Xsub of s ðt Þ) by maximizing the squared product–moment correlation coefficient. In the presented study, this weighting function is used for the analysis of the subseasonal component only (Sect. 4.5). 3.3 Significance Testing and Correlation Analysis 3.3.1 Linear Trends A common nonparametric test for detecting monotonic trends in hydro-meteorological time series is the Mann–Kendall rank-based test. However, serial correlation (autocorrelation) in time series has been shown to heavily influence the power of this test (Yue et al. 2002), and several methods have been proposed to address this issue (Hamed and Rao 1998; Yue and Wang 2004; Hamed 2009). Here, we use the modified Mann–Kendall trend test described by Yue and Wang (2004) on deseasonalized GRACE time series (Xtot - Xseas). In this test, the autocorrelation estimated from the deseasonalized and detrended time series is used to compute an effective sample size, which is then used to correct the Mann–Kendall statistic. In addition, as the trend test is performed locally (i.e. at each grid cell) and due to the high spatial autocorrelation of the GRACE data, there is an increased probability that the null hypothesis is falsely rejected (Wilks 2011). Hence, we additionally control this false discovery rate (FDR) using the approach described by Benjamini and Hochberg (1995), which has shown good performance when applied to climate data (Ventura et al. 2004; Wilks 2006; Gudmundsson and Seneviratne 2015). The trends are considered statistically significant when the p value falls below a critical value (p \ 0.01). 3.3.2 Inter-Annual Anomalies Regarding the inter-annual anomalies (Xinter), the degree of linear association between GRACE and the atmospheric forcing is quantified with the product–moment correlation coefficient. As the inter-annual anomalies correspond to the low-frequency component of the GRACE signal, they exhibit important serial correlation, which prevents the use of conventional hypothesis testing techniques (e.g., t test). Here, we use moving block bootstrapping in order to estimate the null distribution of the correlation coefficient at each grid point (Mudelsee 2014). Although there is no standard recommendation on the selection of an optimal block length, a good starting point is to use a block length larger than the decorrelation time (i.e. the number of time steps after which the serial correlation is not significant anymore). Based on this criterion, we find that a block length of 20 months is sufficient for our application. We perform 10,000 bootstrap replications at each grid point and estimate the 95 % confidence intervals from this null distribution. A correlation coefficient is declared significant when it does not belong to the range of the local confidence interval. 3.3.3 Seasonal Cycle Previous studies have shown that there is often a temporal lag between the seasonal cycle of precipitation and terrestrial water storage (e.g., Papa et al. 2008; Ahmed et al. 2011; Frappart et al. 2013). It is also known that water storage and surface temperature are related through evapotranspiration and snow melt; however, differences in the phasing of GRACE versus these atmospheric variables were, to our knowledge, never surveyed at a global scale. We define the phase shift as the lag (in months) minimizing the residual sum of squares between the standardized seasonal cycles of both GRACE and the atmospheric forcing. When these paired seasonal cycles strongly differ in shape, this procedure can sometimes lead to meaningless lag values. A t-test of the Pearson product–moment correlation between the time–lagged seasonal cycles is used to filter out these potentially misleading values (p \ 0.01). We additionally control the FDR following Benjamini and Hochberg (1995). 3.3.4 Subseasonal Residuals Similarly as for the inter-annual anomalies, the product–moment correlation coefficient is used to quantify the degree of linear association between GRACE and the atmospheric forcing data. The subseasonal residuals (Xsub) correspond to the high-frequency component and are thus the least affected by serial correlation. However, we found that these time series still contain minor but significant serial correlation (not shown). For consistency, we thus use an identical significance testing setting as for the inter-annual anomalies (i.e. a moving block bootstrapping). 3.4 Identifying Droughts in the GRACE Record Here we investigate the average storage deficit during drought events identified using an approach based on Thomas et al. (2014). This approach defines 1) storage deficit as a negative departure (in mm) from the seasonal cycle and 2) drought duration as the number of months with continuous deficits. The average storage deficit simply corresponds to the arithmetic mean of the storage deficit observed during a given drought event and is used as a measure of average drought intensity. Here two differences compared to Thomas et al. (2014) are introduced. First, we remove the linear trends from the time series prior to drought identification. The reason is that strong linear trends can result in one end of the time series being systematically above/below the seasonal cycle. In such a case, the proposed method would have a tendency to underestimate/overestimate the magnitude of dry events. Hence for the purpose of this study, linear trends are removed prior to the analysis and potential decadal drying trends are discussed in a separate section. Our analysis is thus based on the sum of the inter-annual and subseasonal components only (Xinter ? Xsub, also see Fig. 2). Occasionally, drought events occurring at the end or the beginning of the time series can be large enough to influence the trend estimate itself, even when using the Theil–Senn slope to reduce this effect. Hence, it is important to note that, in some cases, removing the linear trends may cause an underestimation of the drought intensity. Second, the minimum duration for considering a drought event is defined as a period of three consecutive months with water storage deficit. Unlike Thomas et al. (2014), we apply this criterion only to the inter-annual component Xinter (see Fig. 2) and not to the sum of the inter-annual and subseasonal components (Xinter ? Xsub). The reason is that, compared to the basin-scale assessment of Thomas et al. (2014), subseasonal variability (Xsub) is larger at the grid level and including it would otherwise considerably reduce the probability of observing long periods with consecutive deficits. 4 Global Hydrological Variability in the GRACE Data 4.1 Distribution of GRACE Variance Among Temporal Components The relative magnitude of the three components extracted from the STL procedure (Xlong, Xseas and Xsub) can be evaluated by comparing each component’s variance to that of the total signal. As shown in Fig. 6, the relative magnitude of each of the different components is subject to high spatial variability across the world. To ease the interpretation, Fig. 6 can also be compared to the standard deviation of the total signal in Fig. 1. As already identified in early studies (Wahr et al. 2004), the seasonal cycle is dominant in many tropical regions like the Amazon basin, Central Africa and India. A notable exception is the IndoAustralian archipelago where the GRACE signal is heavily perturbed by signal leakage from the ocean as well as gravity anomalies consecutive to the 2004 Sumatra earthquake. The seasonal cycle is also dominant at higher latitudes, particularly in Siberia and in northwestern America, although these regions do not have the largest variance in absolute terms. Subseasonal variability (Xsub) is dominant in regions where the GRACE signal has already a relatively low variance (Fig. 1) and is most likely dominated by noise such as in the Sahara desert. Although we do not further investigate this matter, it is interesting to note that Arctic coastal regions such as the coasts of Northeast Siberia and Canada seem to be mostly affected by subseasonal variability. We also observe that many regions of the world are dominated by inter-annual variability (Xlong). The signal found in Greenland and Antarctica, parts of Alaska and the Hudson Bay is the result of the interplay between ice mass loss, other water storage changes and glacial isostatic adjustment. As a result, these regions require a specific treatment before conclusions can be drawn concerning the dominant features of hydrological variability (Velicogna et al. 2014). Other regions particularly dominated by longterm variability include the south-western Central USA as well as the Middle East, some of which are already documented in the literature as being influenced by decadal droughts and long-term trends in groundwater storage (Long et al. 2013; Voss et al. 2013; Forootan et al. 2014b). Other interesting features include the Lake Victoria and the Aral Sea where longterm surface water variations can be related to both human activities and climate variability (Swenson and Wahr 2009; Singh et al. 2012). Finally, some regions in the southern hemisphere like Australia and Argentina were also shown to exhibit an important interannual variability that can be related to the El-Nin˜ o Southern Oscillation (ENSO) (Garc´ıaGarc´ıa et al. 2011; Abelen et al. 2015). 4.2 Linear Trends In this section, we will assess in further detail the relative importance of linear trends (Xlin) versus nonlinear inter-annual variability (Xinter) by looking at the magnitude of each of these two components in the long-term variability (Xlong): 2 rlin Rlin=long ¼ rlong 2 ¼ 1 2 rinter 2 rlong ð10Þ This formulation is also equivalent to the coefficient of determination of the linear trend as estimated by the Theil–Sen slope with respect to the long-term component. The colour scale in Fig. 7a represents the ratio of the linear trend variance to that of the whole longterm component (Eq. 10) and shows how the long-term component (Xlong) variance is partitioned between linear (Xlin) and nonlinear trends (Xinter). The sign, magnitude and significance (p \ 0.01) of the linear trends are shown in Fig. 7b. Note the truncated colour scale for negative trends beyond -30 mm/year. We observe that the long-term variability in large areas of Greenland and Antarctica is dominated by a linear trend. This can be related to ice mass loss once glacial isostatic adjustment has been accounted for (Ramillien et al. 2006; Velicogna and Wahr 2006; Chen et al. 2006b; Wouters et al. 2008; Velicogna 2009; Velicogna and Wahr 2013). Examples Partition of the long-term component (Xlong) into linear trends (Xlin) and inter-annual variability (Xinter) non-linear (b) 0% 25% 50% 75% Linear trends of other regions of the cryosphere concerned with documented linear trends include Alaskan mountain glacier melting (Chen et al. 2006a; Arendt et al. 2013; Larsen et al. 2015) and icefield melting in Patagonia (Chen et al. 2007b; Ivins et al. 2011). Strong linear trends located close to the Hudson Bay have been related to glacial isostatic adjustment (Tamisiea et al. 2007), and recent studies focusing on Arctic regions showed that both isostatic and hydrological trends contribute to the observed signals (Frappart et al. 2011a; Wang et al. 2013). Many pronounced negative trends can also be observed in other regions of the world. One of the most recognized drying trends occurs in north-west India and is related to groundwater depletion (Rodell et al. 2009; Chen et al. 2014). Most of the long-term signal that is dominating GRACE variability over the Middle East, the Caspian Sea and the Aral Sea can be attributed to a drying trend partly due to anthropogenic water abstraction (Singh et al. 2012; Voss et al. 2013; Joodaki et al. 2014; Forootan et al. 2014a). On the contrary, the region of Lake Victoria is dominated by nonlinear variations in the long-term component, which have been related to hydropower dam operations, precipitation anomalies and ENSO (Awange et al. 2008; Swenson and Wahr 2009; Becker et al. 2010; Hassan and Jin 2014). Another important drying trend can be found in Argentina, especially in the southern part of the La Plata basin (Chen et al. 2010b; Abelen et al. 2015). Finally, the drying trend documented by Crowley et al. (2006) in the Congo river basin for the period 2002–2006 is found to be insignificant based on the current analysis. The significance analysis also identified regions with trends which have not been identified yet or have only drawn little attention so far. For instance, the extended drying trends located to the North of both the Black Sea and the Caspian Sea could be potentially investigated in more detail. Interestingly, small but statistically significant drying trends (-3 mm/year) can also be found over large portions of the Sahara desert. So far, little attention has been devoted to GRACE variability in this region as most of the signal is contaminated by noise. Nevertheless, these drying trends have been partly documented (Ahmed et al. 2014; Ramillien et al. 2014) and to some extent attributed to groundwater extraction from fossil aquifers in the Sahara region. Significant positive trends can also be found in southern Africa, near the Upper Zambezi and Okavango river basins as well as in the Sahel and the Niger basin, and these trends have already been well documented (Ramillien et al. 2014). In a comparison with rainfall observations from different sources, Ahmed et al. (2014) have found that the increasing trend in the Niger basin could be related to an increase in precipitation; however, this was not the case for the Upper Zambezi and the Okavango basins. Although Ahmed et al. (2014) suggest that this could be due to longer residence times in these watersheds, we feel that more investigation is still required to explain the very strong positive trend in this region. A positive trend is also found in the Amazon basin and has been described, for instance, in Chen et al. (2010a) and could, to a certain extent, be related to precipitation anomalies based on an analysis of the period 2003–2008 by Xavier et al. (2010). The linear trends derived over north-western China raise some concerns with respect to a possibly spurious origin. The alternating bipolar patterns of positive and negative trends, oriented along the meridian 100 E, could be due to some accidental disturbance originating in the processing of the GRACE data. A very similar pattern could already be found in Fig. 8 of Frappart et al. (2011b), which compared trends derived after different postprocessing methods for the period 2003–2008. It is possible that these trends found over China are specific to the destriping algorithm of Swenson and Wahr (2006) since they are not reproduced by the other postprocessing methods investigated by Frappart et al. (2011b). However, Feng et al. (2013) were also able to relate drying trends in the Beijing region to groundwater observations. Consequently, special care should be taken when interpreting trends from the current GRACE Tellus estimates in that region. 4.3 Inter-Annual Anomalies The inter-annual anomalies correspond to the nonlinear part of the long-term component (Eq. 1). In this section, we assess the degree to which the inter-annual anomalies derived from the GRACE time series can be correlated with the inter-annual anomalies of the atmospheric forcing. Figure 8a depicts the correlation between the inter-annual water storage and precipitation anomalies. We observe moderately high positive correlations Xinter: GRACE vs Precipitation Xinter: GRACE vs Temperature ([0.6) between these two components for parts of Europe, Russia and North America, which indicate a positive effect of precipitation on terrestrial water storage. Correlations are more erratic over subtropical and equatorial regions, possibly resulting from deficiencies in the ERA-Interim precipitation forcing, which are known to be more pronounced, for instance, over Africa and South America (Simmons et al. 2010). In these regions, other precipitation datasets based either on ground measurements or satellite observations may give different results. For example, we find relatively low correlations between inter-annual water storage and precipitation over the region of the African Great Lakes; however, Becker et al. (2010) found a seemingly good agreement with GRACE when using monthly precipitation data from the Global Precipitation Climatology Project (GPCP). For Africa and South America, Morishita and Heki (2008) found that mass changes from GRACE could be related to precipitation anomaly patterns extracted from the Climate Prediction Center Merged Analysis of Precipitation (CMAP). Over the Amazon, Chen et al. (2010a) related inter-annual water storage change to precipitation anomalies (from GPCP) and ENSO indices, while Frappart et al. (2013) found that the inter-annual variability of water storage was in reasonable agreement with precipitation from the Tropical Rainfall Measuring Mission (TRMM). Although correlation coefficients are mostly positive, significant negative correlations can occasionally be found but seem either accidentally caused by perturbations of the long-term gravimetric signal by large earthquakes (e.g., Sumatra region) or would need to be confirmed in a regional investigation (Caspian Sea area). The same analysis was performed with ERA-Interim 2 m temperature (Fig. 8b). We find that the correlation between inter-annual water storage and temperature is negative in most cases. Correlations are moderately strong over parts of North America, South America, southern Africa, India and Australia. Such a negative relationship is usually expected not only since temperature is an important driver for evaporative demand but also because low moisture availability can result in higher temperatures (Seneviratne et al. 2010; Mueller and Seneviratne 2012). However, correlations found in the southern hemisphere could also be related to confounding factors such as atmospheric circulation patterns (e.g., ENSO) and hence do not necessarily indicate a direct link between temperature and terrestrial water storage. In addition, significant positive correlations between long-term temperature and water storage anomalies can be found over parts of Siberia. Long-term water storage anomalies in this region may be related to interactions with permafrost although such relationships are still unclear at this stage (Velicogna et al. 2012; Vey et al. 2013). Non-significant correlations can be due either to other unidentified sources of long-term variability in the GRACE data or errors and biases in the long-term variability of the ERAInterim atmospheric forcing. However, the absence of correlation with either precipitation or temperature in some regions could also indicate that long-term variability in the atmospheric forcing is not controlling or controlled by terrestrial water storage, i.e. that there is no obvious relationship between these variables at the inter-annual time scale. In addition, the literature covered in the section on linear trends already showed that anthropogenic groundwater withdrawal and dam operations could be responsible for longterm changes in the terrestrial water storage variations retrieved from GRACE. Finally, we note that a very large number of locations exhibit moderate correlations (between 0.3 and 0.5 in absolute value), which are actually not significant once serial correlation is taken into account in hypothesis testing (using bootstrapping). This is also an indication that analyses of the inter-annual variability of water storage would greatly benefit from the longer record length that may be provided in the future by the GRACE Follow-On mission scheduled for launch in 2017. 4.4 Seasonal Cycle The STL decomposition provides a data-driven way of estimating the seasonal cycle which, in contrast to the common practice, does not rely on harmonic models (fitted sines and cosines, e.g., Wahr et al. 2004; Hinderer et al. 2006; Schmidt et al. 2008b). Here we characterize GRACE seasonality by mapping the months with the maximum and the minimum of the seasonal cycle of water storage and show that they generally follow latitudinal bands (Fig. 9a, b). In the Northern Hemisphere, the peak in terrestrial water storage generally occurs in spring for the cold and temperate regions and in autumn for the subtropical regions (and vice versa in the Southern Hemisphere). The minimum water storage occurs in autumn for the cold and temperate regions and in spring for the subtropical regions (and oppositely in the Southern Hemisphere). In subarctic regions, there is a clear latitudinal trend towards a later maximum, likely corresponding to the delayed response of snow melt to temperature at higher latitudes. Interestingly, the seasonal Maximum Minimum 1 2 3 4 5 8 9 maximum appears to be also delayed near large inland reservoirs (e.g., the Great Lakes and the Caspian Sea), which potentially reflects the influence of run-off and storage processes and could be subject to further investigations. In most regions, the months with maximum and minimum terrestrial water storage are spaced by an interval of 6 ± 1 months. However, this is not always the case: in northern India, the maximum terrestrial water storage occurs in September and the minimum in May, which is consistent with the effect of the June–September monsoon. These maps can, for instance, be directly compared with hydrological models (see Fig. 6 in Gu¨ ntner et al. 2007b for an example with a closely resembling colour scale). In addition, it is worth mentioning that the seasonal cycle of water storage over African regions located close to the equator (e.g., the Congo basin, Lake Victoria) exhibits a strong secondary peak, which likely corresponds to the oscillation of the inter-tropical convergence zone (not shown). This secondary peak is also present— although less pronounced—over Ecuador, southern India and Indonesia but is completely absent over the Amazon basin (not shown). The phase shift between GRACE and the seasonal cycles of precipitation and temperature is shown in Fig. 10a, b. Over most tropical and subtropical regions, the seasonal Phase shift: GRACE vs Precipitation Phase shift: GRACE vs Temperature -11 -10 -9 -8 -7 -6 (+1) (+2) (+3) (+4) (+5) (+6) -5 -4 -3 -2 -1 0 in phase out of phase in phase months GRACE Atmospheric Fig. 10 Phase shift (in months) between the seasonal cycle of water storage and the seasonal cycle of a precipitation and b temperature. Small phase shifts (-1, 0 or -11(?1) months) indicate that the atmospheric forcing is nearly in phase with water storage, whereas large phase shifts (-7, -6 or -5 months) indicate that they are out of phase. Stippling indicates regions where the correlation between optimally phase shifted seasonal cycles is not significant (p \ 0.05) peak of precipitation typically occurs 1–2 months earlier than the peak in water storage, likely due to the effect of storage processes. Very similar lags have been found for the Amazon subbasins (see Table 3 in Frappart et al. 2013), for selected regions over central Africa (see Fig. 3 in Ahmed et al. 2011) as well as by Rieser et al. (2010) over Australia (all three studies used satellite precipitation data from TRMM). On the contrary, subarctic and inland temperate regions experience the highest precipitation during the warmer summer months, approximately 3–5 months later than the spring maximum in water storage. For coastal subarctic areas, the precipitation maximum tends to occur in autumn due to greater temperature differences between the ocean and land, resulting in a 5- to 7-month phase shift between water storage and precipitation (e.g., Alaska, British Columbia and Scandinavia). More details concerning the phasing of GRACE with snow storage and discharge measurements can be found in Frappart et al. (2011a). The seasonal cycle of temperature is generally out of phase with respect to the seasonal cycle of water storage (Fig. 10b). In most temperate and subarctic regions, the peak temperature typically occurs in summer, 2–3 months earlier than the autumn minimum in water storage. Over tropical regions, the seasonal cycle of temperature is completely opposed to the water storage cycle, with corresponding phase shifts of 4–7 months. This anti-phasing between water storage and temperature is likely related to the effects of both temperature and radiation on evapotranspiration. Over equatorial regions, the seasonality of temperature is much less pronounced but still lagging the water storage cycle by 3–4 months. The southern part of China exhibits a very specific pattern, with maximum temperatures occurring in summer and seasonal water storage peaking in late summer, resulting in an almost perfect phasing between water storage and temperature. 4.5 Subseasonal Residuals Figure 11a shows the correlation between the high-frequency components of GRACE and precipitation averaged with the new averaging scheme presented in Sect. 3.2. Significant positive correlations are found over many regions of the world, indicating that a large fraction of high-frequency GRACE variability can be statistically related to short-term anomalies of the precipitation forcing. Interestingly, significant correlations can also be found over large portions of Indonesia, although the GRACE signal in this region is usually believed to be strongly deteriorated by signal leakage from the ocean. A possible explanation might be that short-term precipitation variability in this tropical monsoon region is large enough to overcome the higher errors associated with coastal and insular regions. A notable exception to the global pattern is the Congo river basin where no significant correlations can be found. This area corresponds to a major convective region for the global climate system which is still poorly represented by atmospheric reanalyses in comparison with other regions (Washington et al. 2013). Many extremely arid regions also show non-significant correlations (Sahara, Atacama, Taklamakan and Gobi deserts), confirming the view that high-frequency GRACE variability in these regions is dominated by noise. In order to assess the influence of the new averaging method, we can visually compare the correlations shown in Fig. 11a to the correlations obtained with a simple monthly arithmetic mean of the daily residuals (Fig. 11b). The exponential decay approach enhances the correlations over most regions of the world, with an average increase of ?0.3 (excluding regions which exhibit non-significant correlations in Fig. 11a). Figure 12a enables an even more direct comparison between the distributions of the correlations shown in Fig. 11a, b. This illustrates the value of the proposed weighting scheme and reveals that using monthly arithmetic averages of precipitation may have resulted in underestimating the relation with water storage on the subseasonal time scale. This finding is of particular interest for studies comparing GRACE data to monthly precipitation time series (e.g., Forootan et al. 2014a) which typically make use of monthly precipitation averages. Values of the calibrated decay parameter s used to compute the monthly averages from the daily precipitation data are shown in Fig. 11c. Overall, decay time scales exhibit systematic spatial variability that could potentially be related to many different factors, including climatic conditions as well as soil and vegetation characteristics. The probability Xsub: GRACE vs Precipitation exponential decay weighting Xsub: GRACE vs Precipitation arithmetic mean 8000 7000 s t in 6000 o p id 5000 r fg 4000 o r eb 3000 uNm2000 1000 Arithm. mean Exp. decay Exp. decay, p<0.05 -1 -0.8 -0.6 -0.4 -0.2 density distribution of this parameter is also shown in Fig. 12b, and we find that significant values generally range between 10 and 200 days with a median value of approximately 50 days. Based on the weighting function (Eq. 9), it can be calculated that for the median value of s = 50 days, the precipitation residuals of the first 100 days preceding the beginning of a GRACE month account for 65 % of the monthly average. On the contrary, days covered by the time interval of a given GRACE month account for only 25 % of the monthly average. This shows that, on the subseasonal time scale, precipitation preceding a GRACE month usually has a higher impact on correlations with terrestrial water storage than the precipitation of the coinciding month. This is due to the fact that the influence of high-frequency precipitation anomalies on regional hydrology tends to decay with time so that precipitation events occurring just before or at the very beginning of a GRACE month have a higher impact on the average water storage of a given month. Conversely, precipitation anomalies occurring during or at the very end of a GRACE month have a lower impact on the water storage anomalies, or may even occur after the latest GRACE overpass. Figure 13a shows the results of the same analysis performed with the temperature data. It reveals regions where the integrated effect of antecedent temperatures can be statistically related to water storage anomalies. Temperature is one of the main controls for evaporative Xsub: GRACE vs Temperature exponential decay weighting Xsub: GRACE vs Temperature arithmetic mean demand, and hence, negative correlations are expected and can indeed be found over many regions of the world, especially over South America, South Africa, the region of the African Great Lakes, India, Indonesia and northern Australia. As for precipitation, the use of the exponential decay approach leads to enhanced correlations when compared to the arithmetic mean (Fig. 13b). However, improvements are less important than for precipitation and are often concentrated in regions where a significant relationship can already be found with the simple mean. Decay time scales over these regions (Fig. 13c) generally fall between 1 to 3 months, yielding a data-driven first-order estimate of how long temperature anomalies can significantly impact the subsequent state of terrestrial water storage. Positive and significant correlations can occasionally be found over some areas, notably over Siberia, Scandinavia and Antarctica. For these regions, we hypothesize that these positive correlations could reflect the tendency of warm winters to be more humid in comparison with cold winters. On the other hand, warm summers are also expected to increase snow melt so that we cannot come to a definitive conclusion for these regions based on the presented results. 4.6 Droughts In Fig. 14a, we show the maximum average storage deficit (see Sect. 3.4) ever observed for all drought events identified in the GRACE record. The year corresponding to this maximum is depicted in Fig. 14b for events with a magnitude larger than 30 mm. This threshold was chosen to mask out smaller features which are difficult to interpret in a global assessment but may still be relevant in a regional context. Many droughts that have previously been documented in the GRACE literature can be identified, notably the 2010 Amazon drought, which is additionally illustrated in Fig. 15a. Drought events in the Amazon basin were shown to be related to precipitation deficits and ENSO (Davidson et al. 2012; Frappart et al. 2013). A multi-year drought is also found for the period 2004–2008 (Fig. 15a) and likely corresponds to the multiple consecutive dry years identified, for instance, by Frappart et al. (2012) and Thomas et al. (2014). Note that this is related to the chosen drought metric, which might not capture all relevant aspects. The ongoing drought in the Central Valley of California is also identified in Fig. 14b, and the time series of the average storage deficit (Fig. 15b) shows that this region also suffered from multiple dry episodes in previous years. This was already identified in previous studies which related these recurrent drought events to severe groundwater depletion (Famiglietti et al. 2011; Famiglietti 2014; Chen et al. 2015). Other documented events identified in Fig. 14b include the 2008–2009 drought in the La Plata basin (Abelen et al. 2015), the 2010–2013 drought in Texas (Long et al. 2013), the 2007–2009 drought in the south-eastern USA (Houborg et al. 2012) and the 2012–2015 North American drought (Chew and Small 2014; Hoerling et al. 2013). We also identify the 2006–2007 dry conditions over Lake Victoria (Swenson and Wahr 2009) and the African Great lakes (Becker et al. 2010) and the 2006–2008 drought in the Zambezi basin (Thomas et al. 2014), which are, in this analysis, captured together as a large-scale and spatially contiguous event. Drought conditions can also be found in northern India for the period 2009–2010 even though the linear trend due to groundwater depletion has been removed from the data prior to drought identification. The year 2009 was indeed shown to be the driest year of the decade for this region in terms of precipitation (Chen et al. 2014) and resulted in higher groundwater abstraction rates. Our analysis shows that the average storage deficit was consecutively maximal in 2010. Maximum average storage deficit Year of the maximum C [mm] A 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 Fig. 14 a Maximum value of the average storage deficit observed in the period April 2002–August 2015, expressed in mm of equivalent water height. b Year corresponding to the maximum value of the average storage deficit, showing only regions with a deficit larger than 30 mm. Letters A–C correspond to the location of the time series in Fig. 15 The Sumatra region also exhibits an important ‘‘deficit’’ which, as confirmed by a local investigation, is probably an artefact caused by the 2004 earthquake. In Australia, multiyear droughts have been related to precipitation deficits (Garc´ıa-Garc´ıa et al. 2011). However, due to the long duration of these decadal drought events, the average storage deficit is lower. Our results also reveal undocumented features found in the GRACE record, such as a dry event from 2012 to 2014 over south-eastern Europe (Figs. 14b, 15c) as well as a severe drought in the Sao Paulo region and a moderate drought over North Korea in 2015 (both still ongoing at the time of writing). A dry period can also be identified during 2010–2011 to the North of the Caspian Sea and is likely associated with the 2010 Russian heatwave. Many other events can also be found over central Russia and were, to our knowledge, never identified using GRACE data. Amazon - 63.5°W/4.5°S Average storage deficit 0 eag ] -10 r m so [m-20 t eag itic -30 re fe -40 vA d -50 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 In this study, we have decomposed the GRACE time series into (1) linear trends, (2) nonlinear inter-annual anomalies, (3) seasonal cycles and (4) subseasonal residuals. The relative importance of each of these components with respect to the original GRACE signal has been evaluated, allowing for a global assessment of the dominant features of temporal variability in terrestrial water storage. In most cases, the GRACE signal is dominated by seasonal or/and long-term variability, while subseasonal variability generally accounts for a small fraction of the total signal variance. Partitioning the longterm variability into linear trends and nonlinear components reveals that some regions are dominated by linear trends, while nonlinear inter-annual variability is prevalent in others. The magnitudes of the linear trends have been quantified using the robust Theil–Sen estimator, reproducing many already documented trends but also revealing some features that had not been identified previously. In addition, the significance of the trends was evaluated using statistical hypothesis testing techniques which take serial correlation (autocorrelation) into account, contrasting the common practice in the GRACE literature. In a more detailed analysis, each component of temporal variability (except linear trends) has been compared with equivalent components extracted from daily precipitation and temperature time series of the atmospheric reanalysis ERA-Interim: • Inter-annual variability We have found that the inter-annual variability of GRACE can be only partly related to anomalies in precipitation and temperature, confirming the results of previous regional studies. Although limitations of the considered atmospheric reanalysis may alter the results at the regional level, this suggests that the inter-annual variability of GRACE is only partly related to the investigated atmospheric drivers, potentially highlighting the role of human water use as additional driver. Seasonal variability We have provided a comprehensive overview on the seasonality of terrestrial water storage and related it to the seasonal cycles of both precipitation and temperature. In tropical and equatorial regions, the seasonal cycle of precipitation was generally found to precede terrestrial water storage with a temporal lag of one to 2 months, while the seasonal cycle of temperature would typically be phase shifted by 6 months with respect to water storage. However, this was clearly not the case in temperate and cold regions, which is probably due to the more complex interplay between precipitation, storage processes, snow dynamics and temperature. Subseasonal variability We have shown that high-frequency variability of the GRACE record can be reconstructed from precipitation anomalies once an adequate averaging filter is applied to the daily precipitation forcing. This filter was designed to explicitly take the effect of earlier precipitation into account when comparing daily precipitation series with monthly GRACE data. This new method yields substantially better results compared to the classical approach based on monthly arithmetic means, providing a new perspective on the hydrological value of subseasonal (month to month) fluctuations of the GRACE signal, which have partly been interpreted as noise in previous studies. Droughts Finally, we have surveyed extreme dry events in the GRACE time series. The most important anomalies in terms of water storage deficits were documented on a global scale and related to droughts already described in the existing literature. Undocumented features were also identified using this global approach and will be subject to further investigation. In summary, we have surveyed key features of temporal variability in the GRACE record and related them to the dominant atmospheric drivers, in contrast to the common practice of comparing GRACE terrestrial water storage to hydrological model simulations. We have related our results to physical interpretations from the rich body of regional GRACE studies, resulting in a comprehensive overview which will both contribute to a general understanding of terrestrial water storage and provide a global observation-based reference for hydrologists and climate scientists. As a novelty, we have shown that highfrequency (month to month) fluctuations of the GRACE signal contain a meaningful hydrological signal, which can be reconstructed from daily precipitation forcing. These findings have important implications for the assessment of GRACE uncertainties as well as for comparisons with hydrological model simulations. Acknowledgments This research was funded by the ERC DROUGHT-HEAT project (Contract No. 617518). GRACE land mass grids were processed by the Jet Propulsion Laboratory (JPL) GRACE Team and are available at http://grace.jpl.nasa.gov, supported by the NASA MEaSUREs Program. The European Centre for Medium-Range Weather Forecasts (ECMWF) and ERA-Interim teams are gratefully acknowledged for providing the atmospheric reanalyses. We also thank the two anonymous reviewers whose comments and suggestions helped improve and clarify this manuscript. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Appendix 1: STL for Unevenly Spaced Time Series STL was introduced by Cleveland et al. (1990) for evenly spaced time series. In this appendix, we present a summary of a simplified STL procedure that includes a modification of the original algorithm for unevenly spaced time series. An important difference with the original STL authored by Cleveland et al. is that the seasonal cycle is defined as constant, i.e. it is not allowed to vary in shape or amplitude over time. We also introduce a modification in the computation of robustness weights that accounts for time series exhibiting large differences in variance (heteroskedasticity) between seasons. Our focus is to decompose the total signal of a given time series X ¼ fx1; . . .; xi; . . .; xng associated with a time vector T ¼ ft1; . . .; ti; . . .; tng into a constantly repeating seasonal component Xseas, a long-term trend component Xlong and the remaining subseasonal residuals Xsub. The STL procedure consists of an inner loop and an outer loop. In the inner loop, seasonal and trend components are estimated from the time series using several passes of smoothing filters. In the outer loop, robustness weights are estimated to reduce the influence of outliers in the time series and are used in the next iteration of the inner loop. This procedure stops once some user-defined stability criterion has been reached. Locally Weighted Regression (Loess) Since the algorithm involves multiple passes of a smoothing filter based on locally weighted regression (Loess), we first present a generic formulation for this smoothing filter. The Loess estimator of x at time tz is denoted gðtzÞ and is given by a weighted polynomial fit of degree d to the values of X that are in the vicinity of tz and are given some weights w. This vicinity is restricted in time by a maximum time-lag parameter k. The weights wi associated with each values of X that are in the vicinity of tz are given by the tricube function (Eq. 11). 8 >< 0"; if jDðti; tzÞj [ k # 3 3 wi ¼ >: 1 jDðti; tzÞj k jDðti; tzÞj k where Dðti; tzÞ is a function of the distance in time between the two points. The values of X closest to tz will thus have the highest weights. Hence, for each prediction point tz, X is associated with a vector of weights W that depends on tz, k and the nature of the distance function D: Wðtz; k; DÞ ¼ fwi; . . .; wng The degree of the polynomial fit gðtzÞ that is fitted at each tz has to be chosen by the user depending on the application (usually a degree 1 or 2 polynomial is sufficient). Once this is performed for every time step in T, we obtain G = fgðt1Þ; . . .; gðtnÞg, the so-called local regression of the original time series. We introduce two different metrics for the distance Dðti; tzÞ between two points in time. The first one is simply an absolute distance in time and is given by: ð11Þ ð12Þ Dabsðti; tzÞ ¼ ti tz ð13Þ The second one is the periodic distance and depends on the length p of the chosen periodicity, for instance, suppose the time vector T is in units of days, then p will be equal to 365.25 days for a seasonal periodicity. The periodic distance is given by: h i p p Dperðti; tzÞ ¼ ti tz þ 2 mod p 2 ð14Þ where mod is the modulo operator. Thus, time points spaced from tz by multiples of exactly p days are assigned a null distance, and then, the distance linearly increases until the time points spaced by multiples of p/2 days are assigned the maximum distance. Note that these distance metrics are used not only for the calculation of weights but also as the input vector of the polynomial fit gðtzÞ to the values of X. Inner Loop The inner loop of the STL procedure aims at estimating the seasonal and long-term components of the time series X through the following steps. Step 1 Compute the detrended time series X* On the very first pass of the inner loop, the long-term trend Xlong will not be known yet so that X* is simply the original time series X. Otherwise, Xlong is subtracted from X, yielding the detrended time series X*. Step 2 Compute the seasonal cycle Xseas The seasonal cycle is estimated at all points in T from X* using Loess with Dper and a free parameter kper that defines the vicinity in terms of periodic distance. The weight vector used in the local polynomial regression is denoted WperðtzÞ ¼ W tz; kper; Dper Note: after the first iteration of the inner and outer loops, Wper is used in combination with the robustness weights Wrobust. The weights simply correspond to the product of Wper and Wper: Wrobust. ð15Þ ð16Þ Step 3 Compute deseasonalized time series XD A deseasonalized time series XD is obtained by subtracting the seasonal cycle Xseas from the original time series X. Step 4 Compute the long-term trend Xlong The long-term trend is estimated after applying Loess to the deseasonalized time series XD with the distance function Dabs and parameter klong which defines the associated weight vector Wlong. WlongðtzÞ ¼ W tz; klong; Dabs As in step 2, these weights are used in combination with the robustness weights Wrobust after the first iteration of the inner and the outer loop. Outer Loop Once an initial run of the inner loop has been carried out, the original time series can be decomposed into X = Xlong ? Xseas ? Xresid, where Xresid corresponds to the residuals. Following Cleveland et al. (1990), extremely large residuals are assumed to correspond to outliers and are assigned a small or a zero weight. These weights are defined using the bisquare function: where h is: 8 >< 0"; wrobust ¼ >: 1 # 2 2 jxresidj [ h jxresidj h h ¼ 6 medianðjXresidjÞ However, a problem arises when time series exhibit seasonal heteroskedasticity because the value of h would change when different seasons of the time series are considered. A typical case of seasonal heteroskedasticity is when precipitation totals are very low during the dry season but exhibit high variability during the wet season. If we followed the approach of Cleveland et al., relatively small outliers in the dry season would be unlikely detected, whereas large but still realistic variations during the wet season would be more likely detected as outliers. This problem can be avoided by introducing seasonally varying estimates of h. This is done by calculating h at each prediction point tz using a weighted median with the weights Wper computed from the periodic distance (step 2): hðtzÞ ¼ 6 median jXresidj; WperðtzÞ Choosing the Parameters The following parameters need to be defined: p the length of each cycle of the seasonal component; d the degree of the weighted polynomial regression; ninner the number of passes through the inner loop; nouter the number of iterations of the outer loop; kper the maximum time lag for the computation of the seasonal component; klong the maximum time lag for the computation of the long-term component Parameter p is obviously determined by the nature of the investigated time series (here, p = 365.25). We chose a polynomial of degree d = 2 in order to hinder smoothing of the low and high peaks of the seasonal cycle. For the trend component, a polynomial of degree d = 1 is sufficient. The number of passes and iterations was chosen so that the resulting decomposition reaches stability. For the number of passes of the inner loop, Cleveland et al. recommend ninner = 2. Regarding nouter, we experimentally determined that nouter = 3 was sufficient for our application. The parameter kper determines the smoothness of the seasonal cycle, with larger values, resulting in a smoother estimate of the seasonal cycle. On the other hand, smaller values of kper will reduce the number of points actually used in the local regression so that the resulting seasonal cycle is more likely to be affected by outliers or sudden changes arising from the uneven spacing of the time series. Hence, the choice of kper is a balanced consideration between accuracy and robustness in the representation of the seasonal cycle. In this paper, a good compromise was experimentally found with kper ¼ 60 days. The parameter klong controls the degree of leakage of the long-term component into the residuals. Larger values of the parameter will result in a smoother estimate of the trend but also cause some of the long-term signal to be incorporated into the residuals. Vice versa, smaller values of this parameter will make the long-term component more sensitive to high-frequency variability. In this application, we followed the recommendations from Cleveland et al. who showed that klong = 1.5 9 p provides a good compromise in most cases. ð17Þ ð18Þ ð19Þ Appendix 2 Analytical Integration of the Weighting Function Integrating Eq. 6 over the interval aj; bj associated with the monthly interval tj must be done with care since the function wðt; tiÞ is discontinuous. In total, three cases can be considered: (1) the continuous case where aj ti (Fig. 4a), (2) the discontinuous case where aj\ti bj (Fig. 4b) and (3) the continuous case where bj\ti. For convenience, we note aj ¼ a and bj ¼ b: The normalized version of W tj ; ti is given by (equivalently to Eq. 8): The denominator has to be decomposed into three continuous parts: wðt; tiÞdt Zb 8a > > > > > > > > < Rti ¼ >>>> a > > > > : 8 > < ¼ >: Za 1 ð20Þ ð21Þ ð23Þ R b R b wðt; tiÞdt R b ti wðt; tiÞdt þ wðt; tiÞdt wðt; tiÞdt a se 1sðb tiÞ þ se 1sða tiÞ se 1sðb tiÞ þ s 0 if if a ti a\ti b\ti a ti b\ti a\ti W^ tj ; ti ¼ R Zb Zþ1 W tj ; ti dti ¼ W tj ; ti dti þ W tj ; ti dti þ W tj ; ti dti ð22Þ W tj ; ti dti ¼ s2e 1sðb aÞ s2e 1sða aÞ s2e 1sðbþ1Þ s2e 1sðaþ1Þ ¼ s2e 1sðb aÞ s2 Za 1 Zþ1 1 The first part yields: The second part yields: And the third part is: Zþ1 1 Zb a Thus after combining equations 23, 24 and 25, equation 22 can be rewritten as: W tj ; ti dti ¼ s2e 1sðb aÞ s2 þ s2 þ sb sa þ 0 ¼ sðb aÞ which is then injected in Eq. 21, yielding the normalized weighting function (Eq. 9): 8 ¼ R e1sðti bÞ b e1sðti bÞ 0 ¼ >> > > > : >>> e 1sðb tiÞ þ e 1sða tiÞ > > < e 1sðb tiÞ þ 1 b 0 a if a ti a\ti if b\ti if a ti a\ti if b\ti s2e 1sðb bÞ þ sb s2e 1sðb aÞ þ sa ¼ s2 þ sb sa W tj ; ti dti ¼ 0 Abelen S, Seitz F (2013) Relating satellite gravimetry data to global soil moisture products via data harmonization and correlation analysis. Remote Sens Environ 136:89–98. doi:10.1016/j.rse.2013.04. 012 Abelen S, Seitz F, Abarca-del-Rio R, Gu¨ntner A (2015) Droughts and floods in the La Plata basin in soil moisture data and GRACE. Remote Sens 7:7324–7349. doi:10.3390/rs70607324 Ahmed M et al (2011) Integration of GRACE (gravity recovery and climate experiment) data with traditional data sets for a better understanding of the time-dependent water partitioning in African watersheds. Geology 39:479–482. doi:10.1130/G31812.1 Ahmed M, Sultan M, Wahr J, Yan E (2014) The use of GRACE data to monitor natural and anthropogenic induced variations in water availability across Africa. Earth-Sci Rev 136:289–300. doi:10.1016/j. earscirev.2014.05.009 Andersen OB, Seneviratne SI, Hinderer J, Viterbo P (2005) GRACE-derived terrestrial water storage depletion associated with the 2003 European heat wave. Geophys Res Lett 32:L18405. doi:10.1029/ 2005GL023574 Arendt A, Luthcke S, Gardner A, O’Neel S, Hill D, Moholdt G, Abdalati W (2013) Analysis of a GRACE global mascon solution for Gulf of Alaska glaciers. J Glaciol 59:913–924. doi:10.3189/2013JoG12J197 ð25Þ ð26Þ ð27Þ Dufresne JL et al (2013) Climate change projections using the IPSL-CM5 earth system model: from CMIP3 to CMIP5. Climate Dyn 40:2123–2165. doi:10.1007/s00382-012-1636-1 Eicker A, Schumacher M, Kusche J, Doll P, Muller Schmied H (2014) Calibration/data assimilation approach for integrating GRACE data into the water GAP global hydrology model (WGHM) using an ensemble Kalman filter: first results. Surv Geophys 35:1285–1309. doi:10.1007/s10712-014-9309-8 Famiglietti JS (2014) The global groundwater crisis. Nat Climate Change 4:945–948. doi:10.1038/ nclimate2425 Famiglietti JS et al (2011) Satellites measure recent rates of groundwater depletion in California’s Central Valley. Geophys Res Lett 38:L03403. doi:10.1029/2010GL046442 Feng W, Zhong M, Lemoine J-M, Biancale R, Hsu H-T, Xia J (2013) Evaluation of groundwater depletion in North China using the gravity recovery and climate experiment (GRACE) data and ground-based measurements. Water Resour Res 49:2110–2118. doi:10.1002/wrcr.20192 Forootan E, Kusche J (2012) Separation of global time-variable gravity signals into maximally independent components. J Geodesy 86:477–497. doi:10.1007/s00190-011-0532-5 Forootan E et al (2014a) Multivariate prediction of total water storage changes over West Africa from multisatellite data. Surv Geophys 35:913–940. doi:10.1007/s10712-014-9292-0 Forootan E et al (2014b) Separation of large scale water storage patterns over Iran using GRACE, altimetry and hydrological data. Remote Sens Environ 140:580–595. doi:10.1016/j.rse.2013.09.025 Frappart F, Ramillien G (2012) Contribution of GRACE satellite gravimetry in global and regional hydrology, and in ice sheets mass balance, water resources management and modeling. In: Purna N (ed) Water resources management and modeling. InTech, p 322. doi:10.5772/34212 Frappart F, Ramillien G, Biancamaria S, Mognard NM, Cazenave A (2006) Evolution of high-latitude snow mass derived from the GRACE gravimetry mission (2002–2004). Geophys Res Lett 33:L02501. doi:10.1029/2005GL024778 Frappart F, Ramillien G, Famiglietti JS (2011a) Water balance of the Arctic drainage system using GRACE gravimetry products. Int J Remote Sens 32:431–453. doi:10.1080/01431160903474954 Frappart F, Ramillien G, Leblanc M, Tweed SO, Bonnet M-P, Maisongrande P (2011b) An independent component analysis filtering approach for estimating continental hydrology in the GRACE gravity data. Remote Sens Environ 115:187–204. doi:10.1016/j.rse.2010.08.017 Frappart F, Papa F, da Silva JS, Ramillien G, Prigent C, Seyler F, Calmant S (2012) Surface freshwater storage and dynamics in the Amazon basin during the 2005 exceptional drought. Environ Res Lett 7:044010. doi:10.1088/1748-9326/7/4/044010 Frappart F, Ramillien G, Ronchail J (2013) Changes in terrestrial water storage versus rainfall and discharges in the Amazon basin. Int J Climatol 33:3029–3046. doi:10.1002/joc.3647 Garc´ıa-Garc´ıa D, Ummenhofer CC, Zlotnicki V (2011) Australian water mass variations from GRACE data linked to Indo-Pacific climate variability. Remote Sens Environ 115:2175–2183. doi:10.1016/j.rse.2011.04.007 Gardner A et al (2013) A reconciled estimate of glacier contributions to Sea level rise: 2003 to 2009. Science 340:852–857. doi:10.1126/science.1234532 Geruo A, Wahr J, Zhong S (2013) Computations of the viscoelastic response of a 3-D compressible Earth to surface loading: an application to glacial isostatic adjustment in Antarctica and Canada. Geophys J Int 192:557–572. doi:10.1093/gji/ggs030 Gudmundsson L, Seneviratne SI (2015) European drought trends. Proc Int Assoc Hydrol Sci 369:75–79. doi:10.5194/piahs-369-75-2015 Gudmundsson L, Tallaksen LM, Stahl K, Fleig AK (2011) Low-frequency variability of European runoff. Hydrol Earth Syst Sci 15:2853–2869. doi:10.5194/hess-15-2853-2011 Gu¨ntner A (2008) Improvement of global hydrological models using GRACE data. Surv Geophys 29:375–397. doi:10.1007/s10712-008-9038-y Gu¨ntner A, Schmidt R, Do¨ll P (2007a) Supporting large-scale hydrogeological monitoring and modelling by time-variable gravity data. Hydrogeol J 15:167–170. doi:10.1007/s10040-006-0089-1 Gu¨ntner A, Stuck J, Werth S, Do¨ll P, Verzano K, Merz B (2007b) A global analysis of temporal and spatial variations in continental water storage. Water Resour Res 43:W05416. doi:10.1029/2006WR005247 Hamed KH (2009) Exact distribution of the Mann-Kendall trend test statistic for persistent data. J Hydrol 365:86–94. doi:10.1016/j.jhydrol.2008.11.024 Hamed KH, Rao AR (1998) A modified Mann–Kendall trend test for autocorrelated data. J Hydrol 204:182–196. doi:10.1016/S0022-1694(97)00125-X Han S-C, Shum CK, Bevis M, Ji C, Kuo C-Y (2006) Crustal dilatation observed by GRACE after the 2004 Sumatra–Andaman earthquake. Science 313:658–662. doi:10.1126/science.1128661 Han S-C, Sauber J, Riva R (2011) Contribution of satellite gravimetry to understanding seismic source processes of the 2011 Tohoku-Oki earthquake. Geophys Res Lett 38:L24312. doi:10.1029/2011GL049975 Ramillien G, Frappart F, Seoane L (2014) Application of the regional water mass variations from GRACE satellite gravimetry to large-scale water management in Africa. Remote Sens 6:7379–7405. doi:10. 3390/rs6087379 Rangelova E, van der Wal W, Braun A, Sideris MG, Wu P (2007) Analysis of gravity recovery and climate experiment time-variable mass redistribution signals over North America by means of principal component analysis. J Geophys Res-Earth 112:F03002. doi:10.1029/2006JF000615 Rangelova E, van der Wal W, Sideris MG, Wu P (2010) Spatiotemporal analysis of the GRACE-derived mass variations in North America by means of multi-channel singular spectrum analysis. In: Mertikas SP (ed) Gravity, geoid and earth observation, vol 135. International association of geodesy symposia. Springer Berlin Heidelberg, pp 539–546. doi:10.1007/978-3-642-10634-7_72 Reager JT, Thomas BF, Famiglietti JS (2014) River basin flood potential inferred using GRACE gravity observations at several months lead time. Nat Geosci 7:588–592. doi:10.1038/ngeo2203 Richey AS, Thomas BF, Lo M-H, Famiglietti JS, Swenson S, Rodell M (2015a) Uncertainty in global groundwater storage estimates in a total groundwater stress framework. Water Resour Res 51:5198–5216. doi:10.1002/2015WR017351 Richey AS et al (2015b) Quantifying renewable groundwater stress with GRACE. Water Resour Res 51:5217–5238. doi:10.1002/2015WR017349 Rieser D, Kuhn M, Pail R, Anjasmara IM, Awange J (2010) Relation between GRACE-derived surface mass variations and precipitation over Australia. Aust J Earth Sci 57:887–900. doi:10.1080/08120099.2010. 512645 Rodell M, Famiglietti JS (2002) The potential for satellite-based monitoring of groundwater storage changes using GRACE: the High Plains aquifer, Central US. J Hydrol 263:245–256. doi:10.1016/S00221694(02)00060-4 Rodell M, Famiglietti JS, Chen J, Seneviratne SI, Viterbo P, Holl S, Wilson CR (2004) Basin scale estimates of evapotranspiration using GRACE and other observations. Geophys Res Lett 31:L20504. doi:10. 1029/2004GL020873 Rodell M, Chen J, Kato H, Famiglietti JS, Nigro J, Wilson CR (2007) Estimating groundwater storage changes in the Mississippi River basin (USA) using GRACE. Hydrogeol J 15:159–166. doi:10.1007/ s10040-006-0103-7 Rodell M, Velicogna I, Famiglietti JS (2009) Satellite-based estimates of groundwater depletion in India. Nature 460:999–1002. doi:10.1038/nature08238 Sakumura C, Bettadpur S, Bruinsma S (2014) Ensemble prediction and intercomparison analysis of GRACE time-variable gravity field models. Geophys Res Lett 41:1389–1397. doi:10.1002/2013GL058632 Sasgen I et al (2012) Timing and origin of recent regional ice-mass loss in Greenland. Earth Planet Sci Lett 333–334:293–303. doi:10.1016/j.epsl.2012.03.033 Schmeer M, Schmidt M, Bosch W, Seitz F (2012) Separation of mass signals within GRACE monthly gravity field models by means of empirical orthogonal functions. J Geodyn 59–60:124–132. doi:10. 1016/j.jog.2012.03.001 Schmidt R et al (2006) GRACE observations of changes in continental water storage. Global Planet Change 50:112–126. doi:10.1016/j.gloplacha.2004.11.018 Schmidt R, Flechtner F, Meyer U, Neumayer KH, Dahle C, Ko¨nig R, Kusche J (2008a) Hydrological signals observed by the GRACE satellites. Surv Geophys 29:319–334. doi:10.1007/s10712-008-9033-3 Schmidt R, Petrovic S, Gu¨ntner A, Barthelmes F, Wu¨nsch J, Kusche J (2008b) Periodic components of water storage changes from GRACE and global hydrology models. J Geophys Res-Solid Earth 113:B08419. doi:10.1029/2007JB005363 Schrama EJO, Wouters B, Lavalle´e DA (2007) Signal and noise in gravity recovery and climate experiment (GRACE) observed surface mass variations. J Geophys Res-Solid Earth 112:B08407. doi:10.1029/ 2006JB004882 Seitz F, Schmidt M, Shum CK (2008) Signals of extreme weather conditions in Central Europe in GRACE 4-D hydrological mass variations. Earth Planet Sci Lett 268:165–170. doi:10.1016/j.epsl.2008.01.001 Sen PK (1968) Estimates of the regression coefficient based on Kendall’s Tau. J Am Stat Assoc 63:1379–1389. doi:10.2307/2285891 Seneviratne SI et al (2010) Investigating soil moisture–climate interactions in a changing climate: a review. Earth-Sci Rev 99:125–161. doi:10.1016/j.earscirev.2010.02.004 Seo KW, Wilson CR, Famiglietti JS, Chen JL, Rodell M (2006) Terrestrial water mass load changes from gravity recovery and climate experiment (GRACE). Water Resour Res 42:W05417. doi:10.1029/ 2005WR004255 Simmons AJ, Willett KM, Jones PD, Thorne PW, Dee DP (2010) Low-frequency variations in surface atmospheric humidity, temperature, and precipitation: inferences from reanalyses and monthly gridded observational data sets. J Geophys Res-Atmos 115:D01110. doi:10.1029/2009JD012442 Wouters B, Bonin JA, Chambers DP, Riva REM, Sasgen I, Wahr J (2014) GRACE, time-varying gravity, Earth system dynamics and climate change. Rep Prog Phys. doi:10.1088/0034-4885/77/11/116801 Wu X et al (2010) Simultaneous estimation of global present-day water transport and glacial isostatic adjustment. Nat Geosci 3:642–646. doi:10.1038/ngeo938 Xavier L, Becker M, Cazenave A, Longuevergne L, Llovel W, Filho OCR (2010) Interannual variability in water storage over 2003–2008 in the Amazon Basin from GRACE space gravimetry, in situ river level and precipitation data. Remote Sens Environ 114:1629–1637. doi:10.1016/j.rse.2010.02.005 Yue S, Wang C (2004) The Mann–Kendall test modified by effective sample size to detect trend in serially correlated hydrological series. Water Resour Manag 18:201–218. doi:10.1023/B:WARM.0000043140. 61082.60 Yue S, Pilon P, Phinney B, Cavadias G (2002) The influence of autocorrelation on the ability to detect trend in hydrological series. Hydrol Process 16:1807–1829. doi:10.1002/hyp.1095 Zaitchik B, Rodell M, Reichle R (2008) Assimilation of GRACE terrestrial water storage data into a land surface model: results for the Mississippi River basin. J Hydrometeorol 9:535–548. doi:10.1175/ 2007JHM951.1


This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1007%2Fs10712-016-9367-1.pdf

Vincent Humphrey, Lukas Gudmundsson, Sonia I. Seneviratne. Assessing Global Water Storage Variability from GRACE: Trends, Seasonal Cycle, Subseasonal Anomalies and Extremes, Surveys in Geophysics, 2016, 357-395, DOI: 10.1007/s10712-016-9367-1