Alaskan Permafrost Groundwater Storage Changes Derived from GRACE and Ground Measurements

The Arctic is in transition from climate-driven thawing of permafrost. We investigate satellite-derived water equivalent mass changes, snow water equivalent with in situ measurements of runoff and ground-survey derived geoid models from 1999 through 2009. The Alaskan Arctic coastal plain groundwater storage (including wetland bog, thaw pond and lake) is increasing by 1.15 ± 0.65 km 3 /a (area-average 1.10 ± 0.62 cm/a), and Yukon River watershed groundwater storage is decreasing by 7.44 ± 3.76 km 3 /a (area-average 0.79 ± 0.40 cm/a). Geoid changes show increases within the Arctic coastal region and decreases within the Yukon River watershed. We hypothesize these changes are linked to the development of new predominately closedand possibly open-talik in the continuous permafrost zone under large thaw lakes with increases of lakes and new predominately open-talik and reduction of permafrost extent in the discontinuous and sporadic zones with decreases of thaw lakes.

Other important permafrost watersheds in Eurasia (Lena, Yenisei and Ob') and North America (Mackenzie) are shown [18]. The seasonal maximum and seasonal minimum water equivalent mass change occur in May and September, respectively.  The Alaskan Arctic coastal plain (Figure 1), an area of about 104,000 krn 2 , is noted for its near-flat topography of tundra-wetland ecosystems, numerous thaw lakes and thick, up to 600 m, permafrost in the continuous permafrost zone [19][20][21]. The Barrow Peninsula area is noted for its large and numerous thaw lakes with their distinctive elliptic-shape and orientation [22]. Patterned ground formed by seasonal freeze-thaw and ice-wedge growing processes are common throughout this region. The Arctic coastal foothills ( Figure I), an area of about 120,000 krn 2 forms a low-relief rolling topography incised by numerous rivers. Some of the rivers crossing the central and eastern foothills have a glacial melt-water component in their summer-time discharge [23 ].
The Brooks Range forms the meteorological-climatologic boundary of the Arctic in Alaska [24].
Permafrost in this area is much thinner and may be less than 200 m in many locations [19].
Precipitation within the Alaskan Arctic coastal plain and foothills is on the order of 100 mm in summer. Over the twentieth century, the number of days below freezing has been on the order of 300.
Winter snowfall has been on the order of 80 cm with snow depths on the order of 20 to 50 cm [25,26].
Annual air temperatures recorded at Barrow, Alaska, from 1949 through 2005 had increases on the order of 2 °C [24]. Seasonal air temperature had increases of 0. 8 °C in autumn and 3.2 °C in winter.
Annual precipitation had decreases on the order of 2 cm. Seasonal precipitation had decreases of 1. 3 cm in winter with summer unchanged, over the same period. Interannual and decadal variability affects the long-term trends, causing short interim periods of rapid changes or periods of neither increasing nor decreasing trend. On a monthly basis, peak precipitation of about 2 cm occurs in August.
The Yukon watershed ( Figure 1) covers an area of 853 ,300 km 2 [27]. It is noted for its lowland wetland boreal ecosystems and high-elevation Alpine tundra, and is underlain by the discontinuous and sporadic permafrost zones [21]. Permafrost, where present, thickness ranges from about 15 m and less in the lowland western portions to about 180 m in the northeastern portions [28]. Headwaters of the tributaries of the Yukon River originate in the high-elevations of the Alaska and Brooks Ranges, in Alaska, and the Coast Mountains in Yukon Territory, Canada. Many are fed by summertime glacial-melt and wintertime base flow from groundwater [23 , 27]. It discharges its freshwater, suspended sediment, and dissolved carbon load at the Yukon Delta on the shores of the Bering Sea [29,30] Interannual and decadal variability strongly affects the long-term trends with the largest increase in precipitation occurred during 1976177. Before and after this shift, neither increasing nor decreasing trends are apparent. On a monthly basis, peak precipitation of about 3 cm occurs in July. Since the early 1990s the NOAA National Geodetic Survey (NGS) have been working on developments of high-spatial resolution and high-accuracy geoid models as part of the US national mapping and geodetic control network modern ization programs [32,33]. The efforts centered on producing geoid models using ground-based gravimetric geoid observations (within the US and territorial coastal w aters) at the US High-Accuracy Reference Network sites with Global Position System station coordinate determinations with long-wavelength undulations from the Earth Gravity Model 1996 (EGM96) a 360-degree and -order 15-arc-minute model and realization with respect to the International Terrestrial Reference Frame (ITR F) in the non-tide geodetic reference convention [34].
The first model was GEOID 1996 (GEOID96) a 2-arc-minute grid. It was followed in 1999 by the l-arc-minute model GEOID99. Models GEOID03 and GEOID06 followed this. The latest efforts have produced GEOID09 that incorporates EGM2008, a 2,159-degree and -order l-arc-minute model.
GEOID99 and GEOID09 are the principle geoid models we use. The geoid models are relative to the WGS84 reference ellipsoid non-tide system in the convention of the International Terrestrial Reference Frame 2005.

Adjustment from GPS Stations
Adjustment for post-glacial and neotectonic (isostatic) changes were implemented by using Global In our study we used the Release-04 (R4) Level-3 grids provided by the GRACE Science Team centers. The global grids are I-arc-degree water equivalent mass change complete to degree and order 40 (ftp:llpodaac.jpl.nasa.gov/pub/tellus/monthly_mass_gridsl). The GRACE solution to the gravity potential formulated as an equivalent water mass change (length scale) is given by LlClm(t), and LlSlm(t): Normalized time-varying Stokes spherical harmonic geopotential coefficients, ae: Earth mean radius, r: spatial radius, k1: Love numbers, Pe: Earth mean density, pw: water density, t: time, and $, A are latitude an d longitude [38]. (1) Beyond degree (order) 40 to 70, the inherent noise level in the mass change signal becomes significant [40]. Processing includes downward propagation and adjustments to remove the time-variable mass change effects from tides and atmosphere and mean variation (GRACE geoid model). A normalized Gaussian smoother filter removes striping artifacts produced by the orbit non-crossing and control-descent geometry [38][39][40]. Differences in processing (de-aliasing) and error sources and products are attributable to differences in assumed zero-degree and order Stokes harmonics, tide (liquid and solid) models and the modeled atmosphere mass change removal, respectively in decreasing order of magnitude [41,42].
The GRACE R4 Level-3 (300 km smoothing kernel) land and ocean monthly grids from August 2002 through December 2008 were combined to give global coverage. Glacial isostatic adjustment (GIA), solid mass flow within the Earth's mantle following the decay of the Pleistocene Laurentide and Cordilleran ice sheets in North America and Euro-Scandinavia is an adjustment we applied to the GRACE grids. We implemented this adjustment using a GIA grid that uses the ICE-5G (VM2) provided by the GRACE science team [43,44].

Snow Water Equivalent
Global snow water equivalent estimates were derived from the NOAA Defense Meteorological Satellite Program by the Special Scanning Microwavellmager (SSM/I) sensor [45]. We use data provided by the National Snow and Ice Data Center (NSIDC, http: //nsidc.orgldatalnsidc-0271.html).
Snow water equivalent estimates in units of millimeters were derived using the horizontally polarized Under-estimations can occur on extreme elevation mountains of complex geometry, along marine-coastlines where signals from ocean bodies are partly convolved, an d the beginning and ending of the snow season when snow depths are thin and moist [47,48]. Noise from transitory atmosphere conditions is removed using a five-day filter. Monthly composite climatology is processed from daily passed-retrievals by averaging. Processed data are gridded in the equal-area scalable Earth (EASE) projection system at 25-km grid intervals.

Wa tershed Runoff
Measurements of surface water discharge (runoff) at gauge stations in the watersheds are provided Kuparuk and Sagavanirktok discharges were summed to give a single monthly time series for the Arctic coastal plain and foothills. We then derived annual time series of runoff by summing the monthly runoff from August through July periods. August through July was chosen as this period coincided with the GRACE time series instead of the typical water-year calendar convention.

Method
We derive area longitude-latitude extents of the Arctic coastal plain, foothills and the Yukon River watershed. The area-extents are used to extract GRACE and SSMII data to compose regionalized time series of the mean, standard deviation and standard error. We compare these with monthly runoff for timing of winter snow water mass loading and spring surface water unloading. We then derive optimal least squares trends for comparison of trend gradients.

1. Theory a/ Water Mass Changes
On a global basis GR ACE measures variations of the gravity field corresponding to the total change of mass from multiple sources [51,52]. On terrestrial non-glacier regions the mass changes will include those from snow, soil and vegetation moisture, lake storage, ground-ice, runoff and groundwater storage. Previously we found that the regionalized gradients of soil and vegetation moisture did not make signifi cant contributions to regionalized total water equivalent storage change on the Arctic permafrost watersheds [18].
In our formalism the regionalized gradient of groundwater storage change includes second-order regionalized gradients of storage changes in lakes, active layers, hyporheic zones and ground-ice bodies. Therefore, Equation (3) forms our theory of measurement. We estimate uncertainties from the trends and formal errors trough the error propagation equation [53]. In this respect the estimate of the uncertainty of groundwater storage change gradient is the root sum of the squares of the uncertainties of GRACE total water storage change gradient and the area-average runoff gradient.
In our formulation there exist the Hilbert Space n, where the functions T, G and R exist and have sub-spaces T, y, and p [54). T, G and R and the sub-spaces are linearly independent. Equation (3) formalizes a reduction of LJ.T by LJR in n to derive LJ.G. This is not the so-called hydrologic balance. In applications, our formalism is the same philosophy as employed in Geopotential Theory, also in Harmonic Analysis, to reduce the total potential (field) by the strongest component potentials (fi elds) to derive an unknown component potential (field) [55]. As such the derived potential (field) LJ.G may contain potentials (fields) of third and lesser order.
We use geoid models GEOID99 and GEOID09 to compute a geoid-difference model. This model is adjusted by using a trend surface derived from GPS vertical rates to remove isostatic effects in the period from 1999 through 2009. In theory, on short time scales and with removal of atmosphere, tide and solid-Earth isostatic effects, geoid changes should reflect changes of surface and groundwater storage and flows. We note that the models GEOID99 and GEOID09 are not completely independent.
Differencing them will remove common features and present a net change of the Earth's geoid in Alaska and western Canada from 1999 through 2009.

1. Regionalized Time Series and Gradients a/ Water Equivalent Mass Changes
Regionalized time series of GRACE, SSMII, and in situ measured runoff show monthly water equivalent mass variations within Alaskan Arctic coastal plain, foothills and the yukon River watershed regions (Figure 3). Gradient and uncertainties of water equivalent mass change in units of thickness (cm) water equivalent mass change are given in Table 1 and graphically in Figure 4.
Runoff trends showed geographic variation north and south of the Brooks Range ( Figure 4 and Table I In all cases the SSMII trends are discordant with those from GRACE and annual runoff ( Figure 4 and Reduction of the GRACE total water equivalent mass change gradient by the runoff gradient yields an estimate of groundwater storage change (see section 4. 1). The results indicate the Arctic coastal plain is gaining water equivalent mass ( Figure 4 and Table 2). This contrasts with the Arctic foothills and the Yukon River watershed, which are losing water equivalent mass.

Geoid-Diff erence Model Changes
The geoid-difference model (   The active layer, the near-surface soil-earth layer that seasonally freezes and thaws is another factor in permafrost watersheds [59]. Changes in the active layer affect soil moisture and groundwater, ecological process and energy fluxes between the atmosphere and land (3]. Dynamics of permafrost, talik and active layer exert significant influences on surface and groundwater hydrology and geomorphology in the Arctic, and they are spatially non-uniform and nonlinear.

1. Alaskan Arctic Coastal Plain
The Arctic coastal plain of Alaska is a region of numerous and some quite large thennokarst thaw ponds, lakes and drained thaw lake basins [22,62J. These lakes reside in the thick ice-rich continuous permafrost zone. Their development and growth come by way of thermal processes at ice-wedge intersections and permafrost thawing with subsidence of the ground surface and increased thickness of the active layer [63]. Their elliptical shapes develop over time from thenno-erosion and wind-oriented waves that produce circulation cells. Over time, small ponds merge into large lakes, and some will drain when thermo-erosion and bank undercutting cause a breach [22,62). Beneath the lake a talik (unfrozen soil) develops (closed talik) [59]. The talik's lateral and vertical extent depend on the internal structure of the permafrost and on the temperature regime at the ground surface. Over time talik becomes an open talik and serves as an aquifer whereby lake water drains to become subsurface groundwater (increased groundwater storage) within or in some cases below the permafrost (where it is relatively thin on the order of lO s of meters). This indicates surface water is being recruited for subsurface groundwater storage (including ice growth) and groundwater residence time is increasing.
Additional groundwater storage changes (increasing) because of increasing numbers and area of thaw bogs, ponds and lakes. Independent research indicates lakes (including thaw bogs and ponds) have been increasing in the continuous permafrost zone in Eurasia and Alaska [17,63]. On the northern Remote Sens. 2011, 3

390
Alaskan Arctic coastal plain there is similar evidence of thaw lakes increasing, lateral spillage channels (outflow) forming on lake margins, and new pingos also appear (Guido Grosse, Personal Communication, 2009). Recent evidence also point to increased thaw penetration, which increases active layer thickness on the Arctic coastal plain of Alaska [64]. In related research we find that GRACE-derived groundwater storage is increasing in the Lena and Yenisei River watershed regions of Eurasia [18]. These watersheds encompass large areas of the continuous and discontinuous permafrost zones and talik [65).

Yukon River Wa tershed
An independent study of Yukon River watershed runoff indicates increasing groundwater (subsurface storage) contribution to streamflow (i. e. , groundwater storage loss) that currently maintains annual discharge [27). This study pointed, in general, to permafrost thawing during the last century climate warming as having enhanced deeper flowpaths. We hypothesize that permafrost degradation is accompanied by development of taliks (which act as subsurface flowpaths) allowing intra-permafrost groundwater (either held trapped or from thawing of ice rich layers) to enter streamflow thus offsetting a negative annual net precipitation (i.e., P-ET) in annual runoff. Furthermore, independent research indicates decreasing lake numbers and areas in the discontinuous permafrost zone [17,59,66].
Therefore, a component of the groundwater loss can be attributed to thaw lake shrinkage ( drainage) in the Yukon River watershed of Alaska and Yukon Territory. In related research we found that GRACE-derived groundwater storage is decreasing in the Mackenzie River watershed of Canada, and groundwater storage is stable, i.e., unchanged, in the Ob' River watershed of Eurasia [18). The Mackenzie watershed is well noted for its predominately discontinuous permafrost extent and the Ob' River watershed has relative little extent of discontinuous permafrost [67).

Hypothesis
We hypothesize the changes of groundwater storage are associated with changes of surface water mass in thaw bogs, ponds, and lakes, and in subsurface water mass, possible including ice, in permafrost zones. Furthermore, we hypothesize the water mass changes are linked to changes in the character and distribution of permafrost, talik, and the active layer: These physical changes in permafrost, talik and active layer serve to seasonally increase surface water storage on the Barrow Peninsula region and also to remove a portion of runoff for recruitment into subsurface groundwater storage.
Within the Yukon River watershed, the lateral reduction of permafrost, development of new open-taliks beneath thaw bog, pond, lake and river bed allows for the removal of surface water storage (thaw ponds and lakes) into groundwater storage and a portion of the groundwater storage is being recruited into runoff. Our findings agree with observations that thaw lakes have increased in the continuous permafrost zone and have decreased in the discontinuous permafrost zone [17,59,63,66].
The changes of permafrost, talik and active layer will affect groundwater residence times. On the Arctic coastal plain of Alaska groundwater residence time will increase as surface water is recruited into groundwater storage. On the Yukon River watershed groundwater residence time will decrease as groundwater leaves storage and becomes part of the river runoff (i. e., discharge). This effect is in part due to the increased presence of taliks on watersheds and beneath riverbeds and serves to enhance exchan ge (both recruitment and export) of shallow groundwater and river water [67].
The character and changes of permafrost is strongly dependent on surface air temperature and snow cover [3,4,16]. Permafrost warming in Alaska was coincident with warming of air temperatures that began most strongly in 1976177 in a step-like function. The Arctic coastal plain air temperatures warmed by up by 2 to 3 cC, on average, and the Yukon River watershed air temperatures warmed up by 1 CC. Surface thawing of permafrost on tundra and forest sites in Alaska was about 0.1 m/a at some locations. Basal thawing at one site was about 0.04 m/a, and accelerated to 0.09 m/a [16]. Talik provides a path for surface water to enter permafrost from above, and for groundwater to join river runoff as base flow.
On longer time scales, significant changes of permafrost and linked responses to climate and ecosystems (water cycle and carbon stocks) have occurred [6,7,10,68]. Modeling of the regional permafrost environments with general circulation model scenarios suggests significant changes in permafrost distribution are likely over the next century [69,70].

Conclusions
Analysis of GRACE and runoff data allow for estimates of the changes in groundwater storage. Our results indicates from 2002 through 2008: 1. Alaska Arctic coastal plain groundwater storage is increasing at 1.15 ± 0.65 km 3 /a equivalent to an area-average is equivalent thickness rate of 1.10 ± 0.62 cm/a.