Application of GRACE to the evaluation of an ice flow model of the Greenland Ice Sheet

Quantifying the Greenland Ice Sheet’s future contribution to sea level rise is a challenging task that requires accurate estimates of ice flow sensitivity to climate change. Forward models of ice flow dynamics are promising tools for estimating future ice sheet behavior, yet confidence is low because evaluation of historical simulations is so challenging due to the scarcity of highly-resolved (spatially and temporally) continental-wide validation data. Recent advancements in processing of Gravity 5 Recovery and Climate Experiment (GRACE) data using Bayesian-constrained mass concentration ("mascon") functions have led to improvements in spatial resolution and noise reduction of estimated monthly global gravity fields. Specifically, the Jet Propulsion Laboratory’s JPL RL05M GRACE mascon solution (GRACE-JPL) now offers an opportunity for ice sheet model evaluation within independently resolved 300 km mascons. Here, we investigate how Greenland Ice Sheet mass balance captured through observations GRACE-JPL differs from that simulated by the ice flow model the Ice Sheet System Model 10 (ISSM). For the years 2003-2012, ISSM is forced with regional climate model (RCM) surface mass balance (SMB), and resulting mass balance is directly compared against GRACE-JPL within individual mascons. Overall, we find good agreement in the Northeast, Southwest, and the interior of the ice sheet, where mass balance is primarily controlled by SMB. In the Northwest, seasonal amplitudes match well, but trends in ISSM are muted relative to GRACE-JPL. In the Southeast, GRACE-JPL exhibits larger seasonal amplitude than that predicted by SMB while simultaneously having more pronounced trends. These 15 results indicate that discrepancies in the Northwest are controlled by changes in ice dynamics that are not currently modeled by ISSM, i.e. transient processes driven by ice sheet hydrology and ice-ocean interaction, while discrepancies in the Southeast are controlled by a combination of these missing dynamics and errors in modeled SMB. Along the margins, we find that transient dynamics are responsible for consistent intra-annual variations in regional mass balance that ultimately contribute to the steeper negative mass trends observed by GRACE-JPL. Consequently, ice-ocean interactions and hydrologically-driven processes at 20 relatively high (monthly-to-seasonal) temporal resolutions must be considered for improving upon ice flow models. 1 The Cryosphere Discuss., doi:10.5194/tc-2015-224, 2016 Manuscript under review for journal The Cryosphere Published: 18 January 2016 c © Author(s) 2016. CC-BY 3.0 License.


Introduction
The Greenland Ice Sheet is a significant source of sea level rise (Gardner et al., 2013;Sto), and its rate of contribution is expected to accelerate in the coming centuries (Pachauri et al., 2014;Church andWhite, 2006, 2011).The quantification of Greenland's future contribution to sea level rise is a challenging task, which requires accurate estimates of ice flow sensitivity to climate change.Often, current observable trends in ice sheet mass balance (MB) are extrapolated (Shepherd and Wingham, 2007;Velicogna, 2009;Rignot et al., 2011) in order to estimate future changes to sea level.An alternative to this method is to take advantage of numerical ice flow models that have been validated against historical data to run forward simulations of the future.
These transient ice flow models are the most promising tools for estimating future ice sheet behavior.However, they are not yet capable of capturing the observed trends in polar ice sheet mass balance.In fact, confidence in these types of future projections is low, especially because evaluation of model historical runs is so challenging due to the scarcity of continentalwide data for evaluation.Behavior also varies dramatically between models (Bindschadler et al., 2013;Nowicki et al., 2013a, b), as each model has its own spinup procedure and a different treatment of boundary conditions/model forcing, implementation of sliding laws, and choice of ice flow equations.As a consequence, model-based estimates of Greenland's historical sensitivity to climate change are plagued with significant spread and large uncertainties.
Main sources of ice flow model uncertainty include (1) surface mass balance (SMB) forcing, (2) limitations to model spinup, and (3) missing dynamical processes, i.e. basal lubrication due to surface runoff reaching the bed (Schoof, 2010;Tedesco et al., 2012), warming of ice due to runoff refreeze (Fau; Phillips et al., 2010Phillips et al., , 2013)), and ice/ocean interaction (including grounding line retreat (Holland et al., 2008;Rignot et al., 2010) as well as changes to flow resistance at the calving front (Walter et al., 2012)).For the Greenland Ice Sheet, recent studies suggest that SMB variability (particularly surface melt) (Enderlin et al., 2014) is the primary contributor to the observed trend in MB, while changes to ice flow dynamics are secondary.Still, the partitioning of MB amongst these sources remains an ongoing challenge (Howat et al., 2007;van den Broeke et al., 2009;Pritchard et al., 2009;Kjaer et al., 2012;Yan et al., 2013;Khan et al., 2014) and a hindrance to the improvement in modelbased projections of future sea level.
Another pressing issue is the lack of observational data for model evaluation.In response to the need for such data, government agencies have deployed a number of instruments for the purpose of monitoring the MB of the polar ice sheets through satellite altimetry, interferometry, and gravimetry (Pritchard et al., 2009;Zwally et al., 2011;Rignot et al., 2011;Shepherd et al., 2012;Velicogna and Wahr, 2013).For more than a decade, the joint U.S.-German Gravity Recovery and Climate Experiment (GRACE) has continuously acquired time-variable measurements of the Earth's gravity field and has provided unprecedented surveillance of MB of the Polar ice sheets (e.g., Wu et al., 2002;Luthcke et al., 2006;Baur et al., 2009;Velicogna, 2009;Wu et al., 2009;Chen et al., 2011;Schrama and Wouters, 2011;Jacob et al., 2012;Sasgen et al., 2012;Velicogna and Wahr, 2013;Schrama et al., 2014).GRACE data are typically processed by estimating gravity field variations using unconstrained spherical harmonic basis functions.These estimates ultimately suffer from a highly correlated error structure in the form of longitudinal stripes, and a variety of methods has been developed to remove these artifacts.Recent advancements in GRACE placement was arbitrary in the derivation of the JPL RL05M solution and was not optimized for any specific application of the GRACE data (i.e.recovering Greenland mass variations).Note also that the estimate of mass in each mascon can be considered relatively independent of the other mascons, as we do not apply any apriori spatial correlation between mascons (this choice is appropriate for high latitudes due to dense ground track coverage), and the formal posteriori covariance matrix indicates small correlations between adjacent mascons after the inversion.This solution has been shown to agree with previously published results (within formal error bars) regarding the total rate of mass loss from the Greenland Ice Sheet (Watkins et al., 2015).
In this analysis we examine mass changes over the ice sheet with the native resolution of the gravity solution (i.e.individual mascons), constituting the highest spatial resolution analysis of the Greenland Ice Sheet from GRACE data thus far.
The GRACE data have certain known limitations, and as such, we apply standard post-processing procedures to correct for these.The C 20 coefficient (defining the oblateness of the Earth), which is poorly observed by GRACE, has been substituted with an estimate derived from Satellite Laser Ranging (SLR) data (Cheng and Tapley, 2004).GRACE does not observe movement of the center of mass of the Earth, since the satellites orbit this point at all times; as such, we use an estimate of geocenter motion from Swenson et al. (2008).The position of the Earth's mean pole has been corrected using the recommendation of Wahr et al. (2015).Known jumps in the background atmosphere dealiasing products (occurring in 2006 and 2010) have been accounted for by using provided GAE and GAF products.Finally, the solid Earth glacial isostatic adjustment (GIA) signal has been removed from the GRACE data using the model provided by A et al. (2013), in an effort to isolate only surface mass variations for a direct comparison against the ice sheet model.
One post-processing algorithm that is unique to the JPL mascon solution concerns the treatment of leakage errors.Many of the Greenland mascons lie on both land and ocean regions (Fig. 1), and as such, the solution for these mascons will contain the average mass of both the land and ocean.We apply a Coastline Resolution Improvement (CRI) filter to the data (Watkins et al., 2015) to separate the land and ocean portions of each of these mascons.As such, in all analyses presented here, the ocean mass from each mascon has been removed, and we are analyzing only the land component of each mascon.For further details on the JPL RL05M mascon solution and the CRI filter, the reader is referred to Watkins et al. (2015).The data are publicly available at www.grace.jpl.nasa.gov.From hereafter, we refer to the JPL RL05M GRACE mascon solution as GRACE-JPL.

