Introduction
The land-ice component of the cryosphere has important connections to the global climate system. Changes in land-ice cover impact the Earth’s radiation budget through variations in albedo, and significantly impact the hydrological cycle as changes in water storage. The Earth’s large areas of land ice, including the ice sheets, glaciers and ice caps, have the potential to significantly impact eustatic sea level, as they hold ∼77% of the world’s fresh water, which amounts to ∼65 m of global sea-level equivalent (Reference Barry and GanBarry and Gan, 2011). Changes in freshwater input from calving and melting Arctic land and sea ice can increase the haline forcing in the North Atlantic and impact global ocean circulation and climate patterns (Reference Häkkinen and RhinesHäkkinen and Rhines, 2004). Variations in land-ice volume are closely correlated to global temperatures, as paleoclimate records show all the major Greenland ice sheet (GIS) changes have been correlated with temperature changes (Reference AlleyAlley and others, 2010). It is important that we accurately quantify the current spatial and temporal mass changes of the Earth’s land ice to improve our understanding, modeling and prediction of its response and contribution to climate change. Since its launch in March 2002, the Gravity Recovery and Climate Experiment (GRACE) mission has acquired ultra-precise inter-satellite K-band range and range-rate (KBRR) measurements between two co-orbiting satellites in a 450 km altitude polar orbit, ∼220 km apart (Reference Tapley, Bettadpur, Ries, Thompson and WatkinsTapley and others, 2004). These data have vastly improved our knowledge of the Earth’s time-variable gravity field. In particular the GRACE mission has been instrumental in quantifying recent land-ice mass evolution (e.g. Reference LuthckeLuthcke and others, 2006a, Reference Luthcke, Arendt, Rowlands, McCarthy and Larsen2008; Reference VelicognaVelicogna, 2009; Reference Chen, Wilson and TapleyChen and others, 2011; Reference Schrama and WoutersSchrama and Wouters, 2011; Reference Jacob, Wahr, Pfeffer and SwensonJacob and others, 2012; Reference King, Bingham, Moore, Whitehouse, Bentley and MilneKing and others, 2012).
There are numerous challenges in applying GRACE to the study of glacier and ice-sheet mass balances. GRACE senses all sources of mass variability at the Earth’s surface, requiring removal of non-glacier changes through incorporation of independent datasets and models. Mass signals associated with glacial isostatic adjustments (GIA) and terrestrial water storage (TWS) are poorly known in many regions and often dominate GRACE error budgets. Owing to the GRACE mission’s limited longitudinal sampling at monthly time intervals, filtering techniques are often applied to mitigate striping errors and to focus signal into a region of interest, often resulting in signal attenuation. When mass is constrained to a particular spatial extent, signal can either leak into or out of that domain and cause additional error. Even when the choice of filtering and models used to correct atmosphere and ocean mass variations is fixed, differences in GRACE solutions can occur, due to different methods used in processing the fundamental GRACE observations (Reference Pritchard, Luthcke and FlemingPritchard and others, 2010).
Here we present a new study to quantify the ice mass evolution of the Antarctic and Greenland ice sheets (AIS and GIS) and Gulf of Alaska (GOA) glaciers from a GRACE global mascon solution. The mascon technique is a time/space domain approach for describing time-variable gravity signals, in which time is discretized in contiguous, non-overlapping finite intervals and, within a given interval, space is divided into gravity contributions due to loading of mass concentrations of equivalent water (mascons) on the Earth surface. This study is unique and differs in many aspects from previous regional and global mascon solutions, including those of Reference LuthckeLuthcke and others (2006a, Reference Luthcke, Arendt, Rowlands, McCarthy and Larsen2008), Reference RowlandsRowlands and others (2010) and Reference Sabaka, Rowlands, Luthcke and BoySabaka and others (2010). The mascons are estimated globally with 10 day temporal and 1 arcdeg equal-area spatial sampling. Ice-sheet mascons are restricted to the grounded-ice areas contributing to eustatic sea-level change. We apply temporal and spatial constraints anisotropically, through seven distinct constraint regions representing the ice sheets (high-interior and low-elevation margins), GOA glaciers, land and ocean. The constraint equation is an exponential function of correlation time and distance applied to all possible pairs of mascons within a constraint region (Reference Sabaka, Rowlands, Luthcke and BoySabaka and others, 2010). A detailed error analysis is provided, including consideration of leakage, and errors in modeling mass variations of the atmosphere and ocean during the GRACE data reduction (forward modeling). In addition, a controlled comparison is performed in which we use the same set of forward modeling and GRACE data reduction, but compute the mass time series using an averaging kernel filter applied to monthly spherical harmonic fields, as is common practice in the literature (e.g. Reference Velicogna and WahrVelicogna and Wahr, 2005). The mascon solution is iterated until convergence, and the impact of the iterated solution on the reduction of the GRACE KBRR residuals is quantified. Furthermore, this study differs from the mascon solution presented by Reference Jacob, Wahr, Pfeffer and SwensonJacob and others (2012), in that we estimate the global mascons directly from the GRACE KBRR data, taking into account the full noise covariance, and iterate the solution while quantifying the minimization of the GRACE KBRR observation residuals. We present the details of the global mascon solution, an error analysis and a discussion of results focused on recent GIS, AIS and GOA glaciers land-ice evolution. This paper also serves to document our mascon solution used in the international Ice Mass Balance Inter-comparison Exercise (IMBIE; Reference ShepherdShepherd and others, 2012).
Mascon Formulation
The GRACE processing centers provide time-variable gravity observations (Level-2 products) in the form of monthly sets of Stokes coefficients or spherical harmonic fields. These monthly fields are estimated directly from GRACE tracking data, typically without the use of constraints or a priori information. The use of these fields for scientific applications, such as measuring the mass evolution of the GIS or hydrological basins, requires filtering due to north–south striping artifacts, which arise from the limited longitudinal sampling of the GRACE orbits within monthly time periods (Reference Swenson and WahrSwenson and Wahr, 2006; Reference KleesKlees and others, 2008). Many filtering methods do not take into account the noise covariance of the estimated Stokes coefficients, and therefore are not optimal filters (Reference KleesKlees and others, 2008). Loss of signal amplitude and limited spatial and temporal resolution, as well as signal leakage in and out of regions of interest, are particular problems when applying filtering methods.
Instead of applying a filtering technique to the GRACE Level-2 monthly Stokes coefficients, we used a mascon technique, applying geolocatable anisotropic constraints to estimate the global mass change directly from GRACE KBRR data, which accounts for the full Stokes noise covariance.
The formulation for mascon parameters is derived from the fact that a change in the gravitational potential caused by adding a small uniform layer of mass over a region at an epoch t can be represented as a set of (differential) potential coefficients which can be added to the mean field. The delta coefficients can be computed as (Reference Chao, O’Connor, Chang, Hall and FosterChao and others, 1987)
where 
is the loading Love number of degree l, to account for the Earth’s elastic yielding which, in general, counteracts the additional surface density (Reference FarrellFarrell, 1972); λ, φ and R are geographic latitude, geographic longitude and the Earth’s mean semi-major axis, respectively; M is the Earth’s mass; H(t) is equivalent water height in centimeters of the mass of the layer over a unit of surface area at the epoch t; and Ω is the solid angle surface area of the mascon region where H(t) is applied, dΩ = cos λ dλ dφ.
The water height, H(t) (cm), can be equated to a surface mass density, σ(t) (kg m−2), by noting that water has a volume density of 1000 kg m−3, which gives σ(t) = 10 × H(t). The estimated mascon parameter, Hj (t), for each mascon region, j, is a scale factor on the set of differential Stokes coefficients for that mascon region, giving a surface mass change (cmw.e.). In matrix notation, an assemblage over j of the above equations provides a set of partial derivatives, L, of the differential Stokes coefficients, Δs, (complete to degree and order 60) with respect to equivalent water height, Δh, such that
The mascon parameters are part of nonlinear functions of the GRACE KBRR measurements and are therefore estimated using the Gauss–Newton (GN) method (Reference Seber and WildSeber and Wild, 1989) of nonlinear least-squares. This method seeks a series of differential corrections to the parameters via a linearization of the observation equations at the current state. This means that the changes in KBRR observations at the current state due to changes in mascon parameters are needed for the estimation. The partial derivatives of the GRACE KBRR observations with respect to the mascon parameters are computed by the chain rule, using the partial derivatives of the KBRR observations with respect to standard Stokes coefficients and the differential Stokes coefficients computed from the equation above for the mascon region of interest:
where ∂Oi
/∂Hj
(t) is the partial derivative of KBRR observation i with respect to the j mascon parameter, Hj
(t), at time t; 
and 
are the partial derivatives of the KBRR observations with respect to the geopotential Stokes coefficients (computed as part of the KBRR reduction and Level-1 data processing); and 
and 
are the partial derivatives of the delta Stokes coefficients with respect to the mascon parameter in region j at time t.
In matrix notation, assemblages over i and j of the above equation provide a set of partial derivatives of the KBRR observations, o, with respect to the mascon parameters, h, by chaining L with partial derivatives, A, of the KBRR observations with respect to the Stokes coefficients.
Following Reference Sabaka, Rowlands, Luthcke and BoySabaka and others (2010), at step k of the GN estimation let the GRACE KBRR residuals, r, be defined as the difference between the observations, o, and a model prediction, a(x
k
), which is a function of a current state of a general parameter set, x
k
. In addition to the current mascon parameter state, this parameter set represents the effects of orbital arc parameters (initial baseline state orbit and accelerometer calibration parameters) and the mean field and forward models of other geophysical processes of mass redistribution (e.g. atmosphere and ocean tides). The details of the orbit arc parameters and forward models are discussed below, in the data processing section. The correction, 
, to the current state of mascon parameters, 
, is then defined as the unique minimizer of the quadratic cost function
where the first term is the weighted misfit and the second term is a measure of mascon complexity, which in this case is a smoothing function that tries to drive the mascon state to zero. The data weight matrix, W, accounts for both measurement noise and the effects of the orbital arc parameters (back substituted), and P
hh
is a mascon regularization matrix, whose influence is controlled by the damping parameter, μ. The updated mascon state, 
, is given by
assuming a step size of unity.
The cost function of Eqn (5) suggests a major philosophical difference between this study and past studies, in that multiple iterations of GN are actually carried out. For instance, Reference RowlandsRowlands and others (2010) and Reference Sabaka, Rowlands, Luthcke and BoySabaka and others (2010) carry out a single step of GN, in which 
is zero in Eqn (5). In this study, the step k update, xk
+1
, which reflects 
, is fed back into the KBRR processing in order to carry out the next iteration. In this sense, the minimum of J(h) is sought instead of the local quadratic approximation to it at 
, i.e. J(Δh
k
). It will be shown below that this leads to a significant improvement in the recovered mascon solutions, particularly in the ability to restore more signal strength.
There are several ways to rewrite the GN iteration k that illuminate different properties of the estimator. For instance, the filtering nature of the estimator and its error sources can be seen by first expressing Eqn (6) as
and then assuming the residuals can be expressed as a first-order linear perturbation, Ah
k
, about a true mascon state, 
, with additive noise, e, such that r = ALΔh
k
+ e. This gives
where K = (LTATWAL + μP
hh
)−1 LTATW and R = KAL. The estimate, 
, differs from the true state, h
k
+1 = h
k
+ Δh
k
, due to two effects: (1) the degradation of resolving power in the filter, due to the presence of regularization, which is manifested as a deviation of the resolution matrix, R, from an identity matrix in the first term of Eqn (8), and (2) the presence of additive noise from the second term in Eqn (8). The effects of the resolution will be discussed in detail below.
Another way to rewrite the GN iteration k is
This is now in the form of an anisotropic, non-symmetric (ANS) filter (Reference KuscheKusche, 2007) that acts upon an unfiltered Stokes state, 
, to produce a filtered Stokes state, 
,
followed by a least-square co-location (LSC) prediction (Reference MoritzMoritz, 1980) of the mascon state, 
, from the filtered Stokes state
where N−1 and 
are the noise covariance and signal auto-covariance matrices of the Stokes coefficients, respectively, and 
is the signal cross-covariance matrix between the mascons and Stokes coefficients. The ability to process the KBRR tracking data directly allows for access to the actual noise covariance matrix of the Stokes coefficients, and thus makes this an example of an optimal ANS filter (Reference KuscheKusche, 2007). Note that the Stokes signal covariances are related to the mascon signal auto-covariance through simple propagation, such that 
and 
. Currently, the mascon signal auto-covariance matrix reflects a first-difference operation between all distinct pairs of mascons in space and time, weighted as an exponential function of temporal and spatial distance between them, and a conservation of mass constraint that keeps the global mascon sum constant over time. Because of the regional aspect of the spatial weighting, the gravity signal is treated as a non-stationary, anisotropic process based upon geolocatable physical properties. Details of the constraint regions, mascon signal auto-covariance construction and selection of the damping parameter are given in the following section.
Grace Data Processing, and Mascon Solution Details
Data processing and forward models
We estimated global surface mass anomalies, or mascons, directly from the reduction of the GRACE Level-1B data, which includes the inter-satellite K-band range-rate (KBRR) data, as well as orbit position, attitude and accelerometer data for each GRACE satellite. We employed a baseline state parameterization (Reference Rowlands, Ray, Chinn and LemoineRowlands and others, 2002) and accelerometer calibration (Reference Luthcke, Rowlands, Lemoine, Klosko, Chinn and McCarthyLuthcke and others, 2006b), which enables the estimation of surface mass variability from GRACE KBRR data alone. The GPS data are used solely to establish an accurate orbital reference and for calibrating accelerometers (Reference Luthcke, Rowlands, Lemoine, Klosko, Chinn and McCarthyLuthcke and others, 2006b). The GRACE KBRR and Level-1B attitude and accelerometer data (including accelerometer calibration) are processed in daily arcs as outlined by Reference Luthcke, Rowlands, Lemoine, Klosko, Chinn and McCarthyLuthcke and others (2006b). Forward models for both static and time-varying gravity are applied in the KBRR data reduction daily arcs, in order to isolate land-ice signal. Static mean gravity is modeled using the complete GGM02C GRACE gravity model through degree and order (d/o) 150 (Reference Tapley, Bettadpur, Ries, Thompson and WatkinsTapley and others, 2004). Ocean tides are modeled with the GOT4.7 tide model (Reference RayRay, 1999; Reference Ray and PonteRay and Ponte, 2003). Atmospheric mass variations are modeled to d/o 90 at 3 hour intervals using European Centre for Medium-Range Weather Forecasts (ECMWF) operational pressure grids. Ocean mass variations are modeled to d/o 90 from the baroclinic Ocean Model for Circulation and Tides (OMCT). Terrestrial water storage (TWS) variations are modeled using the Global Land Data Assimilation System (GLDAS)/Noah Land Surface Model (Noah) 0.25°, 3 hour data, represented as potential coefficients to d/o 90 (Reference EkEk and others, 2003; Reference RodellRodell and others, 2004). The TWS contribution was set to zero for the major mountain glacier areas of the globe using a 0.25° mask (Reference Raup, Kieffer, Hare and KargelRaup and others, 2000) and also set to zero for the GIS and AIS. Glacial isostatic adjustments are modeled based on ICE-5G (VM2) (Reference PeltierPeltier, 2004) and are represented as potential coefficients to d/o 90. Table 1 provides a summary of the modeling applied in the reduction of the KBRR residuals for each of the mascon solutions noted.
Table 1. Summary of forward modeling applied for each global mascon solution
GIA and geocenter
We modified the mascon solution to restore the forward-modeled global GIA, while more recent GIA models are then used in the computation of the mascon solution results presented. For both the GIS and GOA, a GIA model based on the ICE-5G deglaciation history and an incompressible two-layer approximation to the VM2 viscosity profile are used to correct the GRACE mascon solution (Reference Paulson, Zhong and WahrPaulson and others, 2007, computed and provided by A. Geruo, 2011). For the AIS, the IJ05_R2 regional GIA model is used (Ivins and others, in press). We correct our mascon solutions for the viscous component of post-Little Ice Age (post-LIA) GIA following the collapse of the Glacier Bay Icefield. As in Reference Luthcke, Arendt, Rowlands, McCarthy and LarsenLuthcke and others (2008), we apply a regional post-LIA GIA model that is rigidly constrained by ∼100 precise GPS and relative sea-level (RSL) observations of uplift (Reference Larsen, Motyka, Freymueller, Echelmeyer and IvinsLarsen and others, 2005). We further correct our mascon solutions to reflect surface mass variability in the Earth’s center-of-figure frame, rather than the center-of-mass frame in which the solutions are performed. We apply the geocenter correction used in the IMBIE study (Reference ShepherdShepherd and others, 2012) derived from the degree 1 Stokes coefficients determined from Reference Swenson, Chambers and WahrSwenson and others (2008). The geocenter correction accounts for ∼1% of the GIS and GOA ice trends, and <3% for the West AIS and AIS peninsula ice trends. The largest impact of the geocenter correction is for East AIS, at 17% of the trend.
Global mascon definitions
We estimate a global set of 1 × 1 arcdeg equal-area (∼12 390 km2) mascons at a 10 day temporal sampling. The choice of the spatial sampling is a compromise between approximating regional shapes and boundaries (e.g. land, ocean, land ice) and the total number of parameters that can be reasonably iterated for a global solution. There are a total of 41 168 mascons estimated every 10 days, representing a significantly smaller spatial sampling than the fundamental spatial resolution of GRACE (Reference Luthcke, Rowlands, Lemoine, Klosko, Chinn and McCarthyLuthcke and others, 2006b). Furthermore, the number of mascon parameters is significantly larger than the number of d/o 60 Stokes coefficients forming our basis functions (3718). Therefore, as discussed in detail above, we employ a regularization in the form of the mascon signal auto-covariance matrix (Reference Sabaka, Rowlands, Luthcke and BoySabaka and others, 2010). This matrix is constructed from constraint equations which require all mascon differences to be close to zero in a statistical sense. The constraints are applied in groups that pertain to geographical regions, so if two mascons reside within the same region, the constraint weight is a simple exponential function of the distance (characteristic length) and time (characteristic time) between the two mascons. If two mascons reside in different regions, the weight is zero and the constraint vanishes. The details of the constraints and the construction of the regularization matrix, P hh , are given by Reference Sabaka, Rowlands, Luthcke and BoySabaka and others (2010, their section 2.2). The important implementation parameters are the characteristic length and time used in the exponential taper weighting of the mascon first-difference constraint equations, and the damping parameter, which controls the influence of the regularization relative to the data misfit.
Application of regional constraints
The constraints are then applied anisotropically using the following seven constraint regions: (1) GIS margins defined as <2000 m elevation, (2) GIS >2000 m elevation, (3) AIS <2000 m elevation, (4) AIS >2000 m elevation, (5) GOA glacier area, (6) land (including all other ice regions) and (7) ocean and large bodies of water (including ice shelves). In defining the mascons for land and large bodies of water, a minimum enclosed region area of 80 000 km2 was used. Thus, if a land area enclosed within the ocean region is <80 000 km2, then those land mascons were labeled as ocean (similarly for water regions enclosed within land). Figure 1 shows the GIS, AIS and GOA mascons as well as the constraint region boundaries. Also shown in this figure are ice-sheet drainage system divides based on those of Reference Zwally and GiovinettoZwally and Giovinetto (2011). These divides are used for computing mass change as the sum of the mascon time series within a drainage system, but are not used in the mascon estimation.
Fig. 1. Land-ice mascon configurations and drainage systems. (a) The GIS. The numeric labels describe the primary drainage systems and subregions with the first and second digits respectively, and follow definitions based on those of Reference Zwally and GiovinettoZwally and Giovinetto (2011). The yellow border delineates the 2000 m elevation cut-off used in constraining the mascon solution. (b) The AIS. The drainage systems are sequentially labeled from 1 to 36. The 2000 m elevation line defining the constraint regions is again delineated by a yellow border. The three Antarctic regions of EAIS (drainage systems 3–15, 26–36), WAIS (drainage systems 1, 2, 16–22) and AIS peninsula (drainage systems 23–25) are outlined in red. (Note that the yellow elevation line is concealed by the red line along a portion of the boundary between EAIS and WAIS.) (c) The GOA mascons. This region is not further delineated by drainage systems or constraint regions.
We use a characteristic length of 100 km and a characteristic time of 10 days in computing the exponential constraints for the mascon signal auto-covariance matrix. These values represent approximately one parameter sample in time and space. We use a damping parameter (μ in Eqn (5)) of 30 000−1, which provides a similar solution contribution from the signal and data covariance to that determined from an analysis of the GIS, AIS and GOA mascon matrix diagonals. We investigated the impact of the damping parameter by performing solutions using values that span an order-of-magnitude variation, centered on our selected damping parameter value. We find that the impact on the estimated mass balance and annual amplitude is <8% on average over the GIS, AIS and GOA mascons, while the 1σ error estimates of the 10 day mascon solutions can change by as much as 31% on average over those regions. Therefore, varying the signal covariance damping parameter has a larger impact on the solution noise than the signal.
Resolution and leakage
We further analyze the impact of our signal covariance (including our choice of the damping parameter and characteristic length and time) by performing signal leakage analysis within and across constraint regions using the resolution matrix identified in Eqn (8). First, we look at the influence of the signal covariance on individual mascons. Several GIS and AIS mascons were chosen as single-source mascons. For each single-source test we fill the mascon parameter vector with a time series of mass change for that mascon only (we use the time series from our final solution) and zero all other mascons in space and time. The resolution matrix then multiplies the vector of mascons in time and space, where only one mascon contains a time series of mass change. Therefore, before application of the resolution matrix, only one mascon out of all the global mascons contains a nonzero time series of mass change. However, after application of the resolution matrix, several neighboring mascons will contain nonzero time series, as the signal covariance spreads out the signal from the source mascon. The effects of the signal covariance are quantified by differencing the full mascon vector before and after the application of the resolution matrix. The results of several tests show that the spatial signal from a single mascon is ∼100% accounted for by adding all the mascons within a 600 km radius of the original source mascon. Therefore, a single mascon has negligible influence at distances greater than 600 km. This is equivalent to a Gaussian spatial smoothing filter with a 300 km radius. However, the constraint region boundaries modify the spatial pattern of signal spreading as the signal covariance, by design, limits signal leakage across these boundaries. The spatial leakage pattern of a single mascon’s signal is a function of its location with respect to the constraint region boundaries. An error analysis of signal leakage across constraint boundaries is performed for each region considered. The regional leakage analysis and its results are discussed further below. In summary, the spatial resolution of our solutions is equivalent to a 300 km Gaussian spatial smoothing filter within constraint regions, but with significantly reduced leakage across constraint region boundaries.
Mascon Solution Iteration and KBRR Residual Analysis
As described above, we apply a mascon technique to estimate global mass change directly from the GRACE KBRR observations themselves, taking into account the full Stokes noise covariance and a signal covariance. Therefore, we can use the reduction in the GRACE KBRR residuals as a quantitative objective measure of performance. The KBRR residuals are the difference between the GRACE observed KBRR data and KBRR data computed from our forward models and mascon solution using the NASA Goddard Space Flight Center (GSFC) GEODYN precision orbit determination and geodetic parameter estimation software. Unique to this study, we present the improvement in the KBRR residuals from fully iterating the global mascon solution (as in Eqn (6)). Because there are >11 × 106 mascon parameters, we obtain the solution using iterative methods. Following Reference Sabaka, Rowlands, Luthcke and BoySabaka and others (2010), we use the preconditioned conjugate gradient method. Figure 2 presents the KBRR residual performance for each solution defined in Table 1. This figure shows the difference in the daily root mean square (rms) of the KBRR residuals using the v08 forward model and the KBRR residuals computed from subsequent versions of forward models (which include iterated mascon solutions for v10–v12). Negative differences represent improvements due to the new forward model of mass change. The v09–v08 daily rms difference in Figure 2 demonstrates significant reduction in the KBRR residuals (model improvement) from forward-modeling hydrology and GIA. The residual reduction is seasonal, with the largest reduction in daily residual rms in Northern Hemisphere spring and an increase in residuals during a portion of Northern Hemisphere summer. The residuals show that the GLDAS hydrology model improves the global surface mass modeling during most of the year, except the Northern Hemisphere summer.
Fig. 2. Daily KBRR residual rms solution differences. The daily (dots) and 15 day Gaussian smoothed (curves) of the KBRR residual rms differences between the v08 solution and the v09 (blue), v10 (purple), v11 (green) and v12 (orange) solutions. Negative values represent a reduction in the residuals and an improvement in the forward modeling, including an iteration of the mascon solution.
The global mascon solutions represent a reduction in the daily KBRR residual rms, as shown in Figure 2. The first iteration mascon solution, v10 forward model, represents a factor of 2.4 reduction in residuals over the v09 forward model (v08 plus GLDAS hydrology and GIA), and further iterations represent 5% (v11 forward model) and 2% (v12 forward model) reduction over the preceding iteration. The global mascon solutions show reduction in seasonal and long-term trend signals in the daily KBRR residuals, due to improved modeling of the global mass change, dominated by hydrology and land ice. The spatial distribution of the residual reduction is presented in Figure 3 as the trend and annual amplitude of the GRACE average time-differenced KBRR residuals (KBRR-dot). The significant trends in the KBRR-dot residuals correspond to the GOA, AIS and GIS land-ice regions, and in particular, the largest residual trends are found in the West Antarctic ice sheet (WAIS) Amundsen and Bellingshausen coasts, and the southeast and northwest coasts of the GIS (Fig. 3a). Figure 3a demonstrates the significant reduction in these KBRR-dot residual trends using the v12 forward model, due to the inclusion of the v11 iterated global mascon solution.
Fig. 3. Global maps of the (a) trend and (b) annual amplitude of the time-differenced KBRR residuals (KBRR-dot) for the forward models used in the v08 (row 1), v09 (row 2) and v12 (row 3) mascon solutions (Table 1). The differences between v09 and v08 show the impact of the additional forward models included in the v09 solution, while the differences between v12 and v09 show the large reduction in residuals that is achieved with the iterated mascon solution.
Iterating the solution has an important impact on the regional terrestrial ice mass signals obtained from the global mascon solution (Fig. 4; Table 2). Iteration results in increases in ice mass loss trend, annual and acceleration signal computed from the global mascon solution. It accounts for a −16.5, −21.8 and −1.1 Gt a−1 correction to the GRACE observed GIS, AIS and GOA mascon solution trends; and 10.7, 34.2 and 4.4 Gt corrections to the GIS, AIS and GOA annual amplitude (Table 2). Iteration of the global mascon solution is shown to reduce KBRR residuals and increase signal recovery. Therefore, we regard the increased signal recovery through iteration as an important part of the global mascon solution, without which the regional land-ice mass signals estimated from the mascon solution would be in error by the amounts noted in Table 2 (v12–v09).
Fig. 4. The effect of iteration on the mass time series for (a) GIS, (b) AIS and (c) GOA. The third and final iteration solution is shown (v12: gray), along with the differences to the original (v09: blue), the first iteration (v10: purple) and the second iteration (v11: green) solutions.
Table 2. Effect of solution iteration on mascon regional mass change statistics over the time period 1 December 2003 to 1 December 2010. The reported statistics describe the differences of the final and third iteration solution (v12) to the original global solution (v09), the first iteration solution (v10) and the second iteration solution (v11)
Signal and Error Analysis
For each of the 41 168 mascons in our global solution, we obtain a time series of mass anomalies with 10 day sampling. To compute the time series of mass anomalies for any region of interest we sum the individual time series at each 10 day sample for all mascons contained within the region. The 10 day mascon solution 1σ error is estimated from the rms of the misfit of the mascon solution time series and a 10 day width Gaussian filter applied to the time series, in a similar way to Reference Luthcke, Arendt, Rowlands, McCarthy and LarsenLuthcke and others (2008). Also, as in Reference Luthcke, Arendt, Rowlands, McCarthy and LarsenLuthcke and others (2008), the mascon time-series trend and annual amplitude are determined from a least-squares simultaneous estimate of the trend, annual, semi-annual and 161 day periodic. The 161 day periodic compensates for the GRACE alias of an S 2 tide error (Reference Ray and LuthckeRay and Luthcke, 2006). For the trend, we replace the formal error estimate (variance of the least-squares parameter) with one that considers a first-order autoregression structure for the mascon time series (Reference Lee and LundLee and Lund, 2004). The formal estimates of the other parameters (e.g. annual amplitude) are then scaled by the ratio of the trend error with autoregressive structure considered, and the simple formal trend error. The total error is then computed as the root sum square (rss) of the following error estimate sources: statistical; leakage in and out; atmosphere/ocean forward modeling; and GIA. These errors are discussed below.
Leakage errors
Leakage of signals in and out of regions of interest occurs due to the fundamental spatial resolution limitations of the GRACE data themselves and the properties of our mascon signal auto-covariance. Errors arising from leakage of signals into (‘leakage-in’) a region of interest (e.g. GIS or AIS) from unmodeled surrounding mass variation sources, such as ocean, tides and hydrology, are estimated, as well as land-ice mass signal leakage out (‘leakage-out’) of our regions of interest. It should be noted that the hydrology signal leakage error has been significantly reduced through the forward modeling of a TWS model, as previously discussed, and as in Reference Luthcke, Arendt, Rowlands, McCarthy and LarsenLuthcke and others (2008) and Reference Sabaka, Rowlands, Luthcke and BoySabaka and others (2010). We use a representative global mass anomaly model in space and time to estimate the leakage-in and leakage-out errors. For this model, we use the vector h k (Eqn (8)) from our final perferred v12 global mascon solution. To estimate leakage-in for a region of interest, we zero the region mascon parameters in h k from Eqn (8), but leave the remaining global mascon parameter values to produce h kR (0). The estimate of leakage-in error is then the region’s mascon time series after applying the resolution operator to our model, h kR (0). We then perform this procedure for each region of interest. The leakage-out errors are computed using a similar procedure, but in this case all mascon parameters outside the region of interest are set to zero. The estimate of leakage-out error is then the difference between the region’s model mascon parameters and those after applying the resolution operator. A summary of the leakage errors for the major land-ice regions of interest is given in Figure 5a and Table 3. While the leakage analysis does consider sources (1 arcdeg mascons) significantly smaller than the fundamental spatial resolution of GRACE (approximately equivalent to 300 km radius Gaussian smoother), the analysis did not consider point sources and localized signal smaller in spatial extent than an individual mascon. An assessment of the impact of sub-mascon sources on our leakage estimates will need to be evaluated in future analysis.
Fig. 5. (a) The v12 mascon solution mass time series for GIS, AIS and GOA with leakage-in (black circles) and leakage-out (red diamonds) estimated errors. (b) The v12 mascon solution time series for GIS, AIS and GOA with atmosphere and ocean model signal (ECMWF/OMCT: black circles) and estimated model errors (ECMWF/OMCT – MERRA/IB: red diamonds). (c) The GIS v12 mascon solution time series and the difference to the mass time series derived from the averaging kernel filter technique (red diamonds).
Table 3. Effect of leakage-in, leakage-out and the rss of the leakage-in and -out, on the annual amplitude, acceleration and trend of mass changes for each main region and subregion over the time period 1 December 2003 to 1 December 2010
Atmosphere/ocean forward-modeling errors
Following Reference Han, Jekeli and ShumHan and others (2004), we estimate the errors in forward-modeling atmosphere and ocean mass as the difference between two sets of models: (1) ECMWF and OMCT vs (2) atmospheric mass variations derived from NASA/GSFC Modern Era Retrospective-analysis for Research and Applications (MERRA) model and a simple ocean inverted barometer (IB). A summary of the impact of these models on our land-ice solutions is provided in Figure 5b and Table 4. The two models are not completely independent, so we do not scale the differences by 
(Reference Han, Jekeli and ShumHan and others, 2004). Ocean mass variation and ocean tides forward-modeling error contributions to our land-ice mass solutions are included through the leakage-in analysis discussed above. The mass anomaly model used in the leakage analysis includes the ocean mass and ocean tides modeling error observed by GRACE, as the mascon solution reflects the mass anomalies from the forward models.
Table 4. Annual amplitude, acceleration and trend of mass changes due to the atmospheric and ocean signal and errors for each main region and subregion over the time period 1 December 2003 to 1 December 2010
GIA errors
The land-ice mass trends estimated from our global mascon solutions contain errors due to the uncertainties in the GIA model applied. We estimate the GIA error as the difference between the minimum and maximum GIA correction over a range of plausible GIA models divided by 2 for each region of interest. For the GIS and GOA the models considered include the nominal ICE-5G Paulson incompressible two-layer model discussed previously; an ICE-5G compressible model (Reference Geruo, Wahr and ZhongGeruo and others, 2013); and models based on ‘ANU’ ice deglaciation history (Reference Fleming and LambeckFleming and Lambeck, 2004) applying a suite of Earth models (27) over a range of lithospheric thicknesses, and upper mantle and lower mantle viscosities (Reference Barletta, Sørensen and ForsbergBarletta and others, 2012, computed and provided by V.R. Barletta). For the AIS the models considered include the nominal IJ05_R2 used in this study and discussed previously; models based on IJ05_R2 with varying Earth models; and models based on W12a with the ‘best’-fit Earth model and upper- and lower-bound Earth models that give a good fit to the relative sea-level data (Reference Whitehouse, Bentley, Milne, King and ThomasWhitehouse and others, 2012, computed and provided by P. Whitehouse). This is a similar approach to estimate the GIA uncertainty to that used in the IMBIE study (Reference ShepherdShepherd and others, 2012). The post-LIA GIA signal uncertainty is ±30%, estimated from the model comparisons to GPS vertical motions, raised shoreline records of RSL and tide-gauge rates of RSL (Reference Larsen, Motyka, Freymueller, Echelmeyer and IvinsLarsen and others, 2005).
Solution technique error analysis
As another step in our solution evaluation, we compare our mascon solution with that of a solution made using the more common approach of applying a filter to the monthly spherical harmonic fields. In this test we apply a calibrated regional GIS averaging kernel to our monthly spherical harmonic fields estimated from the same Level-1B processing as our mascon solution (Reference Luthcke, Rowlands, Lemoine, Klosko, Chinn and McCarthyLuthcke and others, 2006b). The averaging kernel is constructed, calibrated and applied in a similar fashion to Reference Velicogna and WahrVelicogna and Wahr (2005). The important aspect of this comparison experiment is that both the mascons and the monthly spherical harmonic fields have been estimated from the same Level-1B processing and KBRR reduction, so the results isolate the differences arising from the technique used to compute the total GIS mass time series. Figure 5c shows (in red) the difference between the time series computed from the GIS mascons and that computed from the averaging kernel; the mascon solution time series (black) is also provided in this figure for comparison. The differences in the annual amplitude, acceleration and trend between the two time series are, respectively, 25 Gt, −1 Gt a−2 and −11 Gt a−1. These potential uncertainties are due solely to the differences in the techniques used for the mascon time-series recovery, and are on the order of the estimated uncertainties for our mascon solution, as summarized in Table 5 and discussed in the next subsection.
Table 5. Summary of the v12 global mascon solution for the land-ice regions studied over the time period 1 December 2003 to 1 December 2010. The reported parameters are computed after applying the GIA and LIA model correction to the mascon time series. The total uncertainties include contributions from statistical errors, mismodeling of the atmosphere, leakage errors, GIA model errors and, for Gulf of Alaska, the LIA model errors
Computation of seasonal balances
We use a linear regression analysis to extract the multi-year average trend, acceleration and annual amplitude from our 10 day mascon time series. One of the strengths of GRACE solutions, and of the mascon solution presented here in particular, is the high temporal resolution at which the evolution of land-ice mass balance can be observed and quantified. Therefore, we seek to compute annual mass balances from our GRACE mascon time series, in order to better quantify the more complex spatial and temporal variability of land ice. Reference Luthcke, Arendt, Rowlands, McCarthy and LarsenLuthcke and others (2008) simply used the extrema of the full mascon time series as the start and end of the balance years from which to compute the annual mass balances. Because the mass time series for the sum of GOA glacier mascons has a very clean annual signal with little additional signal or noise, this is an acceptable approach. However, for the more complex mascon time series from the GIS, WAIS, East Antarctic ice sheet (EAIS) and Antarctic ice sheet peninsula (AISP), an adaptive filter approach is necessary. Therefore, we employ an ensemble empirical mode decomposition (EEMD) adaptive filter to first compute the timing of the balance seasons (accumulation and loss), and then to compute the annual mass balances.
The EEMD is a noise-assisted data analysis improvement of the EMD method, which provides an adaptive decomposition of the time series into components known as intrinsic mode functions (IMFs) (Reference Wu and HuangWu and Huang, 2009). The IMFs are monochromatic at a given point in time, but vary in time with both amplitude and frequency. The IMF adaptive decomposition separates out the frequency content of the time series, which represent different physical modes. An ensemble (many runs) of the EMD is performed, applying a new generation of noise to the original signal on each run within the ensemble. The added finite white noise provides a means to sort out signals with similar scale into their proper IMFs and thus mitigate mode mixing. The mean of the ensemble is treated as the final true result, so the effect of noise is minimized. Unlike typical signal analysis techniques, such as the linear regression analysis applied above, the EEMD is an adaptive filter and does not rely on linear and stationarity assumptions.
We apply the EEMD to adaptively find the extrema of the non-stationary annual signal within the mascon time series of interest. The minima and maxima of the annual IMF represent the timing of the mass-balance seasons (beginning of the mass accumulation and mass loss seasons). This is the most important aspect of our EEMD application. The adaptive filter is used to isolate the annual signal, and then to identify the timing of the balance seasons and length of each balance year as the annual signal is non-stationary and evolving. The original signal is then reconstructed by summing the IMFs, and eliminating noise by discarding the first few IMFs (Fig. 8). We then compute the mass value of the reconstructed signal at the beginning of the mass-balance seasons (annual IMF minima). We use 100 runs within each EEMD ensemble, where the added white noise is recomputed on each run within the ensemble, with a noise ratio equal to twice the ratio of the standard deviation of the mascon solution noise (difference between 10 day solution time series values and a 30 day Gaussian filter), and the standard deviation of the signal (with trend removed). The EEMD is then run 50 times, where the overall signal trend is removed before each run and is added back after the reconstruction of the signal from the selected IMFs. The mean (computed over the 50 EEMD runs) timing value of the annual IMF minima and corresponding reconstructed signal mass values are used to compute the annual balance (Fig. 8). The annual balance is computed by differencing the reconstructed mass values at the times of two consecutive mean annual IMF minima. The mass difference, or annual balance, is labeled with the year that is mostly encompassed between the two minima. Therefore, for the GIS and GOA land ice, the balance year starts in the late summer to mid-fall of the previous year, while for the AIS the balance year starts in the Northern Hemisphere spring. The standard deviation of these values over the 50 runs is used as an estimate of the noise contribution to the uncertainty. To estimate a potential systematic error in our balance year timing and values, the EEMD analysis results are differenced from an iterative piecewise trend and annual linear regression analysis of the mascon time series, where the annual pieces are constrained to form a continuous signal. The error estimate of balance year extrema timing and mass value is then computed as the rss of the noise and systematic error estimates. The error estimate of the annual balance is then the rss of the error estimates for the two balance year extrema.
Fig. 8. The EEMD v12 mascon solution time series and error analysis for (a) GIS, (b) GOA, (c) EAIS, (d) WAIS and (e) AIS peninsula. Each plot shows the v12 mascon solution time series (black circles), the mean of 50 EEMD runs, where each case sums the IMFs from the one preceding the ‘annual’ IMF to the final IMF (blue curve), the mean values from the 50 EEMD runs at the times of the ‘annual’ IMF extrema (red triangles) and the extrema error ellipses (green ellipses).
Land-Ice Mass Evolution Results
In this section we provide a summary of our iterated global mascon solution, focusing on the results for the five most significant land-ice regions: GIS, WAIS, EAIS, AISP and GOA. Table 5 provides an overall summary of the regional mass time series computed from the mascon solution. While the time series cover a longer time period, the model parameter fit is applied to an integer number of years, 1 December 2003 to 1 December 2010. We use the terminology defined by Reference CogleyCogley and others (2011) to explain glacier and ice-sheet variations in terms of variations in their climatic balances, B clim, and calving fluxes, D, ignoring the contributions from basal balances.
GIS results
The time series for the sum of the GIS mascons (GIS total), and the mascons above and below 2000 m are shown in Figure 6, and a summary of these time series is provided in Table 5. The spatial distribution of the GIS trend and acceleration is shown in Figure 7. The GIS exhibits the largest rates of mass loss relative to both Antarctica and Alaska, with a trend of −230 ± 12 Gt a−1, an acceleration of mass loss of −10 ± 6 Gt a−2 and an average annual amplitude of 132 ± 29 Gt. Both the magnitude of the GIA correction and our estimated uncertainty determined from the model differences are quite small in comparison with the total corrected trend (Table 5). Our estimates agree to within stated error bars with previous estimates spanning at least the period 2003 to October 2009: −219 ± 38 Gt a−1 for April 2002–November 2009 (Reference Chen, Wilson and TapleyChen and others, 2011); −238 ± 29 Gt a−1 for October 2003–October 2009 (Reference SasgenSasgen and others, 2012); −222 ± 9 Gt a−1 for January 2003– December 2010 (Reference Jacob, Wahr, Pfeffer and SwensonJacob and others, 2012). Our estimated acceleration of −10 ± 6 Gt a−2 is in agreement with −15 ± 7 Gt a−2 for August 2001–August 2010 (Reference SasgenSasgen and others, 2012) and −17 ± 8 Gt a−2 for April 2002–June 2010 (Reference Rignot, Velicogna, Van den Broeke, Monaghan and LenaertsRignot and others, 2011). The GIS mass evolution is dominated by the signal from the coastal margins of the ice sheet, as defined in our solution as mascons <2000 m. Large coastal mass loss rates result from high rates of D associated with the presence of fast-flowing marine-terminating glaciers along the northwest, southwest and southeast regions of the GIS (Reference Moon, Joughin, Smith and HowatMoon and others, 2012), and due to the higher surface melt rates that occur at low elevations on the ice sheet. Above 2000 m there is little statistically significant signal.
Fig. 6. The v12 10 day mascon solution mass time series (circles) and 10 day Gaussian smoothed solution (curves) for GIS total (blue), GIS <2000 m elevation (purple) and GIS >2000 m elevation (green). The red lines show the best-fit trend plus acceleration.
Fig. 7. The GIS v12 mascon solution mass changes computed over the time period 1 December 2003 to 1 December 2010. (a) Trend corrected with ICE5G_Paulson GIA model and (b) acceleration.
The GIS annual balance results from the application of the EEMD filter are shown in Figures 8 and 9. The GIS annual balance anomalies (difference from the mean) show large interannual variability (Fig. 9). The GIS mean annual balance over the balance years 2004–10 is −248 ± 12 Gt a−1. It is important to note that our mean annual balance is slightly different from the trend reported above, primarily due to the differing time periods, where integer balance years are needed for the annual balance computation. The balance years 2005–09 have anomaly magnitudes <20% of the mean annual balance and alternate between positive and negative values. In contrast, the GIS balance years 2004 and 2010 have large positive and negative anomalies, ∼50% of the mean annual balance, giving rise to the statistically significant acceleration determined from the linear regression analysis (Table 5). Therefore our reported GIS mass loss acceleration is, in part, an artifact of the GRACE sampling interval, rather than an indication of a persistent multi-year signal, and longer time series will be needed before making inferences about future evolution of the ice sheet.
Fig. 9. The annual mass balance mean (2004–10) and annual anomalies about the mean for GIS (blue), GOA (green), EAIS (purple), WAIS (orange) and AIS peninsula (cyan) from the EEMD analysis of the v12 mascon solution. The values and corresponding error bars are determined from the mean extrema and extrema error ellipses shown in Figure 8.
In the northwest GIS (drainage systems 81 and 82) we observe an increase in mass losses after 2005 (Figs 10 and 11), attributed to a reduction in snow accumulation (Reference SasgenSasgen and others, 2012) that combines with increased D after 2007 (Reference Moon, Joughin, Smith and HowatMoon and others, 2012) to further enhance mass losses. The timing and extent of the northward progression of mass loss is confirmed by GPS observations, showing enhanced rates of uplift late in 2005 at Thule Station in northwest Greenland (Reference Khan, Wahr, Bevis, Velicogna and KendrickKhan and others, 2010). In the southwest GIS (drainage systems 61, 62 and 72) we observe persistent mass losses, attributed to high D that commenced in the late 1990s, that are dominated by the large calving losses from Jakobshavn Isbræ, the largest marine-terminating outlet on the GIS (Reference JoughinJoughin and others, 2012).
Fig. 10. The GIS annual mass balances and mean annual mass balance determined from the EEMD analysis of the v12 mascon solution. The timing of the seasons is determined by the mean EEMD extrema in Figure 8a, and the mass values are determined by the result of the EEMD method (summing the IMFs from the one preceding the ‘annual’ IMF to the final IMF) for each individual mascon time series at the time of the total GIS extrema.
Fig. 11. The v12 10 day mascon solution time series (circles) and 10 day Gaussian smoothed solution (curves) for selected drainage systems within the GIS and AIS, and mascons within the GOA land-ice region. The colors of the plotted curves correspond to the drainage system and mascon colors in Figure 1.
Large interannual mass variability is observed in drainage system 32, the region of the Geikie Peninsula in central East Greenland that is largely composed of ice not receiving flow contributions from the GIS (Fig. 10) (Reference Jiskoot, Juhlin, St Pierre and CitterioJiskoot and others, 2012). The mass time series for this region shows an increase in mass loss beginning in 2005, followed by mass gains from 2007 to 09 (Fig. 11). The glaciers in this area are smaller and steeper than ice-sheet outlet glaciers and exist in a region with high precipitation rates (Reference EttemaEttema and others, 2009), so they likely respond rapidly to climate variations in a fashion similar to other maritime glacier systems. Approximately 90% of the glaciers in this region drain through tidewater glacier outlets (Reference Jiskoot, Juhlin, St Pierre and CitterioJiskoot and others, 2012), and those draining to the south into the Blosseville Kyst may be particularly sensitive to changes in climate, due to their southerly aspect and concentration of ice at lower elevations.
Several studies have reported synchronous retreat of marine-terminating glaciers in southeast Greenland during 2000–06 (Reference Howat, Joughin, Fahnestock, Smith and ScambosHowat and others, 2008), followed by advance and thickening of Helheim and Kangerdlugssuaq Glaciers (Reference Howat, Joughin and ScambosHowat and others, 2007; Reference JoughinJoughin and others, 2008). The reduction in mass loss associated with these advances has been assumed to be synchronous across the southeast coast (Reference MurrayMurray and others, 2010), although evidence for this was based, in part, on terminus advance rates that did not measure glacier mass or volume changes directly (Reference Howat, Joughin, Fahnestock, Smith and ScambosHowat and others, 2008). Our results show a deceleration in mass loss during the entire 2003–10 period (Fig. 7b), but analysis of the individual mass-balance seasons shows that much of the deceleration occurred due to mass gains in the central east GIS in 2008, and a synchronous reduction in mass losses in 2009 (Fig. 10). These findings are consistent with those of Reference Chen, Wilson and TapleyChen and others (2011), who used GRACE observations to infer a significant reduction in the mass loss of southeast GIS glaciers during September 2007–November 2009; however, our results pinpoint the timing of this reduction in mass loss to the 2009 balance year. The 2009 slowdown of mass loss in the southeast GIS likely occurred due to a complex array of factors, including a deceleration and advance of several southeast GIS glaciers (Reference MurrayMurray and others, 2010), as well as 2008 and 2009 positive accumulation anomalies (Reference SasgenSasgen and others, 2012). We note that observed deceleration and advance of southeast GIS glaciers occurred 1–2 years prior to our observed 2009 reduction in mass losses. This reflects the complex and asynchronous nature of ice-sheet dynamics in the southeast GIS (Reference Moon, Joughin, Smith and HowatMoon and others, 2012).
Our observations show evidence of the response of the GIS to variations in summer air temperatures and the length of the summer melt season. We observe high rates of mass loss across the western and northwestern GIS during 2008, associated with extreme summer snowmelt that summer, detected by passive microwave imaging sensors (Reference Tedesco, Fettweis, Van den Broeke, Van de Wal and SmeetsTedesco and others, 2008). Coastal mass balances during 2009 were less negative than all other years in the GRACE record, and this was followed in 2010 with the most negative mass balances below 2000 m (Fig. 10). These patterns agree well with observed and simulated melt extent records showing low/high melt extents in the 2009/10 summers (Reference Mernild, Mote and ListonMernild and others, 2011).
AIS results
The mascon solution for the AIS regions is summarized in Table 5 and Figures 12 and 13. Unlike the GIS, both the magnitude and uncertainties of the GIA corrections for the AIS are significant with respect to the total trend. The overall AIS trend for the December 2003–December 2010 time period studied is −81 ± 26 Gt a−1. Our trend estimate agrees, within stated uncertainties, with −69 ± 18 Gt a−1 obtained by Reference King, Bingham, Moore, Whitehouse, Bentley and MilneKing and others (2012) for August 2002–December 2010, using a comparable recent regional GIA model. The AIS trend is dominated by significant mass loss from the WAIS, with a trend of −106 ± 16 Gt a−1. A trend of −38 ± 14 Gt a−1 is found for the AIS peninsula, while the EAIS trend is 63 ± 28 Gt a−1. The largest trends are found in the WAIS along the Amundsen and Bellingshausen Sea coasts (drainage systems 19–22), concentrated at the Pine Island embayment (convergence of drainage systems 19–22), and at the northern tip of the peninsula (drainage system 24). These results are in general agreement with −132 ± 60 Gt a−1 for the WAIS and −60 ± 46 Gt a−1 for the AIS peninsula, determined from radar interferometry and regional climate models for the year 2006 (Reference RignotRignot and others, 2008). In addition, our AIS peninsula trend is in good agreement with −42 ± 9 Gt a−1 for January 2003–March 2009, determined from a GRACE spherical cap mascon solution (Reference Ivins, Watkins, Yuan, Dietrich, Casassa and RülkeIvins and others, 2011).
Fig. 12. The v12 10 day mascon solution time series (circles) and 10 day Gaussian smoothed solution (curves) for the AIS total (blue), AIS <2000 m elevation (purple), AIS >2000 m elevation (green), EAIS (orange), WAIS (cyan) and AIS peninsula (gold). The red curves show the best-fit trend and acceleration.
Fig. 13. The AIS v12 mascon solution mass changes computed over the time period 1 December 2003 to 1 December 2010. (a) Trend corrected with the IJ05_R2 GIA model and (b) acceleration. The EAIS, WAIS and AIS peninsula regions are demarcated by the cyan borders.
The most striking signal, with the greatest potential consequence for sea-level rise, is the large accelerated mass loss exhibited by the WAIS at −46 ± 6 Gt a−2. While we do observe interannual variability of the WAIS annual balances, there is a consistent negative trend in these data, representing a persistent multi-year acceleration of mass loss (Figs 8d and 9). The largest accelerated loss is concentrated at Pine Island embayment and along the Amundsen and Bellingshausen Sea coasts with significant mass loss extending to the Zumberge coast (drainage system 1) during the 2009 balance year (Fig. 14). These changes are consistent with observations of grounding-line retreat (Reference RignotRignot and others, 2008), acceleration (Reference Joughin, Smith, Howat, Scambos and MoonJoughin and others, 2010) and drawdown (Reference Wingham, Wallis and ShepherdWingham and others, 2009) of the Pine Island Glacier driven by basal melting of its large ice shelf (Reference Pritchard, Ligtenberg, Fricker, Vaughan, Van den Broeke and PadmanPritchard and others, 2012).
Fig. 14. The AIS annual mass balances and mean annual mass balance determined from the EEMD analysis of the v12 mascon solution. The timing of the seasons is determined by the mean EEMD extrema in Figure 8c, d and e, and the mass values are determined by the result of the EEMD method (summing the IMFs from the one preceding the ‘annual’ IMF to the final IMF) for each individual mascon time series at the time of the extrema for the region in which the mascon is located (east, west or peninsula). The EAIS, WAIS and AIS peninsula regions are demarcated by the cyan borders. The mean annual mass balance is shown for the v12 solution corrected with the IJ05_R2 GIA model.
The EAIS trend (63 ± 28 Gt a−1) and accelerated mass gain (22 ± 22 Gt a−2) (Table 5; Fig. 8c) are not persistent multi-year signals but are due to an exceptional positive surface mass-balance anomaly during the 2009 accumulation season, concentrated primarily along the coast of Dronning Maud Land (drainage systems 5–7), but also along the coast of Wilkes Land (drainage systems 11 and 12) (Figs 8c and 14). These observations have similar patterns to melting days anomalies derived from passive microwave observations which, for the entire AIS, reached a new historical minimum for the 1980–2009 period (Reference Tedesco and MonaghanTedesco and Monaghan, 2009). In addition, our GRACE observations agree quite well with the timing and spatial location of a significant 2009 accumulation anomaly observed in the surface mass-balance record derived from the regional atmospheric climate model RACMO2 (Reference Lenaerts, Van den Broeke, Van de Berg, Van Meijgaard and Kuipers MunnekeLenaerts and others, 2012; Reference ShepherdShepherd and others, 2012).
The northern tip of the AIS peninsula exhibits consistent negative annual balance, with the exception of the 2008 balance year (Fig. 14), which gives rise to the shallow positive acceleration shown in Figure 13. The southern area of the peninsula shows small positive and negative annual balances, with the exception of large negative annual balances for 2007 and 2009, giving rise to the sharp losses shown in the AIS peninsula time series for these years (Fig. 8e). These large negative annual balance anomalies give rise to the acceleration signal in the southern peninsula shown in Figure 13, and the overall peninsula acceleration of −9 ± 5 Gt a−2. Large mass losses on the AIS peninsula result from tributary glaciers that feed into ice shelves that have recently disintegrated (Reference Shuman, Berthier and ScambosShuman and others, 2011), and increasingly negative mass balances resulting from strong climatic warming in the region. Our observations agree with satellite observations of surface lowering (Reference Fricker and PadmanFricker and Padman, 2012) and model estimates of increasing mass losses from peripheral glaciers and ice caps in the region (Reference Hock, De Woul and RadićHock and others, 2009).
GOA results
The mascon solution for the GOA glacier region is summarized in Table 5 and Figures 8b and 15. The overall trend of the GOA glaciers for the time period studied is −69 ± 11 Gt a−1, and −77 ± 11 Gt a−1 if we take the mean of the 2004–10 balance years (Fig. 9). Most of the GIA correction in this region results from uplift corrections due to LIA mass loss in the Glacier Bay region of Alaska. Our estimates very nearly agree to within stated error bars with −46 ± 7 Gt a−1 reported by Reference Jacob, Wahr, Pfeffer and SwensonJacob and others (2012). Differences in the best estimates between these two studies likely results from different treatment of TWS, and slightly different solution domains.
Fig. 15. The GOA mascon trends computed from the v12 mascon solution corrected with the ICE5G Paulson GIA model and the Larsen LIA model over the time period 1 December 2003 to 1 December 2010.
GOA glaciers exhibited large seasonal and interannual variability during the 2003–10 measurement period (Figs 8b, 11 and 16). We observed large mass losses during summer 2004, in response to record high temperatures across Alaska that year (Reference Truffer, Harrison and MarchTruffer and others, 2005). Even larger mass losses occurred in 2009, despite the absence of a strong temperature forcing, and we attribute these losses to the 31 March 2009 eruption of Mount Redoubt that spread volcanic ash across much of the GOA region (Reference SchaeferSchaefer and others, 2012) and likely enhanced melt rates through a reduction in surface albedo. Large 2007/08 winter accumulation was followed by a cool 2008 summer season, resulting in near-balance conditions during 2008. This offset the strong 2004 and 2009 melt seasons and resulted in no significant acceleration in GOA mass balance during the GRACE period.
Fig. 16. The GOA land-ice region annual mass balances and mean annual mass balance determined from the EEMD analysis of the v12 mascon solution. The timing of the seasons is determined by the mean EEMD extrema in Figure 8b, and the mass values are determined by the result of the EEMD method (summing the IMFs from the one preceding the ‘annual’ IMF to the final IMF) for each individual mascon time series at the time of the total GOA extrema.
The spatial distribution of mass balances largely follows the distribution of glacierization in the GOA region, with largest mass losses occurring in the St Elias and Glacier Bay regions. The highly variable distribution of glaciers in the GOA region, with large signals occurring in the center of the domain and much smaller signals at the northwestern and southeastern edges, likely results in smearing and leakage of signal between individual mascons. Therefore, although the total GOA signal is not affected by these problems, additional processing steps are required to further constrain mass changes over drainage systems smaller than our current solution domain.
Concluding Remarks
We have developed a new global mascon solution specifically tuned to observe high-latitude land-ice mass evolution estimated directly from the GRACE inter-satellite range-rate data. The study is unique in many ways including: (1) a detailed error analysis of a constrained global mascon solution with estimates of signal leakage in and out of the regions of interest; (2) a controlled comparison to solutions derived from monthly spherical harmonic solutions using the same GRACE data processing and forward models, providing estimates of uncertainty due to solution technique; (3) fully iterating a global mascon solution to convergence and quantifying the signal recovery improvement and GRACE inter-satellite range-rate residual reduction; and (4) applying the EEMD to compute land-ice annual balances. We estimate the land-ice mass trends for the time period 1 December 2003 to 1 December 2010 to be −230 ± 12 Gt a−1 for the GIS, −81 ± 26 Gt a−1 for the AIS and −69 ± 11 Gt a−1 for the GOA glaciers. For the subregions of the AIS we estimate 63 ± 28 Gt a−1 for the EAIS, −106 ± 16 Gt a−1 for the WAIS and −38 ± 14 Gt a−1 for the AIS peninsula. The total over the land-ice regions studied is −380 ± 31 Gt a−1, equivalent to −1.05 ± 0.09 mm a−1 sea-level rise. Over the same time period we estimate the mass acceleration to be −10 ± 6 Gt a−2 for the GIS, −33 ± 26 Gt a−2 for the AIS and 2 ± 7 Gt a−2 for the GOA glaciers, giving a total of −41 ± 27 Gt a−2, equivalent to −0.11 ± 0.08 mm a−2 sea-level rise. For the subregions of the AIS we estimate the accelerations to be 22 ± 22 Gt a−2 for the EAIS, −46 ± 6 Gt a−2 for the WAIS and −9 ± 5 Gt a−2 for the AIS peninsula. The trends and accelerations determined over our study period have been shown to be highly dependent on significant seasonal and annual balance anomalies, making it difficult to predict future land-ice mass balance and its contribution to sea level. One possible exception is the WAIS acceleration signal, which shows a persistent negative trend in annual balance. The iterated global mascon solution presented here provides the spatial and temporal monitoring of land-ice evolution to aid in improved understanding and modeling.
Acknowledgements
Support for this work was provided by NASA under the Interdisciplinary Research in Earth Science (NNH09ZDA001N-IDS), the GRACE science team (NNH06ZDA001N-GRACE and NNH10ZDA001N-GRACE) and the Cryospheric Sciences (NNH07ZDA001N-CRYO) programs. We thank NASA GSFC NCCS computer services for computational resources used in this research. We gratefully acknowledge the quality of the GRACE Level-1B products produced by our colleagues at the Jet Propulsion Laboratory, California Institute of Technology, Pasadena. We thank J.P. Boy and R.D. Ray for contributions to the forward models applied in our GRACE data reduction and analysis. We acknowledge the software development and maintenance support of T.A. Pennington and D.E. Pavlis. We especially acknowledge the numerous contributions of D.D. Rowlands in developing the foundation of algorithms and software necessary to carry out this research. In addition, we thank D.D. Rowlands for the many technical discussions and support. We also thank the reviewers for thoughtful comments and suggestions.
