Are GRACE-era Terrestrial Water Trends Driven by Anthropogenic Climate Change ?

To provide context for observed trends in terrestrial water storage (TWS) during GRACE (2003–2014), trends and variability in the CESM1-CAM5 Large Ensemble (LE) are examined. Motivated in part by the anomalous nature of climate variability during GRACE, the characteristics of both forced change and internal modes are quantified and their influences on observations are estimated. Trends during the GRACE era in the LE are dominated by internal variability rather than by the forced response, with TWS anomalies in much of the Americas, eastern Australia, Africa, and southwestern Eurasia largely attributable to the negative phases of the PacificDecadalOscillation (PDO) andAtlanticMultidecadalOscillation (AMO).While similarities between observed trends and themodel-inferred forced response also exist, it is inappropriate to attribute such trendsmainly to anthropogenic forcing. For several key river basins, trends in the mean state and interannual variability and the time at which the forced response exceeds background variability are also estimated while aspects of global mean TWS, including changes in its annual amplitude and decadal trends, are quantified. The findings highlight the challenge of detecting anthropogenic climate change in temporally finite satellite datasets and underscore the benefit of utilizing models in the interpretation of the observed record.


Introduction
Changes in water availability are a fundamental driver of impacts related to both climate variability and change, yet quantifying and understanding these variations remain a major challenge, in part due to the complexities inherent to their observation, simulation, and attribution on regional and global scales.While depletion of ground water can arise from a range of factors, including groundwater pumping, irrigation, and other agricultural and water management practices (e.g., [1]), changes in rainfall and temperature also exert a fundamental influence.As Earth warms and the atmosphere moistens, changes in the water cycle are anticipated including the spatial redistribution of precipitation favoring tropical and high latitudes, with potential deficits in many subtropical and midlatitude regions (e.g., [2]).An increasingly variable water cycle is also anticipated, with enhanced frequencies of flooding and drought (e.g., [3,4]).From a practical perspective, these changes pose significant risks for water resources management and for ensuring future water availability [3,5].
Motivated by the importance of these issues, the Gravity Recover and Climate Experiment (GRACE) satellite constellation was launched in March 2002 to monitor changes in near-surface mass, which over land relate primarily to the accumulation or depletion of terrestrial water storage (TWS).Providing the first global estimates at monthly resolution, GRACE revealed major changes across a range of timescales including strong seasonal shifts, interannual variations, and apparent trends [6].While it is clear that both atmospheric and terrestrial hydrologic processes play central roles in these variations, many questions remain, chief amongst which is the potential role of anthropogenic climate change on GRACE-era trends.
Here, the roles of both internal variability (i.e., arising from interactions and chaotic variability within the climate system) and the anthropogenic forced response (i.e., driven primarily by human induced changes in atmospheric composition) are therefore investigated.This analysis is timely in that the GRACE record is likely nearing its end of life, with weakening of the satellite's batteries currently requiring 2 Advances in Meteorology the instrument to be shut off for extended periods.Major questions nonetheless remain unanswered, such as whether the GRACE era is representative of recent decades generally.There are reasons to question whether this is the case given the strong negative phase of the Pacific Decadal Oscillation (PDO) that prevailed through much of the 2000s, accompanying a slowdown in global mean surface warming (e.g., [7]).The role of this and other modes such as the El Niño/Southern Oscillation (ENSO) and Atlantic Multidecadal Oscillation (AMO), in influencing regional TWS, remains similarly unexplored.A range of questions related to climate change can also be asked.Is an intensifying water cycle likely to result in a detectable increased annual cycle in TWS due to enhanced winter snow accumulation, advanced spring melt, and increased summer aridity?Changes in precipitation frequency and seasonality may also occur (e.g., [4]), if so, what is the expected magnitude of associated trends and how do they compare with internal variability?What is the capacity of trends in TWS to be sustained for decades and mitigate sea level rise, a question of particular relevance to the interpretation of the altimeter record [8]?
Addressing these questions solely through observations remains challenging due to the limited duration and accuracy of the global water cycle record and so here a hybrid approach is adopted in which model simulations are validated against available observations and then used to estimate the spatial patterns and magnitude of internal and forced variability on longer timescales.Section 2 describes the methods and data used, including the rationales for their selection, while Section 3 documents observed trends and modes of climate variability.Section 4 addresses questions raised by observed variability and documents patterns associated with the El Niño/Southern Oscillation (ENSO), the Pacific Decadal Oscillation (PDO), the Atlantic Multidecadal Oscillation (AMO), and the forced response under anthropogenic climate change.Section 5 documents simulated variability in key regions while Section 6 addresses characteristics of present-day and future global budgets, followed by a discussion of key findings and conclusions in Section 7.

