Introduction

Recent observations show an accelerated ice loss of the Greenland ice sheet (Reference VelicognaVelicogna, 2009; Reference Rignot, Velicogna, Van den Broeke, Monaghan and LenaertsRignot and others, 2011) during the past few decades, attributed to both changes in surface mass balance and ice dynamics (Reference Howat, Joughin and ScambosHowat and others, 2007; Reference Rignot, Velicogna, Van den Broeke, Monaghan and LenaertsRignot and others, 2011). Retreat of outlet glaciers is observed all around the periphery of the Greenland ice sheet; however, the acceleration has a high temporal and spatial variability (Reference McFadden, Howat, Joughin, Smith and AhnMcFadden and others, 2011; Reference Moon, Joughin, Smith and HowatMoon and others, 2012). Potential triggering mechanisms include, among others, increased calving rate (Reference Benn, Warren and MottramBenn and others, 2007), intrusion of warm water in the fjords (Reference Holland, Thomas, de Young, Ribergaard and LyberthHolland and others, 2008), enhanced basal lubrication (Reference SchoofSchoof, 2010) and cryo-hydrologic warming of ice (Reference Phillips, Rajaram and SteffenPhillips and others, 2010).

To better understand the behavior of the Greenland ice sheet, several large-scale ice-flow numerical models have been employed to simulate the evolution of the ice sheet over the coming centuries (Reference GreveGreve, 1997a; Reference Ritz, Rommelaere and DumasRitz and others, 2001; Reference Greve, Saito and Abe-OuchiGreve and others, 2011; Reference Gillet-ChauletGillet-Chaulet and others, 2012; Reference Seddik, Greve, Zwinger, Gillet-Chaulet and GagliardiniSeddik and others, 2012). These models rely on different ice-flow approximations, parameterizations of physical processes and initialization procedures. Most ice-flow models rely on simplified first-order approximations of the momentum balance equations, such as the shallow-ice approximation (SIA; Reference HutterHutter, 1982), the shelfy-stream approximation (SSA; Reference MacAyealMacAyeal, 1989; Reference Greve, Saito and Abe-OuchiGreve and others, 2011; Reference Aschwanden, Aðalsgeirsdóttir and KhroulevAschwanden and others, 2012a) or a combination of SIA and SSA (Reference Bueler and BrownBueler and Brown, 2009; Reference Pollard and DeContoPollard and DeConto, 2009). Initialization of these numerical models generally consists of running paleoclimate spin-ups over at least the last glacial cycle, in order to obtain a suitable present-day configuration (Reference Greve, Saito and Abe-OuchiGreve and others, 2011; Reference Aschwanden, Aðalsgeirsdóttir and KhroulevAschwanden and others, 2012a). With this method, long-term memory of the ice-sheet evolution is included in several model variables (e.g. ice temperature). The spatial resolution of these models usually varies between 5 and 40 km, due to intrinsic limitations of SIA and computational requirements. Previous studies (Reference Ritz, Rommelaere and DumasRitz and others, 2001; Reference DeConto and PollardDeConto and Pollard, 2003; Reference Pollard and DeContoPollard and DeConto, 2009; Reference Applegate, Kirchner, Stone, Keller and GreveApplegate and others, 2012; Reference RogozhinaRogozhina and others, 2012) illustrated the ability of these simplified models to improve our understanding of ice-sheet paleoclimate, as well as to evaluate the effect of different model parameters.

One major shortcoming of these models is that they operate at a relatively coarse spatial resolution compared with the size of outlet glaciers that control the ice discharge, which affect the total ice-sheet mass balance. Moreover, these models fail to reproduce observed dynamical changes of Greenland outlet glaciers where the model assumptions break down and where glacier size is comparable to the model spatial resolution (Reference Gillet-ChauletGillet-Chaulet and others, 2012).

