Interactive comment on “ Constraining GRACE-derived cryosphere-attributed signal to irregularly shaped ice-covered areas ”

Abstract. We use a Monte Carlo approach to invert a spherical harmonic representation of cryosphere-attributed mass change in order to infer the most likely underlying mass changes within irregularly shaped ice-covered areas at nominal 26 km resolution. By inverting a spherical harmonic representation through the incorporation of additional fractional ice coverage information, this approach seeks to eliminate signal leakage between non-ice-covered and ice-covered areas. The spherical harmonic representation suggests a Greenland mass loss of 251 ± 25 Gt a−1 over the December 2003 to December 2010 period. The inversion suggests 218 ± 20 Gt a−1 was due to the ice sheet proper, and 34 ± 5 Gt a−1 (or ~14%) was due to Greenland peripheral glaciers and ice caps (GrPGICs). This mass loss from GrPGICs exceeds that inferred from all ice masses on both Ellesmere and Devon islands combined. This partition therefore highlights that GRACE-derived "Greenland" mass loss cannot be taken as synonymous with "Greenland ice sheet" mass loss when making comparisons with estimates of ice sheet mass balance derived from techniques that sample only the ice sheet proper.


Introduction
The Gravity Recovery And Climate Experiment (GRACE) satellite constellation measures anomalies in Earth's gravity field.Changes in ice sheet mass balance can be quantified through repeated observations of these gravitational anomalies.Unlike an altimetry approach (e.g., Zwally et al., 2011), in which an ice sheet volume change is converted into a mass change, a gravimetry approach does not require the assumption of an effective density of change or an explicit treatment of the refreezing versus runoff fraction of surface melt.Unlike an input-output approach (e.g., Rignot et al., 2008), in which interferometric synthetic aperture radar ice discharge estimates are combined with modeled surface mass balance estimates, a gravimetry approach does not require precise knowledge of ice geometry or vertical velocity profiles near grounding lines of outlet glaciers (Alley et al., 2007).A gravimetry approach to ice sheet mass balance, however, is sensitive to the models used to isolate the cryospheric mass change signal from other signals, such as mass changes due to land hydrology, ocean and atmospheric mass exchange and solid earth processes.Several GRACE studies have documented the increasingly negative mass balance of the Greenland ice sheet, from an initial estimate of −76 ± 26 Gt a −1 during the May 2002 to July 2004 period (Velicogna and Wahr, 2005) to a more recent estimate of −263 ± 30 Gt a −1 during the January 2005 to December 2010 period (Shepherd et al., 2012).
GRACE-derived spherical harmonic solutions of rate of mass change have coarse spatial resolution, and consequently signal leakage from within defined areas (Velicogna and Wahr, 2006).While GRACE-derived mass change estimates have been assessed at basin-scale resolution over the Greenland ice sheet (Luthcke et al., 2006;Velicogna and Wahr, 2006;Sasgen et al., 2012;Barletta et al., 2013), signal leakage has prevented the development of a GRACE-derived mass change field that is completely constrained to within the irregularly shaped ice-covered areas of Greenland.Deconstructing spherical harmonic solutions, in a non-iterative  (Luthcke et al., 2013).(B) The equivalent spherical harmonic representation.fashion using separate land and ocean filters, has been demonstrated to reduce land-ocean signal leakage (Guo et al., 2010).The use of local mass concentrations ("mascons") offers an alternative approach to improve spatial resolution and reduce signal leakage (Luthcke et al., 2006;Jacob et al., 2012).No interpretation of satellite gravimetry, however, has yet been capable of partitioning mass change due to Greenland peripheral glaciers and ice caps (GrPGICs) from that due to the Greenland ice sheet proper.It is desirable to isolate peripheral glacier mass change from ice sheet mass change, as it is believed that smaller peripheral glaciers and ice caps should have shorter response times than the larger ice sheet to contemporary climate change (Nye, 1960;Jóhannesson et al., 1989).

Method
The fundamental problem addressed in this work is how to extract robust mass variations over spatially limited regions from GRACE data (e.g., Simons et al., 2006).We use a Monte Carlo inversion to infer the most likely 26 km resolution rate of mass change field that, when smoothed with a Gaussian filter, closely reproduces a GRACE-observed spherical harmonic representation of cryosphere-attributed rate of mass change.This inversion of GRACE-observed rate of mass change incorporates new information in the form of observed fractional ice coverage.By constraining the inverted cryosphere-attributed rate of mass change to icecovered areas, this approach seeks to eliminate signal leakage between non-ice-covered and ice-covered areas.We do not, however, purport to resolve spatial heterogeneity in rate of mass change at the scale of individual glaciers and ice caps.Distinguishing sharp spatial heterogeneities in mass changes between adjacent ice-containing nodes would require the introduction of further new information to the inversion.
The variable notation used in the methodology description is summarized in the Appendix Table A1.

