Articles | Volume 18, issue 1
Research article
26 Jan 2024
Research article |  | 26 Jan 2024

Evaluating different geothermal heat-flow maps as basal boundary conditions during spin-up of the Greenland ice sheet

Tong Zhang, William Colgan, Agnes Wansing, Anja Løkkegaard, Gunter Leguy, William H. Lipscomb, and Cunde Xiao

There is currently poor scientific agreement on whether the ice–bed interface is frozen or thawed beneath approximately one third of the Greenland ice sheet. This disagreement in basal thermal state results, at least partly, from differences in the subglacial geothermal heat-flow basal boundary condition used in different ice-flow models. Here, we employ seven widely used Greenland geothermal heat-flow maps in 10 000-year spin-ups of the Community Ice Sheet Model (CISM). We perform two spin-ups: one nudged toward thickness observations and the other unconstrained. Across the seven heat-flow maps, and regardless of unconstrained or nudged spin-up, the spread in basal ice temperatures exceeds 10 C over large areas of the ice–bed interface. For a given heat-flow map, the thawed-bed ice-sheet area is consistently larger under unconstrained spin-ups than nudged spin-ups. Under the unconstrained spin-up, thawed-bed area ranges from 33.5 % to 60.0 % across the seven heat-flow maps. Perhaps counterintuitively, the highest iceberg calving fluxes are associated with the lowest heat flows (and vice versa) for both unconstrained and nudged spin-ups. These results highlight the direct, and non-trivial, influence of the heat-flow boundary condition on the simulated equilibrium thermal state of the ice sheet. We suggest that future ice-flow model intercomparisons should employ a range of basal heat-flow maps, and limit direct intercomparisons with simulations using a common heat-flow map.

1 Introduction

There is presently a tremendous diversity of opinion regarding the geothermal heat flow beneath the Greenland ice sheet due to a paucity of direct heat-flow measurements in the ice-sheet interior. While many subaerial, submarine, and shallow subglacial measurements have been made around the ice-sheet periphery, deep subglacial measurements have only been made at six deep ice coring sites in the interior (Camp Century, DYE-3, GRIP, GISP2, NGRIP, and NEEM). Consequently, the magnitude and spatial distribution of Greenland's subglacial geothermal heat flow remains poorly constrained across the seven unique Greenland heat-flow models presently in widespread use (Fig. 1) (Shapiro and Ritzwoller, 2004; Rezvanbehbahani et al., 2017; Martos et al., 2018; Greve, 2019; Lucazeau, 2019; Artemieva, 2019; Colgan et al., 2022). These individual geothermal heat-flow models are derived from multiple techniques that interpret a variety of geophysical variables (Table 1). We briefly discuss differences in the methodology and geophysical input variables of these existing heat-flow maps.

Figure 1(a–g) The seven geothermal heat-flow maps considered as basal thermal boundary conditions. Color bars saturate about 10 and 100 mW m−2. (h) Ensemble mean. Units for all plots, mW m−2.

Table 1Characteristics of the seven geothermal heat-flow models we explore as basal thermal boundary conditions: methodology used to derive each model, number of geophysical datasets employed by each model, number of in situ heat-flow observations considered by each model, average heat flow (± standard deviation) within a common CISM Greenland ice-sheet area, and the domain coverage of each model. Adapted from Colgan et al. (2022) and arranged from lowest to highest average geothermal heat flow beneath the ice sheet.

Download Print Version | Download XLSX

The Rezvanbehbahani et al. (2017) (R2017), Lucazeau (2019) (L2019), and Colgan et al. (2022) (C2022) heat-flow maps are perhaps methodologically most similar. These three maps use machine learning or geostatistics to predict heat flow as a function of diverse geophysical variables such as topography, tectonic setting, observed gravity, and magnetic field. They differ not only in the applied method but also in the set of geophysical variables utilized and their domains. Whereas Rezvanbehbahani et al. (2017) and Lucazeau (2019) only used global data, Colgan et al. (2022) substituted global datasets with Greenland-specific local data. By contrast, the Shapiro and Ritzwoller (2004) (SR2004), Martos et al. (2018) (M2018), and Artemieva (2019) (A2019) heat-flow maps all employ lithospheric models of varying complexity and more specific geophysical variables to infer heat flow. Shapiro and Ritzwoller (2004) correlate the seismic shear wave velocities of the upper 300 km with heat-flow observations and use this connection to predict heat flow from tomography data in areas without heat flow observations. Martos et al. (2018) use magnetic data to infer the Curie temperature depth, and the steady-state one-dimensional heat conduction equation, with lateral constant values for thermal conductivity and radiogenic heat production. Artemieva (2019) assumes an isostatic equilibrium and translates the corresponding topographic residuals to temperature anomalies which are then converted to a lithosphere–asthenosphere boundary undulation, and then uses individual reference geotherms for the different tectonic settings to derive the geothermal heat flow from lithosphere–asthenosphere boundary topography. Both methods then infer heat flow from the respective isotherms by applying a thermal model. The Greve (2019) (G2019) heat-flow map is unique in using paleoclimatic forcing of an ice-flow model to infer heat flow with a minimum of geophysical variables.

In North Greenland, there is especially poor agreement among the present generation of geothermal heat-flow models. Some models infer a widespread North Greenland high heat-flow anomaly (e.g., Greve, 2019), and some do not (e.g., Lucazeau, 2019). Other models offer products with and without this high heat-flow anomaly (e.g., Rezvanbehbahani et al., 2017). There are numerous secondary disagreements as well, including whether a model (1) infers traces of the Iceland Hotspot Track transiting from West to East Greenland (Martos et al., 2018), (2) infers elevated heat flow in East Greenland in closer proximity to the Mid-Atlantic Ridge (Artemieva, 2019), or (3) infers a low heat-flow anomaly associated with the North Atlantic Craton in South Greenland (Colgan et al., 2022).