Ice flow model
The Ice Sheet System Model (ISSM) is a thermo-mechanical finite-element ice flow model.It relies upon the conservation laws of momentum, mass, and energy, combined with constitutive material laws and boundary conditions.The implementation of these laws and treatment of model boundary conditions are described by Larour et al. (2012).
Here, we simulate the Greenland Ice Sheet with a two-dimensional (2D) Shelfy-Stream Approximation (SSA) (MacAyeal, 1989), implemented within ISSM (Larour et al., 2012).The SSA is based on the Stokes equations, and assumes that vertical shear and bridging effects are negligible.These assumptions reduce the momentum balance equation to the desired 2D sys-The Cryosphere Discuss., doi:10.5194/tc-2015Discuss., doi:10.5194/tc- -224, 2016 Manuscript under review for journal The Cryosphere Published: 18 January 2016 c Author(s) 2016.CC-BY 3.0 License.tem equations.The application and limitations of these equations to the Greenland Ice Sheet using ISSM are discussed by Larour et al. (2012), Schlegel et al. (2013), andSchlegel et al. (2015).

Initialization and Relaxation
The strategy for ISSM Greenland continental initialization, relaxation, and spinup is described in detail by Schlegel et al. (2013Schlegel et al. ( , 2015)).For this study, the anisotropic mesh is made up of 91,490 elements, refined using observed surface (Scambos and Haran, 2002) and velocity (Rignot and Mouginot, 2012) fields (Fig. 2A).Mesh resolution are set to a minimum of 1 km in steep areas with high velocity gradients and to a maximum of 15 km at the ice divides.We initialize the bedrock geometry with 150 m gridded mass-conserving bedrock (Morlighem et al., 2014a) and ice surface data (Howat et al., 2014).Three dimensional ice temperature and ice viscosity are derived from a steady-state higher-order (Blatter, 1995;Pattyn, 2003) three-dimensional formulation of the thermal regime, using observed velocities, surface temperatures, and geothermal heat flux (Larour et al., 2012;Seroussi et al., 2013).Surface temperatures are from Ettema et al. (2009) and geothermal heat flux estimates are from Shapiro and Ritzwoller (2004).We determine the spatially-varying basal drag coefficient using inverse methods (MacAyeal, 1993), following Morlighem et al. (2010), in order to best match the modeled ice surface velocities with InSAR surface velocities (Rignot and Mouginot, 2012).We hold ice viscosity and basal drag constant during the transient simulation.
For relaxation of the model geometry, we run the model to mass balance steady state after Schlegel et al. (2013Schlegel et al. ( , 2015)).To accomplish this, we force the model with 1979-1988 average SMB estimates (SM B) from Box (2013) (Appendix Sect.1), for a total of 30,000 years.This modern period is chosen to serve as our reference climatology because, during this period, the ice sheet was believed to have a mass balance of zero (Rignot et al., 2008).Also during this reference period, all three historical SMB products are defined (Sect.3.3).SMB is interpolated onto the ISSM Greenland mesh and is imposed through a one-way coupling scheme (Schlegel et al., 2013).During relaxation, the total volume of the ice sheet is reduced by 0.5%.The ice surface elevation in the southern ice sheet (below 71 degrees N latitude) increases by a mean of 124 m and decreases by 95 m in the northern ice sheet (above 71 degrees N latitude).The mean change in ice surface elevation is -8 m.

Historical SMB forcing
After relaxation, we spinup the forward model for 173 years, from 1840-2012, using reconstructed monthly SMB after Box (2013).This product, hereafter referred to as BOX, is described in more detail in Appendix Sect. 1.
The 173-year historic run constitutes our base simulation.We plot the total yearly BOX SMB forcing for this run in gray in Fig. 3. On the same figure, in red, we plot the monthly evolution of total ice sheet mass during spinup.From 1840-1900, the model maintains a relative steady-state.After 1900, the overall trend in ice sheet mass balance is dominated by mass loss until 1970, when accumulation over the ice sheet increases.During the following decades -which include the period used as climatological reference period for SM B -through the end of the 1990's, we find that the simulated ice sheet re-achieves a stable condition, growing only slightly over a period of 25 years.
The state of the ice sheet (dictated by ice thickness, bedrock elevation, and ice velocities), at the end of year 1978, is the initial condition for two additional historic simulations.For these simulations, we force the 1978 model state with SMB anomalies The Cryosphere Discuss., doi:10.5194/tc-2015Discuss., doi:10.5194/tc- -224, 2016 Manuscript under review for journal The Cryosphere Published: 18 January 2016 c Author(s) 2016.CC-BY 3.0 License.