Data
Our inversion requires two distinct input data: (i) a spherical harmonic representation of cryosphere-attributed rate of mass change, and (ii) glacier area/extent information used to derive fractional ice coverage at a given node.The spherical harmonic representation we invert characterizes the cryosphere-attributed rate of mass change over a given period ultimately derived from GRACE ( ṀG ), not individual spherical harmonic coefficients.We therefore invert a GRACEderived ṀG field in Cartesian space, rather than in spherical harmonic space.Our time period of interest, over which the ṀG field is determined, is the 1 December 2003 to 1 December 2010 GRACE inter-comparison period used by the Ice Sheet Mass Balance Inter-comparison Exercise (IMBIE; Shepherd et al., 2012).The ṀG field is derived by converting time series of cryosphere-attributed mascons into equivalent spherical harmonic solutions (Fig. 1).
These mascons are derived from a NASA Goddard Space Flight Center (GSFC) data product in which GRACE level 1 (L1) K band inter-satellite range rate (KBRR) data are reduced via forward modeling into a series of iterated monthly mascons (Luthcke et al., 2013).During this L1 data reduction, forward modeling and seasonal detrending removes mass changes associated with terrestrial hydrology, the ocean and atmosphere, as well as glacial isostatic adjustment and ocean tides.Once these constraints have been applied in mascon space to limit leakage and isolate cryospheric mass change signal significantly, the linear mass change trend and 1σ trend error are calculated for each mascon time series.Mascon-derived rates of mass change and associated error are then converted into equivalent spherical harmonics of degree and order 60 (Fig. 2).
Mascons and spherical harmonics have been previously demonstrated to be interchangeable, essentially by representing mascons as a set of differential potential coefficients ("delta coefficients") added to the mean GRACE level 2 (L2) field (Chao et al., 1987;Luthcke et al., 2013).By virtue of this conversion, the spatial distribution of mass change in the resultant spherical harmonic representation appears slightly different than that of the original mascons.Nonzero ocean values acknowledge cryospheric signal leakage from adjacent terrestrial mascons into ocean mascons.This leakage is most pronounced in northern Baffin Bay, where terrestrial mascons surround ocean mascons on three sides.Luthcke et al. (2013) constrain oceanic mass changes via forward modeling in the iterative reduction of L1 data, but acknowledge that the oceanic forward model cannot reproduce heretofore unidentified oceanic mass change signals.This signal leakage, however, is formally incorporated in the error estimate associated with the trend in cryospheric mass change used in this inversion.While this degree-60 spherical harmonic representation of cryosphere-attributed mascons is global in extent, we only utilize the portion that covers an ∼ 11× 10 6 km 2 region of interest centered over Greenland and the Canadian High Arctic.GRACE L2 products include time-variable (monthly) gravity fields derived from spherical harmonics.These fields are often normalized over a given period and described as mass anomalies, typically in units of water equivalent mass.Applying a temporal trend through a suite of monthly mass anomaly fields can yield a spatial distribution of rate of mass change, typically in units of water equivalent mass per time.The method we present can invert such an L2-derived trend product.Our preference for a mascon-derived spherical harmonic stems from the notion that mascons provide a better framework for isolating gravimetry signals associated with a specific surface process than spherical harmonics (cf.Velicogna andWahr, 2006, andJacob et al., 2012).The statistical inversion we present, however, requires the relatively smooth spatial gradients in rate of mass change char-acteristic of spherical harmonics, rather than relatively sharp contrasts in mass change at node boundaries characteristic of mascons.Thus, we adopt a mascon-derived spherical harmonic representation as our input data while acknowledging that the absolute accuracy of mascons makes them a more desirable inversion target.We note that, unlike the approach we present below, which applies an isotropic Gaussian inversion filter to all nodes, a mascon-based inversion requires an anisotropic geometric filter, the precise shape and symmetry of which varies by node.Thus, a mascon inversion represents a non-trivial extension of a spherical harmonic inversion.
The second distinct type of input data we utilize is knowledge of glacier area/extent within our region of interest centered over Greenland and the Canadian High Arctic.We calculate the fractional ice coverage of all glaciers external to Greenland from the Randolph Glacier Inventory (RGI) version 2.0 (Arendt et al., 2012).This ice coverage includes the north and south Canadian Arctic, Iceland and Svalbard (i.e.,RGI regions 3,4,6 and 7), where the RGI polygon accuracy is better than two Landsat pixels (or 30 m).We interpolate fractional Greenland ice coverage from a Greenland glacier inventory that has a polygon accuracy of 10 m (Citterio and Ahlstrøm, 2013).We group the ice coverage within our domain into four types: (i) Greenland ice sheet (GrIS), (ii) Greenland peripheral glacier and ice cap (GrPGIC), (iii) Ellesmere and Devon islands (ED), and (iv) Baffin Island (BA; Fig. 3).The Citterio and Ahlstrøm (2013) glacier inventory classifies GrPGIC ice fraction as the ice coverage associated with glaciers and ice caps demonstrating "no" and "weak" connectivity with the ice sheet proper, while the ice coverage associated with glaciers demonstrating "strong" connectivity with the ice sheet proper is classified as GrIS ice fraction (cf.Rastner et al., 2012).