Geothermal heat flow is a critical basal thermal boundary condition in Greenland ice-sheet models. It can significantly influence basal ice temperature and rheology, which in turn influences basal meltwater production and friction (Karlsson et al., 2021). Given the non-linear relation between ice temperature and rheology, and that most ice deformation occurs in the deepest ice layers, relatively small changes in basal ice temperature can result in large changes in ice velocity (Hooke, 2019). In extreme cases, diminished geothermal heat flow along subglacial ridges may contribute to the formation of massive refrozen basal ice masses (Colgan et al., 2021), or sharply enhanced geothermal heat flow may contribute to the onset of major ice-flow features (Smith-Johnsen et al., 2020).

Despite the clear links between geothermal heat flow and ice dynamics, a standardized geothermal heat flow as the basal thermal boundary condition was not prescribed in the Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6) (Goelzer et al., 2020). Of the 21 Greenland model submissions in ISMIP6, 12 prescribed geothermal heat flow according to Shapiro and Ritzwoller (2004), five prescribed it according to Greve (2019), two prescribed it as a hybrid assimilation of four older geothermal heat-flow models (Pollack et al., 1993; Tarasov and Peltier, 2003; Fox Maule et al., 2009; Rogozhina et al., 2016), and one prescribed a spatially uniform geothermal heat flow.

For Greenland, the ISMIP6 ensemble suggests that  40 % of the ice-sheet bed is frozen, meaning basal ice temperatures below the pressure-melting-point temperature, and  33 % of the ice-sheet bed is thawed, meaning basal ice temperatures at the pressure-melting-point (MacGregor et al., 2022). The ISMIP6 ensemble disagrees on whether the basal thermal state is frozen or thawed beneath the remaining  28 % of the ice sheet. It is unclear what portion of this disagreement is associated with the use of differing geothermal heat-flow boundary conditions across ensemble members. The geothermal heat-flow boundary condition can significantly influence the basal ice temperature and thus change the ice-flow rheology. For example, basal ice that is 1 C below pressure-melting-point temperature Tpmp deforms approximately 10 times more than ice 10 C below Tpmp at the same driving stress (Hooke, 2019).

In preparation for ISMIP7, there is a clear motivation to explore the choice of geothermal heat-flow boundary condition on modeled basal ice temperatures. Here, we spin up an ice-flow model with seven different geothermal heat-flow boundary conditions. This allows us to isolate the influence of the geothermal heat-flow boundary condition on the simulated thermal state and ice flow. We also discuss the pros and cons of these seven Greenland geothermal heat-flow products in the specific context of utility for future Greenland ice-flow simulations.

2 Methods

We use the Community Ice Sheet Model (CISM) (Lipscomb et al., 2019) as configured to spin up the Greenland ice sheet for ISMIP6 simulations (Goelzer et al., 2020). We run CISM on a regular 4 km grid with 10 vertical layers, using a higher-order velocity solver with a depth-integrated viscosity approximation (DIVA) based on Goldberg (2011). Note that the higher-order solver includes both membrane and vertical shear stress. All floating ice is assumed to calve immediately; thus we do not simulate Greenland's small floating ice shelves and ice tongues. For partly grounded cells at the marine margin, basal shear stress is proportionally weighted to the grounded fraction of the cell using a sub-grid grounding-line parameterization (Leguy et al., 2021).

We perform two types of ice-sheet spin-ups that we denote Case 1 and Case 2. The Case 1 spin-up iteratively nudges the friction coefficients in the basal-sliding power law to minimize misfit against observed present-day ice thickness. The nudging method is similar to that of Pollard and DeConto (2012) and was applied to the Antarctic ice sheet by Lipscomb et al. (2021). In this spin-up, we use a classic Weertman-type non-linear basal friction law (Weertman, 1979):

(1) τ b = C u b 1 / m - 1 u b ,

where τb is the basal traction, ub is the basal velocity, and m is a dimensionless constant that we set to 3. C is the spatially varying friction coefficient, in units of Pa yr m−1, that is nudged during spin-up at each basal velocity point. The nudged C is capped at a maximum value of 105 (implying high resistance to basal sliding) and a minimum of 10 (implying little resistance to sliding).

The spun-up ice thickness, by design, is close to observations. In most of the ice sheet, the thickness and velocity fields are in approximate balance, and thus the spun-up ice velocity also agrees well with observations, even though velocity is not a nudging target. The exceptions are regions where ice velocity has recently changed and ice thickness has not had time to adjust. The main drawback of the Case 1 spin-up method is that there is no dependence of basal sliding on basal temperature or water pressure. Thus, the method is not very physically based, and arguably overfits the ice thickness observations.

By contrast, the Case 2 spin-up is unconstrained, meaning that basal friction coefficients are not nudged to match the observed present-day ice thickness. We use a pseudo-plastic sliding law (Aschwanden et al., 2016):

(2) τ b = - τ c u b | u b | 1 - q u 0 q

where τc is the transient yield stress in Pa, q=0.5 is a dimensionless pseudo-plastic exponent, and u0=100 m yr−1 is a threshold speed. The yield stress is computed as τc=Ntan, where N is the effective pressure and is a friction angle. The friction angle varies linearly as a function of bed elevation b between 40 at b=700 m and 5 at b=-700 m. The effective pressure is computed from a local till model (Bueler and Van Pelt, 2015) and is sensitive to the thermal state of the bed; N is equal to the full overburden pressure ρigH (where ρi is the ice density, g is gravitational acceleration, and H is ice thickness) when the bed is frozen, but decreases to a small fraction (0.02) of overburden on thawed beds as the basal water depth rises to a capped value of 2 m. Lipscomb et al. (2019) provide more details.

While the Case 1 spin-up ice geometry closely matches present-day observations, there can be appreciable ice thickness biases for the non-nudged Case 2 spin-up. For ISMIP6, most participating models used data assimilation or nudged spin-ups to obtain a more accurate initial state. It is foreseeable, however, that future ISMIP protocols will encourage unconstrained spin-ups as a complement to nudged spin-ups, especially for simulations over multiple centuries during which basal conditions are likely to evolve. Unconstrained spin-ups are more physically based than nudged spin-ups in that the basal shear stress is closely tied to the modeled bed state (e.g., basal temperature, geology, and/or hydrology). It is more challenging, however, to reproduce a specific (i.e., present-day) ice-sheet configuration in an unconstrained spin-up.

