A global water cycle reanalysis ( 2003 – 2012 ) merging satellite gravimetry and altimetry observations with a hydrological multi-model ensemble

We present a global water cycle reanalysis that merges water balance estimates derived from the Gravity Recovery And Climate Experiment (GRACE) satellite mission, satellite water level altimetry and off-line estimates from several hydrological models. Error estimates for the sequential data assimilation scheme were derived from available uncertainty information and the triple collocation technique. Errors in four GRACE storage products were estimated to be 11–12 mm over land areas, while errors in monthly storage changes derived from five global hydrological models were estimated to be 17–28 mm. Prior and posterior water storage estimates were evaluated against independent observations of river water level and discharge, snow water storage and glacier mass loss. Data assimilation improved or maintained agreement overall, although results varied regionally. Uncertainties were greatest in regions where glacier mass loss and subsurface storage decline are both plausible but poorly constrained. We calculated a global water budget for 2003– 2012. The main changes were a net loss of polar ice caps (−342 Gt yr−1) and mountain glaciers ( −230 Gt yr−1), with an additional decrease in seasonal snowpack ( −18 Gt yr−1). Storage increased due to new impoundments ( +16 Gt yr−1), but this was compensated by decreases in other surface water bodies ( −10 Gt yr−1). If the effect of groundwater depletion (−92 Gt yr−1) is considered separately, subsurface water storage increased by +202 Gt yr−1 due particularly to increased wetness in northern temperate regions and in the seasonally wet tropics of South America and southern Africa. The reanalysis results are publicly available via www.wenfo.org/wald/ .


Introduction
More accurate global water balance estimates are needed, to better understand interactions between the global climate system and water cycle (Sheffield et al., 2012), the causes of observed sea level rise (Boening et al., 2012;Fasullo et al., 2013;Cazenave et al., 2009;Leuliette and Miller, 2009) and human impacts on water resources (Wada et al., 2010(Wada et al., , 2013)), and to improve hydrological models (van Dijk et al., 2011) and initialise water resources forecasts (Van Dijk et al., 2013).The current generation of global hydrological models have large uncertainties arising from a combination of data deficiencies (e.g.precipitation in sparsely gauged regions; poorly known soil, aquifer and vegetation properties) and overly simplistic descriptions of important water cycle processes (e.g.groundwater dynamics, human water resources extraction and use, wetland hydrology and glacier dynamics).Data assimilation is used routinely to overcome data and model limitations in atmospheric reconstructions or "reanalysis".In hydrological applications, there has been an (over-)emphasis on parameter calibration (Van Dijk, 2011), with data assimilation approaches largely limited to flood forecasting.New applications are being developed, however (Liu et al., 2012a), including promising developments towards large-scale water balance reanalyses, alternatively referred to as monitoring, assessment or estimation (van Dijk and Renzullo, 2011).
Here, we undertake a global water cycle reanalysis for the period 2003-2012.Specifically, we attempt to merge global water balance estimates from different model sources with an ensemble of total water storage (TWS) estimates derived from the Gravity Recovery And Climate Experiment (GRACE) satellite mission (Tapley et al., 2004).Various alternative approaches can be conceptualised to achieve this integration, and the most appropriate among these is not obvious.Our approach was to use water balance estimates generated by five global hydrological models along with several ancillary data sources to generate an ensemble of prior estimates of monthly water storage changes.Errors in the different model estimates and GRACE products were estimated spatially through triple collocation (Stoffelen, 1998).Subsequently, a data assimilation scheme was designed to sequentially merge the model ensemble and GRACE observations.The reanalysis results were evaluated with independent global streamflow records, remote sensing of river water level and snow water equivalent (SWE), and independent glacier mass balance estimates.
2 Methods and data sources

Overall approach
We conceptualise TWS (S, in mm) as the sum of five different water stores (s in mm), i.e. water stored in snow and ice (s snow ), below the surface in soil and groundwater (s sub ) and in rivers (s riv ), lakes (s lake ) and seas and oceans (s sea ).We ignore atmospheric water storage changes, which are removed from the signal during the GRACE TWS retrieval process (e.g.Wahr et al., 2006), and vegetation mass changes, which are assumed negligible.The GRACE TWS estimates are denoted by y and have the same units as S but are distinct in their much smoother spatial character.
To date, data assimilation schemes developed for largescale water cycle analysis typically use Kalman filter approaches (Liu et al., 2012a).This requires calculation of covariance matrices and, presumably because of complexity and computational burden, has only been applied for single models and limited regions (e.g.Zaitchik et al., 2008).We aimed to develop a data assimilation scheme that made it possible to use water balance estimates derived "off-line" (i.e. in the absence of data assimilation) so we could use an ensemble of already available model outputs.In the data assimilation terminology of Bouttier and Courtier (1999), our scheme could be described as sequential and near-continuous with a spatially variable but temporally stable gain factor.The characteristics of the data assimilation problem to be addressed in this application were as follows: 1. alternative GRACE TWS estimates (y o ) were available from different processing centres and error estimates were required for each; 2. alternative estimates for some of the stores, s, were available from different hydrological models, with higher definition than y o ; 3. error estimates were required for each store and data source; 4. a method was required to spatially transform between s and y as part of the assimilation.