In order to better simulate the ice dynamics, several next-generation numerical models have been developed (Reference Price, Payne, Howat and SmithPrice and others, 2011; Reference Gillet-ChauletGillet-Chaulet and others, 2012; Reference Larour, Seroussi, Morlighem and RignotLarour and others, 2012a; Reference Seddik, Greve, Zwinger, Gillet-Chaulet and GagliardiniSeddik and others, 2012). These models include improved three-dimensional (3-D) models (e.g. higher-order (HO; Reference BlatterBlatter, 1995; Reference PattynPattyn, 2003) and full-Stokes (FS) models), and use anisotropic mesh refinement (Reference Morlighem, Rignot, Seroussi, Larour, Ben Dhia and AubryMorlighem and others, 2010) to capture narrow outlet glaciers at the ice-sheet periphery. Contrary to SIA- or SSA-based models, HO and FS models are too computationally intensive to perform interglacial spin-ups. In order to perform realistic simulations and start from a suitable representation of present-day conditions, these models either rely on spinups of simpler models (Reference Seddik, Greve, Zwinger, Gillet-Chaulet and GagliardiniSeddik and others, 2012) or use data assimilation to reproduce present-day conditions (Reference Morlighem, Rignot, Seroussi, Larour, Ben Dhia and AubryMorlighem and others, 2010; Reference Larour, Seroussi, Morlighem and RignotLarour and others, 2012a). The first method captures past changes of the ice sheet, but leads to discrepancies between the thermal regime and the stress regime. Indeed, the spin-up model and the model that is used for future scenarios have different spatial resolutions and rely on different ice-flow approximations, leading to an initial shock that influences the results (Reference Seddik, Greve, Zwinger, Gillet-Chaulet and GagliardiniSeddik and others, 2012). The stress regime, and therefore the viscous and basal heating, is not consistent with the initial temperature. The second method consists of using present-day geometry of the ice sheet (ice thickness, surface elevation, ice front position) and inferring unknown parameters (e.g. basal friction) to reproduce observed velocities with data assimilation (Mac-Ayeal, 1993a; Reference Morlighem, Rignot, Seroussi, Larour, Ben Dhia and AubryMorlighem and others, 2010). As this method relies solely on present-day conditions, past climate variability is not accounted for. Temperatures are recovered using a thermal steady state (Reference Morlighem, Rignot, Seroussi, Larour, Ben Dhia and AubryMorlighem and others, 2010; Reference Larour, Seroussi, Morlighem and RignotLarour and others, 2012a), even though the ice is not in thermal equilibrium. The ice velocity and geometry computed using such a method are consistent with present-day observations; however, inconsistencies between velocity and thickness datasets lead to ice flux divergence anomalies and models start by artificially redistributing the glacier mass, not as a realistic projection but to reconcile the inconsistencies (Reference SeroussiSeroussi and others, 2011).

The ice-sheet thermal regime has a large effect on the ice flow, as it affects basal melting and therefore the sliding of ice over its bed, and it also controls ice deformation. Ice-sheet thermal modeling is difficult due to the uncertainties in basal heat flux (Reference GreveGreve, 2005). Few deep ice cores are available to validate the models (Reference AlleyAlley and others, 1993; Reference DansgaardDansgaard and others, 1993; Greenland Ice Core Project (GRIP) members, 1993; Reference MeeseMeese and others, 1994; Reference Dahl-JensenDahl-Jensen and others, 1998; NorthGRIP members, 2004). These ice cores provide temperature measurements through the ice column, but only in slow-moving areas of the Greenland ice sheet, near ice divides. Initializing ice-flow models with an appropriate thermal regime therefore remains challenging, especially for HO or FS models that cannot run interglacial spin-ups.

Here we assess the effect of the initial thermal regime on century-scale simulations of the Greenland ice sheet. To address this question, we (1) compare steady-state modeled ice temperatures against measurements from three deep ice cores, (2) analyze the sensitivity of inferred properties at the ice/bedrock interface to changes in ice temperature, (3) assess the effect of initial temperatures on century-scale projections and (4) compare the ice mass changes due to different initial temperatures against those obtained by changing external forcings, such as basal sliding or atmospheric conditions. We conclude by discussing the influence of the method employed in calculating the thermal regime on century-scale projections of ice-sheet evolution.

Model and Methods

Model description

To model the Greenland ice sheet, we use the ice-sheet system model (ISSM; Reference Morlighem, Rignot, Seroussi, Larour, Ben Dhia and AubryMorlighem and others, 2010; Reference SeroussiSeroussi and others, 2011; Reference Larour, Seroussi, Morlighem and RignotLarour and others, 2012a), an open-source finite-element software that simulates the thermodynamics of ice sheets at continental scale. Here we summarize the main characteristics and equations of the models (a complete description is given by Reference Larour, Seroussi, Morlighem and RignotLarour and others, 2012a).