Algorithm
We invert a 1000 simulation ensemble of GRACE-derived rate of cryospheric mass change.In each simulation within the ensemble, the best estimate ṀG field is randomly perturbed within its associated error (δ ṀG ) to yield a unique ṀG + δ ṀG input field.An iterative approach is then used to update the higher resolution (26 km) rate of mass change field ( ṁ) and its corresponding lower resolution (Gaussian smoothed) rate of mass change field ( Ṁ), with a randomly perturbed variant of the difference between input and Ṁ fields.This continues until the difference becomes sufficiently small between iterations for a system of three key equations to be considered converged (Supplement Animation 1).Thus, we adopt a Monte Carlo approach to perturb the input data across simulations, as well as to perturb the iterative convergence within a given simulation (Fig. 4).Below we describe this iterative inversion algorithm in detail.
A given Monte Carlo simulation is initialized with a spatially variable ṁij field (where i and j are node indices in Cartesian coordinates).This initial ṁij field is comprised of   an array of random numbers uniformly distributed between −100 and +100 kg m −2 a −1 that has been multiplied by fractional ice area (F ij ).The initial ṁij field varies over the same order of magnitude as the anticipated final ṁij field.Proportionately weighing the initial ṁij field by F ij allows a reasonable representation of uncertainty in the interior of the Greenland ice sheet (where F ij = 1) while also ensuring that ṁij → 0 where F ij → 0. Given the relatively limited spatial extent of the inversion domain, and that non-ice-containing nodes constrained by boundary conditions completely surround the transient ice-containing nodes within the inversion domain (Sect.2.3), the final ensemble mean inversion is relatively insensitive to the choice of initial conditions over the range ±0 to ±100 kg m −2 a −1 (Sect.4.1).
In each iteration, the Gaussian smoothed rate of mass change field ( Ṁij ) is calculated by applying a fixed parameter Gaussian filter function (f G ) to the inverted rate of mass change field ( ṁij ): where k denotes a given iteration, and σ is the prescribed characteristic scaling length (or standard deviation) of the isotropic Gaussian filter.While other filters, such as data adaptive cosine windows, conserve more spherical harmonic low-degree energy, thereby potentially allowing a more accurate description of spatial variability, we do not explore alternatives to the conventional Gaussian filter in this study (e.g., Longuevergne et al., 2010).
Relative to a given node, the spatial distribution of Gaussian filter weight across the inversion domain (W ij ; dimensionless) is described by where d ij represents the distances between the given node and all nodes within the inversion domain.The numerator units in the above equation are implicitly the same as those of characteristic scaling length and inter-node distance (i.e., "1 km" when σ and d are in km).As d ij and W ij are inherently unique at each node, they are functionally threedimensional arrays of the form d ijp and W ijp , where p is a unique node coordinate ranging from 1 to i • j .Thus, applying a Gaussian filter to ṁij in order to update Ṁij requires looping through the p coordinate, from 1 to i • j , which is a computationally expensive task.We determine the optimum value of σ through a sensitivity analysis (Figs. 5 and 6).We invert ṁ fields using characteristic Gaussian length scales of 150, 200 and 250 km.We find that the root mean square (RMS) of the difference field between the GRACE-derived ṀG field and the Gaussian smoothed inverted field ( Ṁ) reaches a minimum when σ = 200 km, independent of prescribed boundary conditions described in Sect.2.3 (Fig. 7).We therefore prescribe a characteristic Gaussian length scale of 200 km in our inversion.The combination of a degree-60 spherical harmonic representation and a characteristic Gaussian length scale of 200 km should preserve maximum information of the magnitude and spatial distribution of mass changes through the inversion process while honoring the fundamental spatial resolution of the GRACE satellites.In a given iteration, the difference between the GRACEderived rate of mass change and the Gaussian smoothed inverted rate of mass change fields, denoted k ij , is determined as where δ ṀG ij R represents a perturbation of the input GRACE solution ( ṀG ij ) within its associated 1σ error field, and Ṁk ij represents the Gaussian smoothed ṁk ij field in a given iteration.In each simulation R is a random scalar derived from a normal distribution centered on zero with a standard deviation of one.R varies across simulations, but is constant throughout the iterations of a given simulation.Over an ensemble of simulations, any given region (or ice type) is thus inverted under relatively high and low initial rates of mass change.
The preceding difference field is used to inform the ṁij field of the subsequent iteration according to where R k ij is a spatially variable array of random values between 0 and 1, and F ij is observed fractional ice coverage at a given node.In contrast to the scalar R used to perturb Eq. ( 3), the array R k ij is continually being repopulated by random numbers throughout the iterations of a given simulation (i.e., R ij is not constant throughout the iterations of a given simulation).Over a large number of iterations, a given node is thus perturbed with relatively high and low values during convergence towards a final ṁ value.
Local mascon rates of mass change derived from GRACE are known to be sensitive to the a priori imposition of predefined patterns of mass change (Horwath and Dietrich, 2009).This is not an issue with the mascon solution used here, however, as Luthcke et al. (2013) derive mass changes from the rigorous reduction of the KBRR residuals.In the context of inverting the spherical harmonic representation, employing a random number field in each iteration not only ensures that inferred rates of mass change are not required to be spatially correlated (i.e., subject to a prescribed covariance matrix), but actually enhances the algorithm ability to explore the infinite number of possible solutions efficiently (e.g., Colgan et al., 2012).
Introducing fractional ice coverage information in the fashion of Eq. (4) forces cryospheric rates of mass change to be proportional with the fractional ice coverage at a given node, rather than restricting inferred rates of mass change to a binary field of non-ice-covered and ice-covered nodes (e.g., Barletta et al., 2013).For example, 10 times more mass change is attributed to a node with F = 1.0 than a neighboring node with F = 0.1.While in reality specific mass loss (i.e., mass loss per unit ice-covered area) typically increases to a maximum at the peripheral nodes of ice masses, incor- porating an additional piece of new information would be required in order to distinguish sharp contrasts in mass change between adjacent ice-containing nodes.
We iterate the ṁij field by cycling through Eqs. ( 1)-( 4) until the following condition is satisfied across all GrIS and GrPGIC nodes (Animation 1): where A ij is node area.We therefore deem the system of equations as converged when the total rate of mass change over all Greenland ice varies by less than 0.1 Gt a −1 between iterations.This typically takes between 75 to 100 iterations per simulation, depending on prescribed σ value (Fig. 8).We perform 1000 simulations in order to compute a robust ensemble mean ṁij field.

