Millennial-scale migration of the frozen/melted basal boundary, western Greenland ice sheet

Abstract

The geometry and thermal structure of western Greenland ice sheet are known to have undergone relatively substantial change over the Holocene. Evolution of the frozen and melted fractions of the bed associated with the ice-sheet retreat over this time frame remains unclear. We address this question using a thermo-mechanically coupled flowline model to simulate a 11 ka period of ice-sheet retreat in west central Greenland. Results indicate an episode of ~100 km of terminus retreat corresponded to ~16 km of upstream frozen/melted basal boundary migration. The majority of migration of the frozen area is associated with the enhancement of the frictional and strain heating fields, which are accentuated toward the retreating ice margin. The thermally active bedrock layer acts as a heat sink, tending to slow contraction of frozen-bed conditions. Since the bedrock heat flux in our region is relatively low compared to other regions of the ice sheet, the frozen region is relatively greater and therefore more susceptible to marginward changes in the frictional and strain heating fields. Migration of melted regions thus depends on both geometric changes and the antecedent thermal state of the bedrock and ice, both of which vary considerably around the ice sheet.

1. Introduction

As ice flows from the cold interior of the Greenland ice sheet (GrIS) to the margins, it exchanges energy with surrounding ice and the underlying bedrock, and it undergoes heating as it deforms and generates friction at the basal boundary. At the bed of the ice sheet, these processes result in a thermal pattern where the central core of GrIS is mainly frozen, and the outer flanks are at the pressure melting point (MacGregor and others, 2016). The onset of continuously melted bed conditions along a flowline has been termed the frozen/melted boundary (Brinkerhoff and others, 2011).

The frozen/melted boundary signifies the initiation of melting basal ice, and therefore a transition point on the flowline to basal processes involving liquid water and decoupling of the ice and bedrock. Melted bed conditions are necessary for recharge of groundwater systems (e.g. Flowers, 2015) such that the areal distribution of melted ice-sheet bed governs regional hydrologic processes in the terrestrial and marine areas adjacent to the ice sheet (DeFoor and others, 2011). The production of basal meltwater plays a role in reducing basal traction and enhancing the rates of basal sliding (Zwally and others, 2002; Parizek and Alley, 2004). Whatever the mechanism, regions of the ice sheet with melted bed generally have higher velocities (MacGregor and others, 2016) and therefore greater advection of ice to the margins of the ice sheet, where mass is lost via marine calving or melting in warmer temperatures.

Both climate and Greenland's ice-sheet geometry have varied substantially since the Last Glacial Maximum (Lecavalier and others, 2014). In the Kangerlussuaq sector, the GrIS was advanced ~100 km west of its current boundary, reflecting a colder historical climate (Lecavalier and others, 2014; Briner and others, 2016; Young and others, 2020). While the ice formed at the divide over the Holocene period has yet to flow across the profile from divide-to-margin (Dahl-Jensen and others, 1998), the ice thickness and surface slope of the ice sheet have changed commensurately with the 25% reduction in extent that has taken place (Lecavalier and others, 2014). However, the contemporaneous evolution of the frozen and melted fractions of the bed (i.e. migration of the frozen/melted boundary) remains unclear.

Modeling is the most applicable tool for investigating the evolution of the basal thermal state of western GrIS during Holocene retreat. At the short century timescale, modeling suggests that partitioning between frozen/melted fractions of the bed is highly stable under climate and geometric perturbations (Seroussi and others, 2013). However, models run over ice age cycles, which necessarily employ simplified model physics and coarse resolution, suggest a considerable change in the basal thermal state over this longer timescale (Huybrechts, 1996; Marshall and Clark, 2002). The transient behavior of the bed temperature results from complicated interaction of surface mass balance (which impacts the flow field) and historical surface temperatures. Results are also highly sensitive to prescribed geothermal heat flux (Greve and Hutter, 1995; Brinkerhoff and others, 2011).

Here we explore the migration of frozen/melted conditions in western GrIS over a 11 ka period of Holocene ice-sheet retreat in western GrIS. Prior work has shown that simulations of Holocene retreat across western GrIS are strongly influenced by model resolution in places with rougher and more complex bed topography (Cuzzone and others, 2019). Because the outer flanks of western GrIS contain numerous deep and steep-sided subglacial troughs, we focus our investigation on a single flowline transect so that we can employ high time/space resolution in our model. Our transient model includes higher order stresses, a well-informed paleoclimate (Buizert and others, 2018; Downs and others, 2020), thermally active bedrock, and honors previously published records of terminus retreat (Briner and others, 2016; Lesnek and others, 2020; Young and others, 2020). We track the constitutive components of the basal heat budget during the transient simulation to isolate the energy sources and sinks associated with the migration of frozen/melted bed conditions. The aim of this profile-specific study is to elucidate the changes in the thermal state of the bed during an important period of ice-sheet retreat; changes with implications for hydrologic and ice dynamics processes.

