assessment of key model parametric uncertainties in projections of Greenland Ice Sheet behavior.

. Lack of knowledge about the values of ice sheet model input parameters introduces substantial uncertainty into projections of Greenland Ice Sheet contributions to future sea level rise. Computer models of ice sheet behavior provide one of several means of estimating future sea level rise due to mass loss from ice sheets. Such models have many input parameters whose values are not well known. Recent studies have investigated the effects of these parameters on model output, but the range of potential future sea level increases due to model parametric uncertainty has not been characterized. Here, we demonstrate that this range is large, using a 100-member perturbed-physics ensemble with the SICOPOLIS ice sheet model. Each model run is spun up over 125 000 yr using geological forcings and subsequently driven into the future using an asymptotically increasing air temperature anomaly curve. All modeled ice sheets lose mass after 2005 AD. Parameters controlling surface melt dominate the model response to temperature change. After culling the ensemble to include only members that give reasonable ice volumes in 2005 AD, the range of projected sea level rise values in 2100 AD is ∼ 40 % or more of the median. Data on past ice sheet behavior can help reduce this uncertainty, but none of our ensemble members produces a reasonable ice volume change during the mid-Holocene, relative to the present. This problem suggests that the model’s exponential relation between temperature and precipitation does not hold during the Holocene, or that the central-Greenland temperature forcing curve used to drive the model is not representative of conditions around the ice margin at this time (among other possibilities). Our simulations also lack certain observed physical processes that may tend to enhance the real ice sheet’s response. Regardless, this work has implications for other studies that use ice sheet models to project or hindcast the behavior of the Greenland Ice Sheet.


Abstract.
Lack of knowledge about the values of ice sheet model input parameters introduces substantial uncertainty into projections of Greenland Ice Sheet contributions to future sea level rise. Computer models of ice sheet behavior provide one of several means of estimating future sea level rise due to mass loss from ice sheets. Such models have many input parameters whose values are not well known. Recent studies have investigated the effects of these parameters on model output, but the range of potential future sea level increases due to model parametric uncertainty has not been characterized. Here, we demonstrate that this range is large, using a 100-member perturbed-physics ensemble with the SICOPOLIS ice sheet model. Each model run is spun up over 125 000 yr using geological forcings and subsequently driven into the future using an asymptotically increasing air temperature anomaly curve. All modeled ice sheets lose mass after 2005 AD. Parameters controlling surface melt dominate the model response to temperature change. After culling the ensemble to include only members that give reasonable ice volumes in 2005 AD, the range of projected sea level rise values in 2100 AD is ∼40 % or more of the median. Data on past ice sheet behavior can help reduce this uncertainty, but none of our ensemble members produces a reasonable ice volume change during the mid-Holocene, relative to the present. This problem suggests that the model's exponential relation between temperature and precipitation does not hold during the Holocene, or that the central-Greenland temperature forcing curve used to drive the model is not representative of conditions around the ice margin at this time (among other possibilities). Our simulations also lack certain observed physical processes that may tend to enhance the real ice sheet's response. Regardless, this work has implications for other studies that use ice sheet models to project or hindcast the behavior of the Greenland Ice Sheet.

P. J. Applegate et al.: Parametric uncertainty in Greenland Ice Sheet projections
Ice sheet models provide an additional way of assessing future sea level change due to Greenland ice loss (e.g. Huybrechts and de Wolde, 1999;Greve, 2000;Gregory and Huybrechts, 2006;Price et al., 2011). These models typically include simplified treatments of ice flow, basal sliding, snowfall, and surface melting. The ice sheet modeling community has developed advanced treatments of all these processes, plus enhanced basal sliding due to surface melting, ice shelf growth, calving, and sub-shelf melt (e.g. Parizek and Alley, 2004;Alley et al., 2008;Pollard and DeConto, 2009;Walker et al., 2009;Bueler and Brown, 2009;Robinson et al., 2010;Price et al., 2011;Larour et al., 2012). These new treatments are not implemented in all models at the present time. However, the models show remarkable success in simulating many aspects of ice sheet behavior over millennial time scales and longer (e.g. van Tatenhove et al., 1995van Tatenhove et al., , 1996Greve, 1997;Simpson et al., 2009;Pollard and DeConto, 2009).
Remaining challenges in assessing future Greenland Ice Sheet changes include (1) characterizing model response to parameter choices, (2) establishing an initial state for prognostic simulations, and (3) matching data on the ice sheet's past behavior (van der Veen, 2002;Heimbach and Bugnion, 2008;Aschwanden et al., 2009;Stone et al., 2010;Greve et al., 2011). Ice sheet models have many uncertain parameters, and the choice of parameter values has a strong influence on modeled behavior (Stone et al., 2010;Greve et al., 2011). Because the thermal field within the ice sheet is mostly unknown (cf. Greve, 2005), ice sheet models are "spun up" to the present using reconstructed former surface temperatures and sea levels. Most models are tuned to produce an acceptable match to the modern geometry of the ice sheet (e.g. Ritz et al., 1997;Stone et al., 2010;Greve et al., 2011; for fits to paleo-data, see Tarasov and Peltier, 2003;Lhomme et al., 2005;Simpson et al., 2009). Achieving a good fit between the modeled and observed ice thickness distributions at the end of model spinup is challenging (Aschwanden et al., 2009;Greve et al., 2011), and simulated ice volumes at the end of spinup runs are generally larger than expected (e.g. Heimbach et al., 2008;Stone et al., 2010;Robinson et al., 2010;Vizcaino et al., 2010;Greve et al., 2011;cf. Bamber et al., 2001). Finally, data on past ice sheet variations (e.g. Alley et al., 2010, and references therein) provide a check on ice sheet models. If a model reproduces past changes well, then we can have more confidence in its projections of future changes (cf. Oreskes et al., 1994).
Perturbed physics tuning exercises may help address these challenges. In a perturbed physics ensemble, the model is run many times with different parameter combinations to identify a group of runs that provide a reasonable fit to observations. These "good" ensemble members are likely more reliable estimators of future behavior than the ensemble as a whole (cf. Weigel et al., 2010). This approach is well established in climate modeling (e.g. climateprediction.net; Stainforth et al., 2005;Piani et al., 2005), and a small but growing number of ice sheet modeling studies use ensemble methods (e.g. Tarasov and Peltier, 2004;Napieralski et al., 2007;Hebeler et al., 2008;Stone et al., 2010;Robinson et al., 2011;Born and Nisancioglu, 2011). Other techniques for assessing model sensitivity also exist (Heimbach and Bugnion, 2008;Fitzgerald et al., 2012;Sect. 5.5, below).
Here, we present results from a small perturbed-physics ensemble with the ice sheet model SICOPOLIS (Greve, 1997;Greve et al., 2011;sicopolis.greveweb.net). Our approach builds on existing work (Stone et al., 2010) by using a spinup procedure that takes past climate variability into account (seaRISE partners, 2008; Greve et al., 2011). The results indicate that our present uncertainty about the best values of model parameters translates to a large spread among model-based projections of future Greenland Ice Sheet behavior.
The paper proceeds as follows. Section 2 provides a brief description of ice sheet models and the processes they represent. In Section 3, Methods, we describe the specific ice sheet model that we use (SICOPOLIS; Greve, 1997;Greve et al., 2011), as well as our ensemble design, climate forcing time series, and ensemble culling method. Section 4, Results, describes similarities and differences among the ensemble members during different parts of the model spinup period, and discusses the effects of ensemble culling on the range of model-projected future sea level increases from the Greenland Ice Sheet. Section 5, Discussion, treats the success of our ensemble in addressing the three modeling challenges identified above. Section 6, Conclusions, emphasizes the main outcomes of the study and provides some caveats that should be borne in mind when interpreting our model output.