Ice is modeled as a viscous incompressible material. We use the 3-D HO model (Reference BlatterBlatter, 1995; Reference PattynPattyn, 2003) for the momentum balance equations. This model is derived from the FS model by making two assumptions: (1) the horizontal gradients of vertical velocities are negligible compared with vertical gradients of horizontal velocity and (2) bridging effects are negligible. This leads to a 3-D model where horizontal and vertical velocities are decoupled. The ice viscosity, μ, is assumed to be isotropic and to follow Glen’s flow law (Reference GlenGlen, 1955):

(1)

where is the effective strain rate, n is Glen’s law coefficient (taken as n = 3) and B is the ice hardness. B is mainly temperature-dependent and we rely on the relationship provided by Reference Cuffey and PatersonCuffey and Paterson (2010, p.75). The horizontal velocity is a solution of

(2)

with

(3)

where (u, v, w) are the three components of velocity in a Cartesian coordinate system (x, y, z), with z the vertical axis, ρ is the ice density, g is the norm of the acceleration due to gravity and s is the ice upper surface elevation. Vertical velocity is recovered using the equation of incompressibility.

The basal friction follows a viscous law (Reference Cuffey and PatersonCuffey and Paterson, 2010); the basal condition between the ice and the underlying bedrock is

(4)

where σ is the stress tensor, n is the outward-pointing normal vector at the base, α ² is the friction coefficient, N is the effective pressure, v is the basal velocity and (·)_∥ is a projection operator on the tangent to the bedrock. Effective pressure is N = ρgh, with h the height of the ice sheet above hydrostatic equilibrium as a first-order approximation; a hydrological model would be necessary to obtain a better estimate of basal pressure. The effective pressure varies with time as the ice-sheet geometry evolves. Other boundary conditions include water pressure at the ice front and air pressure for land-terminating glaciers, non-penetrability of ice into the underlying bedrock at the ice/bedrock interface and a free surface at the ice/air interface.

The evolution of the surface elevation is dictated by mass conservation. The mass transport equation reads

(5)

where H is the ice thickness, is the depth-averaged horizontal velocity, is the surface mass balance and is the basal melting. This hyperbolic equation is stabilized using a streamline upwinding finite-element method (Reference Brooks and HughesBrooks and Hughes, 1982).

For the thermal model, we use an enthalpy formulation (as described in Reference Aschwanden, Bueler, Khroulev and BlatterAschwanden and others, 2012b). This method allows both temperate and cold ice to be included in an energy-conserving framework easily, as there is no need to track the cold/temperate transition surface used in polythermal models (Reference GreveGreve, 1997a,Reference Greveb). The enthalpy equation valid for both cold and temperate ice is

(6)

where E is the enthalpy, E _s is the enthalpy of pure ice, E _l is the enthalpy of pure liquid water at the pressure-melting point, T _pmp, Ki is the ice diffusivity coefficient, C _i is the ice heat capacity, k = (1 − ω)k _i + ωk _w is the mixture thermal conductivity (with ω the water fraction and k _i and k _w the thermal conductivity of pure ice and liquid water), k ₀ is a small positive constant and Φ is the internal deformation heat. The case of pure water, where enthalpy is above the enthalpy of pure liquid water, E _l, is not considered here.

The temperature, T, and water fraction, ω, are recovered as follows:

For cold ice, E < E _s:

For temperate ice, E _s < E < E _l:

Contrary to the description of Reference Aschwanden, Bueler, Khroulev and BlatterAschwanden and others (2012b), our model does not include a layer of the Earth’s lithosphere, so the geothermal heat flux is directly applied at the ice base; ice in contact with the bedrock is heated by geothermal and frictional heat flux (Reference Larour, Seroussi, Morlighem and RignotLarour and others, 2012a). Air temperature is imposed on the upper surface.

Data assimilation of surface velocities is used to infer the basal friction coefficient as we try to minimize the misfit between modeled and measured surface velocities. We define the following cost function, :

(7)