2. Methods

2.1. Study transect

We focus our work on an ice flowline located in the land terminating portion of western GrIS (Fig. 1). This flowline allows us to leverage important constraints on our simulations, including reconstructions of past regional climate (Buizert and others, 2018; Downs and others, 2020), detailed basal topography (Morlighem and others, 2017) and observational knowledge of bedrock heat flux (Meierbachtol and others, 2015). The latter is rare for Greenland, offered here by temperature measurements in a bedrock borehole partly extending under the current ice margin (Claesson Liljedahl and others, 2016).

Fig. 1. Ice sheet and bedrock model domains. Initial and final ice-sheet geometries are shown in blue. Dashed lines demarcate all mesh nodes, and solid vertical lines indicate every tenth mesh node. Flowline location in western Greenland is shown in the inset. The modern surface elevation of this flowline is shown with a dashed purple line.

Detailed chronologies of Holocene ice margin retreat have been established for this sector of the ice sheet based on glacial geologic evidence with ¹⁰Be and ¹⁴C age dating (Briner and others, 2016; Lesnek and others, 2020; Young and others, 2020). The ice margin retreated most rapidly during the first half of the last 11 ka. Glacial geologic evidence suggests the ice retreated an unknown distance beyond its modern extent, and then readvanced (van Tatenhove and others, 1996). Models run at the ice sheet and ice age scales suggest the retreat was many tens of km beyond the current margin in this area of GrIS (Simpson and others, 2009; Lecavalier and others, 2014). Downs and others (2020) conducted a focused study of this sector, incorporating the glacial record and ice core climate data into an inverse model to solve for climate forcing and ice margin changes. This work shows the minimum ice terminus position during the last 10 ka likely <~9 km behind the modern position. Other recent studies have shown similar results in western Greenland (Cuzzone and others, 2019; Kajanto and others, 2020). However, the model inversion in Downs and others (2020) is based on the most detailed measurements available for this sector of GrIS, and there are no constraints on the terminus position between ~7 ka and present. As a result, it is unknown to what extent the ice sheet may have retreated beyond its modern position during this time period, and the predictions from the inverse model are limited. Since our approach here is to closely honor the climate and terminus records from Downs' and others (2020) study, this limitation is carried through into this work.

2.2. Ice sheet and bedrock model

We use a thermomechanically-coupled flowline model with a 2-D incompressible Stokes' ice flow solution. The surface of the flowline is dynamically updated according to a surface kinematic boundary condition. Coupled solutions for ice flow, temperature and ice surface elevation are performed using the finite element software Elmer/Ice (Gagliardini and others, 2013). A list of physical parameters used in the model can be found in Table 1.

Table 1. List of physical parameters used in the model

T refers to the temperature of the ice in Kelvin.

A flowline model allows for 1 km horizontal resolution over the entire domain, with manageable computational costs. The domain is split into two subdomains representing ice and thermally active bedrock (Fig. 1). The thermally active bedrock includes five static horizontal layers that extend from the ice-sheet bed to 3 km below sea level. The depth of the bedrock domain was chosen such that temperature changes at the surface are not felt at the base of the bedrock at the timescale of the simulations. The ice domain has ten dynamically rescaled horizontal layers between the ice-sheet bed and surface. Simulations are run at 1-year time steps for a duration of 11.4 ka.

2.2.1. Basal sliding

We adopt a standard viscous linear sliding law, commonly employed in modeling studies (e.g. Seroussi and others, 2013; Larour and others, 2014; Cuzzone and others, 2019) and used for this flowline by Downs and others (2020):

In this relation, the basal shear stress τ _b is proportional to the sliding speed u _b via N, the effective pressure, and β ², a constant basal traction parameter. The basal water pressure is

(1)$$\tau _{\rm b} = \beta ^2Nu_{\rm b}$$

assumed to be a constant fraction of ice overburden pressure P ₀, taken to be 85% of overburden based on observational constraints of basal water pressure (Wright and others, 2016). Effective pressure is then given by β ² and was found by minimizing the misfit between modeled and observed surface velocities (Downs and others, 2020)

(2)$$N = P_0-\;0.85P_0$$

of the present-day flowline data from Mouginot and others (2017). As in Downs and others (2020) we assume the field is unchanging in time along the transect; however, dynamics in basal traction are included through the effective pressure, which is linearly proportional to ice thickness.

2.2.2. Ice surface evolution

The change in the surface of the ice sheet follows

(3)$$\displaystyle{{\partial z_{\rm s}} \over {\partial t}} + u_{\rm x}\displaystyle{{\partial z_{\rm s}} \over {\partial x}}-u_{\rm z} = a_{\rm s}, \;$$