Overview of ice sheet processes and models
In this section, we provide a brief, qualitative overview of selected processes that influence ice sheet behavior and how ice sheet models represent these processes. Alley et al. (2010) also give a succinct description of ice sheet processes. Fuller treatments can be found in Hooke (2005), Blatter (2009), andRutt et al. (2009). Kirchner et al. (2011) provide an accessible description of the capabilities and limitations of many ice sheet models.

Ice sheet processes
Ice sheets collect snow on their high, cold, central areas and translate this material as ice to their lower, warmer margins, where mass is lost. Land ice margins lose most mass through melting. Where ice margins enter water, the ice may go afloat to form ice shelves. Marine ice margins lose mass through calving and melting induced by warm water, and may gain some mass through freeze-on (Walker et al., 2009). The translation of mass from ice sheet accumulation areas to the margins takes place through ice flow and basal sliding. Both these processes are temperature-dependent; ice deforms more easily at higher temperatures (Paterson and Budd, 1982;Greve and Blatter, 2009), and basal sliding operates only where the ice-bed interface is close to the pressure melting point (Hindmarsh and Le Meur, 2001). Basal sliding can be much faster than internal ice deformation, especially if the bed is lubricated by water or sediment. Both the Greenland and Antarctic ice sheets contain ice streams, concentrated zones in which ice velocities are many orders of magnitude faster than on the adjacent, non-streaming ice (e.g. Fahnestock et al., 2001). Finally, the surface mass balance of ice sheets may influence basal conditions through the penetration of surface meltwater to the bed (Zwally et al., 2002;Parizek and Alley, 2004).
The relative importance of the processes listed above is poorly known and varies over time, reinforcing the need for model sensitivity testing and comparison to data on Greenland Ice Sheet's past behavior. At present, mass loss from the Greenland Ice Sheet is divided about equally between surface melting and calving, with large uncertainties (Alley et al., 2010, and references therein;Robinson et al., 2011). Surface melting likely dominated calving during the Eemian, when the retracted ice margin would have had less contact with the ocean. First-order calculations suggest that increased ice discharge will be more important than surface mass balance changes in determining the ice sheet's behavior over the next century (Pfeffer et al., 2008), but this hypothesis needs further testing.