For both Case 1 and 2 spin-ups, the ice sheet was initialized with present-day ice thickness and bed topography (Morlighem et al., 2017). The surface mass balance (SMB) and surface air temperature (Tair) are prescribed from a 1980–1999 climatology provided by the Modèle Atmosphérique Régional, a regional climate model (Fettweis et al., 2017). The initial englacial temperature was initialized to an idealized vertical profile. Where the prescribed SMB is negative, the initial temperature profile in each column is linear, with T= min(Tair, 0) at the surface and T=Tpmp-5 at the bed. Where the SMB is positive, the temperature is initialized to an analytic profile based on a balance between vertical conduction and cold advection (Cuffey and Paterson, 2010, Sect. 9.5.1). The ice sheet was then spun up for 10 000 years with a time step of 1/6 year. The englacial temperature evolves under vertical conduction, horizontal and vertical advection, and deformational heating. Where the bed is frozen (Tb<Tpmp), the basal temperature is computed by prescribing a balance of geothermal heat flow, vertical conductive flux, and frictional fluxes at the ice–bed interface. Where the bed is thawed, meaning the resulting temperature would exceed Tpmp, we set Tb=Tpmp and use the excess energy to melt ice. By the end of the spin-up, the ice sheet is close to equilibrium (the relative change rate of ice volume is below 0.001 %), with transient englacial ice temperatures no longer influenced by the initial temperature profile.

We repeat the Case 1 and Case 2 spin-ups seven times each with the same configuration and execution, only varying the prescribed geothermal heat flow serving as the basal boundary condition (Table 1, Fig. 1). Each of the seven heat-flow maps is re-gridded from its native grid to the CISM grid using bilinear interpolation. For heat-flow maps that are only available onshore, meaning they omit offshore, or submarine, areas of the CISM domain, we similarly infill fjord heat-flow values using bilinear interpolation. These seven maps provide a diverse representation of the magnitude and spatial distribution of Greenland heat flow, with the mean heat flow within the CISM ice-sheet domain ranging from  42 mW m−2 in the C2022 map to  64 mW m−2 in the L2019 map. For R2017 we use the middle range scenario of NGRIP = 135 mW m−2. For A2019, we use the “model 1” scenario, which adopts a deeper continental Moho depth than the “model 2” scenario. For C2022 we use their recommended “without NGRIP” scenario.

Of the seven heat-flow maps, only two are global maps (Shapiro and Ritzwoller, 2004; Lucazeau, 2019), and the remaining five are Greenland-specific. Of these five Greenland-specific maps, all but C2022 are limited to the onshore domain. The seven maps are evaluated against differing numbers of in situ heat-flow observations within a Greenland domain defined as <500 km from Greenlandic shores. The R2017, M2018, and G2019 heat-flow maps employed ≤9 primarily subglacial in situ observations from deep boreholes in the ice-sheet interior. The remaining four maps used significantly more in situ heat-flow observations (≥278), including more subaerial, submarine, and shallow subglacial measurements, associated with progressively improving versions of the International Heat Flow Database (Jessop et al., 1976; Fuchs et al., 2021).

Figure 2Case 1: (a–g) Ice-bed temperature relative to pressure melting point at transient equilibrium using the seven geothermal heat-flow maps. (i) Ensemble mean ice-bed temperature. Units in all plots C below pressure-melting-point temperature. (Compare against Case 2 in Fig. 9).

Figure 3Case 1: (a–g) Relative anomaly from ensemble mean in ice-bed temperature relative to pressure melting point at transient equilibrium using the seven geothermal heat-flow maps. (h) Ensemble mean ice-bed temperature relative to pressure melting point. Units in all plots are C.

3 Results

3.1 Case 1 spin-up

Figure 2 shows the ice-bed temperature Tb relative to Tpmp at the end of each Case 1 spin-up. The C2022 heat-flow map, which has the lowest mean geothermal heat-flow of all seven products, yields the smallest area of thawed basal temperatures (21.8 %) and the lowest basal temperature anomaly relative to the ensemble mean (Fig. 3; Table 2). Conversely, the relatively high M2018 heat-flow map, which has the third highest mean heat flow of all seven products, yields twice the area of thawed basal temperatures (54.4 %) and one of the highest basal temperature anomalies relative to the ensemble mean. Across the seven-member ensemble, however, there is considerable variation in the magnitude and spatial distribution of the ensemble spread in basal ice temperatures (Fig. 4). The seven heat-flow maps yield broadly similar modeled basal ice temperature RMSEs of between 1.0 and 2.8 C in comparison to observed basal ice temperatures at 27 Greenland ice-sheet boreholes (Fig. 5) (Løkkegaard et al., 2023).

Figure 4(a) and (b) Ensemble agreement in basal thermal state (frozen or thawed) across the seven heat-flow maps (a: Case 1, b: Case 2). Units are the fraction of simulations that suggest thawed bed. (c) and (d) Ensemble spread (the difference between maximum and minimum values for different experiments) in basal ice temperature across the seven heat-flow maps (c: Case 1, d: Case 2). Units are C.

Figure 5Modeled ice-bed temperature across the seven heat-flow maps versus observed ice-bed temperature at 27 Greenland ice-sheet boreholes where ice temperatures have been observed. (a–g) Modeled versus observed comparison across the seven geothermal heat-flow maps. The numbers show the deep ice core locations in (h). Case 1 spin-ups shown in blue. Case 2 spin-ups shown in red.

Table 2Thawed-bedded ice-sheet area associated with Case 1 (nudged) and Case 2 (unconstrained) spin-ups of 10 000-year duration for the seven geothermal heat-flow datasets.

Download Print Version | Download XLSX