Boundary conditions
Here we describe the boundary conditions imposed at nonice-containing nodes.By employing a spherical harmonic representation that isolates the mass change associated with terrestrial ice, the mass changes associated with non-icecontaining nodes (where F = 0) are theoretically negligible (i.e., ṁ = 0 kg m −2 a −1 ) .In practice, the mass changes at these non-ice-containing nodes are not truly zero, but rather within uncertainty of zero.We therefore perform a sensitivity study in which we allow ṁij values at non-icecontaining nodes to vary below prescribed absolute threshold values of 0, 15 or 30 kg m −2 a −1 (Figs. 5 and 6).A threshold of 0 kg m −2 a −1 , for example, corresponds to the theoretical expectation of no error in, or perfect isolation of, the cryospheric mass change signal.In comparison to determining the optimal characteristic Gaussian length scale (σ ) via an analogous sensitivity analysis, selecting which non-ice-containing node ṁij threshold to implement is more subjective.The RMS of the difference field between the GRACE-derived ṀG field and the Gaussian smoothed inverted field ( Ṁ) decreases as the non-ice-containing node ṁij threshold increases.RMS would ultimately go to zero when the ṁij values permitted at non-ice-containing nodes are indistinguishable from those at ice-containing nodes (Fig. 7).
We therefore arbitrarily prescribe an absolute threshold ( ṁmax ) of 15 kg m −2 a −1 at non-ice-containing nodes.This boundary condition acknowledges a level of uncertainty in the rate of mass change at non-ice-containing nodes that is representative of the uncertainty typically assessed for GRACE-derived cryosphere-attributed spherical harmonic solutions (Velicogna and Wahr, 2005;Longuevergne et al., 2010), and an order of magnitude less than the ṁij values inferred by our Monte Carlo inversion approach at adjacent icecontaining nodes.This boundary condition is implemented at non-ice-containing nodes according to This modified version of Eq. ( 4) is invoked at non-icecontaining nodes by a Heaviside, or unit step, function (H ij ) of the following form (e.g., Colgan et al., 2012):

Domain
Here we describe the grid spacing and extent of our inversion domain.In the US National Snow and Ice Data Center (NSIDC) polar stereographic projection, our inversion domain extends from −1625 km in the west to 1300 km in the east, and from −125 km in the north to −3800 km in the south.This places the domain boundaries at least one characteristic Gaussian length scale from all major ice masses in Greenland and the Canadian Arctic Archipelago.There is an inherent trade-off between computational burden and horizontal grid resolution, as the former exponentially increases as the latter linearly decreases.We prescribe a uniform grid spacing of 26 km, which results in 113 computational nodes along the easting axis and 142 computational nodes along the northing axis, for a total of 16 046 computational nodes within the model domain.
While grid spacing is a uniform 26 km throughout our domain, the polar stereographic projection inherently introduces increasing distortion away from its central meridian (45 • W) and parallel (70 • N), which influences area calculations.To compensate for area distortion, we use the true area of each individual node across the domain (A ij ) in our calculations.Thus, while our nominal node area is 26 2 km 2 , the true areas of ice-containing nodes vary between 24.5 2 and 27.5 2 km 2 over the domain.26 km is the minimum resolution