where z _s is the ice surface elevation, u _x and u _z are the horizontal and vertical velocity at the surface of the ice sheet from the solution of Stokes' equation, and a _s is the net accumulation/ablation at the surface. A positive degree day model is used to determine a _s, which is described in more detail in Section 2.2.5.

After solving Eqn (3) at each time step, the vertical coordinates of the ice sheet are projected to the new ice-sheet surface. However, we also include a thermally active bedrock layer in our model domain (Section 2.2.4) and this mapping does not occur in that layer. Consequently, a mask was applied to the thermally active bedrock layer while projecting the vertical coordinate of the ice in order to achieve model convergence in Elmer/Ice. Isostatic adjustments to ice and bedrock elevations were not included in this effort, but their effects on the retreat history in this region on this timescale were not found to be significant in Downs and others (2020).

2.2.3. Ice-sheet temperature

Ice-sheet temperatures are solved using the standard heat equation (Gagliardini and others, 2013). This accounts for both advection and conduction of energy within the ice sheet. Additional heat sources from strain heating and frictional heating are also present in the model. The bottom of the ice is thermally coupled to the bedrock (see Section 2.2.4) and temperatures at the surface of the ice sheet are described by

(4)$$T_{{\rm surf}}( t ) = T_{\rm m} + \Delta T( t ) + \alpha \;( {S( t ) -S_{\rm m}} ) .$$

T _m is the modern surface temperature along the flowline based on a 30-year average of data from Box (2013). ΔT(t) is the local temperature anomaly along the flowline as a function of time. Temperature anomaly data are obtained from monthly climate reconstructions from Buizert and others (2018), and a lapse rate, α, of 5 K km⁻¹ (Abe-Ouchi and others, 2007) is used. The surface elevation at the current time step of the model is S(t), and the modern surface elevation of the ice sheet is S _m.

2.2.4. Thermally active bedrock

We used high-resolution topography of the ice-sheet bed from Morlighem and others (2017), and a bedrock layer that is 3000 m thick. Temperature evolution within the bedrock is described using the conductive heat equation. The top of the bedrock is thermally coupled to the base of the ice sheet, and the heat flow in the bedrock at this interface is referred to here as the ‘bedrock heat flux’. A spatially and temporally constant ‘deep geothermal heat flux’ is prescribed at the bottom of the thermally active bedrock layer. Observational constraints near the flowline terminus show the bedrock heat flux to be ~27 mW m⁻² (Claesson Liljedahl and others, 2016), but bedrock heat flux values are potentially higher in the ice-sheet interior (Meierbachtol and others, 2015). These values represent paleoclimatically uncorrected near surface values for the flowline, and the true value for the deep geothermal heatflux is likely ~10mW m⁻² higher than these values.

We recognize that as the bedrock heat flux evolves with the ice sheet, inconsistencies arise between the bedrock heat flux at present-day produced by the model, and present-day measurements in southwestern Greenland. Lacking a constrained record of the long-term thermal evolution of the bedrock layer, starting the simulation with an ice/bedrock interface similar to recent observations was the favorable choice. Further, the common alternative approach of using data from maps (e.g. Shapiro and Ritzwoller, 2004), which tend to overestimate deep geothermal heat flux in this region (Rogozhina and others, 2012), and do not necessarily represent the initial bedrock heat flux at the ice interface, would lead to much larger inconsistencies between the current and measured bedrock heat flux (see Section 4.1). For these reasons, we choose to apply a deep geothermal flux value of 30 mW m⁻² in our simulations that more closely reflect the average bedrock heat flux along the flowline, and provide a lower bound on the deep geothermal heat flux in the region.

To address the uncertainty in the deep geothermal heat flux and inform upscaling of our results to other regions of the ice sheet, we performed an analysis of the sensitivity of the model to different values of the deep geothermal heat flux.

2.2.5. Surface mass balance

Surface mass balance is calculated using a positive degree-day model adapted from Downs and others (2020) and is based on Jóhannesson and others (1995). The precipitation at a given time, P(t), is found using

(5)$$P( t ) = P_{\rm m}\exp ( {\lambda_{\rm p}( {T_{{\rm surf}}( t ) -T_{\rm m}} ) } ) + \Delta P( t ) , \;$$

where P _m is the modern precipitation from a 30-year modern average along the flowline (Box, 2013). A precipitation parameter, λ _p, translates to a change in precipitation of 7% per degree K. T _surf(t) is the surface temperature at time t calculated using Eqn (4), and ΔP(t) is a local, time-dependent precipitation anomaly.