Generally, the ensemble spread in modeled ice-bed temperature approaches zero in the ablation area, especially in Central West Greenland, where the basal thermal state is thawed regardless of the heat-flow map. The ensemble spread is generally largest along the main flow divide of the ice sheet. At South Dome, the ensemble spread exceeds 10 C over an area of  105 km2. This highlights that the heat-flow map has a substantial influence on the simulated basal thermal state over the North Atlantic Craton. While the Northeast Greenland Ice Stream is thawed regardless of the heat-flow map, there is also an area of  105 km2 in Central East Greenland where the ensemble spread exceeds 10 C. Finally, the choice of heat-flow map appears to influence whether the North Greenland ablation area is thawed or frozen.

The Case 1 spin-up nudges the ice-flow model toward present-day ice thickness by iteratively adjusting basal friction coefficients C at each basal velocity point. The ensemble differences in C generally reach a maximum where ice velocities reach a minimum (Fig. 6). Perhaps counterintuitively, the highest surface ice velocities are associated with the lowest geothermal heat-flows (Fig. 7). For example, the high and low heat-flow end members of the L2019 and C2022 maps yield, respectively, low and high ice-velocity end members. Similarly, within the R2017 simulation, the low heat-flow anomaly in southeast Greenland yields a high ice-velocity anomaly. Accordingly, iceberg calving is highest in the lowest heat-flow simulations (Fig. 8). The relatively narrow ensemble spread in iceberg calving ( 1 %; 2 Gt yr−1 ensemble range against 322 Gt yr−1 ensemble mean) is ultimately constrained by the surface mass balance forcing at transient equilibrium when the ice sheet is adjusting during the spin-ups.

Figure 6Case 1: (a–g) The basal friction coefficient at transient equilibrium using the seven geothermal heat-flow maps, expressed as anomalies from the ensemble mean. Units are % and color bars saturate at ±150 %. (h) Ensemble mean basal friction coefficient at transient equilibrium. Units are Pa yr m−1, with the color bar saturating at 106 Pa yr m−1.

Figure 7Case 1: (a–g) Surface ice velocity at transient equilibrium using the seven geothermal heat-flow maps, expressed as anomalies from their ensemble mean. Units are % and color bars saturate at ±150 %. (h) Ensemble mean surface ice velocity at transient equilibrium. Units are m yr−1.

Figure 8Total Greenland ice-sheet calving flux over the 10 000-year spin-up using the seven geothermal heat-flow maps for Case 1 (a) and Case 2 (b). Units are Gt yr−1. The first 500 years of the simulations are not shown due to artifacts associated with model initialization.


3.2 Case 2 spin-up

Similar to the Case 1 spin-up, the Case 2 spin-up yields the smallest area of thawed basal temperatures (33.5 %) with the C2022 lowest mean geothermal heat-flow map and the largest area of thawed basal temperatures (60.0 %) with the M2018 relatively high mean geothermal heat-flow map (Fig. 9). Critically, the thawed-bed area for a given heat-flow map is consistently larger under the Case 2 (unconstrained) spin-up than the Case 1 (nudged) spin-up (Table 2). Basal ice temperatures are accordingly warmer under Case 2 than Case 1 (Fig. 10). As ice-sheet sensitivity generally increases with the thawed-bed area over which basal movement and subglacial hydrology can occur, this may suggest that unconstrained ice-sheet spin-ups are more sensitive than nudged ones. The apparent ice-temperature warming effect of an unconstrained spin-up appears to increase with decreasing heat flow. The shift toward warmer basal temperatures under Case 2 is most apparent for the C2022 heat-flow map, where the temperature difference is >5C beneath a large portion of Central Greenland. All heat-flow maps yield large differences in basal ice temperature between Case 1 and Case 2 spin-ups in regions of fast ice flow around the ice-sheet periphery.

Figure 9Case 2: (a–g) Ice-bed temperature relative to pressure melting point at transient equilibrium using the seven geothermal heat-flow maps. (h) Ensemble mean ice-bed temperature. Units in all plots C below pressure-melting-point temperature.

Figure 10Case 2: (a–g) The difference of ice-bed temperature at transient equilibrium between Case 1 and 2 (TCase2TCase1) using the seven geothermal heat-flow maps. (h) Ensemble mean for (a)(g). Units in all plots are C.

The spatial pattern of the Case 2 ensemble agreement broadly follows that of Case 1 (Fig. 4), although the Case 2 agreement is generally poorer. This is attributable to the unconstrained nature of the Case 2 spin-up. The magnitude and spatial distribution of the ensemble spread in basal ice temperatures under Case 2 is largely similar to that of Case 1. The Case 2 ensemble spread is smaller in Central East Greenland and larger for peripheral ice caps, especially Flade Isblink in Northeast Greenland (Fig. 4). The Case 2 spin-up reproduces the observed basal ice temperatures at 27 Greenland ice-sheet boreholes with an RMSE of between 1.5 and 2.8 C (Fig. 5) (Løkkegaard et al., 2023). This is not significantly different from the RMSE range of the Case 1 spin-up. Basal ice temperatures are better resolved by the Case 1 spin-up for three heat-flow maps (C2022, G2019 and R2017), and better resolved by Case 2 for two maps (A2019 and L2019) with the remaining two maps yielding the same RMSE under both spin-ups (SR2004 and M2018). Empirical temperature observations therefore justify neither the Case 1 nor Case 2 spin-up approach.

In comparison with Case 1, the Case 2 spin-ups generally result in thicker ice in East Greenland and thinner ice in West Greenland (Fig. 11). These substantial differences in ice thickness (i.e., ±100 m) are clearly attributable to the unconstrained nature of Case 2 spin-up. Case 2 spin-ups with different heat-flow maps can yield very different ice thicknesses. For example, the SR2004 and C2022 maps yield substantially thicker than observed ice in North Greenland, where the G2019 and L2019 maps yield substantially thinner than observed ice. Similarly, the ice thickness at South Dome (South Greenland) varies considerably across the seven simulations. The magnitude of ice thickness differences associated with heat-flow maps is non-trivial, and the spatial distribution is complex.

Figure 11Case 2: (a–g) Anomaly in ice thickness at Case 2 transient spin-up, in comparison with Case 1 nudged spin-up, using the seven geothermal heat flow maps. Units in all plots are m and expressed as Case 2 minus Case 1. (h) Ensemble mean of ice thickness anomaly. The color bars saturate at ±150 m.