Data sources
The data used include those needed to derive prior estimates for each of the water cycle stores, the GRACE retrievals to be assimilated and independent observations to evaluate the quality of the reanalysis.All are listed in Table 1 and described below.
Monthly water balance components from four global land surface model estimates at 1 • resolution were obtained from NASA's Global Data Assimilation System (GLDAS) (Rodell et al., 2004).The four models include CLM (Community Land Model), Mosaic, NOAH and VIC (Variable Infiltration Capacity) which, for 2003-2012, were forced with "a combination of NOAA/GDAS atmospheric analysis fields, spatially and temporally disaggregated NOAA Climate Prediction Center Merged Analysis of Precipitation (CMAP) fields, and observation-based radiation fields derived using the method of the Air Force Weather Agency's AGRicultural METeorological modelling system" (Rui, 2011).The models are described in Rodell et al. (2004).From the model outputs we used (i) SWE depth, (ii) total soil moisture storage over a soil depth that varies between models and (iii) generated streamflow, calculated as the sum of surface runoff and subsurface drainage.In addition to GLDAS, we used global water balance estimates generated by the W3RA (World-Wide Water Resources Assessment) model (Van Dijk et al., 2013) in the configuration used in the Asia-Pacific Water Monitor (http://www.wenfo.org/apwm/).For 2003-2008, the model was forced with the "Princeton" merged precipitation, down-welling short-wave radiation, minimum and maximum daily temperature and air pressure data produced by Sheffield et al. (2006).From 2009 onwards, the model primarily uses "ERA-Interim" weather forecast model reanalysis data from the European Centre for Medium-Range Weather Forecasts.For low latitudes, these are combined with near-real-time TRMM Multi-sensor Precipitation Analysis data (TMPA code 3B42 RT) (Huffman et al., 2007) to improve estimates of convective rainfall (Peña-Arancibia et al., 2013).Both were bias-corrected with reference to the Princeton data to ensure homogeneity.W3RA model estimates were conceptually similar to those from GLDAS, except that the model includes deep soil and groundwater stores and sub-grid surface and groundwater routing.
The five hydrological models do not provide estimates of groundwater depletion and storage in rivers, lakes and impoundments; these were therefore derived separately.Wada et al. (2012).The time series were calculated as the net difference between estimated groundwater extraction and recharge.National groundwater extraction data compiled by the International Groundwater Resources Assessment Centre (IGRAC) were disaggregated using estimates of water use intensity and surface water availability at 0.5 • resolution from a hydrological model (PCR-GLOBWB; see Wada et al., 2012, for details).The model also estimated recharge including return flow from irrigation.Groundwater depletion uncertainty estimates were generated through 10 000 Monte Carlo simulations, with 100 realisations of both extraction and recharge (Wada et al., 2010).This method tends to overestimate reported depletion in non-arid regions, where groundwater pumping can enhance recharge from surface water.Wada et al. (2012) used a universal multiplicative correction of 0.75 to account for this.Here, the correction was calculated per climate region rather than worldwide, reflecting the dependency of uncertainty on recharge estimates and their errors, and resulting in values of 0.6 to 0.9.Depletion estimates for 2011-2012 were not available; these were estimated using monthly average depletion and uncertainty values for the preceding 2003-2010 period.Given the regular pattern of depletion in the preceding years this by itself is unlikely to have affected the analysis noticeably.
River water storage was estimated by propagating runoff fields from each of the five models through a global routing scheme.In a previous study, we compared these runoff fields with streamflow records from 6192 small (<10 000 km 2 ) catchments worldwide and found that observed runoff was 1.28 to 1.77 times greater than predicted by the different models (Van Dijk et al., 2013).The respective ratios were used to uniformly bias-correct the runoff fields.Next, we used a global 0.5 • resolution flow direction grid (Oki et al., 1999;Oki and Sud, 1998) to parameterise a cell-to-cell river routing scheme.We used a linear reservoir kinematic wave approximation (Vörösmarty and Moore, 1991), similar to that used in several large-scale hydrology models (see recent review by Gong et al., 2011).The monthly 1 • runoff fields from each of the five models were oversampled to 0.5 • and daily time step before routing, and the river water storage estimates (in mm) were aggregated back to monthly 1 was an inverse linear function of the distance between network nodes and a transfer (or routing) coefficient.For each model, a globally uniform optimal transfer coefficient was found by testing values of 0.3 to 0.9 day −1 in 0.1 day −1 increments and finding the value that produced best overall agreement with seasonal flow patterns observed in 586 large rivers worldwide.These 586 were a subset of 925 oceanreaching rivers for which streamflow records were compiled by Dai et al. (2009) from various sources.We excluded locations where streamflow records were available for less than 10 years since 1980 or less than 6 months of the year.
The resulting river flow estimates do not account for the impact of river water use (i.e. the evaporation of water extracted from rivers, mainly for irrigation).We addressed this using global monthly surface water use estimates that were derived in a way similar to that used for groundwater depletion estimates (full details in Wada et al., 2013).For each grid cell, mean water use rates for 2002-2010 were subtracted from mean runoff estimates for the same period, and the remaining runoff was routed downstream.The resulting mean net river flow estimates were divided by the original estimates to derive a scaling factor, which was subsequently applied at each time step.Lack of additional global information on river hydrology meant that three simplifications needed to be made: (i) our approach implies that for a particular grid cell, monthly river water use is assumed proportional to river flow for that month; (ii) the influence of lakes, wetlands and water storages on downstream flows (e.g. through dam operation) is not accounted for, even though their actual storage changes are (see further on); and (iii) our approach does not account for losses associated with permanent or ephemeral wetlands, channel leakage and net evaporation from the river channel.At least in theory, data assimilation may correct mass errors resulting from these assumptions.
Variations in lake water storage were not modelled, but water level data for 62 lakes worldwide were obtained from the Crop Explorer website (Table 1) and include most of the world's largest lakes and reservoirs, including the Caspian Sea.The water level data for these lakes were derived from satellite altimetry and converted to mm water storage.Measurements were typically available every 10 days.The mean and standard deviation (SD) of measurements in each month were used as respectively best estimate and estimation error for that month.Storage in water bodies without altimetry data was necessarily assumed negligible.This includes many small lakes and dams, but also some larger lakes affected by snow and ice cover (e.g. the Great Bear and Great Slave lakes in Canada) and ephemeral, distributed or otherwise complex water bodies (e.g. the Okavango Delta in Botswana and Lake Eyre in Australia, each of which contains > 10 km 3 of water when full).
New river impoundments lead to persistent water storage increases.A list of dams was collated by Lehner et al. (2011) and was updated with large dams constructed in more recent years with the ICOLD data base (Table 1).For the pe-riod 1998-2012, a total 198 georeferenced dams with a combined storage capacity of 418 km 3 were identified.For the Three Gorges Dam (39 km 3 ), reservoir water level time series (http://www.ctg.com.cn/inc/sqsk.php)were converted to storage volume following Wang et al. (2011).For the remaining dams, we assumed a gradual increase to storage capacity over the first 5 years after construction and assumed a relative estimation error of 20 %.The combined annual storage increase amounted to 21 km 3 yr −1 on average.
Global merged mean sea level anomalies were obtained from the Aviso website (Table 1).The monthly data were reprojected from the native 1/3 • Mercator grid to regular 1 • grids.An estimate of uncertainty was derived by calculating the spatial standard deviation in sea level values within a 4 • by 4 • region around each grid cell during re-projection.When sea level data were missing because of sea ice, we assumed sea level did not change and assigned an uncertainty of 5 mm.Following the recent global sea level budget study by Chen et al. (2013), we assumed that 75 % of the observed sea level change was due to mass increase, and multiplied altimetry sea level anomalies with this factor.
We did not have spatial global time series of glacier mass changes.The five hydrological models have poor representation of ice dynamics, and therefore large uncertainties and errors can be expected for glaciated regions.To account for this, we used the "GGHYDRO" global glacier extent mapping by Cogley (2003) to calculate the percentage glacier area for each grid cell, and assumed a proportional error in monthly glacier mass change estimates corresponding to 300 mm per unit glacier area.This value was chosen somewhat arbitrarily but ensures that a substantial fraction of the regional analysis increment is assigned to glaciers.
Three alternative GRACE TWS retrieval products were downloaded from the Tellus website.The three products (coded CSR, JPL and GFZ; release 05) each had a nominal 1 • and monthly resolution.The land and ocean mass retrievals (Chambers and Bonin, 2012) were combined.The land retrievals had been "de-striped" and smoothed with a 200 km half-width spherical Gaussian filter (Swenson et al., 2008;Swenson and Wahr, 2006), whereas the ocean retrievals had been smoothed with a 500 km filter (Chambers and Bonin, 2012).The data assimilation method we employed is designed to deal with the signal "leakage" caused by the smoothing process, and therefore we did not use the scaling factors provided by the algorithm developers.In addition, gravity fields produced by the Centre National d'Etudes Spatiales (CNES) Groupe de Recherches de Géodésie Spatiale (GRGS) (Bruinsma et al., 2010) at 1 • resolution for 10-day periods were used.The three Tellus data sources had been corrected for glacial isostatic adjustment (GIA); we corrected the GRGS data using the same GIA estimates of Geruo et al. (2013).Initial data assimilation experiments produced unexpectedly strong mass trends around the Gulf of Thailand.Inspection demonstrated that all products, to different degrees, contained a mass redistribution signal associated with the December 2004 Sumatra-Andaman earthquake.To account for this, we first calculated a time series of seasonally adjusted monthly anomalies (i.e. the average seasonal cycle was removed) for the region [5 • N-15 • , 80-110 • E].Next, we adjusted values after December 2004 by the difference in the mean adjusted anomalies for the year before and after the earthquake.

Data assimilation scheme
For each update cycle, the data assimilation scheme proceeds through the steps illustrated in Fig. 1 and described below.1. Deriving the prior estimate for each store.The way to calculate the prior (or background) estimate of storage s b t varied between stores.A systematic and accumulating bias (or "drift") was considered plausible for the deep soil and groundwater components of modelderived subsurface storage due to slow groundwater dynamics (including extraction) and ice storage in permanent glaciers and ice sheets, which may be progressively melting or accumulating.In these cases, the modelestimated change in storage was assumed more reliable than the model-estimated storage itself, and estimates from the five models were used to calculate stor-age change, s b t , for store i (i = 1, . . ., N ) as where x l t is the estimate of storage change from model l (l = 1, . . ., L) between time t − 1 and t, and w l the relative weight of model l in the ensemble, computed as where σ l is the estimated local error for model l based on triple collocation (see Sect. 2.4).Subsequently, s b t was calculated as where s a * t−1 is the posterior (or analysis) estimate from the previous time step.This approach was not suitable for model-estimated seasonal snowpack and river storage, where the ephemeral nature of the storage means that long-term drift is not an issue and Eq. ( 2) could in fact lead to unrealistic negative storage values.For these cases, s b t was computed as where s l t is the storage estimate from model l.The glacier extent map was used to identify whether Eqs. ( 3) or (4) should be used for s snow .Similarly, no drift was expected in the ocean and lake storage data, and these were used directly as estimates of s b t .

Deriving the prior estimate of GRACE-like TWS (y b ).
This estimate was derived by summing all stores s b t as and subsequently applying a convolution operator to transform S b t to a "GRACE-like" TWS y b .The operator was a Gaussian smoother (cf.Jekeli, 1981), written here as where j 1 and j 2 in principle should encompass all existing grid cell coordinates.In practice, was applied as a moving Gaussian kernel with a size of 6 • × 6 • and a half-width of 300 km (see further on).

A. I. J. M. van Dijk et al.: A global water cycle reanalysis (2003-2012)
3. Updating the GRACE-like TWS.The updated GRACElike TWS, y a t , was calculated from the prior (Eq.6) and GRACE observations y o t for time t as (cf.Fig. 1a-d): where δy t is the analysis increment and k a temporally static gain factor derived by combining the error variances of modelled and observed y as follows: where w y,l and w y,m are the weights applied to each of the five GRACE-like TWS estimates and four GRACE data sources, respectively, calculated from their respective error variances σ 2 y,l and σ 2 y,m analogous to Eq. ( 2).
4. Spatially disaggregating the analysis increment to the different stores.The observation model was inverted and combined with the store error estimates in order to spatially redistribute the analysis increment δy t , as follows (cf.Fig. 1e-g): where the redistribution operator can be written as (cf.Fig. 1g) To implement this, spatial error estimates are required for each store.For lakes and seas, the errors were estimated from the observations (see Sect. 2.2).For the model-based estimates, the error was calculated for each time step and store as The resulting error estimates are spatially and temporally dynamic and respond to the magnitude of the differences between the different model estimates.For s sub and s snow we combined the error estimates derived by Eq. ( 11) with the estimated errors in groundwater depletion and glacier mass change, respectively (see Sect. 2.2), calculating total error as the quadratic sum of the composite errors.
5. Updating the stores.In the final step, the state of each store is updated: Subsequently, the procedure is repeated for the next time step.

Error estimation
Spatial error fields are required for all data sets to calculate the gain factor k.Where necessary these were estimated using the triple collocation technique (Stoffelen, 1998).This technique infers errors in three independent time series by analysing the covariance structure.The approach has been applied widely to estimate errors in, among others, satellitederived surface soil moisture (Dorigo et al., 2010;Scipal et al., 2009), evapotranspiration (Miralles et al., 2011) and vegetation leaf area (Fang et al., 2012).A useful description of the technique, the assumptions underlying it and an extension of the theory to more than three time series is provided by Zwieback et al. (2012).Application requires three (or more) estimates of the same quantity.This was achieved by convolving the model-derived storage estimates into large-scale, smoothed TWS estimates equivalent to those derived from GRACE measurements using Eqs.( 5) and ( 6).Inspection of the original Tellus data made clear that the 200 km filter that was already applied as part of the land retrieval had only removed part of the spurious aliasing in the data sets, and propagated these artefacts into the error estimates and reanalysis.Therefore a smoother, 300 km filter was applied to the Tellus TWS data sets.Because conceptual consistency is required for triple collocation, the same filter was applied to the GRGS and model-derived TWS estimates.Several alternative Tellus and model time series were available, and therefore the triple collocation technique could be used to produce alternative error estimates from multiple triplet combinations (i.e.five for Tellus TWS, three for model TWS and 5 × 3 = 15 for GRGS TWS).The agreement between these alternative estimates was calculated as a measure of uncertainty in the estimated errors.Important assumptions of the collocation technique are that (1) all data sets are free of bias relative to each other, (2) errors do not vary over time, (3) there is no temporal autocorrelation in the errors and (4) there is no correlation between the errors in the respective time series (Zwieback et al., 2012).Each of these assumptions is difficult to ascertain, but some interpretative points can be made.First, errors in the GRACE products vary somewhat from month to month depending on data availability, and overall decreased after June 2003.Therefore assumption (2) is a simplification.
Assumption (3) is also unlikely to hold fully for the TWS estimates themselves: there will almost certainly be systematic errors and biases that cause temporal correlation in the errors in the modelled TWS (e.g.due to poorly represented processes causing secular trends such as groundwater extraction or glacier melt).We were able to avoid this assumption by applying the triple collocation to monthly storage changes rather than the actual value of storage, although temporal correlation in storage change errors remains a possibility.Temporal correlation in the GRACE errors is unlikely, however.Therefore, the error in individual monthly mass estimates Hydrol.Earth Syst.Sci., 18, 2955-2973, 2014 www.hydrol-earth-syst-sci.net/18/2955/2014/ was calculated following conventional error propagation theory by dividing the estimated error in mass changes by √ 2. Assumption (4) will not be fully met where estimates are partially based on the same principle or measurement.In this study, arguably the most uncertain assumption is that the GRGS and Tellus errors are to a large extent uncorrelated.The basis for this assumption is that most of the error is likely to derive from the TWS retrieval method rather than the primary measurements (Sakumura et al., 2014).The GRGS time series was selected as the third triple collocation member because the four Tellus products are retrieved by methods that are comparatively more similar than the GRGS method, which uses ancillary observations from the Laser Geodynamics Satellites (Tregoning et al., 2012).Correspondingly, global average temporal correlation among the Tellus TWS time series was stronger (0.61-0.73) than between GRGS and any of the Tellus time series (0.49-0.58).Nonetheless, there may well have been a residual covariance between errors in the GRGS and Tellus products.In triple collocation and subsequent data assimilation, this would cause some part of the differences to be wrongly attributed to the prior estimates rather than the observation products.Therefore, we conservatively inflated the calculated value by including an additional error of 5 mm through quadratic summation before calculating the gain factor (Eq. 8).
Uncertainty in the derived error estimates also arises from sample size, i.e. the number of collocated observations (N = 111).Previous studies have suggested that 100 samples are sufficient to produce a reasonable estimate (Dorigo et al., 2010), although Zwieback et al. (2012) calculate that the relative uncertainty in the estimated errors for N = 111 can be expected to be of the order of 20 %.An uncertainty of this magnitude will not have a strong impact on the reanalysis results.

Evaluation against observations
Evaluation of the reanalysis results for subsurface storage was a challenge: ground observations are not widely available at the global scale, are not conceptually equivalent to the reanalysis terms, require tenuous scaling assumptions for comparison at 1 • grid cell resolution, and many existing data sets contain few or no records during 2003-2012.For example, comparison with in situ soil moisture measurements or groundwater bore data is beset by such problems (Tregoning et al., 2012).Similarly, an initial comparison with nearsurface (< 5 cm depth) soil moisture estimates from passive and active microwave remote sensing (Liu et al., 2012b(Liu et al., , 2011) ) showed that the conceptual difference between the two quantities was too great for any meaningful comparison.
We were able to evaluate the reanalysis for storage in rivers, seasonal snowpack and glaciers, however.Firstly, a total of 1264 water level time series for several large rivers worldwide were obtained from the Laboratoire d'Etudes en Geodésie et Océanographie Spatiales (LEGOS) HYDROWEB website (Table 1).The river levels were retrieved from ENVISAT and Jason-2 satellite altimetry (Crétaux et al., 2011) and included uncertainty information for each data period.From each time series, we removed data points with an estimated error of more than 25 % of the temporal SD.Another 165 altimetry time series were obtained from the European Space Agency (ESA) River&Lake website (Berry, 2009).These were selected to increase measurement period and sample size for the available locations, as well as extending coverage to additional rivers.The ESA time series did not include error estimates; instead data plots were judged visually to assess the likelihood of measurement noise; seemingly affected time series and outlier data points (> 3SD) were excluded.The total 1429 time series were merged for individual 1 • grid cells.In each case, the longest time series was chosen as a reference.Overlapping time periods were used to remove (typically small) systematic biases in water surface elevation between time series; where there was no overlap the time series were normalised by the median water level.The ESA data were used where or when HYDROWEB data were not available, and merged time series with fewer than 24 data points in total were excluded.The resulting data set contained time series for 442 grid cells with an average 61 (maximum 115) data points during 2003-2012.The relationship between river water level and river discharge (i.e. the discharge rating curve) was unknown, and therefore a direct comparison could not be made.The relationship is typically non-linear, and therefore we calculated Spearman's rank correlation coefficient (ρ) between estimated discharge and observed water level.
Secondly, we used the already mentioned discharge data for 586 ocean-reaching rivers worldwide (Dai et al., 2009).From these, we selected 430 basins for which the reported drainage area was within 20 % of the area derived from the 0.5 • routing network.The ratio between reported and modelderived drainage area was used to adjust the reanalysis estimates, and these were compared with recorded mean streamflow.The recorded mean annual discharge values are not for 2003-2012, but we assume that the differences are not systematic and, therefore, that any large change in agreement may still be a useful indicator of reanalysis quality.
Third, snow storage estimates were evaluated with the ESA GlobSnow product (Luojus et al., 2010).This data set contains monthly 0.25 • resolution estimates of SWE (in mm) for low-relief regions with seasonal snow cover north of 55 • N during 2003-2011.The SWE estimates are derived through a combination of AMSR-E (Advanced Microwave Scanning Radiometer-Earth Observing Satellite) passive microwave remote sensing and weather station data (Pulliainen, 2006;Takala et al., 2009).The GlobSnow data were aggregated to 1 • resolution.The root mean square error (RMSE) and the coefficient of correlation (r 2 ) were calculated as measures of agreement.
Finally, we compared the estimated trends in storage in different glacier regions to trends for mountain glaciers 3 Results

Error estimation
The mean errors derived by the triple collocation technique were of similar magnitude for the GRACE and model estimates (Table 2; note that the numbers listed are for storage change rather than storage per se and were not yet adjusted for GRACE error covariance; cf.Sect.2.4).The relatively low values for the coefficient of variation suggest that the error estimates are reasonably robust.
The spatial error in merged GRACE and model storage change estimates were calculated analogous to Eq. ( 8).The resulting GRACE error surface was relatively homogeneous with an estimated error of around 5-20 mm for most regions, but increasing to 20-40 mm over parts of the Amazon and the Arctic (Fig. 2a).The combined model error surface suggest that errors are smaller than those in the GRACE data for arid regions (<10 mm) but higher elsewhere, increasing beyond 80 mm in the Amazon region (Fig. 2b).The mean errors over non-glaciated land areas were similar, at 18.1 mm for the combined model and 13.5 mm for the combined GRACE data.Assuming no temporal correlation and allowing for error covariance among GRACE products reduces the GRACE error estimates to 10.8 mm (i.e.13.5 2 /2 + 5 2 ).

Analysis increments
Inspection of the analysis increments and the overall difference between prior and posterior estimates provides insights into the functioning of the assimilation scheme (Fig. 3).The spatial pattern in root mean squared (rms) TWS increments ( δS 2 ) emphasises the important role of the world's largest rivers in explaining mismatches between expected and observed mass changes, particularly in tropical humid regions (Fig. 3a).Large increments also occurred over Greenland (mainly due to updated ice storage changes) and the seasonally wet regions of Brazil, Angola and south Asia (subsurface storage).When considering the root mean square difference between prior and posterior estimates of actual TWS (as opposed to monthly changes; Fig. 3b) a similar pattern emerges, but with more emphasis on the smaller but accumulating difference in estimated storage over Greenland, Alaska and part of Antarctica (due to updated ice mass changes) and northwest India (groundwater depletion).

Mass balance and trends
At the global scale, the trend and monthly fluctuations (expressed in SD) in mean total water mass should be close to zero, allowing for small changes in atmospheric water content.This provides a test of internal consistency.Among the original GRACE TWS data, the GRGS data showed the smallest temporal SD (0.04 mm) and linear trend (0.007 ± 0.001 SD mm yr −1 ) in global water mass.The three Tellus retrievals showed larger temporal SD (4.7-6.4 mm) and trends (−0.

Regional storage trends
The spatial pattern in linear trends in the merged GRACE TWS (y 0 ) and the reanalysis signal (y b ) agree well (Fig. 4bc), suggesting that the assimilation scheme is able to merge the prior estimates of storage changes and observed storage as intended.Seasonally adjusted anomalies were calculated for the prior and posterior estimates of the different water cycle components by subtracting the mean seasonal pattern.The 2003-2012 linear trends in these adjusted anomalies (Fig. 5) show that the analysis has (i) increased spatial variability in subsurface water storage trends, with amplified increasing and decreasing trends (Fig. 5a, b); (ii) drastically changed trends in snow and ice storage and typically made them more negative (Fig. 5c, d); and (iii) reversed river water storage trends in the lower Amazon and Congo rivers (Fig. 5e, f).
The reanalysis shows a complex pattern of strongly decreasing and increasing subsurface water storage trends in northwest India (Fig. 5b).This may be an artefact from incorrectly specified errors in the groundwater depletion estimates (see Sect. 4.2).Less visible is that the analysis often reduced negative storage trends in other regions with groundwater depletion, that is, decreased the magnitude of estimated depletion.Because all subsurface storage terms were com-bined, an alternative estimate of groundwater depletion cannot calculated directly, but it can be estimated: for all grid cells with significant prior groundwater depletion estimates (> 0.5 mm yr −1 , representing 99 % of total global groundwater depletion) the 2003-2012 trend in subsurface storage change was estimated a priori at −168 ± 3 (SD) km 3 yr −1 , of which 157 km 3 (94 %) was due to groundwater depletion and the remaining −11 km 3 due to climate variability.Analysis reduced the total trend for these grid cells to −103 ± 3 km 3 yr −1 , from which an alternative groundwater extraction estimate of ca.92 km 3 can be derived.
From the seasonally adjusted anomalies, time series and trends of global storage in different water cycle components were calculated.We calculated snow and ice mass change separately for regions with seasonal snow cover, highlatitude (> 55 • ) glaciers and remaining glaciers (Fig. 6).Gt yr −1 ).Some of the effects of the assimilation were to (i) remove the decreasing trend in prior global terrestrial subsurface water storage estimates (Fig. 6a), (ii) change the poor prior estimates of polar ice cap mass considerably (Fig. 6f, g), (iii) reduce the estimated rate of ocean mass increase from 1.84 ± 0.06 (SD) mm to 1.45 ± 0.05 mm (Table 3) and (iv) achieve mass balance closure between net terrestrial and ocean storage changes (cf.Sect.3.3).

Evaluation against river level remote sensing
The rank correlation (ρ) between river water level and estimated discharge for the 445 grid cells with altimetry time series are shown in Fig. 7. Overall there was no significant change in agreement between the prior (ρ = 0.63 ± 0.27 SD) and posterior (ρ = 0.63 ± 0.26) estimates, with an average change of +0.01 ± 0.12.However, ρ did improve for more locations than it deteriorated (286 vs. 159).There are some spatial patterns in the influence of assimilation (Fig. 7c): strong improvements in the northern Amazon and Orinoco basins and most African rivers, except for some stations along the Congo and middle Nile rivers, and reduced agreement for rivers in China (where prior estimates agreed well) and most stations in the Paraná and Uruguay basins (where they did not).In most remaining rivers, agreement did not change much -in some cases because it was already very good (e.g. the Ganges-Brahmaputra and remainder of the Amazon Basin).Altimetry and estimated discharge time series are shown in Fig. 8 for grid cells with the most data points in three large river systems.In these cases, there is reasonably clear improvement in agreement.

Evaluation against historic river discharge observations
The prior estimate of discharge (i.e. the error-weighted average of the four bias-corrected models) provided estimates that were already considerably better than any of the indi-vidual members (Table 4, Fig. 9).Assimilation led to small improvements in RMSE, from 47 to 44 km 3 yr −1 , and a slight deterioration in the median absolute percentage difference from 40 to 41 %.Combined recorded discharge from the 430 selected basins was 20 909 km 3 yr −1 , representing 90 % of estimated total discharge to the world's oceans according to Dai et al. (2009).Assimilation improved the agreement with this number from −11 to −4 %, of which about half (5 %) is due to a closer estimate of Amazon River discharge.However, modelled and observed discharge values relate to different time periods, and so it is not clear whether this should be considered evidence for improvement or merely reflects multi-annual variability.

Evaluation against snow water equivalent remote sensing
The spatial RMSE and correlation between the prior and posterior SWE estimates and the GlobSnow retrievals are shown in Fig. 10.Although RMSE deteriorated in the majority (57 %) of grid cells, correlation remained unchanged at R 2 = 0.79 and average RMSE improved slightly from 23.2 to 22.3 mm.Assimilation appeared most successful for grid cells with large prior RMSE in northern Canada (Fig. 10a-c).

Evaluation against glacier mass balance estimates
Glacier mass changes reported in the literature (Gardner et al., 2013;Jacob et al., 2012) are listed in Table 5 and compared to regional mass trends associated with glaciers and other components of the terrestrial water derived from the analysis.In the polar regions (e.g.Antarctica, Greenland, Iceland, Svalbard and the Russian Arctic) a large part of the gravity signal is necessarily from glacier mass change.Published trends for most of these regions also heavily rely on GRACE data, and hence our estimates are generally in good agreement.Remaining differences can be attributed to the products, product versions and post-processing methods  3.
Table 4. Evaluation of alternative estimates of mean basin discharge using observations collated by Dai et al. (2009).Listed is the agreement for the ensemble models (without bias correction), the merged prior estimate and the posterior estimates resulting from reanalysis.

CLM MOS NOAH VIC W3RA Prior Posterior
Combined discharge (km 3 yr used, without providing insight into the accuracy of our analysis estimates.In the other regions, the glaciated areas are smaller and surrounded by ice-free terrain, which strongly increases the potential for incorrect distribution of analysis increments, as evidenced by the high trend ratios (> 47 %, last column Table 5).As a consequence, glacier mass trends are not well constrained by GRACE data alone and ancillary observations are required.The agreement with independently derived trend estimates varies.For the Canadian Arctic Archipelago, Alaska and adjoining North America, the assimilation scheme assigns only 55 % (68 Gt yr −1 ) of the total regional negative mass trend (−124 Gt yr −1 ) to glacier mass changes, with most of the remainder (40 % or 50 Gt yr −1 ) assigned to subsurface water storage changes.
Excluding regions for which independent storage change estimates are not available (Greenland, Antarctica and Patagonia), our estimate of total storage change in the world's glaciers (−114 km 3 yr −1 ) was 101 km 3 yr −1 less than the estimate of Gardner et al. (2013) (−215 km 3 yr −1 ).

Estimated errors
The triple collocation method produced estimates of errors in month-to-month changes in GRACE TWS estimates of 12.8-14.3mm over non-glaciated land areas.From these, GRACE TWS errors of 10.4-12.0mm can be estimated (cf.Sect.3.1).By comparison, reported uncertainty estimates based on formal error propagation are larger, usually of the order of 20-25 mm (e.g.Landerer and Swenson, 2012;Tregoning et al., 2012;Wahr et al., 2006).One possible explanation is that the 5 mm we assumed to correct for potential covariance in errors between the GRACE products is too low; another is that the formal uncertainty estimates are conservative.Inflating the GRACE error estimates by 10 mm instead of 5 mm reduced the gain by 18 % on average.The resulting uncertainty in the analysis is modest (see next section).Formal error analyses predict that the retrieval errors decrease towards the poles  (Gardner et al., 2013;Jacob et al., 2012) compared to estimates from reanalysis.Uncertainties are given at the 95 % (2 standard deviation) interval.Also listed are regional trends attributed to other parts of the hydrological cycle, and the ratio of the relative magnitude of those residual trends over estimated glacier mass change.

Reported This study
Trend Glacier trend Other components ratio Superscripts refer to estimates derived from GRACE ( g ) or independent methods ( i ) due to the closer spacing of satellite overpasses (Wahr et al., 2006), but we did not find such a latitudinal pattern.The mean errors in monthly changes in prior TWS for the different models were 16.5-27.9mm.We do not have independent estimates of errors in modelled large-scale TWS with which to compare, but the estimates would seem plausible and perhaps less than anticipated.From a theoretical perspective, violation of the assumptions underpinning triple collocation is likely to have produced overestimates of model error, if anything.The calculated error in the prior estimates over oceans and very stable regions such as Mongolia and the Sahara are around 5 mm (Fig. 2).This provides some further evidence to suggest that the 5 mm GRACE error inflation we applied may have been reasonable.The largest errors in the merged model estimates (> 40 mm) were found for humid tropical regions and high latitudes.The former may be attributed to the combination of large storage variations and often uncertain rainfall estimates.Precipitation measurements are also fewer at high latitudes, while poor prediction of snow and ice dynamics and melt water river hydrology are also likely factors.

Assimilation scheme performance
The spatial pattern in analysis increments emphasises the importance of water stores other than the soil in explaining discrepancies between model and GRACE TWS estimates (Fig. 3).Adjustments to storage changes in large rivers, groundwater depletion, mass changes in high-latitude ice caps and glaciers (e.g.Greenland, Alaska and Antarctica) and lake water levels (e.g. the Caspian Sea and the North American Great Lakes) were all considerable within their region, absorbing monthly analysis increments, long-term trend discrepancies or both.
Uncertainty in error estimates for the different data sources affects the analysis in different ways.Incorrect estimation of GRACE and model-derived TWS errors by the triple collocation method primarily affect (i) the weighting of the ensemble members and (ii) the gain matrix.Appropriate weighting only requires that the relative magnitude of errors among ensemble members is estimated correctly (cf.Eq. 2).The average errors for the different GRACE TWS estimates were all within 14 % of the ensemble average (Table 2) and did not have strong spatial patterns, and therefore the analysis would likely have been very similar if equal weighting had been applied (cf.Sakumura et al., 2014).Estimated model errors showed greater differences (up to 52 % greater than the ensemble mean, Table 2) as well as regional patterns.However, the relative rankings and their spatial pattern were robust to the choice of GRACE TWS members in triple collocation, as evidenced by a low coefficient of variation in error estimates (Table 2).This suggests that the errors were correctly specified in a relative sense.For the gain matrix, the relative magnitude of errors in GRACE versus model TWS ensemble means needed to be estimated correctly (cf.Eq. 8).The estimated GRACE TWS ensemble errors are reasonably homogeneous in space (Fig. 1a), which increases our confidence in their validity.The uncertainty due to the correction for assumed correlation between the GRGS and Tellus TWS (see previous section) is further mitigated by the design of the data assimilation scheme: the gain factor determines how rapidly the analysis converges towards the GRACE observations and therefore is important for month-to-month variations, but long-term trends in TWS will always approach those in the GRACE observations (cf.Fig. 4b and c).
The main sources of uncertainty in long-term trends in the individual water balance terms are (i) the removal of nonhydrological mass trends in the GRACE TWS time series and (ii) accurate specification of relative errors in the individual water balance terms, which is needed for correct redistribution of the integrated TWS analysis increments.For example, the analysis results illustrate the insufficiently constrained problem of separating gravity signals due to mass changes in mountain glaciers from nearby subsurface water storage changes.This was particularly evident around the Gulf of Alaska and northwest India, where decreases can be expected not only in glacier mass but also in subsurface storage due to, respectively, a regional drying trend and high groundwater extraction rates (Fig. 5a).We suspect that unexpectedly strong increasing storage trends in parts of northwest India may be because the prior groundwater depletion estimates were too high and the assigned errors too low, causing the analysis update to distribute increments incorrectly.We could have addressed this by inflating the local groundwater depletion estimation errors, but more research is needed to understand the underlying causes.Plausible causes are that groundwater extraction is overestimated, or that extraction is compensated by induced groundwater recharge (e.g. from connected rivers) (see Wada et al., 2010, for further discussion).
Mass balance closure was not enforced and hence provides a useful diagnostic of reanalysis quality.The GRGS product achieved approximate global mass balance closure at all timescales, but the three Tellus products showed a seasonal cycle and long-term negative trend in global water mass.Accounting for atmospheric water vapour mass changes (from ERA-Interim reanalysis and the NVAP-M satellite product, data not shown) could not explain the trends and in fact slightly increased the seasonal cycle in global water mass.Data assimilation reduced the seasonal cycle and entirely removed the trend in total water mass, thanks to the prior estimates of sea mass increase.For comparison, we calculated average ocean mass increases by an alternative, more conventional method, which involved avoiding areas likely to be affected by nearby land water storage changes.Excluding a 1000 km buffer zone produced a 2003-2012 mass trend of +0.58 to +0.72 mm yr −1 for the three Tellus retrievals, +1.12 mm yr −1 for the GRGS retrieval and +0.75 mm yr −1 for the merged GRACE data.Data assimilation produced a stronger trend of +1.22 mm yr −1 due to the influence of the prior estimate of +1.67 mm yr −1 .Our prior estimate followed Chen et al. (2013), who used an iterative modelling approach to attribute 75 % of altimetryobserved SLR to mass increase.Chen et al. (2013) argue that the conventional method produces underestimates of ocean mass increase.Indeed, the trends we calculated for the "buffered" ocean regions are lower than for the entire oceans (+1.22 vs. +1.45mm yr −1 for the reanalysis and +1.67 vs. +1.84mm yr −1 for the prior estimates; Table 3).Nonetheless, the reduction in sea mass change of 0.39 mm yr −1 from prior to analysis does appear to reopen the problem of reconciling mass and temperature observations with the altimetry derived mean sea level rise of +2.45 ± 0.08 mm yr −1 (cf.Chen et al., 2013).

Evaluation against observations
The reanalysis generally did not have much impact on the agreement with river and snow storage observations, with small improvements for some locations and small degradations for others.While a robust increase in the agreement would have been desirable, the fact that agreement was not degraded overall is encouraging.The data assimilation procedure applied has the important benefit of bringing the estimates into agreement with GRACE observations.Moreover, performance improvements with respect to river discharge and level data did occur in the Amazon, where they make an important contribution to TWS changes.Similarly, snow water equivalent estimates were improved in the North-American Arctic, where errors in the prior estimates were largest.This demonstrates that GRACE data can indeed be successfully to constrain water balance estimates, although further development may be needed to avoid some of the undesired performance degradation for water balance components that do not contribute much to the TWS signal.
The models used for our prior estimates provided poorly constrained estimates of ice mass balance changes, and our reanalysis ice mass loss estimates should not be assumed more accurate than estimates based on more direct methods (Table 5).Our analysis is unique when compared to previous estimates based on GRACE, in that data assimilation allowed some of the observed mass changes to be attributed to other water balance components within the same region, depending on relative uncertainties in the prior estimates.Comparison against independent estimates of glacier mass balance changes also demonstrated the challenge of correct attribution, however.Glacier mass balance estimates were in good agreement for several regions, but estimates for North Amer- ican glaciers in particular were questionable: their combined mass loss (−68 Gt yr −1 ) was much lower than the estimates derived by independent means (−124 Gt yr −1 ; Table 5).This can be explained by incorrect specification of errors.Two caveats are made: (i) the GIA signal is relatively large for these three regions (+50 Gt yr −1 ), and hence GIA estimation errors may have had an impact, and (ii) a significant change in subsurface water storage is plausible in principle; for example, higher summer temperatures could be expected to enhance permafrost melting and runoff, as well as enhance evaporation.More accurate spatiotemporal observation and modelling of glacier dynamics are needed to reduce this uncertainty.

Contributions to sea level rise
The reanalysis estimate of net terrestrial water storage change of −495 Gt yr −1 (Table 3) appears a plausible estimate of ocean mass change, equivalent to ca. +1.4 mm yr −1 sea level rise.Our results confirmed that mass loss from the polar ice caps is the greatest contributor to net terrestrial water loss, with Antarctica and Greenland together contributing −342 Gt yr −1 .The next largest contribution was from the remaining glaciers.We combine the reanalysis estimate of −129 Gt yr −1 with another −101 Gt yr −1 estimated to be misattributed (cf.Sect.3.8) and obtain an alternative estimate of −230 Gt yr −1 .A small but significant contribution of −18 Gt yr −1 (Table 3) was estimated to originate from reductions in seasonal snow cover (particularly in Quebec and Siberia; Fig. 5cd).Interannual changes in river water storage were not significant.Small contributions of −10 and +16 Gt yr −1 were attributed to storage changes existing lakes and large new dams, respectively, and compensated each other.The largest change in an individual water body was in the Caspian Sea (−27 Gt yr −1 ; cf.Fig. 5), which experiences strong multi-annual water storage variations depending on Volga River inflows.
Finally, the analysis suggested a statistically insignificant change of +9 Gt yr −1 in subsurface storage globally.Adding back the suspected misattribution of 101 Gt yr −1 associated with glaciers produces an alternative estimate of +110 Gt yr −1 (cf.Fig. 6a).Combining this with the −92 Gt yr −1 attributed to groundwater depletion suggests that storage over the remaining land areas increased by 202 Gt yr −1 .Calculating subsurface storage trends by latitude band suggests that most of this terrestrial water "sink" can be found north of 40 • N and between 0 and 30 • S, and is in fact opposite to the prior estimates (Fig. 11).The main tropical regions experiencing increases are in the Okavango and upper Zambezi basins in southern Africa and the Amazon and Orinoco basins in northern South America (Fig. 5b).Storage increases for these regions are also evident from the original GRACE data (Fig. 4a) and cannot be attributed to storage changes in rivers or large lakes.The affected regions contain low-relief, poorly drained areas with (seasonally) high rainfall.In such environments, the storage changes could occur in the soil, groundwater, wetlands or a combination of these.Further attribution is impossible without additional constraining observations (Tregoning et al., 2012;van Dijk et al., 2011).The 10-year analysis period is short, and this cautions against over-interpreting this apparent "tropical water sink".However it is of interest to note that a gradual strengthening of global monsoon rainfall extent and intensity has been observed, and is predicted to continue (Hsu et al., 2012).In any event, the difference between prior and posterior trends in Fig. 11 illustrates that the current generation of hydrological models, even as an ensemble, are probably not a reliable surrogate observation of long-term subsurface groundwater storage changes.GRACE observations proved valuable in improving these estimates.

Conclusions
We presented a global water cycle reanalysis that merges four total water storage retrieval products derived from GRACE observations with water balance estimates derived from an ensemble of five global hydrological models, water level measurements from satellite altimetry and ancillary data.We summarise our main findings as follows: 1.The data assimilation scheme generally behaves as desired, but in hydrologically complex regions the analysis can be affected by poorly constrained prior estimates and error specification.The greatest uncertainties occur in regions where glacier mass loss and subsurface storage declines (may) both occur but are poorly known (e.g.northern India and around North American glaciers).
2. The error in original GRACE TWS data was estimated to be around 11-12 mm over non-glaciated land areas.Errors in the prior estimates of TWS changes are estimated to be 17-28 mm for the five models.
3. Water storage changes in other water cycle components (seasonal snow, ice, lakes and rivers) are often at least as important and uncertain as changes as subsurface water storage in reconciling the various information sources.
4. The analysis results were compared to independent river water level measurements by satellite altimetry, river discharge records, remotely sensed snow water storage and independent estimates of glacier mass loss.In all cases the agreement improved or remained stable compared to the prior estimates, although results varied regionally.Better estimates and error specification of groundwater depletion and mountain glacier mass loss are required.
6.For the period 2003-2012, we estimate glaciers and polar ice caps to have lost around 572 Gt yr −1 , with an additional small contribution from seasonal snow (−18 Gt yr −1 ).The net change in surface water storage in large lakes and rivers was insignificant, with compensating effects from new reservoir impoundments (+16 Gt yr −1 ), lowering water level in the Caspian Sea (−27 Gt yr −1 ) and increases in the other lakes combined (+16 Gt yr −1 ).The net change in subsurface storage was significant when considering a likely misattribution of glacier mass loss, and may be as high as +202 Gt yr −1 when excluding groundwater depletion (−92 Gt yr −1 ).Increases were mainly in northern temperate regions and in the seasonally wet tropics of South America and southern Africa (+87 Gt yr −1 ).Continued observation will help determine if these trends are due to transient climate variability or likely to persist.

Figure 1 .
Figure 1.Illustration of the data assimilation approach followed using data along a transect through the USA for August 2003.Shown are (a) monthly satellite-derived TWS, y 0 t , and the equivalent prior estimate, y b t ; (b) location of the west-east transect on a map of the gain matrix, k; (c) profile of k along the transect (cf.Fig. 2c); (d) calculation of the TWS analysis increment, δ y t , from k and the innovation, (y 0 t −y b t ); (e) the prior error in the change of each of the stores, σ t (i); (f) the prior and posterior estimate of change in each store, s b t (i) and s b t (i) + δs t (i), respectively; and (g) visual illustration of the disaggregation of the TWS analysis increments to the different stores.All units are in millimetres unless indicated otherwise; see text for full explanation of symbols.Stores shown include the subsurface (green), rivers (blue) and sea (dark red; remaining stores not shown for clarity).

Table 2 .
Spatial mean values (non-glaciated land areas only) of the error in monthly mass change estimates for different GRACE and model sources as derived through triple collocation.Also listed is the number of triple collocation estimates derived (N) and the spatial mean of the coefficient of variation (C.V.) in these N estimates.Mean error Mean C.V. et al. (2013) for 2003-2010 and for Greenland and Antarctica by Jacob et al. (2012) for 2003-2009.In several cases these mass balance estimates were based on independent glaciological or ICESat (Ice, Cloud and land Elevation Satellite) observations, and these were the focus of comparison.Other estimates were partially or wholly based on GRACE data, making comparison less insightful.

Figure 2 .
Figure 2. Triple collocation estimated error in storage change from the merged (a) GRACE and (b) prior estimates, and (c) resulting gain matrix.

Figure 3 .
Figure 3.The impact of GRACE data assimilation on total water storage expressed as (a) the root mean square (rms) analysis increment and (b) the rms difference between prior and posterior storage time series.

Figure 4 .
Figure 4. Trends in GRACE total water storage as derived from (a) prior storage estimates, (b) merged satellite retrievals and (c) posterior estimates.

Figure 5 .
Figure 5. Trends in seasonal anomalies of prior (left column) and posterior (right column) estimates of (a-b) subsurface, (c-d) snow and (e-f) surface (i.e.lake and river) water storage.

Figure 6 .
Figure 6.Time series of the prior (grey lines) and posterior (black lines) estimates of seasonally adjusted global equivalent water height anomalies (SAGEWHA) in different water cycle components.Dashed lines show linear trends for 2003-2012 as listed in Table3.

Figure 7 .
Figure 7. Effect of assimilation agreement with satellite altimetry river water levels: Spearman's rank correlation coefficient (ρ) for (a) prior and (b) posterior estimates and (c) difference between the two.

Figure 10 .
Figure 10.Effect of assimilation on agreement with GlobSnow snow water equivalent estimates, showing (a-c) root mean square error (RMSE) and (d-f) the coefficient of correlation (R 2 ).From left to right, agreement for (a, d) prior and (b, e) posterior estimates as well as (c, f) the change in agreement.

Table 1 .
Description and sources of data used in this analysis.Acronyms are explained in the text.

van Dijk et al.: A global water cycle reanalysis (2003-2012)
• grid cell averages before use in assimilation.The routing function A. I. J. M.

Table 3 .
Calculated linear trends in global mean seasonally adjusted anomalies associated with different water cycle components for 2003-2012.The posterior trend estimates are also expressed in equivalent sea level rise (SLR) and volume.Second number is standard deviation.

Table 5 .
Published trends in glacier water storage