Partitioning mass change and uncertainty
Here we describe how the rate of mass change, and corresponding uncertainty, for each of the four ice types defined in Sect.2.1 (GrIS, GrPGIC, ED and BA) is derived from the inverted ṁ field.By constraining inferred mass changes only to occur at ice-containing nodes, mass changes may be partitioned amongst the distinct ice types within the inversion domain.As nodes containing ED or BA ice type are geographically separate from nodes containing other ice types, the total rate of mass change associated with ED or BA ice type may be quantified via a straightforward summation of the individual rates of mass change at all nodes of either ice type.In contrast, GrIS and GrPGIC ice types can coexist at the same node; the areas of many nodes around the periphery of the Greenland ice sheet are partly covered by both ice sheet and peripheral glaciers.When attributing rates of mass change to either GrIS or GrPGIC ice type at nodes shared by both ice types, we weight attributable rates of mass change by the fractional ice coverage of each ice type.For example, the total rate of mass change for GrIS ( ṁGrIS ) is summed as where A ij is again an array of node true areas, and the superscript denotes GrPGIC and GrIS ice types.While the rates of mass change contributions of GrIS and GrPGIC ice types at a given node are indistinguishable to the inversion algorithm, in that they contribute to inferred mass change in the same way, Eq. ( 8) provides a framework for statistically partitioning ice sheet and peripheral glacier mass change at common nodes.Given the historical expectation that smaller peripheral glaciers and ice caps should have shorter response times than the larger ice sheet to contemporary climate change (Nye, 1960;Jóhannesson et al., 1989), this assumption may underestimate the true contribution of GrPGIC ice type to the inferred mass change at a given node.Simply put, the GrPGIC : GrIS mass change at a given node likely exceeds the GrPGIC : GrIS ice coverage.
After partitioning rates of mass change between distinct ice types, it is desirable to quantify the uncertainty associated with the rate of mass change inferred for each ice type.Previous glaciological applications of Monte Carlo inversion, such as inferring basal sliding velocity from input data of surface velocity observations (e.g., Chandler et al., 2006), and inferring past surface temperature history from input data of borehole temperature profiles (e.g., Muto et al., 2011), have used the perturbation of input data within their associated uncertainty subsequently to estimate uncertainty in the inverted field.We similarly take uncertainty in ṁ at any given node ij (δ ṁij ) as one standard deviation of the given node ṁij values across the ensemble of simulations (Fig. 9).As values of ensemble mean δ ṁij are spatially correlated, we do not derive the uncertainty in the rate of mass change of a given ice type from the point uncertainties of that ice type, but rather from the total rate of mass change of the given ice type across the ensemble of inversions.We take uncertainty in the total rate of mass change of each ice type of interest (GrIS, GrPGIC, ED and BA) as one standard deviation in the total rate of mass change of each ice type across the ensemble of 1000 simulations.

Results
The ensemble mean Gaussian smoothed rate of mass change field closely reproduces the spherical harmonic representation of cryosphere-attributed mascons derived from GRACE observations (Fig. 10).Both the observed and ensemble mean simulated fields exhibit similar patterns of Greenland mass loss, extending from the Geikie Plateau in East Greenland around South Greenland to Humboldt Glacier in Northwest Greenland, as well as mass gain focused in north central Greenland.The difference between the final Ṁ field and the input ṀG field ( ) is < 5 % of the absolute magnitude of the ṀG field throughout the vast majority of the inversion domain.The largest discrepancy occurs in northern Baffin Bay, where the pattern of the spherical harmonic representation over the ocean infers more mass loss than is permitted by the inversion parameters employed, resulting in a local maximum in absolute discrepancy.The nominal 26 km by 26 km resolution inverted mass change field ( ṁ) infers mass loss around the entire low-elevation perimeter of the GrIS, with a broad region of mass gain in the high-elevation central and North Greenland ice sheet interior, similar to the spatial pattern of mass change presented by Barletta et al. (2013;Fig. 11).This pattern is consistent with independent observations of concurrent high-elevation ice sheet thickening and low-elevation ice sheet thinning (e.g., Zwally et al., 2011).As the inversion is designed to eliminate signal leakage between irregularly shaped non-ice-covered and ice-covered areas (i.e., redistribute the mass changes underlying a cryosphere-attributed spherical harmonic representation across only ice-containing nodes), it does not resolve spatial heterogeneity in mass changes between ice-containing nodes (i.e., at the scale of individual glaciers and ice caps).In the Canadian Arctic, all ice-containing nodes exhibit mass loss, except for a small portion of the southeast tip of Baffin Island, where the ensemble mean inversion infers a small area of mass gain.Similar to the local maximum in absolute discrepancy in northern Baffin Bay, we interpret this small region of mass gain on Baffin Island as an artifact of the inversion satisfying the spherical harmonic solution over non-ice-containing nodes.As expected, the Canadian Arctic mass loss is concentrated at numerous discrete ice masses, where F ij → 1.The magnitudes of the rate of mass change inferred by the inversion are broadly consistent with independent altimetryderived estimates of Greenland's high-elevation mass gain (up to 200 kg m −2 a −1 ) and low-elevation mass loss (down to −500 kg m −2 a −1 ) at similar spatial resolution (e.g., Zwally et al., 2011).The accompanying spatial uncertainty field (δ ṁij ) takes on values of ∼ 70 kg m −2 a −1 in the interior of ice masses, increasing to a maximum of ∼ 120 kg m −2 a −1 towards the periphery of the Greenland ice sheet (Fig. 9).
The Monte Carlo inversion approach infers a total Greenland mass change of −251 ± 25 Gt a −1 over the December 2003 to December 2010 period.We note that the absolute magnitude of this mass change is dependent on the mascon solution underlying the input spherical harmonic representation, and therefore does not constitute an independent assessment of GrIS mass loss.Our statistical partition based on fractional ice coverage information infers that the GrIS was responsible for −218 ± 20 Gt a −1 of mass change, while GrPGICs were responsible for −34 ± 5 Gt a −1 of mass change (Fig. 12).This suggests that peripheral glaciers accounted for ∼ 14 % of Greenland's mass change over the study period.Using satellite laser altimetry and a firn compaction model, Bolch et al. ( 2013) estimated an October 2003 to March 2008 mass loss of 28 ± 11 Gt a −1 from GrPGICs substantially independent from the ice sheet, which had an ice-covered area of ∼ 89 × 10 3 km 2 (cf. ∼ 88 × 10 3 km 2 used in this study; Citterio and Ahlstrøm, 2013).When Bolch et al. ( 2013) also included GrPG-ICs with "strong connectivity" (∼ 130 × 10 3 km 2 ), they found a mass loss of 41 ± 16 Gt a −1 over the same period.Using a similar satellite laser altimetry-based approach, Gardner et al. ( 2013 (Gardner et al., 2011), and a more recent altimetry and gravimetry reconciled estimate of 60 ± 8 Gt a −1 over the October 2003 to October 2009 period (Gardner et al., 2013).