Data and Methods
2.1.GRACE.GRACE fields are used from several sources including the recently released mass concentration solutions (Mascon, [9]) and various unconstrained spherical harmonic solutions [10][11][12].The currently available GRACE Mascon record extends from April 2002 through November 2014 while spherical harmonic solutions extend through April 2015.The Mascon fields are provided on a 0.5 ∘ grid at quasimonthly resolution while the spherical harmonic solutions are on a 1 ∘ grid, though the resolution of the grids on which these fields are reported is higher than the 300-400 km resolution of the GRACE retrievals.Due to battery management-related instrument outages starting in 2011, a few GRACE months are either missing or partially averaged across multiple months.The fields are thus converted to monthly means to allow for the removal of the climatological annual cycle prior to computing trends.Ice sheet surface mass variations in GRACE are not evaluated.

The CESM1-CAM5 Large Ensemble. The CESM1-CAM5
Large Ensemble (LE, [13]) consists presently of a 40-member set of simulations spanning 1920-2100 (with member 1 beginning in 1850).A round-off error magnitude (order 10 −14 K) adjustment to tropospheric temperature is used to generate the ensemble spread in a series of branch runs that began in 1920.An extended 2000-yr coupled control run is also included as part of the ensemble.The LE is a unique resource for studying the terrestrial water budget's internal variability and distinguishing it from the forced response as, unlike many multimodel archives, it provides all terms necessary to compute the full terrestrial water budget and allows for an estimation of the forced response from the average of its ensemble members (whereas structural contrasts in multimodel archives preclude doing so).Moreover, the CESM consistently scores among the best models in representing the major modes of climate variability (e.g., [14]).Terrestrial water storage (TWS) is calculated from the LE's land component, version 4 of the Community Land Model (CLM4, [15]).TWS is calculated in CLM4 as the combination of soil liquid and ice, canopy water, snow water, river, and ground water fields.The annual cycle of TWS is captured reasonably in CLM4 [16].As one of the limitations faced by GRACE is its inability to distinguish between these near-surface contributors to the mass budget, the CLM4 offers additional guidance as to how the various land water reservoirs contribute to observed variability.CLM4 does not include the effects of groundwater pumping, irrigation, or dams in its simulation of surface hydrology.
As an example of the fidelity of the CESM in representing internal variability, ENSO teleconnections are explored in Figure 1, where rainfall, near-surface winds, and 500 hPa vertical velocity (w500) are regressed against Niño3.4sea surface temperature (SST) for the LE 2000-yr control run and compared with fields from the European Center for Medium Range Weather Forecasts Interim Reanalysis [17].The fields shown are based on the linear regression of 12month means averaged from July through June (to account for the seasonality of ENSO) after removing any longterm trends.The general character of the observed and simulated patterns is very similar to anomalously strong rainfall in the central equatorial Pacific accompanied by deep convective vertical motion (w500 < 0) surrounded by regions of enhanced subsidence (w500 > 0) and suppressed rainfall, particularly over South America and Australia in a manner that is qualitatively consistent with recent observed variability (e.g., [18]).Extratropical rainfall associations are also generally consistent with those observed as positive regressions prevail over the Mediterranean, the southern tier of the North America, and southern South America.The major region of disagreement resides in Africa and the southern Atlantic Ocean, where regressed wind associations are too weak and rainfall regressions are positively biased in the CESM.

Indices of Variability.
Key to understanding reported changes during the GRACE record is quantifying the character of internal climate modes.Here, indices are employed to characterize ENSO, the AMO, and the PDO in a consistent fashion across observations and models using methods employed in NCAR's Climate Variability Diagnostics Package [14].ENSO state is estimated from sea surface temperature in the Niño3.4region (5 ∘ N-5 ∘ S, 190 ∘ E-120 ∘ W).The AMO index is computed by first removing the monthly global mean SST in ice-free regions (60 ∘ S-60 ∘ N) to separate the Atlantic Multidecadal Oscillation (AMO) from global mean forced changes [19].The index is calculated by averaging detrended annual mean SST anomalies over the North Atlantic Ocean (0 ∘ -60 ∘ N, 80 ∘ -0 ∘ W).This approach differs somewhat from earlier approaches in which the global mean was not considered, a distinction that is of particular relevance for the GRACE period as discussed in Section 3. The PDO index is derived from the standardized leading PC for the first EOF of detrended SST anomalies in the North Pacific Ocean (20 ∘ N-70 ∘ N, 110 ∘ E-80 ∘ W, [20,21]).