where Γ_s is the domain upper surface, Γ_b its lower surface, (u _obs, v _obs) the observed horizontal velocity and ε is a minimum velocity used to avoid singularities. The first two terms in this cost function are the misfit, which measures the square of the difference between model and observations, and the square of the logarithmic difference between model and observations. The combination of these two terms allows best-fitting of observations on both fast-flowing and slow-moving areas (Reference Morlighem, Rignot, Seroussi, Larour, Ben Dhia and AubryMorlighem and others, 2010, Reference Morlighem, Seroussi, Larour and Rignot2013). The last term is a Tikhonov regularization term, which penalizes uncontrolled oscillations of α and prevents overfitting the data. γ ₁, γ ₂ and γ _t are non-dimensionalizing constants taken as γ ₁ = 300, γ ₂ = 1.5 and γ _t = 10⁻⁷ in order to balance these three terms. The algorithm relies on an exact adjoint that includes dependency of ice viscosity into the strain rate (Reference Goldberg and SergienkoGoldberg and Sergienko, 2011) for the derived gradient and quasi-Newton (limited-memory BFGS (Broyden–Fletcher– Goldfarb–Shanno); Reference NocedalNocedal, 1980) for the gradient descent.

Results

The model initial configuration matches the ice geometry and velocity well, as expected, since it uses present-day topography and is initialized with data assimilation of observed velocities. The difference between observed and modeled velocities is 11.0 m a⁻¹ on average. Figure 1 shows the initial ice thickness (Fig. 1a) and modeled velocity (Fig.1b), as well as the absolute difference between modeled and observed velocities (Fig. 1d). The largest differences are found along the periphery of the ice sheet, where velocities are highest (almost 1000 m a⁻¹ difference on Peterman Gletscher). Areas with high friction (Fig. 3i) are concentrated on mountainous regions, while very low friction is found under fast ice streams, as expected. Additional inversions were performed with different initial guesses of basal friction and led to similar results (not shown here), demonstrating the limited influence of the initial value.

Fig. 1. Initial Greenland configuration: (a) ice thickness, H (Reference Bamber, Layberry and GogineniBamber and others, 2001), (b) modeled horizontal surface velocity, V, (c) basal pressure-adjusted temperature, T, and (d) absolute misfit between observed and modeled horizontal surface velocity overlaid on a MODIS (Moderate Resolution Imaging Spectroradiometer) Mosaic of Greenland.

Fig. 2. Observed (black curves) and modeled steady-state (green curves) temperature profiles at Dye3, GISP2 and GRIP drilling sites. Initial temperatures for EXP1 (pink dot-dashed curves), EXP2 (red dashed curves) and EXP3 (blue dotted curves).

Fig. 3. (a–d) Surface velocities, (e–h) basal velocities and (i–l) friction coefficient, α ², for the four temperature profiles. First column corresponds to steady-state temperature (a, e, i), second column to EXP1 (b, f, j), third column to EXP2 (c, g, k) and fourth column to EXP3 (d, h, l).

We first compare the modeled temperatures (that were computed assuming thermal steady state) with temperature measurements along three deep ice cores (Fig. 2). The modeled steady-state profiles compare well with the measurements. GRIP and GISP2 profiles are similar, as expected, as these two ice cores are located only 28 km apart. For these two profiles, modeled ice temperature in the upper 2000 m is in good agreement with the measurements. Differences up to 5°C are concentrated in the lower part of the ice sheet. Differences for Dye3 are smaller, with a maximum difference of 4°C at the base. Standard errors between observed and modeled temperatures are 1.3°C, 1.8°C and 1.9°C for Dye3, GISP2 and GRIP ice cores, respectively.

The steady-state temperatures can be compared with modeled temperatures from paleoclimate reconstructions of Reference Aschwanden, Aðalsgeirsdóttir and KhroulevAschwanden and others (2012) and Reference RogozhinaRogozhina and others (2012). The pattern of basal temperatures (Fig. 1c) is similar, with ice at the pressure-melting point over a large portion of the Greenland ice sheet and colder ice mainly along the divides and in the north.

The three different initial thermal regimes (EXP1, EXP2 and EXP3) are then used to assess the sensitivity of inferred basal friction to initial thermal regime (Fig. 3). All models are able to accurately reproduce the surface velocities, with an average misfit of 10.7, 13.5 and 11.0 m a⁻¹ for EXP1, EXP2 and EXP3, respectively. Differences between steady state and EXP1 are limited, with slightly slower basal velocities (Fig. 3e and f) and higher friction coefficients in the interior for EXP1 (Fig. 3i and j). This higher friction compensates for vertical shearing due to higher temperatures next to the base of the ice sheet. EXP2 and EXP3 initial temperatures are quite different from the steady-state temperature field so it is not surprising to see that results of the inversion are also different. EXP2 leads to significantly higher friction coefficients and lower basal velocities (Fig. 3g and k), especially in the interior of the ice sheet. Results for EXP3 show the opposite, i.e. lower friction coefficient and higher basal velocities (Fig. 3h and l), but the effect is more limited in this case than for EXP2. Overall, we notice that basal friction and sliding seem to be more sensitive to initial temperature conditions in the interior of the ice sheet than in the periphery.