Precipitation and positive degree days for the positive degree day model were calculated for each month of the year (Jóhannesson and others, 1995). To account for melting corresponding to positive degree days, ice-covered surfaces melt at a rate of 8 mm w.e. K⁻¹ day⁻¹, and snowpack melts at a rate of 5 mm w.e. K⁻¹ day⁻¹. Snow in the model initially melts and then refreezes forming ice in the snowpack. Once the ice fraction in the snowpack reaches 60%, all additional melt becomes runoff.

In order to produce a retreat that is consistent with the retreat of the flowline, we used a ΔP(t) that was an adjusted version of the precipitation anomaly for the northern flowline in Downs and others (2020). An adjustment of ~+0.1 m w.e. a⁻¹ was added to account for differences in the initialization procedure and slight differences between the bed topography data used in this model and in Downs and others (2020). The adjustment was found by iteratively updating the reference precipitation anomaly if the model started to diverge from the constraints on the terminus position during the retreat. This adjustment allowed for a good match of the retreat pattern and modern terminus position for the flowline, with an observationally-constrained climate.

2.3. Initialization

The model was initialized in three steps in order to estimate Last Glacial Maximum conditions. First, a steady-state simulation was run beginning with a parabolic ice-sheet surface. This provided a reasonable starting point for the interior temperature state and aided convergence in subsequent transient runs. Second, a transient simulation was run using the steady-state model output as the initial condition. The ΔP was held constant at 0.46 m w.e. a⁻¹ and ΔT was fixed at the initial time step. This second initialization step allowed for the ice-sheet surface to relax into a near steady-state flow configuration with a terminus that matched the observed constraints of terminus position at 11.4 ka BP (Young and others, 2020). This relaxation eliminated all large transient trends and was found to require 3000 years. Third, a secondary steady-state temperature simulation initialized with the relaxed ice sheet was run. This adjusted the interior temperature of the ice sheet to as close to steady state as possible.

2.4. Energy balance near frozen/melted boundary

To evaluate the influence of individual components of the energy budget on migration of the frozen/melted boundary, we monitor the heat fluxes due to various processes in and out of an arbitrary but computationally convenient control volume of 1 m³ at two observation points on the bed. One point is located at 178 km from the ice divide, at a location that is overrun by migration of the frozen/melted boundary (Fig. 2). Once the point warms to the pressure melting point, our methodology described below for partitioning the components of the energy budget no longer apply. Thus, we evaluate a second point located 160 km from the ice divide, 11 km closer to the divide than the final position of the frozen/melted boundary. This point is not affected by the discontinuity in the temperature field at the frozen/melted boundary, but is near enough to the frozen/melted boundary that the documented changes in the processes contributing to the basal energy budget are relevant to the establishment of melted conditions over the full simulation period.

Fig. 2. Simulated temperature of ice and bedrock at the end of the 11.4 ka transient run. Light blue lines show intermediate ice surface profiles at 100-year intervals, and dotted and dashed lines indicate the initial and final position of the ice surface, respectively. The solid red line represents the range of frozen/melted boundary positions. The black arrows indicate the location of the temperature observation points. The dotted red line indicates the depth of the pressure melting point at the end of the simulation.

We partition the energy budget into five sources and sinks: vertical exchanges of energy from conduction within the ice, vertical exchanges of energy from conduction within the bedrock, strain heating, frictional heat produced at the bed and fluxes of energy from the advection of ice along the bed.

Bedrock and ice energy conduction can be approximated using Fourier's law of thermal conduction,

(6)$$q_j = -k_j\displaystyle{{\partial T} \over {\partial z}}, \;$$

where q is the heat flux, k is the thermal conductivity, and ∂T/∂z is the vertical temperature gradient. The subscript j represents values in either ice (i) or bedrock (r). At the bed, horizontal temperature gradients are exceedingly small and do not contribute significantly to the energy budget. For this reason, only the vertical heat flux in the ice is computed.

Frictional heating and strain heating at the bed is output directly by the Elmer/Ice model. Strain heating is output as a volume source, and as such is integrated along the vertical axis for consistency in units. Since the reference volume is small, the amount of strain heating at the bed is also small compared to the frictional heating at the bed.

The final flux of energy is from the advection of the ice. Since ice cannot flow into the bed, vertical velocities near the bed are exceedingly small, and thus vertical energy fluxes from advection do not play a significant role in the energy balance. This leaves the horizontal advection of energy at the bed, which can be found using the observed temperature change. Because the time change in temperature at the bed is output from the model, we can approximate the total energy delivered to the reference volume on the bed, Q, between time steps using

(7)$$Q = c_{\rm i}\rho _{\rm i}V\;( {T( {t + \Delta t} ) -T( t ) } ) , \;$$

where c is the specific heat of the ice, V is the volume of the reference volume, ρ is the density of the ice, and (T(t + Δt) − T(t)) is the difference in temperature between the model time step Δt. Q can also be approximated by summing up the different energy sources integrated over time and space. We assume that each of the individual fluxes (n) is constant over a time step, which gives