Limitations
Individual studies of GRACE-derived total Greenland rate of mass change, over comparable time periods, range from −195 ± 30 Gt a −1 (from August 2003 to July 2009; Pritchard et al., 2011) to −234 ± 29 Gt a −1 (from August 2002 to September 2011; Sasgen et al., 2012).This range in total Greenland mass change inferred by various GRACE studies has been attributed to differences in both temporal sampling period and the treatment of signal corrections (Schrama and Wouters, 2011).The IMBIE has recently reconciled apparent differences between individual studies, and demonstrated that the NASA GSFC mascon solution we employ agrees within the error bars of other GRACE-derived mass change estimates over the IMBIE period (Shepherd et al., 2012).We note, however, that fully iterated NASA GSFC mascon solutions, including those employed in our present study, have  been shown to produce greater mass loss values in comparison to non-iterated NASA GSFC mascon solutions (Luthcke et al., 2013).We therefore acknowledge that initializing the Monte Carlo inversion with a different spherical harmonic representation would yield a different inferred mass loss distribution.Recognizing this limitation, we contend that our present study demonstrates the utility of combining spherical harmonics with additional information in the form of fractional ice coverage.By explicitly constraining cryosphereattributed mass changes to irregularly shaped ice-covered areas, we purport to have completely eliminated signal leakage between non-ice-covered and ice-covered areas.GRACE mascons are conventionally expressed in units of water equivalent thickness per time, under the implicit assumption of F = 1 at all nodes containing terrestrial ice mass (e.g., Barletta et al., 2013).We note that the cryosphereattributed mass change per unit area we describe ( ṁij ) is not equivalent to specific ṁij (i.e., cryosphere-attributed mass change per unit ice-covered area).We may derive equivalent specific ṁij by dividing ṁij by F ij (Fig. 13).The relatively smooth gradients in the specific ṁij field are constrained by the fundamental resolution of the GRACE satellites.While true specific mass loss typically increases to a maximum at the peripheral nodes of ice masses, the inversion we present would need a secondary piece of independent information in order to distinguish sharp contrasts in mass change between adjacent ice-containing nodes.Since fractional ice coverage is the only new information we apply to GRACE data, the inversion only constrains cryosphereattributed mass changes to ice-covered areas, and is not capable of distinguishing spatial heterogeneity in mass loss between adjacent ice-containing nodes.Altimetry data would be a logical second piece of independent information to incorporate, in order to distinguish differential rates of mass change between adjacent ice-containing nodes.Finally, the precise role of initial conditions in interpreting inferred mass changes is not entirely clear.As a sensitivity test, we initialized the algorithm with a spatially uniform ṁij field of 0 kg m −2 a −1 , rather than with random numbers uniformly distributed between −100 and +100 kg m −2 a −1 times fractional ice coverage.Under these revised initial conditions, an ensemble of 1000 simulations produces a virtually identical total Greenland mass loss (250 ± 26 Gt a −1 ) and spatial pattern (not shown here, but available in Colgan et al., 2013).While total Greenland mass loss is therefore insensitive to initial conditions over the range ±0 to ±100 kg m −2 a −1 , the spatial uncertainty field does appear to be sensitive to selection of initial conditions.When mass changes are initialized with zeros, local uncertainty approaches zero in the high-elevation ice sheet interior (Fig. 14).Thus, deriving both ice-sheet-wide and local uncertainties from ensemble spread at convergence is not valid under all initial conditions, as ensemble spread is demonstrably dependent on choice of initial conditions.We contend that by imposing initial conditions that vary over the same order of magnitude as the anticipated inversion field, we have minimized the potential influence of initial condition hysteresis.We acknowledge, however, that we do not yet sufficiently understand the behavior of this approach to encourage other researchers to blindly apply analogous inversions in different contexts.Application of Monte Carlo-type inversion to any system requires a deliberate selection of initial and boundary conditions characteristic of the system under consideration.