GRACE period mass estimates
The purpose of this study is to compare ISSM forward simulations against GRACE-JPL in order to highlight the regions of the modeled ice sheet that differ from GRACE.For this comparison, we attempt to capture the Greenland Ice Sheet's surface features, including topography and surface velocities.The approach of this study is to accomplish an ice sheet topography that is in steady-state with regards to thickness and velocity (Schlegel et al., 2013) through relaxation to SM B. Ice sheet mass balance is comprised of two major components, surface mass balance and discharge (D), such that M B = SM B−D; therefore, after relaxation, the ice flow model is in a state such that the ice outflux from the model margins is equal to SM B. When the relaxed model is forced with historical SMB (Fig. 3), modeled D and MB change with time.Here, changes to D are simulated only through ice dynamic response to perturbations from SMB (Schlegel et al., 2013).Differences between simulated mass change and the GRACE-JPL solution represent missing dynamical processes, limitations to our model spinup, and errors in SMB.
For the comparison with GRACE-JPL, we assess model-based mass changes within all mascons that contain portions of the Greenland land mass (Fig. 1).In this analysis, we consider 2003-2012 SMB anomalies (with respect to SM B), ISSM Greenland Ice Sheet thickness changes, and the temporal evolution of mass beyond the ice sheet margin (hereafter referred to as periphery).The periphery consists of areas of bare rock/tundra, as well as glaciers/ice caps that are not physically attached to the ice sheet.A high resolution (1/120 degree) mask, distinguishing the ice sheet, peripheral ice, and land (Gardner et al., 2013) is relied upon to create the original ISSM domain outline of the ice sheet.We then use this mask to categorize all land within the Greenland mascons (Fig. 1).
Cumulative mass change within each Greenland mascon is determined monthly from 2003-2012.First, over the ice sheet area, we sum the SMB anomalies for each of the forcing products: BOX, MAR, and RACMO.This resultant mass signal for each product represents the SMB forcing for the ISSM historical simulations.Next, we sum the evolution of ISSM Greenland ice thickness for the BOX, MAR, and RACMO historic simulations.This mass signal represents the ISSM response to SMB forcing and model estimates of ice sheet mass balance through time.Note that because GRACE does not capture mass change over floating ice, we remove the mass signal from areas classified as ice shelf (Morlighem et al., 2014a) for this analysis.
Finally, we assess the monthly mass change over the peripheral areas.Because each SMB forcing product represents accumulation and melt differently over the tundra, our analysis methods vary depending on the variables available from each product.Over land areas, we are interested in capturing the annual amplitude of the snowpack load, as we expect such a signal to be captured by GRACE-JPL.Since the BOX product is assessed over the entire continent, and SMB is extended over land, The Cryosphere Discuss., doi:10.5194/tc-2015Discuss., doi:10.5194/tc- -224, 2016 Manuscript under review for journal The Cryosphere Published: 18 January 2016 c Author(s) 2016.CC-BY 3.0 License.we use the provided SMB estimates to determine cumulated snowpack over time.For MAR, we estimate the evolving snow mass using SMB calculated explicitly for the portion of each MAR pixel that is designated as tundra instead of permanent ice (Fettweis et al., 2005).From RACMO, snow mass over tundra is provided as a direct product variable of the model.
For the area covered by peripheral ice, we expect the cumulative mass trend to be captured in the GRACE-JPL signal.For the period 2003-2009, Gardner et al. (2013) ) estimates the total signal of Greenland's peripheral ice to be significant, on the order of −40 Gt/yr.Here, we do not consider the peripheral glaciers and ice caps to be evolving dynamically; therefore, we assume that SMB directly drives the cumulative mass trend.It is important to note that the peripheral ice mask for each SMB product differs, largely because the product grids vary spatially in location and resolution.In addition, the 5-11 km resolution of the products, in many cases, do not properly represent the aerial extent or the topographical features of the peripheral ice.In order to better capture the SMB within these complex areas and to more easily compare their mass balance estimates, we use the 150 m gridded GIMP Digital Elevation Model (DEM) (Howat et al., 2014;Morlighem et al., 2014b) to, separately for each mascon, determine a hypsometric curve for the areas masked by Gardner et al. (2013) as peripheral ice.The curve is binned at every 150 m of surface elevation.For every month, and for each mascon, we plot the SMB of each product separately as a function of elevation, and fit a curve (Gardner et al., 2013).Only SMB values over peripheral ice are considered in these curves.Mascons with similar climates are combined in order to refine the fit.The resulting curve is used to determine the mean SMB value within each elevation bin.Finally, the SMB value within each bin is multiplied by the area of each elevation band, and the results for all elevation bands are summed as the total SMB mass contribution for a particular month.
We acknowledge that this analysis does not capture hydrological processes that may alter the timing of the annual mass cycle captured by GRACE-JPL.The simulation of hydrological processes is beyond the scope of this study, therefore it is important to note that the movement and storage of liquid water within the ice or on the land is present in the GRACE-JPL signal, but not considered in the model-based results presented here.

Quantification of errors and uncertainty
Uncertainty in both the GRACE-JPL mass estimates and the SMB-forced ISSM Greenland simulations are considered in this study.Details on our error assessment are provided below.

Uncertainty in GRACE-JPL
We use the formal posteriori covariance matrix from the gravity field inversion to define the uncertainty in each mascon estimate provided by GRACE-JPL.Typically, the formal covariance matrix is regarded to provide an optimistic estimate of uncertainties, as it is uninformed of certain error sources that affect the GRACE-JPL mass estimates, such as temporal aliasing errors.
We find that this is the case for ocean and land-hydrology regions of the world for which apriori information is derived from geophysical models: the posteriori covariance matrix is too optimistic and must be scaled up to accurately reflect uncertainty.
However, for ice-covered regions, such as Greenland, the apriori information is derived from a bootstrapping methodology from which the magnitude of the K-band range-rate data residuals dictate the spatial variations in the apriori covariance ma-The Cryosphere Discuss., doi:10.5194/tc-2015Discuss., doi:10.5194/tc- -224, 2016 Manuscript under review for journal The Cryosphere Published: 18 January 2016 c Author(s) 2016.CC-BY 3.0 License.trix, and the magnitude of these terms is purposefully left large, to be conservative.As such, we find the resulting posteriori covariance matrix to give an adequate estimate of uncertainty in each mascon.This hypothesis was tested by using spatial variance information from the MAR SMB model to derive an apriori covariance matrix that was used to constrain the GRACE-JPL solution, and analyzing how this impacted the results.Differences in this MAR-constrained solution vs. the relatively unconstrained solution presented here are captured by the formal errors.Furthermore, the posteriori covariance matrix shows adjacent mascons to have small correlations (∼0.2) with each other.As such, we assume all mascons to be uncorrelated with their neighbors.A comparison to ICESat altimetry data (Csatho et al., 2014) validates this assumption, as it agrees extremely well with GRACE-JPL in a mascon-by-mascon comparison framework (Fig. S1).
Leakage errors are considered explicitly by evaluating the expected accuracy of the CRI filter used to separate land and ocean mass components of mascons that lie on coastlines.Simulation results show the CRI filter is effective in reducing leakage errors by greater than 50% globally.Thus, we assume that the estimate of ocean mass for each land/ocean mascon has an error of 50%, and this error is added in a root sum of squares (RSS) to the formal covariance for each mascon.
Finally, since we are ultimately interested in surface mass variations (and as such remove solid Earth Glacial Isostatic Adjustment (GIA) signals using a model), the error in the GIA model must be considered.We take the 1-sigma spread of the ensemble mean of four GIA models.The four models used include the ICE-6G_C (VM5a) model (Peltier et al., 2015), a model by A et al. ( 2013) which uses ICE-5G loading history and a VM2 viscosity profile, a model using ICE-5G loading history and a Paulson viscosity profile (Paulson et al., 2007), and a model by Simpson et al. (2009) which uses the Huy1 (Huybrechts, 2002) ice load history and an independently derived viscosity profile.Using this approach, we derive a GIA uncertainty for the entire Greenland Ice Sheet of ±15 Gt/yr, which matches closely with what is reported in Velicogna and Wahr (2013).The derived uncertainty in GIA is added to the RSS of the formal covariance and leakage error discussed above to arrive at an estimate of uncertainty in surface mass variations for each mascon.Note that in all figures, we show the GRACE-JPL uncertainty increasing linearly with time.This is directly due to uncertainty in the GIA model.Typical GRACE-JPL uncertainty curves are not presented in this fashion; however, since we are trying to compare directly surface mass variations, and identify regions that diverge outside of uncertainties, we present the error in GRACE-derived surface mass as one that grows linearly with time.