We then perform four runs of 100 years using these different initial thermal regimes. Under constant present-day conditions and with the initial steady-state temperature, the model yields an increase in ice mass in the next 100 years of 180 Gt a⁻¹, equivalent to an increase of 0.73% of the initial ice mass of the ice sheet (Fig. 4). Results of Reference BindschadlerBindschadler and others (2013) show that mass increase is a limitation common to about half of the ice-sheet models that participated in the SeaRISE experiments, the other half starting from an initial volume that differs from the present-day volume. Differences between the steady state and EXP1 are limited, with a similar mass gain (0.73% and 0.72%, respectively). EXP3 mass change is also close to the steady-state evolution, with a mass gain of 0.75%. Results from EXP2 are different from the other three, with a lower mass increase of only 0.62% in 100 years. These results are further analyzed below.

Fig. 4. Simulated ice mass change of the Greenland ice sheet (in Gt and %) for 100 year simulations with different scenarios: constant climate with steady-state temperature (green), EXP1 temperature (violet), EXP2 temperature (red), EXP3 temperature (dark blue), AR4 atmospheric conditions (yellow) and increased basal lubrication (gray).

Finally, we run scenarios similar to the SeaRISE initiative (Reference BindschadlerBindschadler and others, 2013; Reference NowickiNowicki and others, 2013) with the steady-state temperature field as initial state. Ice mass change for the climate and sliding scenarios is also presented in Figure 4. These two scenarios respectively lead to a mass gain of 0.50% and 0.36%, which is significantly lower than in all the previous runs. Our 100 year simulations are therefore significantly more sensitive to the applied external forcings than the initial temperature.

Discussion

Large-scale high-resolution simulations require intensive computational resources. In this study, we use a 3-D HO model, as a trade-off between first-order SSA/SIA and FS. Our model comprises 1500 000 elements to capture outlet glaciers while running century-scale simulations in a reasonable amount of time (i.e. ∼4 days for each run of 100 years). Previous studies (Reference HindmarshHindmarsh, 2004; Reference GudmundssonGudmundsson, 2008; Reference Morlighem, Rignot, Seroussi, Larour, Ben Dhia and AubryMorlighem and others, 2010) on both idealized geometries and real glaciers showed that differences between HO and FS exist, especially in the grounding line area (Reference Durand, Gagliardini, Zwinger, Le Meur and HindmarshDurand and others, 2009; Reference Morlighem, Rignot, Seroussi, Larour, Ben Dhia and AubryMorlighem and others, 2010) but remain spatially limited, so using HO for our simulations should not affect our conclusions.

The numerical experiments performed here show that the assumption of thermal steady state is a viable approximation for short-term simulations of the Greenland ice sheet. First, modeled steady-state temperatures match the observations along the three deep ice cores reasonably well. The main features of the temperature profiles (e.g. surface temperature, basal temperature, basal temperature gradient) are captured, with a standard error varying between 1.3 and 1.9°C, even though differences with deep ice cores reach 5°C near the glacier base. These errors correspond to an average ice rigidity difference of 3.5% for Dye3 and 5% for GRIP and GISP. In the upper part of the ice, the modeled temperatures reach a constant value, which is not consistent with the observations. This is due to the air temperature history that is not captured here, as opposed to Reference GreveGreve (2005) and Reference RogozhinaRogozhina and others (2012): the temperature decrease caused by the Last Glacial Maximum is not captured by our steady-state method. Paleoclimate reconstructions, which account for temporal changes in surface conditions, yield a better agreement between modeled temperature and measurements along deep ice core (e.g. <3°C in Reference RogozhinaRogozhina and others, 2012). It is important to note that these deep ice-core drillings were collected near ice divides, where conduction is dominant, and might not be representative of areas where advection is dominant. Comparisons in more dynamic areas would better inform the steady-state model’s ability to capture variations, but measurements in these regions are almost non-existent.