(8)$$Q = \mathop \sum \limits^n ( {q_j} ) _nA\Delta t, \;$$

where (q _j)_n is energy flux, and A is the area of the flux surface. Since each (q _j)_n is accounted for except for the flux from horizontal ice advection, we can equate Eqns (7) and (8) in order to quantify this flux.

3. Results

3.1. Geometry and ice flow

Our transient simulation produced a terminus retreat pattern consistent with Downs and others (2020) and also resulted in a good match to the modern day flowline geometry with a terminus position ~300 km from the ice divide (Fig. S1). The majority of the retreat along this flowline occurs between 11.4 and 8.2 ka BP. Subsequently, the flowline experiences minor variation in terminus position of <10 km (Fig. 3).

Fig. 3. Location of the terminus and frozen/melted boundary over the course of the simulation, plotted on the same timescale to highlight the correlation between terminus and frozen/melted boundary movement. Note that the terminus position scale (right y-axis) varies ~10 times more than the frozen/melted boundary position (left y-axis).

Ice thinning across the profile is most substantial during the initial phase of retreat, but thinning continues throughout the simulation period. The region extending from the ice divide to the frozen/melted boundary area thins by 600–700 m (Fig. 2). As the ice margin initially retreats between 11.4 and 8.7 ka BP, the surface slope increases across the profile. Consequently, the driving stress at our 178 km reference point increases by almost 15% during this time, despite thinning ice (Fig. 4). With enhanced driving stress, ice speed similarly increases, but not in a simple 1:1 fashion. The driving stress peaks at ~8.7 ka BP, whereas the speed peaks ~1200 years later. Peak speeds are ~50% greater than the initial values.

Fig. 4. Evolution of the driving stress and ice speed at the 178 km point. Note that initial increases in driving stress are necessarily the result of steepening surface slopes, since the ice thickness is decreasing throughout the simulation.

A detailed analysis of the sensitivity of our model to bed roughness showed that our results are highly altered if the bed is flat, but are rather insensitive to varying levels of bed undulations given consistent flowline scale topography (Figs S2 and S3). Thus, any potential errors in the Morlighem and others (2017) representation of the bed topography do not appear important to our results.

3.2. Thermally active bedrock

The bedrock layer beneath the ice sheet is potentially a large sink for energy, with subfreezing temperatures penetrating up to 2 km into the bedrock (Fig. 2). The impact of bedrock is demonstrated by a comparison between model runs performed with and without thermally active bedrock. At the end of the simulation period, ice in the model domain had warmed in both runs, but the ice between the divide and the frozen/melted boundary is 1–2 K colder when thermally active bedrock is present (Fig. 5a). However, within several tens of km of the divide where near-vertical ice flow brings cold to the bed, the increased flow speed advected more cold to the bed between about 9.4 and 3.4 ka BP (Fig. 6) resulting in temporary cooling in both runs. In this case, the thermally active bedrock acted as a sink for the added cold content (or a heat source), which caused the ice to be 0.3 K warmer than model runs without a thermally active bedrock layer. Thus, the dampening effect of the thermally active bedrock significantly reduces the magnitude of basal temperature change in our model irrespective of sign, and is therefore an essential aspect of the transient behavior of the basal thermal state.

Fig. 5. Initial (dashed) and final (solid) basal temperature with (blue) and without (orange) thermally active bedrock present for locations within 190 km of the ice divide. Red lines indicate initial (dashed) and final (solid) pressure melting points. Inset indicates the frozen/melted boundary movement that occurs when (a) thermally active bedrock is removed, (b) pressure melting point is kept at the initial value, (c) when thermally active bedrock is present, (c) is the most realistic scenario and is used in all other plots.

Fig. 6. Basal temperature (a) and change in temperature from initial basal temperature (b) along the retreating ice-sheet transect. The dotted orange line highlights the terminus position, and the dotted blue line indicates the frozen/melted boundary position.

3.3. Temperature evolution

The bed within ~100 km of the divide cools slightly or has little change in temperature between 11.4 and 6.4 ka BP, and then warms continuously from 6.4 ka BP until present (Fig. 6). Ice temperature evolution in this region reflects changes in the ice flow field and the increasing temperature of the ice surface as climate warms and surface elevation drops. The bed in this region experienced a net temperature increase over the 11 ka of 1–1.5 K. Further from the divide, the remaining frozen bed warms over the entire simulation, with most of the warming between 11.4 and 8.4 ka BP. The net temperature increase of the bed in this region is 1.5–3 K.