Model structure, limitations, and ongoing improvements
Ice sheet models capture the crucial insight that ice sheets contribute to sea level rise when mass gain from snowfall is smaller than mass loss due to surface ablation and calving. Conceptually, ice sheet models are computer representations of a differential equation in which the time rate of change of ice mass is a function of the instantaneous climate boundary conditions and the present state of the ice sheet. The climate boundary conditions include gridded surface temperatures, surface precipitation values, and (more rarely) ocean temperatures; these boundary conditions primarily affect the surface mass balance of the ice sheet. The ice sheet's state includes the thickness and surface elevation of ice in each map-plane grid cell, as well as the distribution of heat within the ice body; these fields primarily affect ice flow and basal sliding. Ice surface elevations also influence surface mass balance because temperatures, and thus melting, decline with increased elevation (e.g. Born and Nisancioglu, 2011; Sect. 5.2). Many ice sheet models use the shallow-ice approximation (e.g. Hutter, 1983;Greve and Blatter, 2009) to solve for the evolution of ice thicknesses over time. The shallow-ice approximation applies over the bulk of terrestrial ice sheets, but fails in regions of fast flow such as ice streams (see Joughin et al., 2010, for Greenland surface velocity maps) and in ice shelves.
A pattern scaling approach provides the climate boundary conditions for most such models. This method adjusts modern-day temperature and precipitation grids according to data from ice cores or climate model simulations (e.g. Pollard and PMIP participating groups, 2000;Kirchner et al., 2010).
Many processes are missing from standard shallow-ice models such as SICOPOLIS (e.g. Little et al., 2007). Besides the enhanced dynamics that are needed to represent ice streams and ice shelves, these missing processes act at the interfaces between the ice sheet and its bed (including the groundwater system), the atmosphere, and the ocean. Lacking models that resolve these extra processes, it is very difficult to say which ones will matter to the future of the ice sheet and which will not.
Given the inherent complexity and high computational costs of a "super-model" that would resolve all relevant effects, the modeling community is instead developing separate treatments of the missing processes. These treatments include improved methods for solving the equations of ice flow (e.g. Pattyn et al., 2008;Hindmarsh, 2009;Price et al., 2011;Larour et al., 2012) and for coupling ice sheet models to climate models (e.g. Fyke et al., 2011). Higher-order and full-Stokes ice flow models resolve stress components not represented in shallow-ice models, potentially allowing more-accurate simulation of ice streams and ice shelves. Full coupling of ice sheet models to climate models allows better accounting of changes in precipitation distribution and the energy available for melt, among other processes. Improved methods for treating calving (e.g. Alley et al., 2008;Nick et al., 2009), subglacial hydrology (e.g. Schoof, 2010), sediment production and deformation (e.g. Pollard and DeConto, 2003;Rathbun et al., 2008), and basal sliding (Alley, 2000;Hindmarsh and Le Meur, 2001) are also open areas of research. Given these ongoing model development efforts, our work with a shallow-ice model may seem misplaced. However, improved ice sheet models will have to meet the same robustness requirements that we impose on a simple model. Thus, our work provides a template for ice sheet model assessment that later studies can apply to next-generation models.

Methods
In our sensitivity tests, we applied the Latin hypercube ensemble methods of Stone et al. (2010) to the SICOPOLIS ice sheet model, as set up by Greve et al. (2011). The Latin hypercube methods provide a quasi-random sampling of parameter space that is more even than that produced by Monte Carlo methods (Bevington and Robinson, 2003;Saltelli et al., 2008) and avoids wasting model evaluations on uninfluential parameters, as can happen with a grid design (Urban and Fricker, 2009). We can thus make a reasonable www.the-cryosphere.net/6/589/2012/ The Cryosphere, 6, 589-606, 2012 exploration of parameter space with a relatively small number of model evaluations.

Model description
As noted above, we carried out our simulations with the SICOPOLIS ice sheet model (development v. 3.0; current as of 16 November 2010). Greve (1997) and Greve et al. (2011) provide a full description of the model, and we refer interested readers to those papers for more information. SICOPO-LIS is broadly comparable to most other large-scale ice sheet models, such as Glimmer (Rutt et al., 2009).
The model setup that we use is specifically intended for projecting future sea level change (seaRISE partners, 2008; Greve et al., 2011). The model grid has a horizontal resolution of 10 km and includes 81 points in the vertical direction scaled to the local ice thickness. These points are concentrated near the base, where the bulk of ice deformation occurs. The model time step is 1 yr. Consistent with Greve et al. (2011), ice extent is restricted to land grid cells. Ice that flows into the ocean is lost, instead of forming ice shelves.
SICOPOLIS' computational efficiency allows us to incorporate many thousands of years of geological data into model spinup (Sect. 3.3, below). Higher-order models provide improved representations of ice flow where the shallow-ice approximation fails (Pattyn et al., 2008;Hindmarsh, 2009), but typically require more computing time than shallow-ice models.

Ensemble design
We vary five model parameters among 100 ensemble members ( Fig. 1; Supplement). The number of model runs was chosen to achieve a reasonable tradeoff between covering parameter space and minimizing computation time. For comparison, our total number of evaluated model years (12.65 × 10 6 ) is slightly larger than that of Stone et al. (2010), who performed a larger number of shorter model runs using the same time step.
The free parameters (and their ranges) are as follows: 1. The ice flow enhancement factor (1-5, dimensionless) corrects for differences between the rheology of ice as it is measured in the laboratory and that observed on an ice sheet scale. These differences are likely due to impurities and anisotropic fabric within the ice (Greve, 1997, his Eqs. 3 and 4; see also Rutt et al., 2009, their Eq. 9). Larger values of the enhancement factor mean that the ice flows more easily. The enhancement factor likely varies throughout the ice body depending on the age of the ice (Greve, 1997), but we hold it constant.
4. The geothermal heat flux (30-70 mW m −2 ) varies over the Earth's surface, but is difficult to measure under the ice sheet (see discussion in Stone et al., 2010). Thus, it is often taken to be constant for purposes of ice sheet modeling (e.g. Ritz, 1997). A large geothermal heat flux value leads to a larger thawed bed area over which basal sliding can take place.
5. Finally, the basal sliding factor (0-20 m yr −1 Pa −1 ) determines how rapidly the ice slides over its bed where the interface is not frozen (Greve and Otsu, 2007, their Eq. 2). As expected, higher basal sliding factors mean larger sliding velocities. Like geothermal heat flux, the basal sliding factor likely varies under the ice sheet, but the spatial distribution of basal sediment is poorly known.
For the first four parameters, the ranges are roughly the same as those investigated by Stone et al. (2010), who took their ranges from data-based studies in the literature (e.g. Dahl-Jensen and Gundestrup, 1987;Braithwaite, 1995). We expanded the ranges for the positive degree-day factors so that the EISMINT-3 preferred values (Huybrechts et al., 1998) lie within the investigated range, instead of at one end, as in Stone et al. (2010). We also expanded the range of the geothermal heat flux parameter; Stone et al. (2010) found that this parameter was relatively uninfluential, and we hypothesized that a larger range might show an effect. The range we investigate is still within previous estimates (Greve, 2005;Buchardt and Dahl-Jensen, 2007;Stone, 2010). The basal sliding parameter ranges from 0 to about double the best value identified by Greve and Otsu (2007). This list of free parameters is somewhat different from that used by Stone et al. (2010), but consistent with Ritz et al. (1997) in that we fix the atmospheric temperature lapse rates (Fausto et al., 2009) and include the basal sliding factor as a free parameter. In effect, we take the surface temperature and precipitation as given, even though these data sets contribute additional uncertainty to projected ice volumes (van der Veen, 2002;Stone et al., 2010).