Implications
The recent global mass loss trend of small glaciers and ice caps external to Greenland and Antarctica has been estimated as between 148 ± 30 Gt a −1 (over January 2003 to December 2010; Jacob et al., 2012) and 215 ± 26 Gt a −1 (over October 2003to October 2009;Gardner et al., 2013).Our estimate of GrPGIC mass change over a similar period is between ∼ 16 and 23 % of these values.Projecting the future sea level rise contribution of peripheral glaciers in Greenland and Antarctica, the 2007 Intergovernmental Panel on Climate Change report that "the global [glaciers and ice caps] sea level contribution [was] increased by a factor of 1.2 to include [peripheral glaciers] in Greenland and Antarctica" (Meehl et al., 2007).Our results suggest a comparable factor (e.g., ∼ 1.2) would be required to account for GrPGICs alone when using this scaling technique.We regard our GrPGIC contribution as a lower bound, however, due to (i) the potential underestimation of GrPGIC ice extent.The data set of Greenland ice sheet and peripheral glaciers we employ (Citterio and Ahlstrøm, 2013) does not classify as a GrPGIC those glaciers demonstrating "strong" connectivity with the ice sheet proper (e.g., Geikie Plateau or Julianehåb Ice Cap; cf.Rastner et al., 2012), and may therefore be regarded as a conservative estimate of GrPGIC extent.(ii) There is an implicit assumption that mass changes are equally weighted between GrIS and GrPGIC ice at common nodes (Eq.8).Given an apparent equilibrium line ascent by over 500 m between 1994 and 2012 in West Greenland (McGrath et al., 2013), and the historical expectation that smaller glaciers have shorter response times to such climatic perturbations than the larger ice sheet (Nye, 1960;Jóhannesson et al., 1989), we suggest that it would be reasonable to weight preferentially the total mass change inferred at common GrPGIC and GrIS nodes towards GrPGIC mass change.With GrPGICs comprising < 5 % of Greenland's ice-covered area but ∼ 14 % of its total mass loss, the specific contribution of GrPGIC to Greenland total mass loss far exceeds the specific contribution of the GrIS.
The substantial contribution of GrPGIC to Greenland's total mass loss highlights that GRACE-derived estimates of total Greenland mass loss (i.e., GrIS + GrPGIC) should only be compared with mass loss estimates derived by volumetric or input-output approaches that similarly sample both the ice sheet proper and peripheral glaciers.By sampling the GrPGICs in addition to the GrIS, GRACE-derived mass loss estimates can be expected to have ∼ 35 Gt a −1 greater mass loss than volumetric or input-output methods that only sample the ice sheet proper (cf.Alley et al., 2007;Pritchard et al., 2010).While GRACE implicitly samples the mass changes associated with GrPGICs, whether or not these other methods include GrPGICs is dependent on the ice mask employed, which can vary tremendously from study to study (Vernon et al., 2013).Our partition between GrIS and GrPGIC mass loss suggests a first-order correction could be made by multiplying GRACE-derived total Greenland mass loss values by 0.86 to estimate the GrIS contribution.Similarly, multiplying GRACE-derived total Greenland mass loss by 0.14 provides a first-order estimate for the GrPGIC contribution.This partition ratio, however, is likely time-variant,