The frozen/melted boundary began migrating at 10.4 ka BP, several hundred years after terminus retreat. The frozen/melted boundary, initially at ~188 km from the divide, migrated upstream 16 km between 11.4 and 6.4 ka BP, and then stabilized between 171 and 172 km for the remainder of the simulation. Whereas the ice margin retreated ~100 km, the upstream migration of the frozen/melted boundary was just 16 km (Fig. 3).

Ice thinning associated with retreat caused an increase of the pressure melting point across the area near the frozen/melted boundary of 0.6–0.7 K, further restricting basal melting and migration of the frozen/melted boundary. The increase of the pressure melting point results in ~3 km less retreat of the frozen/melted boundary (Fig. 5b), which is ~20% of the total movement of the frozen/melted boundary. Thus, adjustment of the pressure melting point is an important factor in frozen/melted boundary migration during an ice-sheet retreat and thinning but not a primary driver of changes in the frozen/melted boundary location.

Sensitivity experiments with different deep geothermal heat flux prescriptions illustrate the changing sensitivity of frozen/melted boundary migration to the initial frozen/melted boundary location (Table 2, Figs S4–S9). The range of deep geothermal heat flux values given illustrate the sensitivity of the system, but the true value for the deep geothermal heat flux in our region is likely between 30 and 50 mW m⁻². Simulations with increased deep geothermal heat flux (relative to our 30 mW m⁻² model runs) indicate a frozen/melted boundary position that is initially located closer to the ice divide. Frozen/melted boundary migration is directly linked to its initial location. Net migration distances decrease when the initial frozen/melted boundary location is closer to the ice divide, and under the highest prescribed deep geothermal heat flux, the frozen/melted boundary actually temporarily advances toward the ice margin. For values >90mW m⁻², our flowline geometry produces no frozen regions, and therefore represents an upper bound on the system.

Table 2. Initial frozen/melted boundary positions and net frozen/melted boundary movement for a large range of geothermal heat flux values (Q _geo).

Note that the deep geothermal heat flux in our region is likely between 30 and 50 mW m⁻²_.

3.4. Bed energy budget

The bed reference point at 178 km is overrun by the upstream-migrating frozen/melted boundary at 8.2 ka BP. Individual heat fluxes at this point demonstrate the evolving heat balance at the bed associated with frozen/melted boundary migration (Figs 7a–d). Nearly all of the warming at the 178 km point, a full 2 K, occurs during the ~1900 years after 10.3 ka BP. The warming begins just after the frozen/melted boundary, initially 10 km down-flow from the observation point, begins to migrate toward the divide. During the bed warming episode, both negative and positive heat fluxes become elevated. Most substantial are the increases in frictional and strain heating, which warm the bed by increasing ~80 and 150%, respectively. Conversely, greater ice flow speed associated with terminus retreat drives enhanced advective fluxes of heat away from the bed. This also increases conductive heat loss to cold overlying ice. And, the bedrock heat flux declines as the bed warms. The changes in bedrock and ice conduction naturally lag the changes in the other heat sources and sinks because these processes are driven by bed temperature. Despite decreasing heat flux from the underlying bedrock as the bed temperature increases, and increasing fluxes of heat away from the bed from advection and ice conduction, the energy balance remains positive and eventually the ice at the observation point begins to melt. Upon melting of the 178 km observation point, our gradient-based methods for computing individual heat fluxes are no longer applicable.

Fig. 7. Heat sinks (a, e), heat sources (b, f), net heat flux (c, g), and temperature at the observation point (red), and frozen/melted boundary position (purple) (d, h) during the simulation. Panels (a–d) refer to bed conditions at a point 178 km from the divide, and panels (e–h) refer to conditions at a point 160 km from the divide.

The 160 km point offers insights into changes in basal heat balance near the frozen/melted boundary over the full simulation period. The 160 km point demonstrates a similar energy balance to that at 178 km for the period of overlap, but the individual fluxes from friction and ice conduction are slightly larger in magnitude at the 178 km location (Figs 7e–h). After migration of the frozen/melted boundary stabilizes ~6.6 ka BP, the 160 km point experiences small and continuous warming as the result of a slightly positive net heat flux, not exceeding 1 × 10⁻⁵ mW m⁻². The larger heat fluxes decrease in magnitude by 22 − 44 mW m⁻² except for the bedrock heat flux, which remains roughly constant. The temperature increase of a few tenths of a degree is balanced by the increasing pressure melting point as ice thins, preventing any further migration of the frozen/melted boundary.

At the profile scale, the bed-warming resulting from driving stress and velocity changes during retreat varies substantially over the profile (Fig. 8). Within ~100 km of the ice divide, the zone with a fully frozen bed, driving stress is small and any changes in surface slope are relatively balanced by thinning ice such that the driving stress is nearly constant over the simulation period. As a result, there are minimal changes to the frictional or strain heat fluxes. Beyond ~100 km, the frictional and strain heat fluxes begin to increase. The driving stress of the profile in this region also increases and, most notably, increases with time as the ice-sheet retreats. Consequently, heat generation at the bed from frictional and strain processes also increases in this region, resulting in frozen/melted boundary migration.