Initial condition and climate forcing time series
All runs were driven by the same surface temperature, precipitation, and sea level forcings. The paleoclimate spinup (Fig. 2) closely resembles that of Greve et al. (2011). We began with the observed modern ice thickness and bedrock elevation grid (Bamber et al., 2001) at −125 ka, during the Eemian interglacial. This initial condition is not ideal; the The Cryosphere, 6, 589-606, 2012 www.the-cryosphere.net/6/589/2012/ Fig. 1. Parameter combinations used in the perturbed-physics ensemble, as projected onto two-dimensional slices through the fivedimensional space (all parameters were varied simultaneously). Dashed lines indicate EISMINT-3 best estimates for most model parameters (Huybrechts, 1998;Stone et al., 2010) except the basal sliding factor, which comes from Greve and Otsu (2007). Blue crosses indicate parameter combinations that are consistent with the modern ice volume after model spinup (within 10 % of the estimated modern ice volume in 2005 AD; Sect. 3.4).
ice sheet contained ∼30-85 % of its present volume during the Eemian (Alley et al., 2010, and references therein). We comment on the potential effects of this initial condition in Sect. 5.4. From −125 ka onward, we drove the model using a temperature anomaly curve based on the GRIP ice core oxygen isotope record (Dansgaard et al., 1993;Johnsen et al., 1997) and background sea levels from the SPECMAP compilation of oxygen isotope measurements in deep-sea sediment cores (Imbrie et al., 1984). The transfer functions for converting oxygen isotope measurements to surface temperatures and sea levels are given in Greve et al. (2011). We assume a constant relation between GRIP ice core oxygen isotope values and surface temperatures (cf. Cuffey and Clow, 1997). Precipitation changes by ∼7 % for each degree of temperature change relative to the present (Greve et al., 2011, their Eq. 6;cf. van der Veen, 2002;. Modern-day surface temperature and precipitation grids are from Fausto et al. (2009) and Ettema et al. (2009); these patterns are scaled in the model according to the calculated temperature and precipitation anomalies.
Near the end of the paleoclimate spinup, we substituted an instrumental record of southwestern Greenland mean annual temperature anomalies (Vinther et al., 2006) for the GRIPbased temperatures (Fig. 2). The Vinther et al. (2006) compilation covers the years 1784-2005 AD, but we chose to begin the instrumental period in 1840 when the record becomes more complete. The Vinther et al. (2006) temperatures help us to capture the interannual variability that could be important in explaining modern mass balance trends (Alley et al., 2010;Zwally et al., 2011;Sect. 5).
After 2005 AD, the surface temperature anomaly increases according to (1)  Greve et al., 2011) and are based on oxygen isotopes in ocean sediment cores (Imbrie et al., 1984) and in central-Greenland ice cores (Dansgaard et al., 1993;Johnsen et al., 1997). 1840-2005 AD temperatures come from southwestern Greenland observations (Vinther et al., 2006). Future temperatures assume an asymptotic increase to ∼5 Greve (2000). The final temperature rise is reasonable, given that Greenland is expected to warm ∼1.5-2 times as much as the global average (Church et al., 2001). We used this temperature forcing curve instead of one produced by a climate model, because our goal is to highlight parametric uncertainties within the ice sheet model. Differences among surface temperature change projections Meehl et al., 2007) represent another layer of uncertainty.
We hold the sea level anomaly constant between 1840 AD and the end of the simulations, consistent with some earlier simulations of future Greenland evolution (Huybrechts, 1998). This assumption is conservative, in that allowing sea level to rise would increase modeled ice loss. However, the real ice sheet is likely much less sensitive to changes in sea level than temperature, even considering the large-amplitude sea level changes that take place over glacial-interglacial cycles (Alley et al., 2010;Sect. 5.4). Future sea level around Greenland is a function of thermal expansion of ocean water and the behavior of ice masses other than the Greenland Ice Sheet, which we do not include in our modeling exercise.
Moreover, large mass loss from an ice sheet tends to lower sea level nearby, due to gravitational effects (Gomez et al., 2010). Thus, sea level is likely a second-order control on future Greenland mass loss, and even the sign of this local sea level change is uncertain (note that global mean sea level will probably rise over the next century; Meehl et al., 2007;Pfeffer et al., 2008).

