Total land water storage change over 2003–2013 estimated from a global mass budget approach

We estimate the total land water storage (LWS) change between 2003 and 2013 using a global water mass budget approach. Hereby we compare the ocean mass change (estimated from GRACE space gravimetry on the one hand, and from the satellite altimetry-based global mean sea level corrected for steric effects on the other hand) to the sum of the main water mass components of the climate system: glaciers, Greenland and Antarctica ice sheets, atmospheric water and LWS (the latter being the unknown quantity to be estimated). For glaciers and ice sheets, we use published estimates of ice mass trends based on various types of observations covering different time spans between 2003 and 2013. From the mass budget equation, we derive a net LWS trend over the study period. The mean trend amounts to +0.30 ± 0.18 mm yr−1 in sea level equivalent. This corresponds to a net decrease of −108 ± 64 km3 yr−1 in LWS over the 2003–2013 decade. We also estimate the rate of change in LWS and find no significant acceleration over the study period. The computed mean global LWS trend over the study period is shown to be explained mainly by direct anthropogenic effects on land hydrology, i.e. the net effect of groundwater depletion and impoundment of water in man-made reservoirs, and to a lesser extent the effect of naturally-forced land hydrology variability. Our results compare well with independent estimates of human-induced changes in global land hydrology.


Introduction
Liquid fresh water on land is stored in various reservoirs: rivers, lakes, man-made reservoirs, wetlands and inundated areas, root zone (upper few meters of the soil) and aquifers (groundwater reservoirs). Terrestrial reservoirs continuously exchange with the atmosphere, oceans and land, through vertical and horizontal mass fluxes (precipitation, evaporation, transpiration of the vegetation, surface runoff and underground flow). Land water storage (LWS) varies with change in mean climate and climate variability. Human activities also directly affect LWS through water extraction from aquifers, building of dams along rivers, urbanization, wetland drainage, land use and land cover changes, and deforestation. All these effects modify the water budget in river basins, and because of water mass conservation in the climate system, cause sea level changes. Studies based on hydrological modeling have not reported any clear long-term trend in global LWS over the past 60 years but only interannual variability (e.g., Ngo-Duc et al 2005). This is unlike human-induced factors such as dam building (Chao et al 2008) and groundwater extraction (Konikow 2011, Pokhrel et al 2012, Wada et al 2012, Wada 2015. Although their contributions to the global mean sea level (GMSL) are of opposite sign (<0 for dams, >0 for groundwater pumping), their net effect is responsible for a significant longterm positive trend at least for the recent decades (other human-induced factors have negligible contributions to the GMSL). Building on the results from Konikow (2011) and Wada et al (2012), Church et al (2013) estimated that the net effect of dams and groundwater depletion (i.e., groundwater abstraction minus recharge; e.g., Wada 2015) on the GMSL amounted 0.38±12 mm yr −1 over 1993-2010. This represents 12% of the observed GMSL rate of rise over this time span, an amount of the same order of magnitude as the Antarctic ice mass loss (see table 13.1 in Church et al 2013). Because of such a significant contribution to sea level, it is worth to examine this component in more detail. In addition, uncertainty of this component has direct impact on our capability to close the sea level budget, thus constrain missing contributions (due to lack of data) such as the deep ocean thermal expansion (see discussions on that topic in Llovel et al 2014, Dieng et al 2015a.
The effect of dams and man-made reservoirs has been estimated by Chao et al (2008). They reconstructed the history of water impoundment in the nearly 30 000 reservoirs built during the twentieth century and estimated the contribution to sea level by dams and artificial reservoirs (including seepage) at −0.55±0.08 mm yr −1 in sea level equivalent (SLE) during the last half-century, with a stabilization in recent years. Estimates of groundwater depletion are based on three methods (see Wada 2015): (1) volumebased method, (2) flux-based method, and (3) satellite observations from the GRACE space gravimetry mission. Each method has strengths and weaknesses. For example, GRACE gives a vertically integrated estimate of the water mass change; thus surface waters and soil moisture must be known and removed to estimate the ground water contribution. In addition, GRACEbased estimates do not yet have full global coverage for the estimation of groundwater depletion (Famiglietti 2014, Wada 2015. The volume-and flux-based methods lack global information and suffer from model uncertainties. Thus estimating groundwater depletion remains very challenging, as is the global dam contribution. In this study, we develop another approach based on the global water mass budget of the climate system to estimate the total LWS change. Focusing on the January 2003-December 2013 time span (for which GRACE data are available), we compare the GRACEbased ocean mass change to the sum of mass components (glaciers, Greenland and Antarctica ice sheets, atmospheric water vapor and LWS). We neglect other mass components such as permafrost because global data are lacking, as well as change in the snow pack, previously shown to give negligible contribution to the GMSL beyond time spans larger than 1 year (Biancamaria et al 2011). In this mass budget approach, we use estimates of each component from different observational data sets, except for the net total LWS, the unknown quantity to be estimated. A mean LWS trend is first estimated over 2003-2013. Then, accounting for increasing rate of change (acceleration) of several components (ocean mass, ice sheet mass balances) over the study time span, we investigate whether the LWS rate varies with time. To validate our results, we perform a similar analysis but instead of using GRACE, we estimate the ocean mass term from the satellite altimetry-based GMSL corrected for steric effects (i.e., effects of ocean temperature and salinity).
All results are expressed in terms of SLE change. Units are given in mm yr −1 and mm yr −2 for trend and acceleration respectively.

Method
To estimate the contribution of LWS change to sea level, we can simply consider the conservation of water mass in the Earth's system (e.g., Llovel et al 2010). Of course, LWS change could be derived from GRACE data over the continents, as done previously in a number of studies. However, considering that the GRACE resolution (∼300-500 km) may be problematic in separating nearby masses (e.g., river basins and glaciers), our objective here is to use a different approach.
On time scales of years to decades, water mass changes inside the solid Earth (e.g., in the crust) can be neglected, so that only changes in land reservoirs, ocean and atmosphere need to be considered, with the mass conservation equation written as follows: Atm .

LWS
where ΔM(t) represents changes with time t of water mass in the different reservoirs: ocean, glaciers (including small ice caps), Greenland ice sheet, Antarctica ice sheet, atmosphere and land water stores. Note that ΔM(t) may be either positive or negative. Using equation (1), we deduce the LWS component by simple rewriting as: Antarct.

Atm
As mentioned above, all contributions are expressed in terms of SLE.

Ocean mass
For estimating the ocean mass component, we apply two approaches : (1) use of GRACE space gravimetry data over the oceans, and (2) estimate of the GMSL corrected for steric effects.
3.1.1. GRACE-based ocean mass Three different data sets of the GRACE Release 05 products have been considered: (1) from the Texas University (CSR RL05), (2) from the German GeoForschungsZentrum (GFZ RL05) (3) from the Jet Propulsion Laboratory (JPL RL05).
To study the ocean mass evolution, a specific processing has been carried out by D. Chambers (described in Johnson and Chambers 2013; geocenter terms included; data available at https://dl. dropboxusercontent.com/u/31563267/ ocean_mass_orig.txt). The data are provided as global mean (averaged over the 90°S-90°N°domain) time series at monthly interval with associated uncertainty. The GIA (Glacial Isostatic Adjustment) effect is corrected for using the GIA correction computed in Chambers et al (2010). In the following, we consider the mean of the three data sets.

Ocean mass estimated from the GMSL corrected for steric effects
Changes in the GMSL result from steric effects plus ocean mass changes. Thus, the ocean mass component can be also derived from the difference 'GMSL minus steric effects'. For that purpose we used the mean of six different satellite altimetry-based GMSL data sets: (1) Validation and Interpretation of Satellite Oceanographic (AVISO; http://aviso. altimetry.fr/en/data/products/ocean-indicatorsproducts/actualitesindicateurs-des-oceansniveaumoyen-des-mersindexhtml.html); (2) University of Colorado (CU Release 5; http://sealevel.colorado. edu/); (3) National Oceanographic and Atmospheric Administration (NOAA; http://star.nesdis.noaa. gov/sod/lsa/SeaLevelRise/LSA_SLR_timeseries_ global.php); (4) Goddard Space Flight Center (GSFC version 2; http://podaac-ftp.jpl.nasa.gov/dataset/ MERGED_TP_J1_OSTM_OST_GMSL_ASCII_ V2); (5) Commonwealth Scientific and Industrial Research Organization (CSIRO; http://cmar.csiro. au/sealevel/sl_data_cmar.html); (6) The European Space Agency/ESA Climate Change Initiative/CCI sea level data (http://esa-sealevel-cci.org/). Details on these data sets can be found in Dieng et al (2015aDieng et al ( , 2015b. For the steric component, instead of using Argo that suffer from gaps in the data coverage (e.g., in the Indonesian region; Dieng et al 2015b), we make use of the ORAS4 reanalysis (Balmaseda et al 2013) that provides ocean temperature and salinity down to 5350 m and global coverage. Note that over their common geographical and depth coverage, Argo-based and ORAS4-based steric sea level are in good agreement (see Dieng et al 2015b for a discussion).
Using the mean of the six GMSL products, we compute the ocean mass component by subtracting the ORAS4 steric component. It is simply called below 'GMSL minus ORAS4'.

Atmospheric water vapor mass
To estimate change in atmospheric water vapor mass, we used the vertically integrated water vapor grids from the ERA Interim reanalysis performed by the European Center for Medium Range Weather Forecast/ECMWF (Dee et al 2011). The data are provided as 1.5°×1.5°grids at monthly interval. We compute a globally averaged water vapor time series and express it in terms of SLE (see Dieng et al 2014 for details).

Greenland and Antarctica mass
For the ice sheet mass balances, we used two approaches: (1) time series given by Velicogna et al (2014) and from the ESA CCI Ice Sheet project  Note that the four glaciers estimates considered in this study do not include Greenland and Antarctica peripheral glaciers. Trend values are given in the SI.

Data analysis
When time series are used to estimate trends and accelerations, the annual & semi annual signals are removed by fitting 12-month and 6-month sinusoids. Figure 1 shows the GRACE-based ocean mass (called GOM) time series over 2003-2013. This time series is an update of that previously used by Dieng et al (2015aDieng et al ( , 2015b to examine the closure of the sea level budget. A mean trend of 1.85±0.1 mm yr −1 is estimated over the study time span. We fitted a degree 2 polynomial to the data, from which we deduce an acceleration of 0.29±0.04 mm yr −2 . The acceleration is defined as 2 times the adjusted coefficient of the polynomial t 2 term. The quoted uncertainties represent 1 sigma errors estimated from the leastsquares fit and accounting for the time series errors. In figure 1 is superimposed the 'GMSL minus ORAS4'ocean mass time series and associated uncertainty (note that the GOM uncertainty is not shown because smaller than the latter). The mean 'GMSL minus ORAS4' trend over 2003-2013 amounts to 2.03±0.11 mm yr −1 . Besides, the acceleration is found to be almost zero over the study time span.
Similarly, figure 2 shows the global atmospheric water vapor time series. The mean trend estimated from the time series is slightly negative (equal to −0.04±0.04 mm yr −1 SLE), indicating a small but not significant increase in atmospheric water vapor content. Dieng et al (2014) considered other water vapor datasets and found little differences in terms of interannual variability and trend. As for the ocean mass data, we fitted a degree 2 polynomial but found zero acceleration.     A treatment similar to that applied for the ice sheets was performed for the glaciers mean trend and acceleration (see We are also aware that the estimated acceleration needs to be used with caution due to the few available glacier observation-based data sets.

LWS trend over 2003-2013
In figure 6, we present a chart of the mean trends over 2003-2013 for ΔM Ocean (from GRACE and from 'GMSL minus ORAS4'), the sum of ΔM (atmospheric water vapor plus glaciers plus ice sheets) and the residuals (ΔM Ocean -sum of ΔM). For the sums (and residuals as well), 3 trend values are considered. Sums a, b and c correspond to: (a) the average trend estimated with the published results (i.e., sum of mean Greenland Residuals 1 and 2 (called res1 and res2) in figure 6 are based on ΔM Ocean from GRACE and from 'GMSL minus ORAS4', to which 'sum a' is subtracted. We also compute the mean of res1 and res2. Corresponding trends amount to 0.21±0.18 mm yr −1 (res1), 0.39±19 mm yr −1 (res2) and 0.30±0.18 mm yr −1 (mean). Residuals 3 and 4 (called res3 and res4) are based on sums b and c respectively, using the mean value of ΔM Ocean from GRACE and 'GMSL minus ORAS4'. These are also plotted in figure 6 as well as LWS trends estimated by Wada et al (2012)-based on the flux method-for 2 cases: (1) only dams and groundwater depletion are accounted for, and (2) in addition to dams and ground waters, account of deforestation and wetland drainage (called Wada1 and Wada2 hereinafter). Trends over 2003-2013 for res3, res4, Wada1 and Wada2 amount to 0.44±24 mm yr −1 , 0.50±0.20 mm yr −1 , 0.39±0.11 mm yr −1 and 0.54±0.12 mm yr −1 , respectively.
All computed residuals appear rather consistent in spite of the quite different data sets used. They also compare rather well with Wada1 & 2. Trend values are gathered in the SI.

LWS acceleration over 2003-2013
As mentioned above, figures 3 and 4 show clear acceleration for the Greenland and Antarctica mass balances (unlike the glacier and water vapor components). The ΔM Ocean based on GRACE also displays important acceleration over the study time span. Using the mass budget equation, we can deduce the acceleration of the LWS residuals for this case. For res1, it amounts to +0.17±0.04 mm yr −2 . However if we consider ΔM Ocean based on 'GMSL minus ORAS4', acceleration of the residual time series (res2) becomes negative (and equal to −0.08±0.05 mm yr −2 ). Besides accelerations of Wada 1 and 2 are very small and non significant (of 0.008±0.010 mm yr −2 ). With the data currently available, it does not seem possible to estimate any reliable LWS acceleration, nor to identify which term of the mass budget equation compensates the ice sheet mass balance acceleration. Our results suggest nevertheless that the LWS acceleration is not significantly different from zero over the 2003-2013 time span.

Interannual variability in LWS trends over 2003-2013
We computed short-term trends of the ΔM Ocean time series based on GRACE and 'GMSL minus ORAS4', over successive 2-year time spans (with 1 year overlap). To these short-term trends, we removed the acceleration term of sum a (i.e., the combined acceleration of glaciers, ice sheets and water vapor). The corresponding curves are shown in figure 7 (labeled res1 and res2). In figure 7 are superimposed four additional LWS 2-year trend curves (also with 1-year overlap) using: (1) LWS determined by Yi et al (2015) using GRACE over continental river basins, (2) (Döll et al 2014a(Döll et al , 2014b. The MERRA dataset is the version 5.2.0 of the GEOS-5 data service. We used the total water storage in land reservoirs product that includes the groundwater component. As for ISBA/TRIP and WGHM, it is available as gridded time series at monthly interval over the 2003-2013 time span. We computed geographical averages, applying a cosine latitude weighting. All six curves (expressed in SLE trends) exhibit large interannual variability, mostly related to El Nino-Southern Oscillation (ENSO) events (note for example the minimum corresponding to the 2011 La Nina). On average, a good correlation at interannual time scale is noticed between these six curves. But we note closer agreement between 2-year trend curves from res1 and Yi et al LWS on the one hand, and ISBA/TRIP and MERRA LWS on the other hand. The Wada 1 and Wada 2 short-term trends are also shown in figure 7 (bottom curves; signal amplified by a factor of 10). It is interesting to note that these also display a minimum in 2011. Although the latter only represent the direct anthropogenic components, the 2011 minimum likely reflects increased groundwater recharge during this La Nina episode.
All trends and accelerations estimates presented above are gathered in the SI.

Discussion
The mean trend in LWS estimated by the global mass budget approach developed in this study is found to be positive in terms of SLE over the 2003-2013 time span. The mean of the two estimates based on two different values of ΔM Ocean is 0.30±0.18 mm yr −1 SLE. This corresponds to an annual decrease in net LWS of −108±64 km 3 yr −1 . This quantity represents the combined effects of natural climate variability, anthropogenic climate change and direct anthropogenic factors. The uncertainty of this estimate directly relies on the ocean mass trend uncertainty. Here we used two independent methods to estimate ΔM Ocean and associated uncertainty, with quite consistent results. Wada et al (2012)'s results for the direct anthropogenic LWS components are only slightly larger (of +0.39±0.11 mm yr −1 SLE or −140±40 km 3 yr −1 LWS trend) for the net effect of dams and ground water depletion. While the rate of reservoir impoundment exceeded groundwater depletion over most of the 20th century, for the recent years, groundwater depletion exceeds impoundment, thus the net effect leads to a positive contribution to the GMSL. In addition to Wada et al (2012) (2014) LWS trends (both based on GRACE data processing on land) are respectively slightly positive and negative (+0.07±0.04 mm yr −1 and −0.06±0.09 mm yr −1 respectively). This dispersion of LWS values and large associated uncertainties based on GRACE is not totally surprising. The GRACE LWS rate estimates are much dependent on the study period, considering the importance of the interannual variability (see figure 7), as discussed in Jensen et al (2013). Besides, the GRACE resolution does not allow unambiguous separation between nearby sources (Longuevergne et al 2010). This is particularly true in the region of the Ganges basin and Himalayan glaciers. Finally, studies considering a limited number of river basins (e.g., Llovel et al 2010) miss part of the signal, and the GIA correction, important in high latitudes, adds significant uncertainty. All together, the GRACE LWS estimates remain uncertain, in particular if study time spans are short (Landerer and Swenson 2012).
In this study, we proposed a different approach to estimate the net LWS contribution to GMSL change over a 10 year time span. Using a large number of different data sets for the mass components, we came up to a positive value for the LWS trend over 2003-2013 (in terms of SLE; i.e., decrease of total water storage on land) that likely reflects the net anthropogenic component, i.e., the dominant contribution of groundwater depletion versus dams, in good agreement with the Wada et al (2012) estimates.