Fig. 8. Driving stress (a), friction (b) and strain heating (c) during the simulation from the ice divide to 2 km past the initial frozen/melted boundary location. Profiles are displayed at 11.4, 9.9, 8.4 and 0 ka BP. The colored diamonds show the frozen/melted boundary location at the respective times during the retreat of the frozen/melted boundary.

4. Discussion

4.1. Model shortcomings

Interpretation of our results requires consideration of the modeling approach and key assumptions we have made. To start, we point to the inherent limitations of a flowline model compared to a 3-D or flowband model. Surface velocity vectors in this region of the ice sheet exhibit relatively simple flow from divide-to-margin, with no strong points of convergence to outlets and no fast-flow areas. While this supports the fidelity of our flowline model, the frictional heating may be imprecisely estimated where ice is forced over bedrock features with no ability for flow around features. Further, we have not represented complex ice flow over roughness features at the sub-1 km scale of our grid. We believe this is most problematic in high friction areas of the melted zone, but our frictional heat term must be interpreted in this context.

The spin-up and relaxation method we follow has eliminated transients in the geometry, temperature and flow fields at the start of the model run. Nevertheless, these fields are assumed to be in steady state at the start of the model run. While work from Lecavalier and others (2014) suggests there is a relatively stable terminus position from ~16 to 12 ka BP, this does not necessarily correspond to steady-state temperature, flow or even geometric conditions. This is an acknowledged shortcoming, necessary because we require initial conditions, and data constraining the terminus position and prior climate history for this location are not currently available for a longer-term spin-up. As a result, the frozen/melted boundary at the start of our run may not be accurately located. We note that our objective here is to examine change in the frozen/melted boundary over time, and we see no reason why the spin-up procedure would alter the driving processes behind frozen/melted boundary migration.

Our model assumes a well-drained bed, meaning that any meltwater produced is instantly drained away such that water effects do not impact the basal heat balance. Our focus here is on delineating where the bed is frozen vs melted, not necessarily capturing the full suite of physical processes of the melted regions. However, at the sub-gridscale transition from frozen to melted, our assumption would fail to capture potential circumstances where ice melts and refreezes. Unfortunately, we lack observational constraint on detailed bed conditions around frozen/melted transition zones. Further downstream, the radiostratigraphy lacks plume structures associated with large-scale and widespread basal freeze (Florentine and others, 2018).

Prescription of deep geothermal heat flux is based on measurements in the uppermost 500 m of the bedrock, but we then prescribe this deep geothermal heat flux at the base of our model domain. In this way, the spin-up achieves the proper bedrock heat flux to the base of the ice. However, the bedrock heat flux at the ice/bed boundary is prone to variation with time as changing ice cold content is conducted into the underlying bedrock. Our simulation results show that, at the location of the bedrock heat flux measurements near the terminus, the simulated bedrock heat flux varies from 27–31 mW m⁻², and so our prescribed deep geothermal heat flux may reduce geothermal heat fluxes by as much as up to 10%. This impacts the certainty with which we report the frozen/melted boundary migration, as our sensitivity analyses show that frozen/melted boundary migration is dependent on the initial location of the frozen–melted transition. However, importantly, our suite of simulations based on variations to the prescribed deep geothermal heat flux illustrate that the low frozen/melted boundary migration distances we present are likely to be even lower if our prescribed deep geothermal heat flux were higher.

The sliding relation used in this work inverts modern observations of ice velocity to determine the value of β ². While the sliding relation is time-dependent due to the effective pressure term, the static value of β ² means that sliding velocities may not generalize for the entire period simulated. However, the lack of a generalized sliding law for glaciology leads us to favor our approach over ad-hoc methods that trigger sliding based on the basal thermal state of the ice sheet, require extensive tuning, and often result in a large misfit between simulated and modern ice thicknesses (Martin and others, 2011). In spite of our sliding relation, the history of ice-sheet extent matches observations, and our ice thicknesses are consistent with that extent.

4.2. Frozen/melted boundary migration