Second, for all four temperature fields the basal drag inversions yield a good agreement between modeled and observed velocities. The temperature field affects the friction coefficient pattern, especially with the extreme cold and warm cases (EXP2 and EXP3), but the basal stress remains similar in all four cases. This result corroborates the conclusions of Reference SchäferSchäfer and others (2012), who showed that the friction distribution of Vestfonna ice cap, Nordaustlandet/Svalbard, does not depend significantly on ice temperature. The surface velocity is optimized to fit observations and should not vary from one simulation to the next. Warmer ice deforms more easily and leads to more vertical shear; this larger vertical shear is compensated for by a lower basal velocity in order to obtain the same surface velocity. This effect is the foundation of the SIA, that is valid for most ice-sheet interiors, but does not hold for fast-flowing outlet glaciers, where sliding dominates ice motion.

The pattern of inferred basal friction based on steady-state temperatures is consistent with the results of Reference Larour, Seroussi, Morlighem and RignotLarour and others (2012a) and Reference Gillet-ChauletGillet-Chaulet and others (2012) for Greenland, and Reference Joughin, Fahnestock, MacAyeal, Bamber and GogineniJoughin and others (2001) and Reference SeroussiSeroussi and others (2011) for the northeast Greenland ice stream, with motion mainly due to sliding on fast-flowing areas and differences from observations concentrated on the fast-flowing outlet glaciers. These models use different approximations of FS and different temperature fields (thermal steady state for Reference Larour, Seroussi, Morlighem and RignotLarour and others, 2012a; an interpolation from a paleoclimate reconstruction for Reference Gillet-ChauletGillet-Chaulet and others, 2012) and confirm that ice temperatures have little impact on the basal friction inferred from surface velocities.

Finally, the projection of the ice-sheet mass balance over the next 100 years shows an increase in mass of ∼0.73%. This result is not consistent with the current trends that show a rapid and accelerated decrease of the Greenland ice sheet mass (e.g. Reference Rignot, Velicogna, Van den Broeke, Monaghan and LenaertsRignot and others, 2011). However, a mass increase is not unusual for ice-sheet models of the Greenland ice sheet. About half the models that participated in the Greenland SeaRISE experiments show a mass increase for the control run that keeps present-day conditions constant (Reference BindschadlerBindschadler and others, 2013), although they all use different ice-flow approximations, spatial resolutions and initialization procedures. The other models show a mass decrease but start from an initial volume that differs from the actual volume. Mass increase often happens for models initialized with inverse methods and is caused by data that have inconsistencies (Reference SeroussiSeroussi and others, 2011), incorrect position of the ice-sheet margin, mesh resolutions too coarse (Reference Gillet-ChauletGillet-Chaulet and others, 2012) and/or missing physical processes, such as ice/ocean interaction and basal hydrology. Models by Reference Seddik, Greve, Zwinger, Gillet-Chaulet and GagliardiniSeddik and others (2012) and Reference Gillet-ChauletGillet-Chaulet and others (2012) show a volume increase of 0.81% and 0.28% in 100 years, respectively, after an initial relaxation of 100 and 50 years, respectively. Such a relaxation allows the reconciliation of inconsistent model inputs that lead to unphysical flux divergences (Reference SeroussiSeroussi and others, 2011) and limit the high rates of thickness change. However, no ice-sheet model is able to replicate the present-day configuration and evolution today. We did not use such a relaxation technique here, which could lead to slightly different initial conditions (ice thickness, initial volume, driving stress, etc.), in order to ensure that the ice-sheet response is due only to changes in the ice thermal regime.

For century-scale projections, our experiments show that the initial thermal regime has only a limited effect on the evolution of the Greenland ice sheet. Ice temperature of models that do not start with a thermal steady state see their temperatures gradually adjusted, but this is a slow process that has a limited effect for short-term simulations. Influence of geothermal heat flux will have a similar impact. Reference Larour, Morlighem, Seroussi, Schiermeier and RignotLarour and others (2012b) showed that large variations of geothermal flux affected the ice rigidity by up to 5% in slow-moving areas, where basal friction and geothermal energy are the same order of magnitude. This is within the range of the initial temperature variations that we tested here. We therefore expect that uncertainties in geothermal heat flux do not significantly affect the behavior of the ice sheet for century-scale simulations. This is consistent with the results of Reference Larour, Morlighem, Seroussi, Schiermeier and RignotLarour and others (2012b), who concluded that uncertainties in geothermal heat flux affect glacier mass balance by <1%. This is not the case for long interglacial runs, where temperature plays a major role in ice-sheet growth; Reference GreveGreve (2005) and Reference RogozhinaRogozhina and others (2012) noticed that small variations in geothermal flux led to differences in ice volume of up to 10%. By contrast, changes in external forcings, such as atmospheric conditions or basal friction, lead to significant variations in ice-sheet mass over 100 years, despite the moderate forcings applied in our forcing experiments. Results from the SeaRISE experiments (Reference BindschadlerBindschadler and others, 2013) and Reference Gillet-ChauletGillet-Chaulet and others (2012) also show the strong influence of these forcings over short timescales.