There are considerable velocity differences across the seven Case 2 spin-up simulations. Generally, these velocity differences appear to be negatively correlated with ice thickness differences (Fig. 12). For example, the SR2004 and C2022 heat-flow maps that yield substantially thicker ice in North Greenland also yield lower ice temperatures there. Similarly, the G2019 and L2019 maps that yield substantially thinner ice in North Greenland also yield faster velocities there. While relative velocity differences in the ice-sheet interior can appear striking in both magnitude and extent, there are also velocity differences around the ice-sheet periphery (Fig. 13), which strongly influences the iceberg calving from tidewater glaciers. Iceberg calving under Case 2 has a greater ensemble spread ( 5 %; 18 Gt yr−1 ensemble range against 365 Gt yr−1 ensemble mean) than under Case 1 (Fig. 8). Similar to Case 1, however, the C2022 lowest heat-flow map again has the highest iceberg calving flux, while the relatively high M2018 and G2019 heat-flow maps have substantially lower iceberg calving fluxes at equilibrium.

Figure 12Case 2: (a–g) Surface ice velocity at transient equilibrium using the seven geothermal heat-flow maps, expressed as anomalies from their ensemble mean. Units are % and color bars saturate at ±150 %. (h) Ensemble mean surface ice velocity at transient equilibrium. Units are m yr−1.

Figure 13Case 2: (a–g) Anomaly in ice surface speed at Case 2 transient spin-up, in comparison with Case 1 nudged spin-up, using the seven geothermal heat-flow maps. Units in all plots are m and expressed as Case 2 minus Case 1. (h) Ensemble mean of ice surface speed anomaly. The color bars saturate at ±100 m yr−1.

4 Discussion

The apparent association of higher ice velocities with lower geothermal heat flows under Case 1 appears to be a clear artifact of nudging the basal friction coefficient during spin-up. This effect has previously been described as the surface velocity paradox, whereby constraining an ice-flow model to match observed ice thickness results in underestimating deformational velocities where basal sliding is present, and overestimating deformational velocities where basal sliding is absent (Ryser et al., 2014). Avoiding this surface velocity paradox is the main motivation for the Case 2 spin-up, in which basal friction coefficients are not nudged. Under Case 2, there is more variation in the geometry, velocity, and thermal state of the ice sheet at the end of the 10 000-year spin-up. Perhaps counterintuitively, however, the highest iceberg calving fluxes remain associated with the lowest heat-flow maps (and vice versa for lowest iceberg calving fluxes). In unconstrained Case 2 simulations, this behavior cannot be attributed to a model artifact from the surface velocity paradox associated with nudging in Case 1. We instead speculate that a substantial portion of this variability simply reflects increased ice thicknesses under decreased heat flow.

The potential influence of anomalously high geothermal heat flow on contemporary local ice-sheet form and flow has been previously highlighted, with suggestions including the following: the onset of the Northeast Greenland ice stream may be associated with elevated geothermal heat flow (Fahnestock et al., 2001); there may be a feedback between deeply incised glaciers and topographic enhancement of local geothermal heat flow (van der Veen et al., 2007); and that the transit of the Iceland hotspot may have deposited anomalous heat into the subglacial lithosphere that influences ice flow today (Alley et al., 2019). Our evaluation suggests that knowledge of where anomalously low geothermal heat flow may be influencing contemporary regional ice-sheet form and flow can help constrain the choice of heat-flow map. For example, the widespread presence of Last Glacial Period ice in the ablation area across North Greenland suggests that heat flow must be sufficiently low to prevent basal melt across the region (MacGregor et al., 2020). This broad condition is only characteristic of a minority of the heat flow maps we evaluate, specifically the SR2004, R2017 and C2022 maps.

South Dome appears to be the most sensitive portion of the ice sheet to the choice of the geothermal heat-flow basal boundary condition. There, the choice of heat-flow map results in an ensemble spread in ice-bed temperature of >10C over an area the size of Iceland. There is currently a poor level of scientific understanding of whether South Dome persisted through the Eemian interglacial, with some ice-sheet reconstructions suggesting persistence of the ice sheet's southern lobe (Quiquet et al., 2013; Stone et al., 2013) and others suggesting local deglaciation (Otto-Bliesner et al., 2006; Helsen et al., 2013). Our evaluation specifically highlights substantial disagreement over geothermal heat flow within the North Atlantic Craton that underlies South Dome. Similar to the contemporary persistence of Last Glacial Period ice in North Greenland, we speculate that paleo-ice-sheet simulations that adopt the low heat flow beneath South Dome characteristic of the R2017 map are more likely to yield an Eemian-persistent South Dome than paleo-ice-sheet simulations that adopt the high heat flow characteristic of the L2019 map. Simply put, the heat-flow map influences not only contemporary simulations of ice-sheet form and flow, but also paleo-ice-sheet simulations.We should note that, despite basal heat flow being a key factor in controlling ice dynamics, some other important physical processes (e.g., subglacial hydrology) are not considered in this study. The influence of different basal heat-flow models may not fully capture the role of enhanced basal meltwater in a warming climate. By holding basal friction coefficients fixed in time, Case 1 ignores the effects of evolving basal hydrology. Case 2 allows the thawed-bed area to change, but using a local till model that ignores subglacial water transport. Thus, Case 2 might be overly sensitive to local temperature changes, whereas more realistic hydrology changes would be spread over larger scales.

Furthermore, some higher-order ice-sheet models use data assimilation approaches (e.g., Hoffman et al., 2018) instead of spin-up, which may result in different model behaviors when applying different basal heat-flow datasets during initialization. Also, since our study focuses on the overall impacts of basal heat flow on Greenland ice-sheet dynamics, a more detailed understanding of the relative importance of thermal model components, such as ice frictional heating, heat advection and diffusion, is still required to improve the thermodynamic knowledge of the deep layers of the Greenland ice sheet.

5 Summary remarks

