Hindcasting to measure ice sheet model sensitivity to initial states

Validation is a critical component of model development, yet notoriously challenging in ice sheet modeling. Here we evaluate how an ice sheet system model responds to a given forcing. We show that hindcasting, i.e. forcing a model with known or closely estimated inputs for past events to see how well the output matches observations, is a viable method of assessing model performance. By simulating the recent past of Greenland, and comparing to observations of ice thickness, ice discharge, surface speeds, mass loss and surface elevation changes for validation, we find that the short term model response is strongly influenced by the initial state. We show that the thermal and dynamical states (i.e. the distribution of internal energy and momentum) can be misrepresented despite a good agreement with some observations, stressing the importance of using multiple observations. In particular we identify rates of change of spatially dense observations as preferred validation metrics. Hindcasting enables a qualitative assessment of model performance relative to observed rates of change. It thereby reduces the number of admissible initial states more rigorously than validation efforts that do not take advantage of observed rates of change.


Introduction
Realistic projections of ice sheet response to a changing climate should be based on physical understanding of the processes involved, rather than trend extrapolation of historical observations (Arthern and Hindmarsh, 2006).Ice sheet system models are multi-physics models incorporating such physical process understanding.The main sub-systems are ice dynamics, surface and basal processes, and thermody-namics.Ideally verification and validation should be done for each model sub-system independently.Verification (i.e. the comparison of results from a numerical approximation to exact solutions of the same continuum model equations) is currently only possible for a few sub-systems (e.g.Bueler et al., 2005;Leng et al., 2013) and in a simple coupled system (Bueler et al., 2007).In engineering, validation is commonly defined as the process of comparing model results to a set of observations adequate to falsify a model (Roache, 1998).Such validation is challenging to apply in ice sheet modeling; nonetheless attempts have been made (Robison et al., 2010;Burton et al., 2012).Direct validation of substantial sub-systems such as basal hydrology, thermodynamics, and ice dynamics is difficult or impossible as most or all observations available for validation are not linked to a single process, but are the consequence of a complex interplay between sub-systems.
Here we do not seek to isolate, e.g., the model of ice dynamics, and evaluate it independently; instead we evaluate how the system responds to a given forcing.In other words we ask the question: "how successful is a state-of-the art ice sheet system model (i.e. the combination of physical models, their numerical approximations and implementations, and particular choices of boundary forcing and initial states) in reproducing observations of quantities such as ice thickness, and their temporal changes?"Even if all sub-systems could be verified and validated independently, testing of the system as a whole is indispensable.Consequently, our definition of "validation" differs from the engineering literature.
The aim of this study is to demonstrate that hindcasting, i.e. forcing a model with known or closely estimated inputs for past events and comparing model results to timedependent observations, is a viable method of assessing the performance of an ice sheet system model.First we initialize an ice sheet system model.Second the initialized model is integrated forward in time for a period where a wealth of data are available for validation.This produces a hindcast of this time period.Finally, the hindcasts are compared to available observations and the performance of the system model is evaluated.In an ideal case, when using appropriate model physics paired with realistic initial states and boundary conditions, a hindcast and observations should agree to within error bars associated with observations.Such an agreement with all available observations may not occur due to several reasons.Disagreement may arise from uncertainties in boundary conditions, unrealistic initial states, incomplete model physics, and inadequate choices in numerical methods, to name a few.Using multiple data sets can help to pin down the source of disagreement.As models evolve, the level of disagreement between observations and model simulation is expected to decrease.
Similar to short-term weather forecasts, a realistic initial state is essential for accurate simulations of the future evolution of an ice sheet (Arthern and Gudmundsson, 2010).In this study we use hindcasting to assess an ice sheet system model's sensitivity to its initial state.We show that it is possible to get good agreement with observations when too few data sets are used.One may wrongly put too much confidence in a given model simulation, not realizing that important properties are misrepresented.In particular, previous studies have often used ice volume (e.g.Stone et al., 2010;Rogozhina et al., 2011;Applegate et al., 2012) for validation of the initial state.Spatially dense observations, especially rates of change, are, however, better metrics to evaluate the quality of initial states (Vaughan and Arthern, 2007).As an example, we show that good agreement between observed and simulated total mass loss can be achieved, and only by comparing observed vs. simulated ice discharge is it possible to identify that the apparent good agreement of the former is due to the wrong reason.
Our manuscript is intended to serve as a guidance for future studies developing better practices in ice sheet system model validation.We do not identify a preferred initialization procedure.Hindcasting is not limited to our particular choices of an ice sheet model, boundary conditions, and initial states; but is transferable to combinations of these.However, conclusions regarding the performance of initialization procedures may not be transferable.Hindcasting may also be helpful to study other aspects of the ice sheet system model's behavior.It is, for example, important to assess the sensitivity to climate forcing; but this is not within the scope of our study.