Summary remarks
Our iterative inversion, in which inferred mass changes are constrained only to occur at ice-containing nodes, produces a ground-level cryosphere-attributed rate of mass change field that has a spatial distribution and magnitude that is fully consistent with an input GRACE-derived spherical harmonic representation, and thus spherical harmonic coefficients, without assuming that rates of mass change are constant within or across pre-defined regions (e.g., delineated drainage systems).While this approach can technically be used to invert any given mass change represented in the form of spherical harmonic coefficients (e.g., GRACE L2-derived trends) within the constraint of ground-level spatial distribution, we caution researchers from blindly applying such an algorithm in different contexts without taking appropriate precautions (e.g., careful selection of initial and boundary conditions).The method we have presented is unconventional, and its ultimate degree of success in cryospheric application remains to be verified through forward modeling of either observed GRACE KBRR data or analogous simulated synthetic data.
In the context of Greenland mass change, the inference that GrPGICs, which comprise < 5 % of Greenland's ice covered area, are contributing to ∼ 14 % of Greenland's total mass loss observed by satellite gravimetry, highlights that GRACE-derived estimates of "Greenland" mass loss cannot reasonably be taken as synonymous with "Greenland ice sheet" mass loss.Comparisons of GRACE-derived mass loss should therefore be limited to other mass balance techniques (i.e., altimetry or input-output) that sample both the ice sheet and its peripheral glaciers, or alternatively GRACE-derived mass loss values should be adjusted to account for GrPGIC mass loss.
While the inferred mass change field we produce offers the potential to compare GRACE-derived estimates of cryospheric mass change with other observational and modeled spatial data sets, such as surface mass balance estimates, at higher spatial resolution than previously possible, we note that 26 km is nominal resolution across the inversion domain.Incorporating further independent information, such as altimetry-derived ice surface elevation changes, offers the potential to distinguish trends in mass change between adjacent ice-containing nodes better.This may permit cryosphere-attributed mass changes to be assessed at 26 km (or better) actual spatial resolution, and thus provide a more confident assessment of the magnitude and spatial distribution of the cryospheric mass change than we present here.

Fig. 2 .
Fig. 2. (A) Spherical harmonic representation of the cryosphereattributed rate of mass change ( ṀG ij ) derived from GRACE over the December 2003 to December 2010 period (Fig. 1b).(B) Accompanying 1σ trend error field (δ ṀG ij ).Outline of the Greenland ice sheet shown for reference.

Fig. 4 .
Fig. 4. Flowchart overview of inverting a given lower (spherical harmonic) resolution GRACE-derived rate of mass change field ( ṀG ij ) into an ensemble of higher (26 km) resolution rate of mass change fields ( ṁij ).An ensemble of 1000 simulations is performed, with each simulation comprised of an iterative inversion to convergence as defined by Eq. (5).

Fig. 8 .
Fig.8.The system of equations in each simulation is deemed converged when the inferred total rate of mass change over all Greenland ice coverage (both GrIS and GrPGIC) varies by less than 0.1 Gt a −1 between iterations (Eq.5).The whisker plot at the right shows one (thick line) and two (thin line) standard deviations from the ensemble mean.

Fig. 9 .
Fig. 9.The uncertainty in rate of mass change at any given icecontaining node ij (δ ṁij ) is taken as one standard deviation of the given node rate of mass change values across the ensemble of 1000 simulations.

Fig. 10 .
Fig. 10.(A) Ensemble mean inferred rate of mass change field ( ṁ) at ice-containing nodes.Grey shading denotes non-ice-containing nodes and the white contour line denotes 0 kg m −2 a −1 .(B) Gaussian smoothing (σ = 200 km) of the ensemble mean inferred rate of mass change field ( Ṁ). (C) Spherical harmonic representation of the trend in GRACE-derived cryosphere-attributed rate of mass change ( ṀG ).(D) Difference between ṀG and Ṁ ( ).Black contour lines denote irregularly shaped ice-containing nodes within the domain.

Fig. 11 .
Fig. 11.Inverted rate of mass change constrained to irregularly shaped ice-covered areas in Greenland and the Canadian High Arctic (Fig. 10a).Color scale saturates at −750 and +150 kg m −2 a −1 .Map extent follows Fig. 1, with grey shading denoting areas beyond the inversion domain.
) assessed a GrPGIC mass loss of 38 ± 7 Gt a −1 over the October 2003 to October 2009 period.Our inversion attributes a total Canadian High Arctic mass change of −48 ± 8 Gt a −1 , comprised of a contribution of −28 ± 4 Gt a −1 from Ellesmere and Devon islands and −20 ± 4 Gt a −1 from Baffin Island.The error bars of this Canadian Arctic Archipelago mass loss overlap with those of a combined altimetry, gravimetry and input-output estimate of 61 ± 7 Gt a −1 over the March 2003 to December 2009 period

Fig. 12 .
Fig. 12. Ensemble probability density functions of rate of mass change inferred by Monte Carlo inversion over the Greenland ice sheet (GrIS), Greenland peripheral glaciers and ice caps (GrPGICs), Ellesmere and Devon islands (ED) and Baffin Island (BA) over the December 2003 to December 2010 period.

Fig. 13 .
Fig. 13.Specific rate of cryospheric mass change per unit ice area ( ṁij divided by F ij ).

Fig. 14 .
Fig. 14.Sensitivity analysis of the spatial distribution of uncertainty in inferred mass changes (δ ṁij ).(A) Uncertainty at a given node when the inversion is initialized with a constant array of 0 kg m −2 a −1 .(B) Uncertainty at a given node when the inversion is initialized with an array of random numbers uniformly distributed between −100 and +100 kg m −2 a −1 times fractional ice coverage (identical to Fig. 9).(C) Difference.

W.
Colgan et al.: Constraining GRACE signal to ice-covered areas and therefore only pertains to the December 2003 to December 2010 study period.