1. Introduction
In chapter 10 (‘Global climate projections’) of the Fourth Assessment Report (AR4) of the United Nations Reference SolomonIntergovernmental Panel on Climate Change (IPCC), an increase of the mean global sea level by 18–59cm during the 21st century (more precisely: 2090–99 relative to 1980–99) is projected for the six SRES (Special Report on Emissions Scenarios) marker scenarios B1, B2, A1B, A1T, A2 and A1FI (Reference Meehl and SolomonMeehl and others, 2007). The main causes are thermal expansion of ocean water and melting of glaciers and small ice caps, and, to a lesser extent, changes of the surface mass balance of the Greenland and Antarctic ice sheets. However, recent observations suggest that ice flow dynamics could lead to significant additional sea-level rise, as stated in the AR4: ‘Dynamical processes related to ice flow not included in current models but suggested by recent observations could increase the vulnerability of the ice sheets to warming, increasing future sea level rise. Understanding of these processes is limited and there is no consensus on their magnitude’ (IPCC, 2007).
The conjectured dynamical processes are basal sliding accelerated by input of surface meltwater to the ice/bed interface and reduced flow buttressing due to a loss of floating ice (ice shelves and ice tongues) (Reference Lemke and SolomonLemke and others, 2007; Reference Meehl and SolomonMeehl and others, 2007). On the observational side, results from satellite altimetry, satellite gravimetry and other methods indicate a mass loss from the Greenland ice sheet (the focus of this study) in the range 100–250 Gt a-1 (~0.3– 0.7 mm SLE a-1; SLE means sea-level equivalent) during recent years (Reference Lemke and SolomonLemke and others, 2007). Furthermore, major outlet glaciers (Jacobshavn Isbræ, Kangerdlugssuaq Gletscher and Helheimgletscher) have sped up dramatically since the 1990s (Reference Rignot and KanagaratnamRignot and Kanagaratnam, 2006; Reference Howat, Joughin and Scambos.Howat and others, 2007; Reference JoughinJoughin and others, 2008).
The scientific community has reacted to the need for improved predictions of sea-level rise from ice-sheet models. Coordinated research projects have been launched, such as the European-led ice2sea programme funded by the European Union Framework-7 scheme (http://www.ice2sea.eu/) and the US-led, community-organized SeaRISE effort (Sea-level Response to Ice Sheet Evolution; http://websrv.cs.umt.edu/isis/index.php/SeaRISE_Assessment,http://oceans11.lanl.gov/trac/CISM/wiki/AssessmentGroup).The Japanese ice-sheet modelling community is committed to contribute to both ice2sea and SeaRISE as part of several funded research projects. In this study, we present first results obtained with the models SICOPOLIS and IcIES for the Greenland ice sheet within the framework of the SeaRISE effort.
2. Ice-Sheet Models Sicopolis and Icies
2.1. Sicopolis
SICOPOLIS (SImulation COde for POLythermal Ice Sheets) simulates the large-scale dynamics and thermodynamics (ice extent, thickness, velocity, temperature, water content and age) of ice sheets three-dimensionally and as a function of time (Reference GreveGreve, 1997; for the latest version, 3.0, used here see http://sicopolis.greveweb.net/). It is based on the shallow-ice approximation (Reference HutterHutter, 1983; Reference MorlandMorland, 1984) and the rheology of an incompressible, heat-conducting, power-law fluid (Glen’s flow law). The thermomechanical coupling is described by the temperature- and water-content-dependent rate factor in the form of Reference Greve, Weis and HutterGreve and others (1998). Isostatic depression and rebound of the lithosphere due to changing ice load is modeled by the elastic-lithosphere– relaxing-asthenosphere (ELRA) approach (see Reference Greve and BlatterGreve and Blatter, 2009, and references therein).
A particular feature of the model thermodynamics is that it distinguishes between cold ice with a temperature below the pressure-melting point and temperate ice with a temperature at the pressure-melting point, the latter being considered as a binary mixture of ice and small amounts of water. The interface that separates cold and temperate ice is tracked through the use of Stefan-type energy-flux and mass-flux matching conditions (this procedure is referred to as the ‘polythermal mode’).
Basal sliding, v b, is described by a Weertman-type sliding law with an allowance for sub-melt sliding (Reference GreveGreve, 2005),
where p and q are sliding exponents, C
b the sliding coefficient, γ the sub-melt-sliding parameter, τ
b the basal drag (shear stress), N
b the basal normal stress (counted positive for compression) and 
the basal temperature relative to pressure melting (in ˚C, always ≤0˚C). In the shallow-ice approximation, the basal normal stress, N
b, is equal to the hydrostatic basal pressure, P
b,
where ρ is the density of ice, g the gravitational acceleration and H the ice thickness. The basal drag, τ b, is equal to the basal pressure times the surface slope,
where h is the surface elevation and grad the gradient operator in the horizontal plane (Reference Greve and BlatterGreve and Blatter, 2009). The minus sign in Equation (1) indicates that the direction of basal sliding is antiparallel to the basal drag. The parameters have the values C b = 11.2ma-1 Pa-1 (sliding coefficient), p = 3, q = 2 (sliding exponents) and γ =1˚C (sub-melt-sliding parameter).
The model domain covers the entire land area of Greenland and the surrounding oceans, projected on a polar stereographic grid with standard parallel 71˚N and central meridian 39˚W. Distortions due to this projection are accounted for as metric coefficients in the model equations. The present geometry (surface and basal topographies, ice thickness, equilibrated bedrock elevation) is derived from the ‘Greenland Developmental Data Set’ (Greenland_5km_dev1.2.nc) provided on the SeaRISE website, resampled to a horizontal resolution of 10 km. In the vertical direction, sigma coordinates are used; the cold ice column, the temperate ice layer (if present) and the thermal boundary layer of the lithosphere are mapped separately to [0,1] intervals. The cold ice column is discretized by 81 gridpoints (concentrated towards the base), and the temperate ice and lithosphere layers are discretized each by 11 equidistant gridpoints. A fixed time-step of 1 year is used for all simulations in this study.
2.2. IcIES
Similar to SICOPOLIS, the model IcIES (Ice sheet model for Integrated Earth system Studies; Reference Saito and Abe-OuchiSaito and Abe-Ouchi, 2005, Reference Saito and Abe-Ouchi2010) is a three-dimensional, large-scale, dynamic/ thermodynamic ice-sheet model with shallow-ice approximation dynamics. The main differences compared to SICOPOLIS are:
Cold ice mode is employed (i.e. The Stefan-type conditions at the cold/temperate transition surface are ignored, computed temperatures above pressure melting are reset to pressure melting with no computation of the englacial water content).
Isostasy: Local-lithosphere–relaxing-asthenosphere (LLRA) approach.
Basal sliding: Different parameters for the Weertman-type sliding law (Equation (1)), namely C b = 1.61 × 10-6 ma-1 Pa-2, p = 3, q = 1 and γ→0˚C (no sub-melt sliding; instead an abrupt switch between warm-based sliding and cold-based no-slip) (Reference Huybrechts and de WoldeHuybrechts and de Wolde, 1999).
Present geometry derived from the ‘Greenland Standard Data Set’ (Greenland_5km_v1.1.nc; this affects only the immediate vicinity of Jacobshavn Isbræ).
Discretization of the vertical direction by sigma coordinates. Ice column: 26 gridpoints (concentrated towards the base); thermal lithosphere: 17 equidistant gridpoints.
Adaptive time-step, with a maximum value of 0.125 years.
The values of the main physical parameters used for the simulations in this study, for both SICOPOLIS and IcIES, are listed in Table 1.
Physical parameters of the ice-sheet models SICOPOLIS and IcIES
3. SeaRISE experiments
3.1. Paleoclimatic spin-up
In order to obtain suitable present-day configurations of the ice sheet that can be used as initial conditions for future-climate experiments, paleoclimatic spin-ups over a full glacial cycle (from 125 ka BP until today) are carried out with SICOPOLIS and IcIES. The forcing follows that specified by the SeaRISE experiments. The surface air temperature is parameterized as a function of surface elevation, h, latitude, ϕ, longitude, λ, and time, t, following Reference Fausto, Ahlstrøm, van As, Bøggild and JohnsenFausto and others (2009):
where T ma and T mj are the mean annual and mean July (summer) surface temperatures, respectively, the temperature constants are d ma = 41.83˚C and d mj = 14.70˚C, the mean slope lapse rates are γ ma = -6.309˚Ckm-1 and γ mj = -5.426˚C km-1, the latitude coefficients are c ma = -0.7189˚C (˚ N)-1 and c mj = -0.1585˚C (˚ N)-1, and the longitude coefficients are ma = 0.0672˚C (˚ W)-1 and kmj = 0.0518˚C (˚ W)-1.
The purely time-dependent anomaly term ΔT(t) describes the deviation from present-day conditions. It is based on the oxygen isotope record (δ18O) from the Greenland Icecore Project (GRIP) ice core (Reference DansgaardDansgaard and others, 1993; Reference JohnsenJohnsen and others, 1997), which was converted to a record of temperature variation from 125 ka BP to the present by the formula
The ΔT/δ18O conversion factor, d = 2.4˚C%-1, is the standard value used by Reference HuybrechtsHuybrechts (2002). The result of Equation (5) is shown in Figure 1a.
(a) Surface temperature anomaly, ΔT(t), derived from the GRIP δ18O record (Reference DansgaardDansgaard and others, 1993; Reference JohnsenJohnsen and others, 1997). (b) Sea level, z sl(t), derived from the spectral-mapping project (SPECMAP) marine δ18O record (Reference Imbrie, Berger, Imbrie, Hays, Kukla and SaltzmanImbrie and others, 1984).
For the present-day mean annual precipitation rate, P ma, present(λ, ϕ), recent data by Reference EttemaEttema and others (2009) are used. For any time, t, between 125 ka BP and today, the data are modified according to
which corresponds to a 7.3% change of precipitation rate for every 1˚C of temperature change (Reference HuybrechtsHuybrechts, 2002).
Conversion from the mean annual precipitation rate, P ma, to the snowfall rate (solid precipitation), not specified by SeaRISE, is done by SICOPOLIS on a monthly basis using the empirical relation (Reference MarsiatMarsiat, 1994)
where S mm is the mean monthly snowfall rate and T mm the mean monthly surface temperature, computed from T ma and T mj assuming a sinusoidal annual cycle. Mean monthly rainfall (liquid precipitation) is obtained as the difference between precipitation and snowfall.
IcIES handles the conversion of precipitation to snowfall rate on the basis of instantaneous rather than mean monthly surface temperature. It is assumed that precipitation falls entirely as rain if the surface temperature is above T rain = 0˚C, or entirely as snow otherwise. Similar to the PDD method (see below), the instantaneous surface temperature is described statistically as the mean monthly temperature plus superimposed Gaussian noise, so that
(cf. Reference Huybrechts and de WoldeHuybrechts and de Wolde, 1999), where the standard deviation σ = 4.5˚C.
Surface melting is not specified by SeaRISE either, and the standard implementations of SICOPOLIS and IcIES are used. For SICOPOLIS, it is parameterized by Reference ReehReeh’s (1991) positive degree-day (PDD) method, supplemented by the semi-analytical solution for the PDD integral of Reference Calov and GreveCalov and Greve (2005). Following Reference Tarasov and PeltierTarasov and Peltier (2002), the PDD factors for ice melt, ice, and snowmelt, snow, depend on the mean July surface temperature, T mj,
and
The limiting temperatures are T
w = 10˚C (warm conditions), T
c = -1˚C (cold conditions), and the limiting PDD factors are 
i.e. d–1 ˚C–1, 
i.e. d–1 ˚C–1, 
i.e. d–1 ˚C–1, 
i.e. d–1 ˚C–1 (where i.e. means ice equivalent). The standard deviation of short-term, statistical air-temperature fluctuations is σ = 5.2˚C, and the saturation factor for the formation of superimposed ice is chosen as P
max = 0.6 (Reference ReehReeh, 1991).
The method for IcIES is similar, except that constant values for the PDD factors are used, namely ice = 8 mm i.e. d–1 ˚C–1, snow = 3mm i.e. d–1 ˚C–1 (Reference Huybrechts and de WoldeHuybrechts and de Wolde, 1999), and the standard deviation of short-term, statistical air-temperature fluctuations is set to the slightly higher value σ = 5.5˚C, while the saturation factor for the formation of superimposed ice is the same as for SICOPOLIS (P max = 0.6).
Sea-level forcing, z sl, which determines the land area available for glaciation, is derived from the spectral-mapping
project (SPECMAP) marine δ18O record (Reference Imbrie, Berger, Imbrie, Hays, Kukla and SaltzmanImbrie and others, 1984) converted to global sea level by
This parameterization produces a Last Glacial Maximum (LGM) sea-level minimum of –130m at 19 ka BP and an Eemian sea-level high of 5.9 m at 122 ka BP (see Fig. 1b).
The geothermal heat flux is that of Reference Shapiro and RitzwollerShapiro and Ritzwoller (2004). Both SICOPOLIS and IcIES apply this heat flux 5 km below the ice base in order to account for the thermal inertia of the lithosphere (Reference RitzRitz, 1987).
3.2. Future-climate experiments
Currently, SeaRISE specifies four different future-climate experiments:
Experiment C1_E0: Constant climate control run; beginning at present (more precisely, the epoch 2004-1-1 0:0:0, corresponding to t = 0) and running for 500 years holding the climate steady to the present climate.
Experiment C1_E1: Like C1_E0 (constant climate forcing), but with increased basal sliding assumed. This is implemented in a simple fashion by doubling basal sliding everywhere it occurs (that is, the value of the sliding parameter, C b, is doubled).
Experiment C2_E0: AR4 climate control run; starts with the same present-day condition, but the climatic forcing (mean annual temperature, mean July temperature, precipitation) is derived from an ensemble average from 18 of the AR4 models, run for the period 2004–98 under the A1B emission scenario; beyond 2098 the climate persists to the end of the run 500 years into the future.
Experiment C2_E1: Like C2_E0 (AR4 climate forcing), but with increased basal sliding.
Computation of the solid precipitation (snowfall) and surface melting rates is the same as described in section 3.1. Further experiments are still under discussion within the SeaRISE community and will be specified later.
4. Results
4.1. Standard paleoclimatic spin-up
The present-day surface topographies of the Greenland ice sheet obtained with the standard set-up described in section 3 are shown in Figure 2. The agreement with the observed topography (also shown) is not satisfactory. In the southern half of Greenland, both SICOPOLIS and IcIES produce an ice sheet that is too large and extensive and leaves very little ice-free land. In the northern half, this is again the case for IcIES, while SICOPOLIS produces vast areas of ice-free land north of ~77˚N which have no counterpart in reality. This difference between the two models is a consequence of the different values of the PDD factors (in particular, ice) employed by the two models; the temperature-dependent parameterization of SICOPOLIS favors melting in the colder areas of North Greenland. In the far north, both models produce glaciation of the actually ice-free Peary Land (in the case of SICOPOLIS as an isolated ice cap), which may be due to either underpredicted melting or large prescribed precipitation rates. These problems point to the necessity of additional tuning in order to achieve a better agreement with observations.
Present-day surface topography of the Greenland ice sheet. (a) Result of the standard paleoclimatic spin-up with SICOPOLIS. (b) Result of the standard paleoclimatic spin-up with IcIES. (c) Observed (according to the dataset Greenland_5km_dev1.2.nc). Contour spacing is 200 m, labels are in km a.m.s.l. Brown areas mark ice-free land.
4.2. Paleoclimatic spin-up with additional tuning
For SICOPOLIS, additional tuning of the paleoclimatic spinup was achieved by modifying the PDD factors as follows:
where the PDD modification factor changes with latitude and is different for West and East Greenland. Using the observed ice margin (Fig. 2c) as a tuning target, values >1 were chosen in regions where the simulated ice sheet shown in Figure 2a has advanced too far, while values <1 were chosen where the simulated ice sheet has retreated too far. A trial-and-error procedure led to the values shown in Table 2, which provide a reasonably good fit between the simulated and observed ice margins (Fig. 3a). The remaining problem, highlighted in Figure 3b, is that the ice is generally too thick (with major exceptions in the southwest, the far southeast and the northwest), so that the simulated ice volume (8.26mSLE) exceeds the observed one (7.2mSLE) by 14.9%. In future work, we will attempt to improve the misfit by tuning the basal sliding and/or the flow enhancement factor to interferometrically measured surface velocities.
Result of the paleoclimatic spin-up with SICOPOLIS and additional tuning of the PDD factors. (a) Present-day surface topography. Contour spacing is 200 m, labels are in km a.m.s.l. Brown areas mark ice-free land. (b) Difference of simulated and observed present-day ice thickness. The thick white lines indicate the simulated (a) and observed (b) ice margins.
Latitude-dependent PDD modification factors for the tuned paleoclimatic spin-up with SICOPOLIS, defined separately for West (W) and East (E) Greenland (west and east of 44˚W, respectively). For arbitrary latitudes, linear interpolation between these values is employed
For IcIES, additional tuning of the paleoclimatic spin-up has not yet been carried out. Therefore, the standard spin-up is used for this study. The resulting present-day ice-sheet volume for IcIES is 8.23m SLE, which is 14.5% more than the observed one, but in close agreement with the present-day ice volume predicted by SICOPOLIS.
4.3. Future-climate experiments
Using the tuned spin-up of SICOPOLIS and the standard spin-up of IcIES as initial conditions, the future-climate experiments described above were carried out. For the experiments with SICOPOLIS the tuned PDD factors according to Equation (12) were used in order to guarantee consistency with the spin-up.
Simulated ice volumes over time are shown in Figure 4. The ice sheet reacts differently to the imposed scenarios; future climate warming results in more rapid ice mass loss than the ice-dynamical scenario of increased basal sliding. For all experiments, the sensitivity of SICOPOLIS is significantly larger than that of IcIES. Relative to the constant climate control run C1_E0,
Ice-volume (V) changes simulated with (a) SICOPOLIS and (b) IcIES for experiments C1_E0 (constant climate control run), C1_E1 (constant climate forcing, doubled basal sliding), C2_E0 (AR4 climate control run) and C2_E1 (AR4 climate forcing, doubled basal sliding).
increased basal sliding (C1_E1) leads to volume losses of 0.22 mSLE for SICOPOLIS and 0.12 mSLE for IcIES,
future climate warming (C2_E0) leads to volume losses of 0.79 mSLE for SICOPOLIS and 0.25 mSLE for IcIES,
future climate warming plus increased basal sliding (C2_E1) leads to volume losses of 1.01 mSLE for SICOPOLIS and 0.38mSLE for IcIES,
after 500 years of model time. It is also interesting to note that the sensitivity of the two increased basal sliding experiments relative to their control runs (that is, C1_E1 – C1_E0 and C2_E1 – C2_E0) is virtually identical for both models:
This supports the contention that the experiment minus control approach is a viable means to examine the sensitivity to ice-dynamical processes and can largely remove the artifacts of spin-up differences between models.
5. Discussion
The experiences gained with the paleoclimatic spin-ups show that bringing an ice-sheet model into a reasonable present-day configuration is a challenge (as already reported by Reference Aschwanden, Khroulev and BuelerAschwanden and others, 2009). Standard implementations tend to produce poor agreement with the observed state of the ice sheet, and careful additional tuning is required to obtain satisfactory results. We have demonstrated this using the observed surface topography as a tuning target and modifying the model PDD factors. Further observations (surface velocities, temperature profiles from deep ice cores, isochrones) can be used as additional constraints in order to tune the basal sliding and/or the flow enhancement factor, and we will attempt to do this in future work in order to optimize the representation of the present-day state of the Greenland ice sheet.
The future-climate experiments show that sensitivity between models can vary greatly. We found that the ice-sheet mass loss varies between SICOPOLIS and IcIES by a factor of ~2 for sliding experiments and a factor of ~3 for climate-warming experiments, which is the result of different parameterizations for basal sliding and surface melting used in the two models. The basal sliding parameterization of SICOPOLIS favors relatively more rapid sliding for thin ice compared to the parameterization in IcIES, because the sliding exponents p = 3 and q = 2, used in SICOPOLIS, yield v b ∝ H |grad h |3, while the sliding exponents p = 3 and q = 1, used in IcIES, yield v b ∝ H 2|grad h|3 (see Equations (1–3)). Since thin ice occurs near the margins where ice flow is generally fast, the impact of doubled basal sliding on ice-sheet decay is more pronounced in SICOPOLIS than in IcIES. The tuned PDD factors employed in SICOPOLIS lead to very high melt rates, especially in the southern half of the ice sheet. By contrast, the standard PDD factors used in IcIES produce significantly smaller melt rates, leading to the different responses of the two models in the future-climate-warming experiments.
Despite the differences between SICOPOLIS and IcIES, both models show a larger sensitivity to future climate warming than to a doubling of basal sliding. of course, the significance of this statement is limited by the fact that the assumed scenario of doubled basal sliding is a mere assumption. However, R. Greve and S. Sugiyama (http://arxiv.org/abs/0905.2027) came up with a more rational, observation-based parameterization of basal sliding accelerated by surface meltwater, and simulated the response of the Greenland ice sheet to the ‘WRE1000’ scenario (stabilization of atmospheric CO2 at 1000 ppm within the next few centuries) with SICOPOLIS. They also found that the sensitivity to the direct climate forcing is larger than that to the ice-dynamical effect, and concluded that basal sliding accelerated by surface meltwater accelerates the decay of the Greenland ice sheet as a whole significantly, but not catastrophically, in the 21st century and beyond.
In any case, the factor of 2–3 differences in sensitivity seen here alert us to the uncertainty in model outputs as a result of uncertainties in model input parameters. Surface mass balance and basal sliding in SICOPOLIS and IcIES are both parameterized within the range of uncertainties in our understanding of these processes, and still produce significantly different predictions for the ice-sheet mass loss. This highlights the need to focus on improving our understanding of these processes.
Acknowledgements
We thank R.A. Bindschadler, S. Nowicki and others for getting and keeping SeaRISE up and running, J.V. Johnson and others for compiling and maintaining the SeaRISE datasets, and K. Takahashi for carrying out the SeaRISE experiments with IcIES at the Japan Agency for Marine–Earth Science and Technology. The comments of the scientific editor, D.R. MacAyeal, and the reviewers, S.F. Price and J.V. Johnson, helped considerably to improve the clarity of the manuscript. This study was supported by a Grant-in-Aid for Scientific Research A (No. 22244058) from the Japan Society for the Promotion of Science.