Given the non-linear dependence of deformational velocity on ice temperature, properly resolving the thermal state of the Greenland ice sheet is critical for generating reliable ice-flow simulations. We have performed both nudged and unconstrained ice-sheet spin-ups of 10 000-year duration with seven geothermal heat-flow models. Under a nudged spin-up, we find that the thawed-bed ice-sheet area ranges from 21.8 % to 54.4 % across these heat flow models. Under an unconstrained spin-up, the thawed-bed ice-sheet area is consistently larger, ranging from 33.5 % to 60.0 %. This suggests that unconstrained spin-ups generally yield a warmer Greenland ice-bed interface than constrained, or nudged, spin-ups. The unconstrained spin-up also yields inter-simulation differences in both ice thickness and velocity that are large in magnitude and extent. This ensemble of simulations highlights that sector-scale ice flow, both peripheral and interior, is at least moderately sensitive to the choice of heat-flow map.

The recent effort to compile all Greenland englacial temperature observations into a standardized database now permits the thermal state of ice-sheet simulations to be evaluated against all empirical data. Here, we evaluate simulated basal temperature against observed basal temperature at 27 selected Greenland boreholes. Despite the fact that the spatial resolutions of several basal heat-flow models are quite coarse in comparison with that of CISM, this evaluation still provides some insights on which heat flow map or spin-up approach is most locally suitable. Rather than quantitative comparisons against point temperature observations, however, there seems to be value in qualitative comparisons between heat-flow map and large-scale ice-sheet features, such as evaluating which heat-flow map can yield a widespread frozen bed in North Greenland under contemporary conditions. Naturally, evaluation of these seven heat-flow maps would be strengthened by using more than a single ice-flow model, as we do here.

Within our simulation ensemble, the unconstrained spin-ups may possibly be regarded as simulating more sensitive ice sheets than the nudged spin-ups, as the unconstrained spin-ups yield greater thawed-bed area and higher iceberg calving flux. While most recent ice-sheet simulations projecting Greenland's future sea-level contribution have focused on nudged spin-ups, our simulation ensemble unsurprisingly suggests that unconstrained spin-up is required to fully resolve the choice of geothermal heat-flow boundary condition on ice-sheet geometry and velocity. Given the strong influence of geothermal heat flow on ice dynamics that we document, it seems prudent to limit the direct intercomparison of ice-sheet simulations with those using a common heat-flow map. Similar to employing a range of commonly prescribed climate forcing scenarios, it would be ideal for future ISMIP ensembles to employ a range of commonly prescribed basal forcing conditions.

Data availability

To help accelerate community efforts toward exploring the influence of geothermal heat flow on ice-sheet simulations, we have deposited a copy of the seven geothermal heat-flow maps that we evaluate here at Zenodo (, Zhang et al., 2023). Interpolated versions of these seven geothermal heat-flow datasets are provided on a common coarse-resolution netCDF grid that conforms with CISM standards.

Author contributions

TZ and WC conceptualized this study and were responsible for formal analysis. AL and AW provided data curation. TZ, CX, WL and GL provided funding, resources, and software. All authors participated in interpreting the data and writing the manuscript.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


Tong Zhang and Cunde Xiao thank the National Key Research and Development Program of China (2023YFF0805200), the Natural Science Foundation of China (42271133), Faculty of Geographical Sciences, Beijing Normal University (2022-GJTD-01), an the State Key Laboratory of Earth Surface Processes and Resource Ecology (2022-ZD-05) for financial support. Anja Løkkegaard and William Colgan thank the Independent Research Fund Denmark (Sapere Aude 8049-00003) and the Novo Nordisk Foundation (Center for Sea-Level and Ice-Sheet Prediction) for financial support. Agnes Wansing thanks the European Space Agency and the German Research Council (DFG) for their financial support through the projects 4D-Greenland and GreenCrust. Gunter Leguy and William Lipscomb were supported by the NSF National Center for Atmospheric Research, which is a major facility sponsored by the U.S. National Science Foundation under Cooperative Agreement No. 1852977. Computing resources were provided by the Climate Simulation Laboratory at NSF NCAR’s Computational and Information Systems Laboratory (CISL). We thank two reviewers, Lu Li and Felicity McCormack, and the editor Alexander Robinson, for their help in improving the manuscript.

Financial support

This research has been supported by the National Key Research and Development Program of China (grant no. 2023YFF0805200), the National Natural Science Foundation of China (grant no. 42271133), the Faculty of Geographical Sciences, Beijing Normal University (grant no. 2022-GJTD-01), the State Key Laboratory of Earth Surface Processes and Resource Ecology (grant no. 2022-ZD-05), the Independent Research Fund Denmark (Sapere Aude (grant no. 8049-00003)), the Novo Nordisk Foundation (Center for Sea-Level and Ice-Sheet Prediction), the projects “4D-Greenland” and “GreenCrust” under the European Space Agency and the German Research Council (DFG), and the National Science Foundation under Cooperative Agreement no. 1852977.

Review statement

This paper was edited by Alexander Robinson and reviewed by Lu Li and Felicity McCormack.


Alley, R., Pollard, D., Parizek, B., Anandakrishnan, S., Pourpoint, M., Stevens, N., MacGregor, J., Christianson, K., Muto, A., and Holschuh, N.: Possible role for tectonics in the evolving stability of the Greenland Ice Sheet, J. Geophys. Res.-Earth, 124, 97–115,, 2019. 

Artemieva, I. M.: Lithosphere thermal thickness and geothermal heat flux in Greenland from a new thermal isostasy method, Earth-Sci. Rev., 188, 469–481,, 2019. 

Aschwanden, A., Fahnestock, M., and Truffer, M.: Complex Greenland outlet glacier flow captured, Nat. Commun., 7, 10524,, 2016. 

Bueler, E. and van Pelt, W.: Mass-conserving subglacial hydrology in the Parallel Ice Sheet Model version 0.6, Geosci. Model Dev., 8, 1613–1635,, 2015. 

Colgan, W., MacGregor, J., Mankoff, K., Haagenson, R., Rajaram, H., Martos, Y., Morlighem, M., Fahnestock, M., and Kjeldsen, K.: Topographic Correction of Geothermal Heat Flux in Greenland and Antarctica, J. Geophys. Res.-Earth, 126, e2020JF005598,, 2021. 