Culling the ensemble
We evaluate the trustworthiness of each run by comparing the simulated total ice volume in 2005 AD to the modern ice volume (∼7.2 m sea level equivalent; Bamber et al., 2001), which is somewhat uncertain. The modern ice volume estimate (Bamber et al., 2001) is based on kriging of geographically distributed ice thickness measurements made from airborne radar units. Adjacent measurements do not always agree exactly, reflecting some aggregate of bed roughness and measurement uncertainty. Moreover, observation density varies over the ice sheet; flight lines lie close to one another near airfields and become sparse farther south. Thus, the true modern ice volume could be either larger or smaller than the central estimate, but how wide the range of possible values might be is difficult to estimate. We used a windowing approach to account for uncertainty in the modern ice volume. All runs that fall within a certain distance of the modern ice volume in 2005 AD (dashed line in Fig. 3) are kept, whereas the others are discarded. We investigate window widths of ±20 %, 10 %, 5 %, and 2.5 % of the modern ice volume. These estimates of uncertainty in modern ice volume are somewhat ad hoc, but as we show later, they have little influence on our projection uncertainty.
This method of confronting model results with data (Hilborn and Mangel, 1997) is highly simplified; we discuss alternative approaches in Sect. 5.5. Total ice volume is a reasonable comparison metric because we are interested in future sea level change, and because inferences about the past state of the ice sheet are usually stated in terms of volume changes relative to the present (e.g. Alley et al., 2010, their Fig. 13). Other metrics for comparing simulated and observed ice sheets include the ice-covered area, maximum ice thickness (Ritz et al., 1997;Stone et al., 2010), the root mean squared ice thickness (e.g. Greve and Otsu, 2007), the distribution of ice surface velocities (e.g. Aschwanden et al., 2009), and the partitioning of mass loss between surface melting and calving . Using an aggregate measure such as total ice volume helps avoid nontrivial statistical issues with autocorrelation, in which adjacent residuals between observations and model predictions are not independent of one another (e.g. Bloomfield and Nychka, 1992).
As noted above, data on the ice sheet's past behavior (Alley et al., 2010, and references therein) also provide constraints on model behavior. Consistent with most earlier Greenland Ice Sheet modeling studies, we neglect this information in our model tuning; instead, we use it to evaluate the reasonableness of the culled ensemble (Sect. 5.3 and 5.5, below).

Results
In this section, we discuss the behavior of the full ensemble (Fig. 3, gray and blue curves) before treating those runs that reproduce the modern ice volume within reasonable limits (Fig. 3, blue curves). Because the initial condition and forcings are identical for all model runs (Fig. 2), variability among runs (Figs. 3, 4) is due solely to parameter choice (Fig. 1).

Eemian through the early glacial period (−125 ka to −75 ka)
As noted above, all simulations start from the modern ice geometry, with an ice volume of 7.2 m sea level equivalent (Fig. 3). Temperatures begin a few degrees below modern values (Fig. 2)  1997). From this maximum, temperature and sea level generally decline. Ice volumes increase over the same period, stabilizing sometime before −75 ka (Fig. 3).
Despite the general resemblance of the model curves to one another (Fig. 3), there is substantial divergence among ensemble members. The spread among model runs is most noticeable during the Eemian warmth; some model realizations produce a nearly ice-free Greenland, whereas in others the ice volume changes only slightly.

Early glacial period through the early Holocene
(−75 ka to −10 ka) Between −75 and −10 ka, simulated ice volumes are remarkably stable and the spread among runs, while substantial, is much smaller than during the Eemian (Fig. 3). The stability of ice volumes seems counterintuitive given the large temperature fluctuations in the GRIP record (up to ∼15 • C; Fig. 2). However, these fluctuations are only a few ka long, and the resulting ice losses due to mass balance changes are small compared to simulated ice volumes.