Uncertainty in ISSM
Regional climate model SMB products are considered to be more mature than ice dynamic models on decadal time scales, which define the length of our comparison study.Therefore, we define the ISSM errors to be the 1-sigma spread between the three different ISSM runs forced with RACMO, MAR, and BOX SMB timeseries.As such, the errors capture uncertainty in the SMB models solely.This approach allows us to explicitly identify regions for which GRACE-JPL and ISSM diverge outside of formal errors, and attribute these differences to likely errors in ISSM (which could be due to spinup, lack of a basal hydrology model, unmodeled ocean-ice interactions, errors in bedrock, errors in the basal drag coefficient, or resolution limitations with the mesh size).Note that since all SMB products are based on output from an RCM that is forced at the boundaries with the ERA-I reanalysis, there could be common mode errors in the SMB products that are not considered here.
For model comparison against the GRACE-JPL timeseries on a mascon-by-mascon basis, we attempt to capture the Greenland Ice Sheet's surface features, including topography and surface velocities, through the initialization, relaxation, and spinup procedure described above.To aid in evaluation of the spinup procedure with respect to model velocities, we plot the change in ice velocity from the beginning of 2003 to the end of 2012 and the simulation mean difference in ice velocity between InSAR observed surface velocities and ISSM Greenland observed velocities at the end of 2012 (Rignot and Mouginot, 2012) in Fig. 2B and Fig. 2C, respectively.It is clear that in the Greenland interior, model velocities are generally too low, especially in large outlet glaciers, including Jakobshavn Isbrae, Petermann Glacier, and the outlets of the Northeast Greenland Ice Stream (Fig. 2C).The smaller fast-flowing outlet glaciers on the Northwest coast also have velocities lower than observed.In the Southeast, marginal ice speeds are greater than observed, but continuously slow down over the ten years of simulation (Fig. 2B).In general, throughout this period, the modeled ice velocities along the margins are slowing, while interior ice velocities are accelerating.Additionally, we plot the simulation mean change in ice thickness from the beginning of 2003 to the end of 2012 and the difference in ice thickness between ISSM Greenland at the end of 2012 and the mass conserving thickness estimated by Morlighem et al. (2014b) in Fig. 2E and Fig. 2F, respectively.The mean 10-year cumulative SMB anomaly forcing is plotted for comparison in Fig. 2D.To reach steady-state, we find that the model thins in the Northern interior and thickens in the Southern interior (Fig. 2F).In general, the ice sheet thickens along the margins, especially along the Southeast coast (Fig. 2F), where the thinning is accompanied by increases in velocity (Fig. 2C).The 10-year study period is, overall, dominated by marginal thinning (Fig. 2E), and this thinning is driven strongly by a decrease in SMB during the GRACE period (Fig. 2D).
For an overall comparison between GRACE-JPL and the ISSM Greenland simulation results, we plot the total cumulative Greenland mass over the 10-year study period in Fig. 4A, and the total interior cumulative mass in Fig. 4B.These plots include the mean total SMB anomaly RCM-derived forcing over the ice sheet (SMB-GrIS), the mean simulation results of the ISSM Greenland historical runs over the ice sheet (ISSM-GrIS), and the mean ISSM-GrIS simulation results over the ice sheet plus the calculated mass change over the Greenland periphery (ISSM-GrIS+P).All timeseries are plotted as cumulative mass through time and are fixed to begin at zero at the start of 2003.Plots showing the timeseries for the three individual runs plus periphery are also provided for reference (Fig. S2).Linear 10-year trends for each individual RCM and means are summarized in Table 1.It is immediately clear that the model estimates of the trend in cumulative total Greenland mass are less negative than captured by GRACE-JPL, and account for only 50% of the total GRACE signal.Inclusion of the periphery accounts for an additional 13% of the total MB trend captured by GRACE-JPL.The seasonal variability, on the other hand, appears to be well captured by the ISSM-GrIS+P estimates.Note that the reported GRACE-JPL trend of -291 Gt/yr includes mascon 33 which includes a portion of Ellesmere Island, so it is not a true estimate of mass change solely over Greenland.
In Fig. 5, we plot the difference in trend spatially, per mascon.The trend is obtained by simultaneously fitting a trend, annual, and semiannual to each timeseries of mass.We find that the majority of the discrepancy occurs in specific regions: mascon 167 (Kangerdlugssuaq), mascon 266 (Southeast glaciers), and mascons 86/87 (Northwest glaciers).The Southwest also contributes, The Cryosphere Discuss., doi:10.5194/tc-2015-224,2016 Manuscript under review for journal The Cryosphere Published: 18 January 2016 c Author(s) 2016.CC-BY 3.0 License.though to a lesser extent, i.e. mascon 165 (Jakobshavn Isbrae) where the models result in a smaller negative trend and mascon 212 where the models result in a larger negative trend, each by about 10 Gt/yr.
In Fig. 6 we plot the mean annual amplitudes for each mascon.The annual amplitude is calculated by first removing a 13-month running mean from each mascon timeseries, and then simultaneously fitting an annual and semiannual to each timeseries of mass.Overall, the annual amplitudes are well captured by the model results, suggesting that the spatial and temporal variability of SMB are well represented by the three forcing products.This is especially the case for the mascons that disagree the most in trend (i.e.86, 87, 167, and 214), suggesting that errors in SMB are not dominantly responsible for the differences between modeled mass trends and GRACE-JPL.More likely, these differences are caused by transient dynamics not captured by ISSM, including the effects of hydrology and ice-ocean interaction.
In the interior, we find that the total signal for GRACE-JPL is positive throughout the study period (Fig. 4B), while the models suggest that mass increases until 2006 and remains neutral for the second half of the simulation.The disagreement between GRACE-JPL and ISSM-GrIS in the ice sheet interior is driven almost exclusively by mascon 58 (Fig. 5C).The other four interior mascons agree well, suggesting that most of the interior of the ice sheet is indeed close to being in a steady state.
In mascon 58, however, GRACE-JPL captures thickening well outside the spread of the model simulations.Overall, we find that the interior signals are small and that the interannual variability is within the inherent error of GRACE-JPL (Fig. 4B).This limitation makes it difficult for studies to use GRACE to assess the variability and trends in SMB on decadal timescales in the Greenland interior.Still, it is important to note that in the Greenland interior, the three RCM-derived SMB products agree on a monthly timescale.Also notable in Fig. 4A is the difference in trend between ISSM-GrIS and the SMB-GrIS.The ISSM-GrIS trend is less steep than that of the SMB forcing anomalies, indicating that model dynamics play a role in decreasing ice discharge (Table 1).In Fig 7, we plot the dynamic effects of the model on trend and amplitude, for each mascon.Spatially, this effect on trend is largest in the south, especially in the Southeast, and also in the Northeast, where the dynamics of the model work to reduce mass loss through time.We also find that along the margins, the model mutes the annual cycle, suggesting that the effect is related to melt.Dynamically, it is the interior of the ice sheet that increases in velocity over the study period (Fig. 2C), yet at same time the margins thin (Fig. 2E), flatten (Fig. S3B), and slow down (Fig. 2B).This results in dynamic thickening (ice thinning at a lesser rate than that predicted by the SMB forcing), especially in the Southeast and in the large outlet glaciers in the north (Fig. S3A).We find that the thinning of the margins works to decrease discharge and that the flattening works to decrease the driving stress along the margins (Fig. S3C).Overall, this behavior leads to the model conserving ice, and we find that this dynamic behavior works to reduce the spread between the three simulations (Fig. 4A and Table 1).The interior mascons are less affected, though there is evidence that minor thinning and resultant dynamic thickening occur in the south (Fig. 2E), i.e. mascons 166 and 124 (Fig. 1 and Fig. S3A).Overall, however, we find a close match between the ISSM-GrIS simulation and the SMB-GrIS in the interior (Fig. 4B and 7).
Another significant component of the mass signal is the contribution from peripheral glaciers.In Fig. 8, we plot the spatial contribution of Greenland's periphery on trend and amplitude.We find that inclusion of the periphery contributes negatively to the trend (Fig. 8A), particularly in the Southwest, in mascon 33 (i.e.Ellesmere Island) and in mascon 167, but it contributes 10 The Cryosphere Discuss., doi:10.5194/tc-2015Discuss., doi:10.5194/tc- -224, 2016 Manuscript under review for journal The Cryosphere Published: 18 January 2016 c Author(s) 2016.CC-BY 3.0 License.positively to the amplitude of the annual signal (Fig. 8B).For the overall ice sheet (Fig. 4A), it is clear that during the summer months the periphery glaciers (included in the red line, ISSM-GrIS+P) are responsible for an increase in mass loss and during the winter months, an increase in mass gain.The winter mass gain in driven largely by seasonal snow load on tundra, while summer melt of peripheral land ice dominates the signal and contributes to the overall negative trend.
We estimate that the peripheral glaciers are responsible for a total trend of −37±25 Gt/yr (Table 1), which agrees well with Gardner et al. (2013).However, it is important to note that this estimate is associated with large uncertainty (Fig. 4A difference between red and blue shading).In fact, while the inclusion of the periphery glaciers allows us to account for a part of the discrepancy between GRACE-JPL and ISSM-GrIS, doing so also increases the estimated uncertainty of the model results (Fig. 4A).Comparison between RACMO, MAR, and BOX simulation trends reveals that indeed, there is a substantial spread between periphery estimates (Fig. S4 and Table 1).The discrepancy between the models is particularly large in the southern mascons.In the Southeast, where slopes and gradients in SMB are large along the ice sheet margins, there is inconsistency between the RCMs, even in the sign of the trend.Analysis of the periphery in RACMO, results in slightly positive trends in the south, particularly in the Southeast.Likely, this is partially due to the lower resolution of the RACMO product (which is not downscaled to a higher resolution in post-processing, like the other two RCM products).The lower resolution leads to difficulty in resolving the ice margin near the coast.Comparison of mass trends in GRACE-JPL with annual altimetry estimates (which do not include periphery) from 2003-2009, offers an observational estimate of peripheral mass trend during the first portion of the GRACE record (Fig. S1).Results suggest that peripheral estimates from MAR have the best overall agreement with the observed in the south, while BOX has the best overall agreement in the north.
Timeseries and mean seasonal cycle of mass change estimated by GRACE-JPL and by ISSM-GrIS+P are compared in individual non-interior mascons in Figs. 9 -12.We also include estimates of error which, for the model estimates, are defined by the spread in the three different simulations (as discussed in Sect.5.2).The mascons are organized geographically, by ice sheet region (i.e.Northeast, Southeast, Southwest, and Northwest).In some cases, mascons contain small fractions of the ice sheet margin.In these cases, mascons are combined with a neighboring mascon.We include timeseries for all individual mascons in the Supplement, Figs.S5 -S9.
In the majority of the mascons, ISSM-GrIS+P captures the seasonal cycle observed by GRACE-JPL.The largest discrepancies between ISSM-GrIS+P and GRACE-JPL occur during the summer.During this time, we also find that there is the largest disagreement between the ISSM-GrIS+P runs.These results suggest, in agreement with Velicogna et al. (2014), that errors in estimates of runoff within the SMB products are largely responsible for driving diverging error bars through time (most notably in the Southwest, Fig. 9).Consistent errors in RCM estimates of runoff may also be partially responsible for the trend differences between ISSM-GrIS+P and GRACE-JPL, particularly in the mascons that agree well during the winter months but differ during the summer.Overall, it is difficult to pinpoint a consistent bias that is associated with a particular region or SMB model.Specifically, when comparing these products at the spatial scale of a mascon (∼300 km), discrimination of the sources of uncertainty becomes more complicated.
In the Northeast, for instance, GRACE-JPL and ISSM-GrIS+P agree well in overall trend for mascons 59, and 89+90, while the annual cycle in 89+90 is consistent with GRACE-JPL estimates (Fig. 10).However, the annual cycle for mascon 59 in The Cryosphere Discuss., doi:10.5194/tc-2015-224,2016 Manuscript under review for journal The Cryosphere Published: 18 January 2016 c Author(s) 2016.CC-BY 3.0 License.ISSM-GrIS+P is more exaggerated than GRACE-JPL, with greater accumulation during the winter months and greater mass loss during the summer months.For mascon 35 and 125+126, we find that GRACE-JPL has a more negative trend than ISSM-GrIS+P.In mascon 35, this is due to an ISSM-GrIS+P underestimation of mass loss in the summer relative to GRACE-JPL, while in 125+126, it is due to an ISSM-GrIS+P overestimation of mass gain during the winter relative to GRACE-JPL (Fig. 10).In the Southwest, such a discrepancy in trend for mascon 165 is due to a combination of the two (Fig. 9), with the spread between ISSM-GrIS+P and GRACE-JPL increasing non-linearly through time (Fig. S6).Mascons 212 and 265 have more negative trends in ISSM-GrIS+P than GRACE-JPL, but agree very well in the seasonal cycle.This area is well covered by observations (including the K-transect), and often RCMs are evaluated in this area.Here, we find that the differences are due to a higher estimate of runoff during the summer months (Fig. 9), predominantly causing a divergence between ISSM-GrIS+P and GRACE-JPL from 2003-2009 (Fig. S6).Finally, for mascon 324, we find that the ISSM-GrIS+P estimates have a large spread in trend, largely due to uncertainty in periphery estimates (Fig. S4).The periphery estimate for mascon 324 also contributes to an exaggeration of the annual amplitude for this mascon (Fig. 8), resulting in a perceived overestimate of mass gain in the winter and mass loss in the summer (Fig. 9).While errors in SMB forcing in these regions likely play the dominant role in differences between GRACE-JPL and ISSM-GrIS+P, there is evidence that missing model processes may also play a factor.Specifically, during the winter, GRACE-JPL captures month-to-month variability that is beyond the spread of the three ISSM-GrIS+P runs (e.g.mascons 165 and 324).In comparison, the ISSM-GrIS+P results are smooth and do not exhibit the same type of variability.
Such variability is also present in the Southeast (Fig. 11) and Northwest (Fig. 12) sectors, where the majority of mascons show a significant discrepancy in trend between GRACE-JPL and ISSM-GrIS+P.Mascons 213 and 33+34 are clear exceptions, and match well in both trend and annual amplitude.However, note that for mascon 213, ISSM-GrIS+P and GRACE-JPL agree prior to 2010 and then continuously diverge through 2012 (Fig. S8).These results agree with observations of velocity for the region, in particular the acceleration of Ikativaq region in 2009 and of Helheim Glacier in 2010 (Joughin et al., 2008;Moon et al., 2012;Khan et al., 2014).The rest of the mascons in the Southeast and the Northwest, where dynamics are believed to play a large role in recent mass changes, have GRACE-JPL signals that show large negative trends.These negative trends are consistently underestimated by the ISSM-GrIS+P runs, even in the locations where the seasonal signal appears to be well captured (e.g.Fig. 11, Mascon 214 and Fig. 12, Mascon 56).In the Northwest, in particular, it is clear that in this region, it is not just the summer season that is responsible for the difference between GRACE-JPL and ISSM-GrIS+P.In Mascons 56, 86+87, and 123, we find a distinct difference between GRACE-JPL and ISSM-GrIS+P during the entire year (Fig. 12).
In these cases during the winter, GRACE-JPL indicates that mascon regions continue to lose mass, while the SMB forcing (represented by the ISSM-GrIS+P runs) remains positive.In fact, for mascons 56 and 86+87, the summer melt appears to be well represented by the models.For Mascon 123, we find that summer runoff is overestimated; yet the negative trend is simultaneously underestimated due to an overestimate of mass gain during the rest of the year.In order to evaluate model estimates of Greenland mass balance, we present regional comparisons between model results and the GRACE-JPL solution.For these comparisons, we focus specifically on differences in the amplitudes of the mean seasonal cycle and on differences in trends over a 10-year GRACE period (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012).Model estimates are presented as the mean of the results from three independent ice sheet flow model forward simulations, where each simulation is forced with different RCMand observationally-based products of monthly SMB.Estimates of snow load over peripheral tundra and estimates of mass balance in peripheral glaciers are also included for direct comparison with the GRACE-JPL solution on a mascon-by-mascon basis.
Results suggest that surface mass balance is the primary driver of simulated mass balance (Fig. 2D,E) (Enderlin et al., 2014).We find that about 63% of mass loss is from SMB and about 37% is from ice dynamics, in agreement with previously reported GRACE-based estimates of Greenland mass balance (Velicogna et al., 2014).With respect to seasonal variability, variations in monthly SMB forcing dominate over the modeled changes in ice sheet dynamics over the short study period (Fig. 6 and Fig. 7).Overall, we find that the model actually dampens the annual SMB signal, particularly during the summer when marginal thinning drives a reduction of ice discharge (Goelzer et al., 2013;Enderlin et al., 2014).It is important to note that this negative feedback between marginal ice thickness and discharge ultimately serves to reduce the spread between the initial SMB products (Fig. 4 and Table 1), which suggests that the model state -common to all three simulations -may play a role in dictating the results of the three different simulations.As a consequence of the dampening of mass loss during the summer, we find that the model serves to decrease the magnitude of the negative mass trend along the margins.Simultaneously, in response to marginal thinning, the model simulates an increase in interior velocities throughout the simulation (Fig. 2B).While this simulated acceleration in ice dynamics is not a trivial component to ISSM-GrIS estimates of mass balance, we find that it is relatively small compared to the contribution from the SMB forcing.Indeed, a total ice sheet comparison against GRACE-JPL (Fig. 4 and Table 1) reveals that the dynamic changes simulated by the ice flow model represent a minor source of uncertainty in the modeled trends.We find, therefore, that other sources of model uncertainty (including dynamical processes not included in the ice sheet model) are responsible for substantial differences between the ISSM-GrIS and GRACE-JPL mass estimates on a total Greenland-scale, accounting for 108 ± 35 Gt/yr of total mass loss (Table 1).These results are consistent with recent studies (e.g.Moon et al. (2014)) which report observed seasonal accelerations in local ice flow of magnitudes far larger (by a factor of 10) than the changes in ice velocity modeled by ISSM-GrIS over the ten-year simulation period (Fig. 2B).
Mascon-scale comparisons between the model mass estimates and the GRACE-JPL solution offer a more detailed view of the spatial distribution in model uncertainty.Regional comparisons suggest that trend differences are concentrated in the Northwest and the Southeast sectors of the ice sheet, and that in the rest of the ice sheet, SMB dominates local trends.Results also suggest that the spatial and temporal variability of SMB are well represented by the three forcing products, as annual amplitudes are well captured by the model results over the majority of the ice sheet (Fig. 6).This is the case even in the Northwest, where the model and GRACE-JPL mass estimates differ significantly in trend (Fig. 12).Additionally, we find that inclusion of the periphery area proves to be an important component of the seasonal cycle captured by GRACE-JPL with 13 The Cryosphere Discuss., doi:10.5194/tc-2015Discuss., doi:10.5194/tc- -224, 2016 Manuscript under review for journal The Cryosphere Published: 18 January 2016 c Author(s) 2016.CC-BY 3.0 License.marginal mascons.Bare rock, peripheral glaciers and ice caps contribute to an overall larger negative trend as well as larger annual amplitudes (Fig. 4 and Fig. 8), resulting in better overall agreement between GRACE-JPL and model-estimated mass balance.However, results suggest that the periphery is also a source of large uncertainty (Table 1).Particularly in periphery areas of steep, complex terrain (i.e. the Southeast) -where validation is difficult, and RCM spatial resolution of 5-11 km is inadequate to capture gradients in SMB -we find that the SMB forcing products do not have good agreement (Fig. S4 and Fig. 11).Future improvements in RCM resolution, snowpack models for tundra regions, and simulation of climate over the peripheral glaciers and ice caps will be essential for future comparisons and validation against seasonal-scale mascon-style GRACE products.
Overall, results from this comparison suggest that largest discrepancies in mass trend between the model and GRACE-JPL are in the Northwest sector of Greenland.Here, such discrepancies are likely due to consistent ocean forcing, hydrology-driven events, errors in modeling the bedrock, or errors in the ice model spinup.Modeled velocities in these areas are lower than measured by InSAR, so differences in model state exist at relaxation, and then persist throughout the simulation (e.g.errors due to insufficient bedrock, lack of physics, or not enough melt in relaxation SMB).Mean annual plots of GRACE-measured mass change (e.g.mascons 86/87, 123, and 56) reveal that the Northwest continues to lose mass throughout the entire year, even during winter months (Fig. 12).Comparison between the mean annual cycles and GRACE-JPL indicate that it is ice dynamics -not captured by the model -that dominate the mass trends here.Indeed, it is in this region where we find SMB plays less of a role in determining mass balance, particularly in areas where modeled mass and GRACE-JPL disagree outside of estimated uncertainty.Since SMB is positive during the winter, and the SMB products have strong agreement in this region during the fall, winter, and spring, dynamics are most likely responsible for the strong discrepancies between GRACE-JPL and ISSM-GrIS+P during non-summer months.During these months, GRACE-JPL exhibits mass loss inconsistent with SMB, which suggests that the total mass in the Northwest is strongly out of balance.This finding is supported by the model behavior during relaxation to steady-state SMB.During relaxation, we find that many glaciers in the Northwest slow down in order to be in balance with the SMB forcing (Fig. 2C).These results suggest that our assumption of historical steady-state is likely invalid for the Northwest region of Greenland.
In the Southeast, it is more difficult to pinpoint a particular factor that drives the differences between GRACE-JPL and the model estimates of mass change.Mean seasonal plots (Fig. 11) of the mascons in this area reveal that the GRACE-JPL signal exhibits larger seasonal variations than estimated by the models.This suggests that discrepancies may be controlled by errors in modeled SMB, including errors in mass contribution from the periphery (i.e.trend from glaciers and annual signal from load on bare rock and tundra).The topography in the Southeast is steep, mountainous, and generally plagued by the largest uncertainties in modeled snowfall estimates, yet we find that the SMB products represented here tend to agree well with GRACE-JPL during the majority of the year.The largest discrepancies with GRACE-JPL occur during the summer months, which also happens to be when the discrepancies between the SMB forcing products are the largest.Such results suggest that SMB errors may contribute to model uncertainty in the Southeast, and RCM estimates of runoff for both the ice sheet and periphery glaciers may not be accurately captured.This may particularly be the case in mascon 266 (Fig. 11), where the steep terrain creates a very narrow ablation zone that is difficult to capture at the resolution of the SMB forcing.In contrast 14 The Cryosphere Discuss., doi:10.5194/tc-2015Discuss., doi:10.5194/tc- -224, 2016 Manuscript under review for journal The Cryosphere Published: 18 January 2016 c Author(s) 2016.CC-BY 3.0 License.with the other Southeast mascons, mascon 266 exhibits poor agreement in annual amplitude.In this mascon, we find that a consistent divergence between GRACE-JPL and the model estimates of mass loss occurs almost exclusively during the summer months (Fig. 11), suggesting that a seasonal phenomenon may be responsible for mass loss in this region.It is important to note that observational evidence indicates that glaciers in this region are characterized by a slowdown during the summer, not acceleration.In addition, there is little evidence of an overall background dynamic acceleration after 2005 (Moon et al., 2014).
Therefore, we find that the most likely contributor to seasonal discrepancies in these Southeast glaciers is common runoff error in all of the SMB forcing products utilized in this study.
There is additional evidence that, in agreement with recent publications, e.g.Csatho et al. (2014); Khan et al. (2014); Moon et al. (2012Moon et al. ( , 2014)); Velicogna et al. (2014), the transient dynamics (not captured by ISSM) play a role in mass loss in this area.For instance, mass balance within the Helheim and Ikativaq region (Mascon 213) is well captured by the model overall, but it is clear that GRACE-JPL and the models diverge between 2005 and 2006 and then again in 2010 (Fig. S8).This trend in dynamics is consistent with observations of high velocities in Helheim in 2005, followed by a slowdown, and then acceleration in Ikativaq in 2009 and Helheim in 2010 (Joughin et al., 2008;Moon et al., 2012;Khan et al., 2014).Similarly, a well-documented shift in dynamics is captured by GRACE-JPL at Kangerdlugssuaq Glacier (mascon 167) in 2005 (Fig. S8).Observational evidence suggests that such changes in sensitive tidewater glaciers are strongly coupled to calving events and the position of the glacier terminus, especially during periods of rapid advancement during the spring and early summer (Joughin et al., 2008).This is consistent with the behavior of the mean GRACE-JPL seasonal signal in most of the Southeast mascons (i.e.167/168, 213, and 214), which appear to have a much noisier signal during the spring than is simulated by the models (Fig. 11), including single months of high mass loss.These results suggest that the GRACE-JPL solution is capable of capturing monthly-scale dynamic changes in large outlet glaciers, and therefore it is possible to quantify dynamic mass loss by removing the ISSM-GrIS+P from the GRACE-JPL signal.However, it is clear that with regards to the seasonal cycle, where model results fall within the GRACE-JPL error bars, we cannot confidently distinguish between errors in SMB and signals of ice dynamics.This is the case in many of the mascons, particularly in the Southeast, with the exception of mascon 266 where (as discussed above) we can confidently conclude that SMB is a significant contributor to the divergence between GRACE-JPL and ISSM-GrIS+P.Overall, in the Southeast, discrepancies are likely caused by a combination of errors including lack of ocean forcing, poor bedrock, inadequate mesh representation of the smaller and steeper glaciers in ISSM, as well as uncertainty in SMB forcing due to the steep terrain and narrow ablation zone.
In the SMB-dominated areas (i.e. the Northeast and Southwest), however, we find overall good agreement between the models and GRACE-JPL in both amplitude and trend (Figs. 5,6,9,and 10).In the Southwest, this result is not surprising, as most of the glaciers are land-terminating and the position of the glacier termini are not affected by ice-ocean interaction (Khan et al., 2015).In the Northeast, the Northeast Greenland Ice Stream (mascon 59) is well captured in trend, though we find that the annual amplitude of the GRACE-JPL signal is highly muted, particularly during the summer.Such a discrepancy could be caused by common mode errors in the SMB forcing, but the match in trend suggests that unmodeled hydrological processes may be responsible for this discrepancy (e.g.storage and delayed release of runoff) (Willis et al., 2015).Reconciling The Cryosphere Discuss., doi:10.5194/tc-2015Discuss., doi:10.5194/tc- -224, 2016 Manuscript under review for journal The Cryosphere Published: 18 January 2016 c Author(s) 2016.CC-BY 3.0 License.these results with observations of ice elevation (i.e.altimetry measurements), when available at a monthly to seasonal temporal resolution, could shed light on the key processes responsible for continued mass loss in this area.
In the Southwest region, our results also capture signals that may be related to ice dynamics.Specifically, we find that in this region, the relationship between SMB and mass change is not consistent through time.For instance, the model and GRACE-JPL diverge for mascons 212 and 265 between 2005 and 2010, but then agree well for the remainder of the simulation (Fig. S6).In contrast, mascon 165 (i.e.Jakobshavn Isbrae) has good agreement between the model and GRACE-JPL at the beginning of the simulation, but begins to diverge around 2008 (Fig. S6).In fact, the mass loss in Jakobshavn Isbrae appears to be accelerating through time.These results are consistent with published observations of minor speedups in velocity for Jakobshavn Isbrae beginning in 2008, as well as a general velocity decrease in the Southwest between 2005 and 2010 (Moon et al., 2012) (corresponding to mascons 212 and 265).We find that the disagreement between the model and GRACE-JPL at the beginning of our study period is responsible for the discrepancies in trend for these regions.The SMB in this area is well validated (i.e.K-transect), and annual amplitudes agree well.Therefore, it is likely that transient dynamics, not captured by ISSM, contribute to the divergence between the models and GRACE-JPL within mascon 265.The monthly-scale variations in regional mass loss evidenced in the GRACE-JPL seasonal cycle (Fig. 9) is most likely driven by dynamics of the few, but active, marineterminating glaciers in the region, including the effects of hydrology and ice-ocean interaction (calving events/position of glacier terminus) (Holland et al., 2008).It is clear that over the course of just a decade, consistent dynamic responses to these types of climate-driven forcing can ultimately perturb regional mass trends, even in the regions where mass loss appears to be dominated by SMB.
In the high altitude interior of the ice sheet, where the mass balance is dominated by accumulation, comparisons show good agreement in trend (Fig. 4b and Fig. 5), with GRACE-JPL showing a slight discrepancy of about 20 Gt throughout the ten year period.Both the models and GRACE-JPL suggest that Greenland was gaining mass within its interior until 2006, and that this was followed by a relatively neutral period through 2012.The variability in the GRACE-JPL signal results in periodic disagreement between GRACE-JPL and models, outside of the assessed uncertainty (Fig. 4b).One possible explanation is that the RCMs -which are commonly forced by ERA-I and agree well in the interior -are not capturing stochastic accumulation events that occur during roughly the same time every year.However, evidence suggests that this is not the case; in particular, analysis reveals that the MAR3.5.2 product forced with NCEP/NCAR Reanalysis 1 (Kalnay et al., 1996) exhibits similar temporal variability (not shown) and annual amplitude (Fig. S10A) to the ERA-I SMB products considered here.A more likely explanation is that the discrepancy is caused by a combination of noise in a locally small signal in GRACE-JPL coupled with modest leakage errors from neighboring coastal mascons with much stronger signals.If so, these results indicate that mass signals in the high altitude interior of Greenland are sufficiently small enough to push the limits of GRACE utility for model evaluation, both temporal and spatially.The use of GRACE in this area is additionally complicated by a GIA correction that is significant in comparison to the magnitude of the GRACE signal.We expect that advances in GRACE mascon processing (including increases in spatial resolution), GIA modeling, and progressions in RCM estimates in SMB (including improved validation in the interior using satellite data and data from a growing network of in-situ stations, and the diversification of RCM forcing products) may, in the near future, help clarify these discrepancies.
The Cryosphere Discuss., doi:10.5194/tc-2015Discuss., doi:10.5194/tc- -224, 2016 Manuscript under review for journal The Cryosphere Published: 18 January 2016 c Author(s) 2016.CC-BY 3.0 License.Despite the complexities in the variability of mass balance within the Greenland interior, it is possible to focus on comparison between the interior mass trends (which agree within the products' respective bounds of uncertainty) in order to evaluate ISSM-GrIS spin-up procedures.Specifically, upon relaxing the model using a historical period of mass balance neutrality, the spin-up procedure adopted for this study assumes that the Greenland Ice Sheet was in near steady-state during the recent past.More specifically, we relax the model to a virtual steady-state, using a mean climate forcing from the 1979-1988 period -a period in which the rate of ice sheet mass loss was negligible compared to the mass loss captured by GRACE during the last decade.
We adopt this procedure in order to remove spurious transients from the model that may manifest due to mismatched input including: bedrock and ice surface elevation maps, surface ice velocities, and SMB.After relaxation, ISSM-GrIS discharge is equal to the mean SMB forcing, and resultant perturbations to ice thickness or velocity that occur in the forward model will be solely in response to anomalies in the transient SMB forcing.Since trends match well between the models and ISSM-GrIS, particularly in the beginning of the timeseries (Fig. 4b), we find that steady-state is indeed a reasonable spin-up assumption.
We acknowledge that this assumption does result in differences between modeled and observed ice thicknesses (Fig. 2), and in turn may cause the model to exhibit second-order deviations in ice velocities.This is especially the case considering thateven though the SMB products have been validated against observations -these observations are sparse, and all SMB products are associated with systematic errors that may impact spinup and propagate into the simulation.However, it is clear that in the absence of a long-term model spin-up (on the order of thousands of years) the assumption of steady-state is adequate for short term (annual-to-decadal scale) simulations.In fact, we find that our results are fairly insensitive to the SMB product chosen as forcing for relaxation (Figs.S10B and S11B).These results suggests that -on the temporal and spatial scales considered here -background trends in mass balance (that occur in response to paleo-climate forcing) may play a minor role in dictating present-day evolution of Greenland mass balance when compared to SMB anomalies and seasonal-to-annual scale variations in ice velocities.Indeed, we find that the regions that are currently in the strongest imbalance are also affected by seasonal-to annual-scale ice dynamics that appear to be associated with transient processes not captured by ISSM-GrIS.
Though it is clear that SMB accounts for the majority of Greenland mass balance, these transient processes are strongly linked to changes in climate forcing.In addition, they may be responsible for positive feedbacks on ice velocity that drive monthly-scale regional mass variations large enough to be captured by the GRACE-JPL signal.The model representation of these processes, most importantly the evolution of the hydrological system and ice-ocean interactions, must be implemented in order to improve ice sheet model skill for decadal-scale simulations.Such model improvements are difficult, because these processes are the not well understood and are associated with large uncertainties.However, in order for the glaciological community to take full advantage of the array of new observational products that are available (and will be made available) for model evaluation in the near future, it is essential that simulations include hydrological effects on the system.The future success of Greenland ice flow model simulations, including hindcasts as well as future projections of ice sheet mass balance, will require high confidence in SMB forcing and physically-based representations of transient processes.Simulation of these transient processes at a monthly-to-seasonal temporal resolution is a key element in model evaluation using observational datasets like GRACE, since results indicate that consistent intra-annual variations combine to ultimately dictate regional trends in mass balance.Multiple station records now contribute to the near surface air temperature for each given year, month and grid cell in the domain, while in Box (2013) data from the single highest correlating station yielded the reconstructed value.The estimation of values is made for a domain that includes land, sea, and ice, which is an expansion to the Box (2013) product that estimates T over ice only.A physically-based meltwater retention scheme (Pfeffer et al., 1990(Pfeffer et al., , 1991) ) replaces the simpler approach used by Box (2013).The RACMO2.3 output have a higher native resolution of 11 km as compared to the 24 km Polar MM5 output used by Box (2013) for air temperatures.In addition, the revised SMB product ends two years later, in year 2012.The annual accumulation rates from ice cores are dispersed into a monthly temporal resolution by weighting the monthly fraction of the annual total for each grid cell in the domain evaluated using 1960-2012 RACMO2.3output.

Methods for defining peripheral SMB
For this study, we define peripheral ice as isolated permanent ice that exists outside of the ISSM Greenland domain (see Fig. 1).A land-ice-ocean mask accompanies all the SMB products considered here.The masks differ for each product; therefore we must interpret them independently, with reference to the specific mask defined for a particular product.For BOX, valid peripheral ice SMB fall within mask values ranging from greater than 0 and less than or equal to 1.For MAR, valid SMB values fall within the mask values ranging from greater than 1 and less than or equal to 2. For RACMO, valid SMB values fall within the mask 'IceSheetMask' values ranging from greater than 0.5 and less than or equal to 1.
For determining snow load outside of areas of permanent ice, we define peripheral tundra as area masked as land on Greenland or Ellesmere Island, within our Greenland mascons (see Fig. 1).Once again, because the product masks differ, we must consider each mask independently.For BOX, valid peripheral tundra falls within mask values greater than 1 and less than 2.
For MAR, valid peripheral tundra falls within mask values greater than 0 and less than 2. For RACMO, valid peripheral tundra falls within 'IceSheetMask' values greater than or equal to 0 and less than 1 and 'LandSeaMask' values greater than 0 and less than or equal to 1. On all grid points that contain only fractional areas of tundra, the snow load is scaled to the percentage of   (Rignot and Mouginot, 2012), peripheral tundra (beige), and peripheral permanent ice (green) (Gardner et al., 2013).Periphery -37±25 -14 -63 -32 Table 1.In the right three columns are the cumulative mass trends (Gt/yr) for individual RCM SMB products (SMB-GrIS), ISSM forced with each individual RCM SMB product (ISSM-GrIS), and ISSM-GrIS plus cumulative mass estimates for the ice sheet periphery (ISSM-GrIS+P).Presented in the left column are the mean cumulative mass trend and trend uncertainty calculated for the timeseries presented in Fig. 4A.Also reported for each of these columns are the trends for the total Ice Dynamics (difference between SMB-GrIS and ISSM-GrIS) and the total Periphery (the total SMB calculated over peripheral permanent ice defined in Fig. 1).Along with the trends, 1-sigma uncertainties (Sect.5) are displayed. 35 The Cryosphere Discuss., doi:10.5194/tc-2015Discuss., doi:10.5194/tc- -224, 2016 Manuscript under review for journal The Cryosphere Published: 18 January 2016 c Author(s) 2016.CC-BY 3.0 License.
The CryosphereDiscuss., doi:10.5194/tc-2015Discuss., doi:10.5194/tc--224, 2016     Manuscript under review for journal The Cryosphere Published: 18 January 2016 c Author(s) 2016.CC-BY 3.0 License.The CryosphereDiscuss., doi:10.5194/tc-2015Discuss., doi:10.5194/tc--224, 2016     Manuscript under review for journal The Cryosphere Published: 18 January 2016 c Author(s) 2016.CC-BY 3.0 License.temperature (T) and SMB components is based on a 53-year overlap period.Note that the overlap period for the calibration of snow accumulation rate is shorter, since ice core data availability drops after 1999.Calibration is made using linear regression coefficients for 5 km grid cells that match the average of the reconstruction to RACMO2.3.The RACMO2.3 output are resampled and reprojected from the native 0.1 deg (∼10 km) grid to a 5 km grid better resolving areas where sharp gradients occur, especially near the ice margin where mass fluxes are largest.To create the BOX SMB forcing used here, several refinements are made to the Box (2013) T and SMB reconstruction. Figures

Figure 2 .Figure 3 .
Figure 2. (A) The ISSM Greenland mesh; (B) change in modeled surface velocities during the 10-year ISSM simulation period (2003-2012); (C) the departure of modeled surface velocities (December 2012) from observations (Rignot and Mouginot, 2012); (D) the total SMB anomaly forced upon the ISSM Greenland model during the simulation period; (E) change in modeled ice thickness during the simulation period; and (F) the departure of modeled ice thicknesses (December 2012) from observationally based data Morlighem et al. (2014b).Model output (B)-(F) is presented as the mean of three different ISSM simulation runs (forced with BOX, MAR, and RACMO).

Figure 5 .
Figure 5. Spatial representation of trend in surface mass over Greenland from (A) GRACE-JPL, (B) ISSM-GrIS+P, and (C) the difference between GRACE-JPL and ISSM-GrIS+P.

Figure 6 .
Figure 6.Spatial representation of annual amplitude of surface mass over Greenland from (A) GRACE-JPL , (B) ISSM-GrIS+P, and (C) the difference between GRACE-JPL and ISSM-GrIS+P.

Figure 7 .
Figure 7. Difference in (A) trend and (B) annual amplitude of (ISSM-GrIS -SMB-GrIS) showing that ISSM mutes both the trend and annual amplitude of surface mass relative to what SMB anomalies would predict.

Figure 8 .
Figure 8. Difference in (A) trend and (B) annual amplitude of (ISSM-GrIS+P -ISSM-GrIS) showing that including estimates of mass from the periphery both increases the magnitude of the annual amplitude and the negative trend of surface mass.