Colgan, W., Wansing, A., Mankoff, K., Lösing, M., Hopper, J., Louden, K., Ebbing, J., Christiansen, F. G., Ingeman-Nielsen, T., Liljedahl, L. C., MacGregor, J. A., Hjartarson, Á., Bernstein, S., Karlsson, N. B., Fuchs, S., Hartikainen, J., Liakka, J., Fausto, R. S., Dahl-Jensen, D., Bjørk, A., Naslund, J.-O., Mørk, F., Martos, Y., Balling, N., Funck, T., Kjeldsen, K. K., Petersen, D., Gregersen, U., Dam, G., Nielsen, T., Khan, S. A., and Løkkegaard, A.: Greenland Geothermal Heat Flow Database and Map (Version 1), Earth Syst. Sci. Data, 14, 2209–2238,, 2022. 

Cuffey, K. M. and Paterson, W. S. B.: The Physics of Glaciers, Elsevier, Burlington, ISBN 978-0-123-69461-4, 2010. 

Fahnestock, M., Abdalati, W., Joughin, I., Brozena, J., and Gogineni, P.: High Geothermal Heat Flow, Basal Melt, and the Origin of Rapid Ice Flow in Central Greenland, Science, 294, 2338–2342,, 2001. 

Fettweis, X., Box, J. E., Agosta, C., Amory, C., Kittel, C., Lang, C., van As, D., Machguth, H., and Gallée, H.: Reconstructions of the 1900–2015 Greenland ice sheet surface mass balance using the regional climate MAR model, The Cryosphere, 11, 1015–1033,, 2017. 

Fox Maule, C., Purucker, M. E., and Olsen, N.: Inferring magnetic crustal thickness and geothermal heat flux from crustal magnetic field models, Danish Climate Centre Report, 9–09, 2009. 

Fuchs, S., Beardsmore, G., Chiozzi, P., Gola, G., Gosnold, W., Harris, R., Jennings, S., Liu, S., Negrete-Aranda, R., Neumann, F., Norden, B., Poort, J., Rajver, D., Ray, L., Richards, M., Smith, J., Tanaka, A., and Verdoya, M.: A new database structure for the IHFC Global Heat Flow Database, International Journal of Terrestrial Heat Flow and Applied Geothermics, 22, 1–4,, 2021. 

Goelzer, H., Nowicki, S., Payne, A., Larour, E., Seroussi, H., Lipscomb, W. H., Gregory, J., Abe-Ouchi, A., Shepherd, A., Simon, E., Agosta, C., Alexander, P., Aschwanden, A., Barthel, A., Calov, R., Chambers, C., Choi, Y., Cuzzone, J., Dumas, C., Edwards, T., Felikson, D., Fettweis, X., Golledge, N. R., Greve, R., Humbert, A., Huybrechts, P., Le clec'h, S., Lee, V., Leguy, G., Little, C., Lowry, D. P., Morlighem, M., Nias, I., Quiquet, A., Rückamp, M., Schlegel, N.-J., Slater, D. A., Smith, R. S., Straneo, F., Tarasov, L., van de Wal, R., and van den Broeke, M.: The future sea-level contribution of the Greenland ice sheet: a multi-model ensemble study of ISMIP6, The Cryosphere, 14, 3071–3096,, 2020. 

Goldberg, D. N.: A variationally derived, depth-integrated approximation to a higher-order glaciological flow model, J. Glaciol., 57, 157–170,, 2011. 

Greve, R.: Geothermal heat flux distribution for the Greenland ice sheet, derived by combining a global representation and information from deep ice cores, Polar Data Journal, 3, 22–36,, 2019. 

Helsen, M. M., van de Berg, W. J., van de Wal, R. S. W., van den Broeke, M. R., and Oerlemans, J.: Coupled regional climate–ice-sheet simulation shows limited Greenland ice loss during the Eemian, Clim. Past, 9, 1773–1788,, 2013. 

Hoffman, M. J., Perego, M., Price, S. F., Lipscomb, W. H., Zhang, T., Jacobsen, D., Tezaur, I., Salinger, A. G., Tuminaro, R., and Bertagna, L.: MPAS-Albany Land Ice (MALI): a variable-resolution ice sheet model for Earth system modeling using Voronoi grids, Geosci. Model Dev., 11, 3747–3780,, 2018. 

Hooke, R.: Principles of Glacier Mechanics, Cambridge University Press, ISBN 978-1108446075, 2019. 

Jessop, A., Hobart, M., and Sclater, J.: The World Heat Flow Data Collection-1975, Geothermal Series Number 5, Geological Survey of Canada, Ottawa, Canada,, 1976. 

Karlsson, N., Solgaard, A., Mankoff, K., Gillet-Chaulet, F., MacGregor, J., Box, J., Citterio, M., Colgan, W., Larsen, S., Kjeldsen, K., Korsgaard, N., Benn, D., Hewitt, I., and Fausto, R.: A first constraint on basal melt-water production of the Greenland ice sheet, Nat. Commun., 12, 3461,, 2021. 

Leguy, G. R., Lipscomb, W. H., and Asay-Davis, X. S.: Marine ice sheet experiments with the Community Ice Sheet Model, The Cryosphere, 15, 3229–3253,, 2021. 

Lipscomb, W. H., Price, S. F., Hoffman, M. J., Leguy, G. R., Bennett, A. R., Bradley, S. L., Evans, K. J., Fyke, J. G., Kennedy, J. H., Perego, M., Ranken, D. M., Sacks, W. J., Salinger, A. G., Vargo, L. J., and Worley, P. H.: Description and evaluation of the Community Ice Sheet Model (CISM) v2.1, Geosci. Model Dev., 12, 387–424,, 2019. 

Lipscomb, W. H., Leguy, G. R., Jourdain, N. C., Asay-Davis, X., Seroussi, H., and Nowicki, S.: ISMIP6-based projections of ocean-forced Antarctic Ice Sheet evolution using the Community Ice Sheet Model, The Cryosphere, 15, 633–661,, 2021. 