As mentioned in the previous section, mass increase is higher for EXP3 than for the steady-state experiment and lower for EXP1 and EXP2. We would expect the opposite: for identical surface velocities, depth-averaged velocities should be smaller for warmer ice, as there is more vertical shear and basal velocity from the inversion shows a lower basal velocity for EXP2 than for the steady-state case. Ice discharge is mainly controlled by fast-flowing outlet glaciers, and as motion due to vertical shear only represents a fraction of the ice velocity in these areas, ice flux is not affected. Conversely, differences between observed and modeled velocities have a stronger effect. Contrary to EXP1 and EXP3, Greenland mass evolution with EXP2 is quite different from the steady-state case. This is due to the initial velocity after inversion, which is higher than observed, especially in the eastern and southeastern part of Greenland (Fig. 5a), where the difference reaches 50 m a⁻¹ over large areas. In all other cases, the modeled velocity is usually smaller than the measurements (e.g. EXP3 in Fig. 5b: the difference between observed and modeled velocities is ∼0 m a⁻¹ in those same areas). This higher velocity leads to a discharge at the ice front that is 10% higher for EXP2 than for the steady state and EXP1 and EXP3 cases. The inversion for EXP2 did not converge as well as the others, even though the absolute misfit is only slightly higher for this experiment than for the others (13.5 m a⁻¹ for EXP2 vs 10.7 and 11.0 m a⁻¹ for EXP1 and EXP3). This suggests that the model is highly sensitive to the initial modeled velocity and that this parameter has a significantly larger impact on simulations than changes in initial temperatures.

Fig. 5. Modeled minus observed velocities for (a) EXP2 and (b) EXP3. Modeled velocities are higher than observations in red areas and lower in blue areas.

The basal friction law employed here is not temperature-dependent, and only the effective pressure changes as the ice-sheet geometry evolves. Projections of the Greenland ice sheet evolution could be more sensitive to the initial temperature field if a temperature-dependent friction law were used, as nonlinear effects could develop. However, most of the ice sheet is already at the pressure-melting point (Fig. 1c), which limits the impact of such a law, and our conclusions should not be significantly affected. A basal hydrological model that models effective pressure and water content (de Fleurian, 2013) is under development and will improve this aspect in the future.

Realistic modeling of the ice thermal regime remains a major objective and more in situ measurements need to be acquired to calibrate and validate models. There are few temperature measurements available in fast-flowing areas, and uncertainties associated with the geothermal heat flux, which is a key control on the thermal regime, remain large. Our results indicate that variations of <5°C in parts of the ice sheet in the initial temperature field have a limited effect on both initializations with data assimilation and ice-sheet evolution for short-term simulations of large-scale models. Other factors and parameters are of greater importance: an accurate representation of small outlet glaciers is necessary to capture the ice discharge; a realistic description of the ice-sheet present-day configuration with improved bedrock topography data is needed to reduce uncertainties, as well as better constrain mass flux and unknown parameters (Reference SeroussiSeroussi and others, 2011; Reference Gillet-ChauletGillet-Chaulet and others, 2012). These aspects are more limiting for ice-sheet numerical modeling and affect sea-level rise projections at the continental scale more significantly than small variations in the ice thermal regime.

Conclusions

In this study, we have presented simulations performed with a 3-D thermomechanical HO model of the Greenland ice sheet. Anisotropic mesh adaptation and data assimilation allow the numerical model to capture outlet glaciers with better accuracy and to match the present-day conditions in a realistic manner. We show that modeled temperatures, assuming an ice sheet in thermal steady state, are in agreement with in situ measurements along three deep ice-core profiles. The impact on ice rigidity is 5% or less, and such changes do not have a significant influence on basal friction inversion and ice-sheet evolution for century-scale projections. The ice-sheet model is far more sensitive to changes in external forcings, namely basal sliding and atmospheric conditions, than its initial temperature. Using a steady-state approach combined with present-day conditions to compute the ice thermal regime is a viable alternative to spin-up when looking at short-term simulations of the Greenland ice sheet at the continental scale.