During a period in which the ice margin retreated ~100 km, our simulations show a far smaller migration of the frozen/melted boundary of just ~16 km. Thus, as the divide-to-margin distance shortened from 400 to 300 km, the fraction of the profile experiencing frozen conditions at the bed actually increased by ~10%. The implication of this slow migration of the frozen/melted boundary during ice-sheet retreat is that the ice-sheet forcing on hydrologic systems does not scale constantly with ice-sheet size. Comparison of this finding with prior modeling efforts is difficult because of dissimilarities in modeling approaches and scales. Huybrechts (1996) modeled the GrIS over ice age cycles, and reported that the GrIS-wide frozen fraction decreased by less than a few percent over the last 10 ka. While their finding of decreased frozen fraction contrasts with our finding of increased frozen fraction, both results agree that a relatively large change in ice-sheet area is accompanied by a small change in frozen fraction. Elsewhere, Marshall and Clark (2002) showed that geometric changes on the Laurentide ice sheet during the build-up to and retreat after the LGM resulted in large increases and decreases in the frozen fraction at the ice bed. The scale of ice-sheet change in these simulations is much larger than the retreat we simulate in western GrIS, but underscores the complexity of heat transfer processes at the ice bed in response to changing ice geometry.

The bedrock acts as an effective heat sink in our transient model, dampening the impact of thermal perturbations from climate as they move across the bed. Although not all ice-sheet modeling studies have employed a thermally active bedrock layer (e.g. Larour and others, 2012; Seddik and others, 2012; Wang and others, 2012), we are not the first to incorporate this feature (Greve and others, 2011; Aschwanden and others, 2012; Rogozhina and others, 2012). In our case, the deep geothermal heat flux is a constant value which continuously warms the ice over the time scales we simulate. However, because the ice is far colder than the bedrock, a slight increase in the ice temperature initially flattens the temperature gradient in the bedrock, reducing the bedrock heat flux at the ice/bed boundary until the bedrock temperature profile warms to the new temperature. Our reference point experiences a 10 mW m⁻² reduction in heat flux from the bedrock as the ice warms over the simulation period. The thermal disturbance in bedrock from this 11 ka long simulation of the deglaciation phase extends around 1 km vertically, and the timescale of full recovery is longer than our simulation period. Thus, much of the heat from the warming climate over our simulation timescale is diverted into warming the bedrock rather than melting more of the bed.

Dissecting the primary controls on the heat balance at our reference point 11 km upstream from the frozen/melted boundary reveals a close balance of competing processes. The conduction away from the bed to cold overlying ice is effective at cooling the bed, nearly balancing the heat generated from friction. A consequence is that despite the fact that the strain heating at the bed is trivially small (1/500th of the frictional term), the sign of the heat balance is heavily influenced by this term. Taken together, changes in the frictional and strain heat sources through time are first-order controls on the basal heat balance. Both of these processes are strongly dependent on the ice-sheet driving stress. Our simulations show that the driving stress sensitivity to ice-sheet geometry (i.e. surface slope and ice thickness) change increases away from the ice divide (Fig. 8). Thus frictional and strain heat processes are also more prone to enhancement in response to ice-sheet surface steepening, resulting in a variable sensitivity of frozen/melted boundary migration to ice-sheet change, depending on its location relative to the ice divide.

The present day bedrock heat flux in our study domain is relatively well constrained by direct observations, and is lower in magnitude than elsewhere in Greenland as represented by commonly referenced maps (Shapiro and Ritzwoller, 2004). The frozen/melted boundary in our simulations is positioned outward from the divide by the lower deep geothermal heat flux of the region, and as noted above, our model assumptions and spin-up procedure may not precisely locate its position. Nevertheless, the processes we observe in our simulations should apply to our transect and elsewhere in Greenland: (1) antecedent cold content in the bedrock acts to buffer migration of the frozen/melted boundary as climate warms; and (2) the further the frozen/melted boundary is from the divide, the more its migration is enhanced by ice geometry changes. Current consensus shows that the location of the frozen/melted boundary, relative to the ice divide, is highly variable around the GrIS (MacGregor and others, 2016). Our findings thus imply that past and future rates of migration of frozen and melted conditions along the ice-sheet bed are likely to vary around the ice sheet and cannot be extrapolated based on a single-study transect or ice sheet-wide integration.

5. Conclusions

Our transient and thermo-mechanically coupled model with thermally active bedrock simulates the last 11.4 ka along a flowline in western GrIS, with constraints provided by prior work on ice margin chronology and climate forcing. Results suggest that during the period, the ice margin retreated about seven times more than the subglacial boundary between frozen and melted bed conditions. Partitioning the transient heat balance into constituent components reveals the processes driving the thermal evolution of the bed as climate warmed and the ice margin retreated ~100 km. The thermally active bedrock layer dampened thermal perturbations irrespective of sign, with the net effect of acting as a heat sink, reducing the net temperature increase in the frozen region of the ice sheet by 1–2 K. The 16 km retreat of the frozen area was driven by enhanced frictional and strain heating fields associated with margin retreat and steeping surface slopes toward the ice margin. Our results from this transect imply the migration rate of frozen/melted boundaries is not directly proportional to change in ice extent, and is a function of the thermal state of the bedrock and changes in geometry, which can be highly variable in time and space. This may be an important consideration for long-term changes in associated hydrologic systems.