Løkkegaard, A., Mankoff, K. D., Zdanowicz, C., Clow, G. D., Lüthi, M. P., Doyle, S. H., Thomsen, H. H., Fisher, D., Harper, J., Aschwanden, A., Vinther, B. M., Dahl-Jensen, D., Zekollari, H., Meierbachtol, T., McDowell, I., Humphrey, N., Solgaard, A., Karlsson, N. B., Khan, S. A., Hills, B., Law, R., Hubbard, B., Christoffersen, P., Jacquemart, M., Seguinot, J., Fausto, R. S., and Colgan, W. T.: Greenland and Canadian Arctic ice temperature profiles database, The Cryosphere, 17, 3829–3845,, 2023. 

Lucazeau, F.: Analysis and Mapping of an Updated Terrestrial Heat Flow Data Set, Geochem. Geophy. Geosy., 20, 4001–4024,, 2019. 

MacGregor, J., Fahnestock, M., Colgan, W., Larsen, N., Kjeldsen, K., and Welker, J.: The age of surface-exposed ice along the northern margin of the Greenland Ice Sheet, J. Glaciol., 66, 667–684,, 2020. 

MacGregor, J. A., Chu, W., Colgan, W. T., Fahnestock, M. A., Felikson, D., Karlsson, N. B., Nowicki, S. M. J., and Studinger, M.: GBaTSv2: a revised synthesis of the likely basal thermal state of the Greenland Ice Sheet, The Cryosphere, 16, 3033–3049,, 2022. 

Martos, Y., Jordan, T., Catalán, M., Jordan, T., Bamber, J., and Vaughan, D.: Geothermal heat flux reveals the Iceland hotspot track underneath Greenland, Geophys. Res. Lett., 45, 8214–8222,, 2018. 

Morlighem, M., Williams, C., Rignot, E., An, L., Arndt, J., Bamber, J., Catania, G., Chauché, N., Dowdeswell, J., Dorschel, B., Fenty, I., Hogan, K., Howat, I., Hubbard, A., Jakobsson, M., Jordan, T., Kjeldsen, K., Millan, R., Mayer, L., Mouginot, J., Noël, B., O'Cofaigh, C., Palmer, S., Rysgaard, S., Seroussi, H., Siegert, M., Slabon, P., Straneo, F., van den Broeke, M., Weinrebe, W., Wood, M., and Zinglersen, K.: BedMachine v3: Complete bed topography and ocean bathymetry mapping of Greenland from multi-beam echo sounding combined with mass conservation, Geophys. Res. Lett., 44, 11051–11061,, 2017. 

Otto-Bliesner, B., Marshall, S., Overpeck, J., Miller, G., and Hu, A.: Simulating Arctic Climate Warmth and Icefield Retreat in the Last Interglaciation, Science, 311, 1751–1753,, 2006. 

Pollack, H. N., Hurter, S. J., and Johnson, J. R.: Heat flow from the Earth's interior: Analysis of the global data set, Rev. Geophys., 31, 267–280,, 1993. 

Pollard, D. and DeConto, R. M.: A simple inverse method for the distribution of basal sliding coefficients under ice sheets, applied to Antarctica, The Cryosphere, 6, 953–971,, 2012. 

Quiquet, A., Ritz, C., Punge, H. J., and Salas y Mélia, D.: Greenland ice sheet contribution to sea level rise during the last interglacial period: a modelling study driven and constrained by ice core data, Clim. Past, 9, 353–366,, 2013. 

Rezvanbehbahani, S., Stearns, L., Kadivar, A., Walker, J., and van der Veen, C.: Predicting the Geothermal Heat Flux in Greenland: A Machine Learning Approach, Geophys. Res. Lett., 44, 12271–12279,, 2017. 

Rogozhina, I., Petrunin, A., Vaughan, A., Steinberger, B., Johnson, J., Kaban, M., Calov, R., Rickers, F., Thomas, M., and Koulakov, I.: Melting at the base of the Greenland ice sheet explained by Iceland hotspot history, Nat. Geosci., 9, 366–369,, 2016. 

Ryser, C., Lüthi, M., Andrews, L., Hoffman, M., Catania, G., Hawley, R., Neumann, T., and Kristensen, S. Sustained high basal motion of the Greenland ice sheet revealed by borehole deformation. J. Glaciol., 60, 647–660,, 2014. 

Shapiro, N. and Ritzwoller, M.: Inferring surface heat flux distributions guided by a global seismic model: particular application to Antarctica, Earth Planet. Sc. Lett., 223, 213–224,, 2004. 

Smith-Johnsen, S., de Fleurian, B., Schlegel, N., Seroussi, H., and Nisancioglu, K.: Exceptionally high heat flux needed to sustain the Northeast Greenland Ice Stream, The Cryosphere, 14, 841–854,, 2020. 

Stone, E. J., Lunt, D. J., Annan, J. D., and Hargreaves, J. C.: Quantification of the Greenland ice sheet contribution to Last Interglacial sea level rise, Clim. Past, 9, 621–639,, 2013. 

Tarasov, L. and Peltier, W. R.: Greenland glacial history, borehole constraints, and Eemian extent, J. Geophys. Res.-Sol. Ea., 108, 2143,, 2003. 

van der Veen, C., Leftwich, T., Von Frese, R., Csatho, B., and Li, J. Subglacial topography and geothermal heat flux: Potential interactions with drainage of the Greenland ice sheet, Geophys. Res. Lett., 34, L12501,, 2007. 

Weertman, J.: The Unsolved General Glacier Sliding Problem, J. Glaciol., 23, 97–115,, 1979.  

Zhang, T., Colgan, W., and Wansing, A.: Greenland Geothermal Heatflow, Zenodo [data set],, 2023. 

Short summary
The geothermal heat flux determines how much heat enters from beneath the ice sheet, and thus impacts the temperature and the flow of the ice sheet. In this study we investigate how much geothermal heat flux impacts the initialization of the Greenland ice sheet. We use the Community Ice Sheet Model with two different initialization methods. We find a non-trivial influence of the choice of heat flow boundary conditions on the ice sheet initializations for further designs of ice sheet modeling.