Observed Trends and Variability
The GRACE-era TWS trends from 2003 to 2014 based on Mascon solutions are shown in Figure 2, with regions of disagreement in the sign of trends across the various GRACE products indicated by stippling.The trends are characterized by large increases in midlatitude North America, the Amazon and Parana basins of South America, western and southern Africa, eastern Australia, and portions of eastern Eurasia.Large decreases are evident across the northern and southern tiers of North America, Patagonia, western Eurasia, and the Himalaya.In general, the various GRACE products are consistent in depicting trends, with disagreement limited to regions of very small Mascon signals (<2 mm yr −1 ) and western Europe where GRACE retrievals are complicated by challenges in retrieving variability along coastlines.Many of the trends in the GRACE record are associated with significant potentially water cycle driven socioeconomic disruptions, including droughts in the southern United States [22], conflict in western Eurasia (e.g., [23]), and major floods in Australia that impacted sea level worldwide (e.g., [18]).Many of the regions in which pronounced trends are evident in GRACE are also strongly influenced by internal variability, particularly in rainfall [21,[24][25][26].Understanding the evolution of such internal variability is therefore critical for interpreting trends and evolution of ENSO (using the Multivariate ENSO Index, MEI, [27]), the PDO, and AMO is shown in Figure 3. Here, the MEI is used instead of Niño3.4SST as it is a more holistic measure of ENSO, being based on sea level pressure, surface wind, SST, surface air temperature, and cloud fraction, though both mean anomalies in the MEI (0.01) and Niño3.4SST (<0.01 K) are negligible for the GRACE era.As MEI is unavailable for model fields, Niño3.4 is used in assessing LE simulations.In contrast, the PDO and AMO are anomalous for much of the period with strong negative anomalies prevailing over much of the middle and latter part of the record.Of note, however, in regard to the AMO index, is the fact that when estimated by earlier techniques that do not account for global change [24] a positive anomaly is reported, whereas, once corrected for such effects, the methodology of Trenberth and Shea [19] shows it to be modestly negative (−0.33 K; −0.35 ).Here, the latter convention is adopted in light of the strong high latitude warming that occurred during the early 2000s on a planetary scale (e.g., Figure 9 of [7]).

Simulated Internal and Forced Variability
A major advantage of the LE is that it incorporates multiple realizations of the 20th and 21st centuries with a single model.It is therefore possible to accurately estimate the forced response by averaging together the various ensemble members to remove internal variability, the influence of which can be viewed as being random at any given point in time.Here, the forced response in TWS is estimated for the GRACE era and for the full span of the LE (Figure 4), with regions for which the sign of trends lacks consistency across at least three-quarters of the ensemble members which are stippled.The forced GRACE-era trend (Figure 4(a)) is very similar to the overall trend from the beginning of the LE to the present (1920 to 2014, not shown) and also exhibits many of the features of the forced response from 1920 to 2100 (Figure 4(b)), with reductions across western North America and the Amazon and Mediterranean regions and increases at high latitudes, across Australia and central Africa.However, there is substantial spread across ensemble members at regional and zonal mean scales, with the sign of trends in each being highly variable across such a brief period.In contrast, when the full LE record is considered, trends in these regions both intensify and become robust across ensemble members.Zonal mean averages across the full record become particularly robust, with increases at high latitudes and northern tropics (albeit with significant spatial structure) and decreases near 40 ∘ N and especially north of 70 ∘ N.
Given the strong role for internal variability suggested in Figure 3, the spatial structures of major modes over variability are explored (Figure 5) by regressing indices of climate variability against TWS in the LE's 2000-yr control simulation.For these internal variations, the CLM simulates the greatest contribution to TWS from liquid soil water and ground water terms.At high latitudes (>50 ∘ ), ice soil water and snow cover can also plan an important role.For ENSO (Figure 5(a)), strong positive associations are evident across the southern tier of North America, southern South America, central Africa, and northern Eurasia, with deficits in northwestern North America, north and eastern South America, west and east Africa, Southeast Asia, and Australia.The TWS pattern is therefore spatially correlated with rainfall teleconnections (Figure 1), suggesting a strong (and perhaps obvious) physical linkage between the two.Many of the features of ENSO are also evident in the PDO, which again may be expected given their strong interconnectedness [21].Regional associations with the AMO however are in many ways distinct from those of ENSO and the PDO, with positive anomalies in northern South America, the Sahel, and southeast Asia, and negative anomalies spanning much of eastern North America and western Africa.Other features, such as regressions in the southwestern United States, central South America, western Eurasia, and south Africa, are qualitatively similar to ENSO and the PDO.The roles of these features contributing to observed changes will be discussed further in Section 7.

Regional Trends and Variability
The complexity of the spatial patterns in both the forced response and modes of internal variability motivates an analysis of the time evolution of regional structures (Figure 6).For global land, there is a positive TWS trend (∼0.3 mm yr −1 ), associated with an intensified and spatially redistributed water cycle (e.g., [2]), and an increase in internal variability (red dotted line), associated with elevated risks of episodic drought and flooding [3,4].The increase in variability between the early 20th and late 21st centuries is substantial (∼50%).
Regional characteristics can deviate considerably from the global mean.In Australia (Figure 6(b)), a stronger increase in storage (∼0.6 mm yr −1 ) and a weaker increase in variability (∼30%) than for global land are evident.A pronounced Amazonian drying trend (Figure 6(c), −1.2 mm yr −1 ) accompanies an increase in variance (∼35%).Trends that are less monotonic are also evident, such as for the Colorado river basin (Figure 6(d)), where an overall negative TWS trend (0.7 mm yr −1 ) and increase in variance (∼25%) are evident.The Mississippi river basin also dries (0.7 mm yr −1 ) while evolution of its variance is in itself highly variable.In California (Figure 6(f)), mean changes are highly variable while variability increases (∼35%).

Global Budgets
Among the most important indicators of climate change and a key driver of its associated impacts is sea level.Its interpretation however relies critically in being able to distinguish between contributions arising from changes in ocean heat content and mass.TWS is a key part of this budget (e.g., [6,18]).Cazenave et al. [8] use a land model forced by estimated observed rainfall to speculate that the deceleration in sea level observed during the altimeter era has been driven by TWS increases over the past decade, requiring a net increase in TWS of almost 20 mm since 2003, relative to the prior decade.The plausibility of such a large increase is explored here.In addition, it has been hypothesized that an amplified annual cycle in TWS may serve as a useful indicator  of global water cycle intensification [28] and that instruments such as GRACE may be able to detect such effects.The plausibility of these hypotheses is explored here using the LE simulations.
To address these issues, the frequency distribution of 10yr global land TWS trends in the LE is explored (Figure 7).In this computation, the forced response estimated from the 40-member mean is removed from each ensemble member and running 10-yr trends are computed for all members from 1920 to 2100.Consideration is given to the states of ENSO, the PDO, and the AMO and the frequency of occurrence by intensity is then normalized.The mean ( ẋ ), standard deviation (), and skewness () of the various distributions are also shown.Overall, decadal trends in TWS tend to be small with  of 3.1 mm dec −1 and are skewed positively (), likely due to skewness in the rainfall distribution ( ẋ = 0 as the forced change is removed).Perhaps unexpectedly given their pronounced influence on regional budgets (e.g., Figure 5), modes of variability do not tightly constrain the mean global TWS trend, with ẋ being indistinguishable from zero for all modes.There is some suggestion however of increased  for the negative phases of the PDO and AMO, which would be consistent with its expected influence in enhancing terrestrial rainfall floods.However, skewness also increases for El Niño events, for which no such strong expectation exists.Thus, while aspects of the TWS frequency distribution remain to  be understood, the LE demonstrates two key points: mean decadal trends are rarely greater than 6 mm dec −1 and the trends themselves are not strongly determined by the internal modes of variability considered here.
Lastly, the capacity for using the annual amplitude of TWS as an indicator of the strength of the water cycle is explored.The mean annual cycles of global land TWS for the early 20th and late 21st centuries are shown in Figure 8 along with a ± 1 range of interannual variability, where the annual mean (Figure 6(a)) has been removed.While the forms of the annual cycles are similar, the mean boreal summer minimum occurs about a month sooner (Aug.versus Sep.) and is about 1 mm lower in general in the late 21st versus early 20th centuries.These differences are statistically significant given the very large number of years included in the composite (1200).Nonetheless, the expectation that such a shift in the satellite record of a decade or so will be detectable is highly questionable, given the large internal variability inherent to the annual cycle and the demands of such detection related to the absolute accuracy of retrievals.These challenges are suggested by the LE despite a robust increase in terrestrial rainfall in the ensemble, increasing from about 2.24 mm d −1 to 2.44 mm d −1 , or about 9%, from 1920 to 2100.

Discussion and Conclusions
The GRACE record represents a remarkable technological achievement, providing significant and unprecedented insight into climate variability and change as manifested in terrestrial hydrology.However, the brevity of the record makes it subject to the same tensions that have long existed between the relatively short record of high-quality global satellite observations and the typically longer period of time necessary to separate the convoluted influences of forced and internal variability.Here, it is shown that models can play a fundamental role in clarifying these roles, albeit with their own inherent uncertainties.Based on the LE, the forced TWS response is generally overwhelmed by internal variability on decadal timescales, with the AMO and PDO playing particularly influential roles in many key regions.Largely by chance, the negative phases of these modes (Figure 2) and related teleconnections (Figure 4) occurred during the GRACE era, leading to TWS increases in South America and Australia and drying in the southern tier of North and South America and western Eurasia.As these variations are also consistent with the forced response, internal variability has provided what can be reasonably viewed as an accelerated realization of many of the changes eventually anticipated under climate change.For example, this finding also suggests that there is little expectation that reported trends during subsequent satellite missions (e.g., the launch of GRACE Follow On is anticipated in 2017) will necessarily be similar to those observed during GRACE on decadal timescales.Nonetheless, despite this obfuscating role of internal variability, the sustained evolution of the forced response over time is likely to emerge into a clear pattern of response and be associated with significant shifts in terrestrial hydrology that can be either obscured or enhanced by internal variability at any given time, portending significant future impacts.
Examining global budgets, a number of findings also appear to be robust.Of particular note is the fact that the LE suggests that decadal trends in TWS rarely exceed 5 mm and are not strongly influenced by the PDO, AMO, or ENSO.Rather significant, apparently random scatter exists in the distribution of decadal trends.Explanations for recent sea level trends based on TWS (e.g., [8]) increases arising from the influence of the PDO therefore appear to be improbable according to the LE and alternative hypotheses should be explored.Similarly, evidence for the assertion that the seasonal cycle of TWS may act as a useful indicator of future water cycle intensification is found to be weak, as it is undermined by the large monthly variability and small signal that are anticipated as the water cycle strengthens.
In addition to robust global-scale changes, variability on regional scales shows an increase in many regions according to the LE.The realization of these regional changes in coming decades is likely to be much more uncertain than many of the findings above, however, due to structural model uncertainties, potential associated error, and the increased influence of internal variability on regional changes, even on multidecadal timescales (e.g., [29]).
Lastly, GRACE has provided a data record that is likely to yield additional benefits far into the future.The fields have provided developers with a new dataset that they can use to scrutinize and develop models [30], allowing for increased fidelity in simulating the movement and exchanges of moisture at the land surface.Given the importance of water in many terrestrial processes (e.g., carbon, vegetation) and its coupling to the atmosphere, the benefits of these improvements are expected to extend well beyond TWS simulation alone.While many of these improved models have yet to be used for large-scale projects such as the LE, it is only a matter of time before the insights gained and associated legacy of GRACE reach their full potential.

Figure 2 :
Figure 2: GRACE-era trends (2003-2014) in surface liquid water equivalent (mm c −1 ) from JPL Mascon fields [9], stippled where the sign of the trend is inconsistent across available GRACE datasets (see text).

Figure 3 :
Figure 3: Evolution of ENSO (based on the MEI), the PDO, and AMO during the GRACE era.

Figure 4 :
Figure 4: CESM forced trends in TWS from (a) 2003-2014 and (b) 1920-2100.Regions in which the ensemble mean trend is less than the standard deviation among ensemble members are stippled.

Figure 5 :
Figure 5: Regressed contributions to TWS of (a) ENSO as estimated from Niño3.4 SST and (b) the PDO and (c) AMO indices from the LE control run.

Figure 6 :Figure 7 :
Figure 6: Time series of TWS anomalies and their interquartile range (shading) and internal variability based on the spread across ensemble members at each time (red dotted) for (a) global land, (b) Australia, (c) the Amazon, (d) Colorado, and (e) Mississippi river basins and (f) California.A 10-year running mean smoothing is applied to all time series.