Methods
Our choice of ice sheet system model comprises the Parallel Ice Sheet Model (PISM; Khroulev and the PISM Au-thors, 2012) unidirectionally coupled to the high-resolution regional climate model HIRHAM5 (Christensen et al., 2006).HIRHAM5 provides climatic mass balance (sum of the surface and internal mass balances, Cogley et al., 2011) and 2 m air temperature for the period 1989-2011 (ERA-interim; (ERA-interim;Dee et al., 2011).HIRHAM5 performs well over Greenland (e.g.Rae et al., 2012) but formal error estimates are not available.Additional information on PISM and HIRHAM5 is given in the Supplement.
We obtain three initial states by using distinct approaches to modeling the 1989 state of the Greenland ice sheet, in each case based on forward modeling and assumptions about past climate.We choose three initialization methods that have been used in the published literature (e.g.Rogozhina et al., 2011;Greve et al., 2011;Price et al., 2011).Below we describe how we obtain the three initial states.
The first initial state, "constant-climate", is obtained by running the model for 125 ka using mean values for 1989-2011 of 2 m air temperature and climatic mass balance from HIRHAM5.The 2 m air temperature serves as the boundary condition for the conservation of energy equation, and a lapse rate of 7.1 • C km −1 (Steffen and Box, 2001) is used to correct for differences between the fixed surface elevation that HIRHAM5 uses and the time evolving modeled elevation.No lapse rate correction is applied to the climatic mass balance."Constant-climate" is expected to be in equilibrium with modeled present-day climate.
The second initial state, "paleo-climate", closely follows the SeaRISE initialization procedures (Bindschadler et al., 2013).To obtain this initial state we make the following modeling choices.First, we use a positive degree-day scheme to compute the climatic mass balance from surface temperature (Fausto et al., 2009) and model-constrained precipitation (Ettema et al., 2009).The degree-day factors are the same as in Huybrechts (1998).Second, we account for paleo-climatic variations by applying a scalar anomaly term derived from the GRIP ice core oxygen isotope record (Dansgaard et al., 1993) to the temperature field (Huybrechts, 2002).Then we adjust mean annual precipitation in proportion to the mean annual air temperature change (Huybrechts, 2002).Finally, sea level forcing, which determines the land area available for glaciation, is derived from the SPECMAP marine δ 18 O record (Imbrie et al., 1984).
The third initial state, "flux-corrected", is similar to "paleo-climate", but the climatic mass balance has been modified during the last 2000 yr of the simulation to obtain a present-day geometry in close agreement with observations (Supplement).This initial state is not in equilibrium with the applied forcing fields (Price et al., 2011).To prevent a shock in climate forcing that could lead to a model drift (Gupta et al., 2012), climate forcing is applied as anomalies at the hindcasting stage (Supplement).Such flux-correction methods have been applied in coupled ocean-atmosphere general circulation models (e.g.Sausen et al., 1988).
Our modeling goals are to produce physically consistent initial states of the Greenland ice sheet in 1989 that mimic the overall dynamic state of the whole ice sheet.One cannot achieve agreement with all available observations except through introduction of many tunable parameters, which we do not do.
Also not all observed processes can be adequately represented in the model.For example, increased ice discharge since the late 1990s has been observed (Joughin et al., 2004;Howat et al., 2007) and attributed to an increase in ocean temperatures (Holland et al., 2008;Murray et al., 2010).Unfortunately, theoretical models to predict such behavior are currently not available (Holland et al., 2008).Therefore no forcing at the ocean boundary is applied.
Properties at the ice sheet bed cannot be directly observed.Ice sheet models may achieve a close fit to observed surface speeds by using data assimilation techniques that invert for a field of basal parameters, for example adjusting the value of the basal traction at each grid location (Morlighem et al., 2010;Price et al., 2011;Gillet-Chaulet et al., 2012).However, the resulting parameter fields may not be applicable for prognostic simulations, as these parameters evolve with time.Moreover, validation cannot rely on observations used in the data assimilation process (van der Veen, 1999).
In this work we adjust only three spatially uniform scalar parameters (Supplement), controlling ice dynamics and basal processes, to achieve a good match with observed surface speeds.Thus observed surface speeds are not a part of the initialization procedure and can be used to validate the model.
At the hindcast stage, all initial states are integrated forward in time from 1989 to 2011 using monthly fields of climatic mass balance and 2 m air temperature computed by HIRHAM5.Note that identical parameter choices are used to produce all initial states, as well as in the 1989-2011 integration, allowing the influence of initialization procedures to be isolated.All hindcasts are carried out on a grid with horizontal resolution of 2 km.
The hindcasts are then compared to observations listed below.Our study does not use all available observations, but the ones used allow drawing robust conclusions about hindcasting as a method of assessing model performance.Additional observations that could be used are mentioned in Sect. 4.
-Ice thickness and ice volume (Bamber et al., 2013).Error estimates in bed elevation range from a minimum of ±6 m to about ±200 m, as a function of distance from an observation and local topographic variability.
-Surface speeds derived from synthetic aperture radar (SAR) data from RADARSAT representing an average of the velocity maps for the period 2006-2009.The mean measurement and processing errors are less than 2 m a −1 in areas with low surface slope and additional slope-dependent error of 3 % in steeper areas (Joughin et al., 2010).
-Time series of mass changes from the Gravity Recovery and Climate Experiment (GRACE) for the Greenland Ice Sheet, derived from the 10-day 1 • × 1 • NASA GSFC mascon solutions for the period January 2004 to December 2010 (Luthcke et al., 2013).Using a linear least-squares fit a mass loss trend of −223 Gt a −1 with an uncertainty of ±14.3 Gt a −1 is obtained.

1989 initial state
Figure 1 shows the simulated surface speeds of the three initialized states together with the SAR-derived speeds.All initial states show well-defined ice divides and fast flowing outlet glaciers, with largest differences between simulated and observed speeds occurring in outlet glaciers.The root mean square error (RMSE) of the surface speeds is 43-46 m a −1 (Table 1), having a similar magnitude as obtained in a data assimilation study (38 m a −1 , Price et al., 2011).
Relative differences in ice thickness and surface speeds are shown in Figs. 2 and 3, respectively; absolute differences are shown in the Supplement.
Two initializations, "constant-climate" and "paleoclimate", overestimate ice thickness near the margin, but underestimate ice thickness in the interior of the ice sheet and in northwest Greenland (Fig. 2).The flux-correction method used to obtain the "flux-corrected" initial state significantly reduces the mismatch between observed and simulated ice thickness except in coastal areas in southeast Greenland.
All three initializations underestimate surface speeds in most major outlet glaciers (Fig. 3).

Hindcasts
Figure 4 shows the time series of mass change since the beginning of the GRACE period (January 2004).The GRACE time series are compared to modeled ice mass changes using a correlation and linear regression analysis for the period January 2004 to December 2011 (Table 1   coefficient, r, is greater than 0.92 in all experiments, indicating that both the amplitude and the phase of the mass change signal are well reproduced by the model.The detrended correlation coefficients, r, are 0.86 ("constant-climate", "paleoclimate") and 0.83 ("flux-corrected"); these are a measure of how well the seasonality is reproduced.Simulated mass loss trends are −139 Gt a −1 for the hindcast based on the "constant-climate" initial state (subsequently shortened to "constant-climate" hindcast), −299 Gt a −1 for the "paleoclimate" hindcast, and −198 Gt a −1 for the "flux-corrected" hindcast, differing by a factor of more than 2. Mass loss trends are consistent with respect to the grid resolution (Supplement Table 1) for horizontal grid spacings ≤ 10 km (Supplement Fig. 2).Modeled mean ice discharge is −511, −581, and −168 Gt a −1 for the "constant-climate", "paleoclimate" and "flux-corrected" hindcasts, respectively, over the period 1989-2011.None of the ice discharge time se-ries shows a significant trend.The modeled mass loss trend is thus due to the evolving climatic mass balance (Fig. 5).
Observed and simulated surface elevation changes between October 2003 and November 2009 are shown in Fig. 6.The "constant-climate" hindcast shows a thinning pattern similar to observations except for the drainage basins of Jakobshavn Isbrae, Helheim, and Kangerdlugssuaq.The "paleo-climate" and "flux-corrected" hindcasts show distinctively different patterns."Paleo-climate" exhibits areas of both strong marginal thinning and thickening; near-marginal thickening occurs almost everywhere, except near the northwest coast.The hindcast based on the "flux-corrected" initialized state reveals strong marginal thickening and interior thinning.(c) "flux-corrected" hindcast.Blue and red colors indicate the model under-and overestimates surface speed, respectively.Areas of observed surface speeds smaller than 100 m a −1 (92 % of observed speeds) are masked out because small absolute speed differences result in large relative speed differences.Also uncertainties in SAR-derived surface speeds are higher in slow-flowing areas.

Discussion
The discussion is split into two parts.First we discuss the initial states and hindcasts relative to observations.Second we address the suitability of hindcasting as a method to evaluate the performance of an ice sheet system model.

Initial states and hindcasts
Despite having comparable ice volumes, the three initial states respond differently to the applied climate forcing, and the difference is significant.Therefore reproducing observed ice volume is not a sufficient metric for evaluating initial states.
Surface speeds of most major outlet glaciers are underestimated by all three initial states relative to the 2006-2009 average.First, simulated surface speeds of initial states should be compared to the observations from 1989, but observations for this year are not available.Between 1989 and 2006 surface velocities of many outlet glaciers have experienced large increases (Joughin et al., 2004;Luckman and Murray, 2005;Howat et al., 2005).Therefore initial states seemingly underestimate speeds in these outlet glaciers.Second, a different explanation is required for the "Northeast Greenland Ice Stream" because the inland part has not undergone significant changes in flow speed (Joughin et al., 2010) and no initial state reproduces the fast-flow pattern.For the onset region of "Northeast Greenland Ice Stream" Fahnestock et al. (2001) estimate a geothermal flux 15 to 30 times higher than continental background.Such a feature is not present in the data set by Shapiro and Ritzwoller (2004)  100 m a −1 , our parameter choice favors slow flowing ice.RMSE ranges from 43 to 48 m a −1 when computed for six different combinations of year of simulation and observation (Supplement Table 4.2).This indicates that the RMSE is not strongly affected by the choice of the combination of year of observation and year of simulation.
The simulations using a "constant-climate" and a "paleoclimate" initialization underestimate and overestimate the mass loss trend, respectively, whereas the "flux-corrected" initialization shows a very good match with observations (i.e.within twice the uncertainty of the GRACE mass loss trend).As mentioned earlier, observations show a rapid increase in ice discharge since the late 1990s, which was attributed to changes at the ocean boundary.Our model does not include ocean forcing that could lead to an increase in simulated ice discharge.Therefore, underestimation of the simulated mass loss trend is expected.
Compared to an observed cumulative mass change of Because we do not attempt to simulate the observed increase in ice discharge, we compare our results to an estimate of 1996, representative of the time before the rapid increase.For the "constant-climate" and "paleo-climate" initializations, the simulated ice discharge compares well with an estimate of −480 Gt a −1 (van den Broeke et al., 2009, Fig. S3).However, the good agreement is a consequence of consistently overestimating ice thickness (Fig. 2) and underestimating surface speeds (Fig. 3).The "flux-corrected" initialization underestimates the ice discharge as a result of underestimating flow speeds in fast-flowing outlet glaciers.We find that simulated ice discharge is most sensitive to errors in ice thickness, which is in line with Larour et al. (2012b).
Surface elevation changes at Jakobshavn Isbrae, Helheim, and Kangerdlugssuaq are mostly attributed to changes in marine outlet dynamics (Howat et al., 2011) and are therefore not reproduced by our model.The "constant-climate" initialization is expected to be in equilibrium with its climate, and, consequently, surface elevation changes should result from the climate forcing, either directly through the applied climatic mass balance or, indirectly, through dynamical adjustment due to changes in climatic mass balance.
For the "constant-climate" and "paleo-climate" initialization, forcing is applied directly, whereas anomalies are used for the "flux-corrected" initialization.For the "paleoclimate" initialization this may result in a shock at the begin- ning of the hindcast.To study this effect we perform a fourth hindcast where the ERA-interim forcing is applied as anomalies relative to the climate forcing at the end of the "paleoclimate" initialization.Results reveal a mass loss trend of −410 Gt a −1 , greater than when the ERA-interim forcing is applied directly, with a surface elevation change pattern very similar to the "flux-corrected" hindcast (not shown).
Generally speaking, a model initialization may generate transients, both wanted (e.g.ongoing adjustments due to forcing history) and unwanted (e.g. from grid refinement), leading to a model drift (Gupta et al., 2012).Even sufficiently long "constant-climate" initializations are not expected to be free of transients; for example, basal hydrology may prevent a steady-state configuration in the mathematical sense (e.g.Kamb et al., 1985).An approach to mitigate the effect of model drift involves calculating model drift and subtracting it from the experiment, implicitly assuming that legacy behavior does not feed back on the experiment (Bindschadler et al., 2013).
We calculate model drift by performing a reference run for each initial state where the mean climate of the last 10 yr of the initialization is applied.For the GRACE epoch, cumulative mass changes of −50 Gt ("constant-climate"), −281 Gt ("paleo-climate"), and −205 Gt ("flux-corrected") are attributed to model drift, corresponding to 5, 13, and 15 % of the total mass change, respectively.After subtracting model drift, simulated mass loss trends are −132 Gt a −1 for the "constant-climate" hindcast, −258 Gt a −1 for the "paleoclimate" hindcast, and −169 Gt a −1 for the "flux-corrected" hindcast (Table 2).
The "constant-climate" hindcast exhibits the smallest drift, only minimally affecting the mass loss trend.In comparison, the drift of the "paleo-climate" and "flux-corrected" hindcasts are ∼ 5.5 × and 4 × higher.Now the mass loss of "fluxcorrected" is of the right order of magnitude, while "paleoclimate" is still too high.Ideally the "flux-corrected" hindcast should exhibit little to no drift at all, as flux-correction methods are supposed to compensate for drift.This is, however, not the case in our simulation.Similarly, Rahmstorf (1995) observes a residual drift in flux-corrected coupled ocean-atmosphere general circulation model (AOGCM) simulations.Improved physics, higher resolutions, and more physically consistent coupling have rendered flux corrections mostly unnecessary in AOGCMs.(Luthcke et al., 2013) and simulated daily values starting from the three initial states.The inset figure shows the regression analysis; trends are given in the legend and correlation coefficients are listed in Table 1.
1 9 9 1 1 9 9 3 1 9 9 5 1 9 9 7 1 9 9 9 2 0 0 1 2 0 0 3 2 0 0 5 2 0 0 7 2 0 0 9 year Model drift is removed by subtracting the surface elevation change time series of the drift experiment from the hindcast.Drift-corrected surface elevation changes are shown in Fig. 7.The surface elevation changes pattern changes little, still exhibiting large transients not present in observations.There is thus a question of how drift-correction affects the system's response.This effect remains unquantified and thus quantitative studies are highly desirable.

Hindcasting as a method of assessing an ice sheet system model
A comparison of observed and simulated rates of change requires a reference period covering both observations and simulations.Hindcasting provides simulated rates of change for this reference period.In other words, it adds a temporal dimension to validation efforts.Hindcasting enables a qualitative assessment of model performance relative to observed rates of change.This reduces the number of admissible initial states more rigorously than validation efforts that do not take advantage of observed rates of change.As an example, hindcasting allowed us to realize that our flux-correction method produces unrealistic surface elevation changes between 2003 and 2009.At present there are both theoretical and practical limitations for using hindcasting as a method of assessing the performance of an ice sheet system model.Theoretically, the appropriate timescale for hindcasting is unknown.Hindcasts are short (decades) compared to the timescale associated with changes in energy (thousands of years).As a consequence, even a hindcast showing good agreement with all available observations may not capture the system's true behavior.Hence hindcasting does not identify the initial state representing the system's true state.Unfortunately the distribution of energy within an ice sheet cannot be measured directly.The age field, however, exhibits similar timescales as energy and may thus serve as a surrogate.Practically, the duration of hindcasts is limited by the length of observational records.In consequence of those limitations, a quantitative assessment of model performance is currently not possible with hindcasting.
We have not used all available observations.Additional observations that can be used for validation are, for example, spatially distributed mass loss estimates from GRACE (e.g.Luthcke et al., 2013) and temporal variation in ice surface velocity (Joughin et al., 2010;Moon et al., 2012).The aforementioned modeled age field could be compared to dated isochrones.Such a comparison provides a strong validation metric because the age field is, first, representative of long energy timescales and, second, a three-dimensional field (in contrast to, for example, surface elevation).For those reasons an accidental agreement between simulated and observed age of the ice is quite unlikely.Fortunately the number of both remotely sensed and in situ observations is constantly growing.Future studies should take advantage of these data sets for validation.

Conclusions
It has become clear that an initial state achieving a good fit with both spatially dense surface speeds and rates of change of total ice mass may still misrepresent important characteristics (ice discharge and ice thickness are two examples).Informally speaking, it is possible to "get the right result for the wrong reason".Our analysis identifies spatially dense time series of observations as preferred validation metrics, thus corroborating the assertion of Vaughan and Arthern (2007).
In this study the thermal and dynamical states (i.e. the distribution of internal energy and momentum) remain misrepresented by the "paleo-climate" and "flux-corrected" initial states, with potential implications for prognostic simulations.Removing the model drift in this comparison changes the mass loss trend, but the spatial patterns are not significantly affected.Conclusions regarding the performance of initialization techniques, however, may only apply to our particular choice of ice sheet system model, and may not be transferable to other combinations.
We have demonstrated that ice sheet system models are capable of transforming climate forcing into simulated mass loss from the ice sheet in a realistic way.Our model successfully reproduces the seasonal signal and decadal trends in mass loss as well as the velocity structure using only a few ice-sheet-wide parameters."constant-climate" and "paleoclimate" show reasonable agreement with the observed ice discharge of 1996.However, the subsequent increase is not captured because of missing model physics at the ocean boundary.Clearly, future work should focus on theoretical models of ice-ocean interaction, and their numerical implementations.
The comparison of initialization strategies shows that the model response to climate forcing is sensitive to the initial state on decadal timescales, emphasizing the importance of well-validated initial states for reliable projections.We recommend thorough validation of initial states used for prognostic simulations, and hindcasting offers a viable validation strategy.Hindcasting is not limited to our particular choice of an ice sheet system model.
While a direct validation in the engineering sense may not be applicable, our view of ice sheet models as ice sheet system models allows validation of another type.Furthermore, hindcasting can be part of a concerted effort to validate ice sheet models.Other parts include formal sensitivity analyses to assess error propagation in forward models as, for example, carried out by Larour et al. (2012a,b).Ultimately, validation may be integrated in statistical frameworks to quantify uncertainties in ice sheet evolution due to different sources of model and observation uncertainty (c.f.Steinschneider et al., 2012, for an example in hydrologic modeling).

Fig. 2 .
Fig. 2. Relative difference in ice thickness (model-observation)/observation.(a) "constant-climate" hindcast.(b) "paleo-climate" hindcast.(c) "flux-corrected" paleo-climate hindcast (only shown for comparison, as ice thickness was part of the initialization).Blue and red colors indicate the model under-and overestimates ice thickness, respectively.Brown color indicates areas where observed ice thickness is zero but the model simulates an ice thickness in excess of 100 m.

Fig. 3 .
Fig. 3. Relative difference in surface speeds (model-observation)/observation.(a) "constant-climate" hindcast.(b) "paleo-climate" hindcast.(c)"flux-corrected" hindcast.Blue and red colors indicate the model under-and overestimates surface speed, respectively.Areas of observed surface speeds smaller than 100 m a −1 (92 % of observed speeds) are masked out because small absolute speed differences result in large relative speed differences.Also uncertainties in SAR-derived surface speeds are higher in slow-flowing areas.
used as the basal boundary condition for the conservation of energy equation.Finally, our choice of the three spatially uniform scalar parameters (Supplement), controlling ice dynamics and basal processes, attempts to minimize the root mean square error.Because 92 % of observed surface speeds are smaller than www.the-cryosphere.net/7

− 1695
Gt for the 2004-2010 period, simulated values are −1005 Gt ("constant-climate"), −2134 Gt ("paleo-climate"), and −1358 Gt ("flux-corrected").van den Broeke et al. (2009) report an equal split between surface processes and ice dynamics for 2000-2008 mass loss.By assuming the same equal split between surface processes and ice dynamics during the 2004-2010 analysis period, we expect a simulated cumulative mass change of ∼ −848 Gt from changes in climatic mass balance, corresponding to a trend of ∼ −113 Gt a −1 .Of the three hindcasts only "constantclimate" agrees reasonably well with this estimate.

Fig. 4 .
Fig.4.Observed and modeled cumulative mass changes.10-day solutions from GRACE observations(Luthcke et al., 2013) and simulated daily values starting from the three initial states.The inset figure shows the regression analysis; trends are given in the legend and correlation coefficients are listed in Table1.

Fig. 5 .
Fig. 5. Simulated mass fluxes from 1990 to 2010.Climatic mass balance (solid line) and ice discharge (dashed line).For comparison the ice discharge estimate for 1996 (van den Broeke et al., 2009) is shown (black dotted line).Time series are smoothed with a 13month running-average filter.

Table 1 .
Initial (Luthcke et al., 2013)rmance statistics.The RMSE is given for the comparison of simulated surface speeds of the initial state (1989) and the observed 2006-2009 average(Joughin et al., 2010); see Supplement Table4.1 for details.Also listed are observed ice volume and GRACE trend(Luthcke et al., 2013), the correlation coefficient, r, and the detrended correlation coefficient, r.Ice volume is calculated by integrating ice thickness in space and multiplying by the density of ice (910 kg m −3 ).* Ice thickness was used in the initialization, hence this metric is not available for validation, and only shown for comparison.

Table 2 .
Simulated mass loss trends before and after drift removal for the three hindcasts.