Early Holocene to the beginning of the preindustrial period (−10 ka to −0.16 ka, or 1840 AD)
Despite relatively stable Holocene temperatures (Fig. 2), simulated ice volumes generally decrease between −10 ka and −0.16 ka (1840 AD; Fig. 3). Many runs show a very slight growth near the end of this period, reflecting the "Little Ice Age" as seen in the central Greenland ice cores (Alley et al., 2010, and references therein). As for the Eemian, the Holocene warmth produces considerable spread among the individual model runs. One run even grows for a few ka before slowly shrinking. This realization (#95, Supplement) has exceptionally low values of both the ice and snow positive degree-day factors, which allow increased precipitation (Sect. 3.3; Greve et al., 2011) to overtake ablation temporarily.

1840 AD to 3000 AD
The Vinther et al. (2006) instrumental temperatures begin the climb out of the "Little Ice Age", and the assumed future climate trajectory builds off the end of their record. Ice volumes decline up to ∼10 cm sle between 1840 AD and 2005 AD, with larger decreases thereafter. As during the Eemian and earlier in the Holocene, the spread among model realizations is large.

Long-term future response
All runs continue to lose mass after 3000 AD. At the end of the simulations in 3500 AD, all modeled ice sheets have negative mass balances except those with very small remaining ice volumes (Fig. 3). Thus, our model runs show no equilibration to imposed climate warming, at least within ∼1000 yr of temperature stabilization (Fig. 2).

The culled ensemble
Culling the ensemble reduces the divergence among model runs during warm periods (Figs. 3, 4). The 27 ensemble members that lie within 10 % of the estimated modern ice volume in 2005 AD change by comparable amounts between warm and cold periods. Curiously, the spread among these runs during cold periods is almost as large as that of the full ensemble.
Depending on the culling window width (Sect. 3.4, above), we have different numbers of ensemble members remaining (Fig. 4). An assumed 20 % uncertainty in modern estimated ice volume leaves 52 ensemble members out of a possible 100. For 10 %, 5 %, and 2.5 % uncertainty in modern ice volume, the ensemble culling leaves 27, 12, and 6 model runs, respectively.

Future sea level change in the culled ensemble
The spread in projected future sea level change among model runs is large, even after ensemble culling (Fig. 4). For 2100 AD and an assumed 10 % uncertainty in modern Greenland ice volume, the median change in global mean sea level due to Greenland mass loss is ∼0.085 m, and the range among the 27 model runs that meet the modern ice volume criterion is ∼0.065 m, for a fractional uncertainty of ∼75 %. Much of this uncertainty persists even if we use stricter culling criteria -for example, assuming that we know the modern ice volume to within 2.5 % gives a fractional uncertainty of ∼40 % in 2100 AD.

Discussion
Our results highlight the challenge presented by parametric uncertainty for projections of future Greenland Ice Sheet behavior, building on previous work in this area (e.g. van der Veen, 2002;Stone et al., 2010). In the Introduction, we identified three additional modeling problems that we hoped to address with our perturbed-physics ensemble. These were (1) characterizing model response to parameter choice, (2) establishing an initial state for prognostic simulations, and (3) matching data on the ice sheet's past behavior. We address the ensemble's success in meeting these goals (Sects. 5.1-5.3) before discussing the model's sensitivity to initial conditions, the relative influence of sea level and temperature (Sect. 5.4), and alternative methods and data sets for assessing model performance (Sect. 5.5). Fig. 4. Histograms of modeled ice volume change relative to 2005 AD (top) and the effects of different assumed uncertainties for the modern ice volume on the median and range of ice volume change projections and hindcasts (bottom). Y -axis scaling is the same for all panels in the top row, but differs among panels in the bottom row. Color coding is the same as in Fig. 2; gray: all model runs; dark blue: model runs that lie within 10 % of the estimated modern ice volume in 2005 AD; red line: "paleoclimate spinup with additional tuning" model run from Greve et al. (2011). The green points with error bars in the top panels indicate assessed changes in the ice sheet, relative to the modern, from Alley et al. (2010, their Fig. 13). None of our model runs produce a smaller-than-today ice volume during the mid-Holocene (−5 ka). For the future time slices (2100, 2300, and 3000 AD), the range of potential future ice volume changes is always at least 40 % of the median, regardless of how strictly the ensemble is culled.

Parameter choice, simulated modern ice volume, and ice sheet sensitivity
Given these model results (Figs. 3, 4), we might ask which parameter values are most consistent with the modern observed ice volume. If there are one or more parameter combinations that match the modern condition well, these values can be used in other modeling experiments. This question was previously posed by Stone et al. (2010), who noted that high values of the ice flow enhancement factor and ice positive degree-day factor yielded the best matches with observed total ice volume; the other parameters played smaller roles. Our results are consistent with Stone et al. (2010) in that the ice PDD factor has a dominant influence on simulated ice volumes at the end of the spinup (Fig. 5; cf. Fig. 7 in Stone et al.). However, we found that treating the basal sliding factor as a free parameter reduces the influence of the ice flow enhancement factor.
We can make few general statements about optimal parameter combinations for ice sheet spinup (Fig. 5). None of our "best" runs, those that fall within 10 % of the modern ice volume in 2005 AD, have an ice positive degree-day factor greater than ∼15 mm day −1 • C −1 . This value is well within the range identified by Braithwaite (1995). Otherwise, the "best" runs span the entire free range of the four remaining parameters, indicating that any value of these parameters is potentially consistent with the modern ice volume. There is some suggestion that large basal sliding factors are most compatible with smaller values of the ice flow enhancement factor and vice versa (Fig. 1), but this apparent tradeoff might disappear given more model runs.
In our ensemble, the ice positive degree-day factor largely determines the near-term future ice sheet response (Fig. 6). As with modern ice volumes (Fig. 5), there is a strong relation between the ice PDD factor and 2005-2100 AD ice volume change, but the other parameters appear to be relatively unimportant (Fig. 6).
www.the-cryosphere.net/6/589/2012/ The Cryosphere, 6, 589-606, 2012 Fig. 5. Relation between input parameter values and simulated ice volumes in 2005 AD for all ensemble members (gray crosses) and runs that lie within 10 % of the estimated modern ice volume (blue crosses). Vertical dashed lines indicate EISMINT-3 best estimates for most model parameters (Huybrechts, 1998;Stone et al., 2010) except the basal sliding factor, which comes from Greve and Otsu (2007). Influential parameters are indicated by points that are tightly arranged about a curve that dips steeply from one side of the plot to the other (Saltelli et al., 2008). In our ensemble, the ice positive degree-day factor has the greatest influence on simulated modern ice volume, with the snow PDD factor and the basal sliding factor taking second and third places. However, ice PDD factors greater than ∼15 mm day −1 • C −1 appear inconsistent with the modern ice volume constraint. Fig. 6. Relation between input parameter values and simulated ice volume changes between 2005 and 2100 AD (Fig. 3). Color coding is the same as in Fig. 5. The ice positive degree-day factor appears to dominate the short-term future behavior of the model. Comparison of the spatial distribution of mean ice thickness among the 27 "best" model runs to the observed modern ice geometry. By "best" runs, we mean those that give ice volumes within 10 % of the estimated modern ice volume in 2005 AD (Bamber et al., 2001;Sects. 3.4, 4.6). 90 % range: difference between 95th and 5th percentiles of ice thickness values within each grid cell. Black line: modern ice margin; gray line: coast (Bamber et al., 2001, as gridded by Greve et al., 2011); white line: contour bounding areas where the differences between the observed and mean modeled ice thicknesses are less than 250 m. Simulated ice thicknesses are generally too large near the margins, but there are large areas of too-thin ice in the northern part of Greenland and upflow from Jakobshavn (about a third of the way northward on the western side).

Simulated modern ice thicknesses
By construction, our culled ensemble matches the modern ice volume. However, problems persist in the modeled ice thicknesses (Fig. 7). Consistent with earlier results using a similar model setup (Greve et al., 2011, their Fig. 2), the ice is generally too extensive in the south and has large gaps in the north. Our tuning exercise did cover part of the falsely icefree area in northern Greenland noted by Greve et al. (2011). Ice thickness errors of 1.5 km or more are present around the edges of the ice sheet. In particular, the ice is too thin upflow from Jakobshavn, about a third of the way northward on the west coast. It is possible that an appropriately-tuned higher-order model would produce a better fit to the observed ice thickness grid than we have yet achieved with SICOPOLIS. However, the large-scale shape of the ice sheet is more strongly controlled by surface mass balance than ice flow. Near the ice margins, thin ice and low slopes lead to small ice fluxes that are easily overwhelmed by negative mass balances (Greve, 1997;Alley et al., 2010;Born and Nisancioglu, 2011; see also Kirchner et al., 2010). We attribute the bulk of the remaining errors in geographically distributed ice thickness values to problems with the modeled mass balance.
These errors in the ice thickness grid have implications for the future evolution of the ensemble members, both in terms of ice flow and surface mass balance. As noted by Alley et al. (2010), ice flow is primarily a function of ice thickness; thus, if the ice is too thin or too thick, flow will be correspondingly too slow or too fast. This dependence is not linear, but goes with some power of the thickness, depending on whether the ice is frozen to the bed. Similarly, a local negative mass balance leads to a positive (self-sustaining) feedback, because surface lowering due to ablation brings the surface to a lower, warmer elevation where ablation progresses more rapidly (Born and Nisancioglu, 2011).

Comparison of modeled results to assessed past volume changes
As noted in the Introduction, the ice sheet's past behavior provides a check on ice sheet model results. Alley et al. (2010) give assessments of ice volume changes, relative to the present, for three time slices covered by our runs -the Eemian, Last Glacial Maximum, and mid-Holocene. These assessments are based on a combination of isostatic rebound studies (Peltier, 2004;Fleming and Lambeck, 2004) and data-constrained ice sheet modeling (Cuffey and Marshall, 2000;Lhomme et al., 2005). Because some of these studies use ice sheet models to translate observed ice core oxygen isotope values into ice sheet changes, these assessed ice volume changes are not truly independent of our results. However, they do provide a first-order check on our model output.
The culled ensemble simulates assessed ice volume changes well during the Eemian (−115 ka; Fig. 4), but has problems during the mid-Holocene and Last Glacial Maximum (−5 ka and −20 ka, respectively). In particular, none of the model runs produce a smaller-than-modern ice sheet during the mid-Holocene (−5 ka), as expected from paleo-data. The overlap between our "good" model runs and the estimated Last Glacial Maximum (−20 ka) ice volume change is also minimal, but this discrepancy could be reduced by a different weighting function for evaluating the model runs. Fig. 8. Effects of initial condition on volume curves, and the relative influence of sea level and temperature on model output. All runs displayed in this figure use a parameter combination that closely matches the estimated modern ice volume (Bamber et al., 2001) using our standard setup (run #29; Supplement). The initial condition exerts some control over simulated ice volumes (especially during the late Eemian), whereas sea level is clearly a second-order control on simulated ice volumes when compared to temperature change.
Moreover, the agreement between our simulated late-Eemian ice volume changes and previously estimated values may be fortuitous; the peak in the GRIP temperature record at −115 ka could be due to flow disturbances in the ice core (Chappellaz et al., 1997;Cuffey and Marshall, 2000).
The apparent disagreement between simulated and assessed ice volume changes in the mid-Holocene (Fig. 4) could be due to the temperature forcing curve used to drive the model, or to the assumed-constant exponential relation between surface temperature anomaly and precipitation (among other possibilities). Air temperatures over some parts of the Greenland Ice Sheet were likely warmer during the mid-Holocene than geographic scaling of the GRIP oxygen isotope record would indicate (Y. Axford, personal communication, September 2011;Dahl-Jensen et al., 1998;Young et al., 2011;see also Vinther et al., 2009). If accounted for in a model simulation, these warmer regional temperatures might bring modeled ice volumes closer to the estimated values. For example, Simpson et al. (2009) simulated a smaller ice volume at this time by imposing an additional up-to-2.5 • C warming on the scaled GRIP oxygen isotope values (cf. van Tatenhove et al., 1995(cf. van Tatenhove et al., , 1996. Further, our model runs assume a constant ∼7 % increase in precipitation per degree Celsius temperature increase (Greve et al., 2011). This relation has long been controversial (see discussion in van der Veen, 2002;, and it may be especially poor for the mid-to late-Holocene (Cuffey and Clow, 1997). Further research is needed to reduce the discrepancies between model-and data-based reconstructions of past Greenland Ice Sheet configurations.

Model sensitivity to initial conditions; relative influence of sea level and temperature change
As noted in Sect. 3.3, we use the modern ice thickness grid as the initial condition and hold the sea level anomaly constant at zero between 1840 AD and the end of the runs (Fig. 2). To test the effects of these design choices on the model output, we performed three additional model runs using parameter values that produce the best match to the modern ice vol-ume in our normal setup (run #29; Supplement). The first of these extra runs begins from an ice-free initial condition with a fully relaxed bedrock surface. The other two runs use the modern ice thickness grid initial condition, but hold either the sea level anomaly or the temperature anomaly constant at zero while the other forcing varies as usual. These tests suggest that our choice of initial condition likely affects the model results, but the lack of sea level forcing in the historical and future part of the runs is a second-order effect. Somewhat surprisingly, the ice-free initial condition leads to a simulated modern ice volume that is ∼16 % larger than when the modern ice thickness grid is used (Fig. 8). Rogozhina et al. (2011) suggest that an ice-free initial condition is a poor choice for simulations that start under colder-than-modern conditions (as we do here; Fig. 2), because the new snowfall will be colder and more resistant to flow than the real ice sheet's basal ice. We estimate the uncertainty in simulated modern-day ice volumes from this source to be ∼10 %. The uncertainty in simulated Eemian ice volumes due to the choice of initial condition is probably much larger (Fig. 8).
If the model is spun up using the ice core-based temperatures and a constant sea level anomaly of zero, the simulated ice volume is ∼2 % smaller than the base case. On the other hand, if the temperature anomaly is held constant at zero while sea level varies, the simulated modern ice volume is ∼24 % smaller. The changes in sea level over the spinup period (>120 m; Fig. 2) far exceed expected sea level rise by 3500 AD. These results confirm that sea level forcing is not a dominant control on the model output.

Alternative sensitivity testing methods and data sets for calibration
In this study, we have used a simple perturbed-physics approach to assess the effects of parametric uncertainty on ice sheet model projections. Two more advanced techniques for assessing model sensitivity that have been applied in a glaciological context are the adjoint method (Heimbach and Bugnion, 2008) and the elementary effects method Fig. 9. Comparison of simulated and observed ice sheet mass balance over the last half-century. Color coding for model results follows that in Fig. 3. The temperature data are from Vinther et al. (2006). The mass balance estimates (Rignot et al., 2008) are shown with 2σ uncertainties for clarity. As expected, simulated mass balance is closely tied to temperature; low temperatures are associated with a more positive mass balance (top and middle panels). The black curves in the bottom two panels are based on Monte Carlo sampling from the individual mass balance estimates, assuming that these estimates are independent of one another. The vertical lines in these panels are based on sampling the model curves for only those years in which mass balance observations are available, allowing for a clean comparison between the model and the data. Lines that fall close to the centers of the black curves indicate runs that agree well with the mass balance data. Because many runs that pass the modern volume test (Fig. 3) do not agree well with the mass balance data, these results suggest that recent mass balance estimates provide a useful additional constraint for ice sheet model tuning. (Fitzgerald et al., 2012). Both of these methods are more computationally efficient than our approach; moreover, the adjoint method allows investigation of spatial parameter dependencies. We chose our approach because it is conceptually simple and because it allows some insight into projection uncertainties.
The windowing approach that we use to constrain the ensemble has some precedent in the climate modeling literature (Knutti et al., 2002). However, the probabilistic insight that this method yields is limited, because model runs are assigned a binary score based on whether or not they fall within some distance of an observational constraint (in this case, the estimated modern ice volume). In contrast, a likelihood function assigns relative weights to the different model runs, but at the price of higher complexity. For examples of the likelihood function approach, see Urban and Keller (2010) and Olson et al. (2012).
Our results indicate that tuning against aggregate metrics of the modern geometry (i.e. the observed ice volume, area, and maximum thickness) provide only very weak constraints on appropriate parameters for ice sheet modeling ( Fig. 5; Stone et al., 2010). To illustrate how other data sets may contribute to ice sheet model tuning, we present a comparison to mass balance data from Rignot et al. (2008;Fig. 9). These data are broadly consistent with other estimates of Greenland Ice Sheet mass balance over the last half-century (Alley et al., 2010, their Fig. 1). The agreement between the mass balance observations and the model output is fairly good (note that the observed mass balances depend on the temperature records synthesized by Vinther et al., 2006, so the model and observations are not fully independent of each other). However, the model tends to overpredict mass losses from the ice sheet, even within the culled ensemble (Fig. 9). This comparison suggests that mass balance serves as a useful additional constraint on ice sheet model simulations: it is an achievable, but not trivial, goal.

Conclusions
Our results (Fig. 4) suggest that parametric uncertainty in ice sheet model-based projections of Greenland Ice Sheet behavior is on the order of 40-70 %, expressed as the range of www.the-cryosphere.net/6/589/2012/ The Cryosphere, 6, 589-606, 2012 plausible model outcomes divided by the median. This outcome is not sensitive to the parameter ranges we investigated, but rather depends on the uncertainty in modern ice volumes. Our results do not provide a probabilistic assessment of future ice sheet changes. That is, we make no statement about the relative plausibility of different model runs within our culled ensemble. Our analysis neglects several sources of uncertainty that will tend to increase the ranges shown in Fig. 4. In particular, we assume that the climate and subglacial topography boundary conditions are well known. Estimating the effects of errors in these data sets on modeled ice volumes is complex, but Stone et al. (2010) found that updating older EISMINT-3 data sets (Huybrechts, 1998) to their modern equivalents increased simulated equilibrium ice volumes by ∼17 % if model input parameters were held constant. It is clear that the best-estimate model parameter values depend strongly on the input data sets. The ice core-inferred paleotemperatures used to spin up ice sheet models also contribute to projection uncertainty (Stone et al., 2010; see also Rogozhina et al., 2011). Finally, future ice volumes depend on uncertain emissions trajectories and the broader climate system's response (Meehl et al., 2007).
Even our most responsive model runs may underestimate mass loss from the real ice sheet. If precipitation remains constant in the future, instead of increasing at ∼7 % • C −1 of warming (Sect. 3.3; Greve, 2011), then the ice sheet will shrink more rapidly than we project. Many mechanisms for rapid ice loss are not included in this ensemble. For example, surface melting may lead to basal lubrication and enhanced transport of ice to the margin (Zwally et al., 2002;Parizek and Alley, 2004;Bartholomew et al., 2010). We neglect this possibility here (see Greve and Otsu, 2007, for model runs with SICOPOLIS that includes this effect). Additionally, ocean warming may contribute to mass loss where the ice is in contact with the water (Straneo et al., 2010;Yin et al., 2011), and the resulting rapid thinning of marine ice margins could then propagate up ice streams to the central parts of the ice sheet. This scenario cannot be captured by shallow-ice models like SICOPOLIS, but is expected to appear in higher-order models if the ocean boundary conditions are correctly represented. Complex models typically have more parameters than simpler ones, so sensitivity experiments with higher-order models (e.g. Price et al., 2011) might lead to a wider range of future Greenland states (cf. Saltelli et al., 2008).
Given these problems, both our uncertainty estimates and our projections of future ice sheet mass loss may be too small. Despite the large variation among individual model runs, all of our modeled ice sheets lose mass from 2005 AD onwards. Thus, our work agrees with the scientific consensus, which says that sea level rise due to enhanced mass loss from the Greenland Ice Sheet in the face of surface temperature increases is very likely (Lemke et al., 2007).
The supplement includes processed model output and Octave scripts for generating the figures. Full output files from the ice sheet model are available from the Bert Bolin Centre for Climate Research data repository, http://www.bbcc.su.se/ data/.