|Home | About | Journals | Submit | Contact Us | Français|
The basal thermal state of an ice sheet (frozen or thawed) is an important control upon its evolution, dynamics and response to external forcings. However, this state can only be observed directly within sparse boreholes or inferred conclusively from the presence of subglacial lakes. Here we synthesize spatially extensive inferences of the basal thermal state of the Greenland Ice Sheet to better constrain this state. Existing inferences include outputs from the eight thermomechanical ice-flow models included in the SeaRISE effort. New remote-sensing inferences of the basal thermal state are derived from Holocene radiostratigraphy, modern surface velocity and MODIS imagery. Both thermomechanical modeling and remote inferences generally agree that the Northeast Greenland Ice Stream and large portions of the southwestern ice-drainage systems are thawed at the bed, whereas the bed beneath the central ice divides, particularly their west-facing slopes, is frozen. Elsewhere, there is poor agreement regarding the basal thermal state. Both models and remote inferences rarely represent the borehole-observed basal thermal state accurately near NorthGRIP and DYE-3. This synthesis identifies a large portion of the Greenland Ice Sheet (about one third by area) where additional observations would most improve knowledge of its overall basal thermal state.
A substantial portion of the ongoing change in the flow pattern of the Greenland Ice Sheet (GrIS) is attributed either directly or indirectly to externally forced changes in subglacial processes [e.g., van de Wal, 2008; Andrews et al., 2014; Moon et al., 2014]. Non-negligible basal sliding and deformation of subglacial till (hereafter “basal motion”) both require a thawed basal state. Hence, assessing the present thermal state of the bed of the GrIS, even at the simple level of a binary distinction between a frozen and a thawed bed, could advance our understanding of GrIS dynamics and the importance of temperature-related englacial and subglacial processes [e.g., Seroussi et al., 2013; Colgan et al., 2015; Poinar et al., 2015]. Such knowledge is essential for improving predictions of the future contribution of the GrIS to sea-level rise [e.g., Alley et al., 2005; Price et al., 2011; Nowicki et al., 2013].
Multiple borehole-temperature studies within the GrIS interior (>100 km from its margin) have generally reported a frozen ice-sheet bed [Weertman, 1968; Gundestrup and Hansen, 1984; Cuffey et al., 1995; Dahl-Jensen et al., 1998] (Figure 1; Table 1). However, these boreholes are often located at or near ice divides, where a frozen bed was predicted, so that a longer ice-core record (in time) could be recovered. Hence, these boreholes have a sampling bias toward frozen beds. NorthGRIP was an unintended but ultimately invaluable exception to this rule [Dahl-Jensen et al., 1997, 2003]. In contrast, borehole studies closer to the ice-sheet margin generally report a thawed bed [Thomsen et al., 1991; Iken et al., 1994; Lüthi et al., 2002, 2015; Harrington et al., 2015].
Because boreholes are inherently sparse compared to the areal extent of an ice sheet, thermomechanical modeling has often been applied to evaluate the thermal state of the GrIS bed [Jenssen, 1977; Funk et al., 1994; Huybrechts, 1994, 1996; Greve and Hutter, 1995; Calov and Hutter, 1996; Greve, 1997, 2005; Tarasov and Peltier, 2003; Marshall, 2005; Heimbach and Bugnion, 2009; Rogozhina et al., 2011, 2012, 2016; Aschwanden et al., 2012; Brinkerhoff and Johnson, 2013; Petrunin et al., 2013; Seroussi et al., 2013; Poinar et al., 2015]. These models generally predict what we term a “scalloped frozen core” for the GrIS, i.e., a central region where the bed is frozen, surrounded by thawed regions that can extend from the ice-sheet margin to hundreds of kilometers inland. These models typically agree with contemporaneous borehole-temperature measurements, but there is less obvious inter-model agreement elsewhere, and their basal temperature outputs have yet to be synthesized. In terms of importance for constraining the basal thermal state, these studies have variously emphasized the multi-millennial memory of ice sheets, the spatiotemporal variability of key boundary conditions (accumulation rate and geothermal flux), conservation of energy in a polythermal ice sheet, model initialization and fidelity to the flow and geometry of the modern GrIS. Inferences of basal motion from surface-velocity patterns also suggest a spatially heterogeneous GrIS basal thermal state, similar to that predicted by thermomechanically coupled models [Rignot and Mouginot, 2012; Sergienko et al., 2014].
Studies of internal and basal reflections recorded by radar sounding have also produced evidence of spatiotemporally variable basal melting, freezing and motion beneath the GrIS [Fahnestock et al., 2001; Dahl-Jensen et al., 2003; Oswald and Gogineni, 2008, 2012; Palmer et al., 2013; Bell et al., 2014; Christianson et al., 2014; Keisling et al., 2014; Wolovick et al., 2014]. With the recent exception of Rogozhina et al. , none of this radar-derived information regarding the basal thermal state has been synthesized with thermomechanical models of the entire GrIS, leaving a gap between contributions to our knowledge of its basal thermal state from radar and modeling.
Given the range of existing estimates of the present basal thermal state of the GrIS, both additional approaches to resolving this state and a synthesis of existing estimates are warranted. In this study, we produce multiple new estimates of the GrIS basal thermal state and synthesize these results with earlier estimates from three-dimensional (3-D) thermomechanical ice-sheet models. Our new estimates are derived from radiostratigraphy-constrained ice-flow modeling and analysis of surface velocity and imagery. Our goal is to map the portions of the GrIS where these methods agree that the bed is frozen or thawed, versus where there is poor agreement between these methods.
To quantitatively evaluate basal temperatures for the entirety of an ice sheet, numerical thermomechanical modeling is required. Such modeling explicitly solves the coupled mass-, momentum- and energy-conservation equations over the entire ice sheet, given certain boundary conditions. We evaluate the basal temperature outputs from the eight 3-D thermomechanical models of the GrIS included in the Sea-Level Response to Ice Sheet Evolution (SeaRISE) effort described by Nowicki et al. . That study also described the models’ varied resolutions, initializations, thermodynamic representations and boundary conditions in detail. More 3-D model representations of the present basal thermal state of the GrIS exist than this ensemble of SeaRISE models alone. We opted to use the SeaRISE ensemble only because it enables a straightforward evaluation of multiple models to which similar modern boundary conditions were applied. Several tradeoffs presently exist between different types of models, such as those that are initialized across glacial–interglacial cycles (e.g., SICOPOLIS; Rogozhina et al. ) versus those that assimilate modern observations (e.g., ISSM; Seroussi et al. ). We specifically set aside the myriad issues facing such comparisons and focus instead on simply evaluating the present agreement between a relatively large number of models.
We consider only the SeaRISE basal temperature fields at the end of their 100-year control runs (“CC” in the nomenclature of Nowicki et al. ), because these fields represent relaxed states characteristic of present day and with no prescribed additional forcings. We correct all modeled basal temperatures (Tbed) to be relative to the pressure-melting point using the ice-thickness field for each model at the end of its control run and a value of 8.7 × 10−4 K m−1 for the rate of decrease in the pressure melting point with increasing ice thickness [Cuffey and Paterson, 2010].
For three of the SeaRISE models (Community Ice Sheet Model, CISM; Ice Sheet System Model, ISSM; Parallel Ice Sheet Model, PISM), we use subsequent revisions to the models described below. All models use an enthalpy formulation for conservation of energy that accounts for the latent heat of liquid water within temperate ice (i.e., ice at the pressure-melting point) [Aschwanden et al., 2012]. For the instance of CISM employed here (v2.0), the momentum balance is still based on the 3-D first-order Stokes-approximation (the “Blatter-Pattyn” approximation), but the equations are now formulated following the variational approach of Dukowicz et al.  and discretized using finite elements on a fixed, regular grid. Model tuning follows that for SeaRISE [Price et al., 2011], but the spin-up has been extended to last 350 ka to better approximate equilibrium englacial temperatures. The instance of ISSM we use here is the same as the steady-state run described by Seroussi et al. . This model also used the SeaRISE datasets for several boundary conditions, but adjusted geothermal flux to better match borehole constraints. This model instance assumes thermal steady state, although its englacial temperatures compare favorably with radar-inferred depth-averaged temperatures [MacGregor et al., 2015b]. The instance of PISM we use here was initialized using the “paleo-climate” method described by Aschwanden et al. , and the prescribed geothermal flux follows Shapiro and Ritzwoller , as for the SeaRISE experiments.
The local rate of basal melting and/or motion necessary to explain radar-observed depth–age relationships (dated radiostratigraphy) can be inferred from ice-flow modeling of varying degrees of complexity [e.g., Fahnestock et al., 2001; Dahl-Jensen et al., 2003; Keisling et al., 2014; Wolovick et al., 2014; Wolovick and Creyts, 2016; Koutnik et al., 2016]. Hence, analysis of dated radiostratigraphy can inform our understanding of the basal thermal state across large regions of an ice sheet.
Here we extend the one-dimensional (1-D; vertical), steady-state ice-flow modeling approach of Fahnestock et al.  to produce estimates of the local basal melt rate , basal shear layer thickness h and shape factor that best match Holocene radiostratigraphy across the GrIS. These three quantities are derived from two analytical 1-D ice-flow models (Nye + melt, Dansgaard–Johnsen) that make distinct assumptions about flow within the ice column. Figure 2 illustrates these two models schematically. To constrain these 1-D models, we use the same GrIS-wide dated radiostratigraphy developed and described by MacGregor et al. [2015a] and further investigated by MacGregor et al. [2015b, 2016]. Figure 1a shows the spatial coverage of the radar-sounding data.
To estimate , we use the “Nye + melt” ice-flow model described by Fahnestock et al. , which predicts the following depth–age relationship:
where a(z) is the age a of the ice at ice-equivalent depth z, H is ice thickness and is the surface-accumulation rate. In this first 1-D ice-flow model, the ice column is assumed to deform in pure shear only (uniform horizontal velocity with depth).
The shape factor is the ratio between the depth-averaged horizontal speed of the ice column ū and its surface speed us [Cuffey and Paterson, 2010]:
Where approaches unity, vertical shear in the ice column is negligible, therefore the rate of basal motion ub must approach us and the ice deforms in pure shear, as in the Nye + melt model. Where < 1 ice flow must be also accommodated by simple shear [Cuffey and Paterson, 2010]. Assuming that the Glen’s flow law exponent n = 3 and a uniform rate factor for the ice column, a lower limit for is approximately (n + 1)/(n + 2) = 0.8. Following MacGregor et al. , the geometric relationship between and h is
We estimate h (and hence ) using a second 1-D ice-flow model, introduced by Dansgaard and Johnsen , which predicts the following depth–age relationship:
In this second 1-D ice-flow model, the ice column is assumed to be deforming in pure shear only where z ≤ H − h (uniform vertical strain rate), but also in simple shear where z > H − h, i.e., within the basal shear layer (linearly decreasing vertical strain rate).
To estimate either or h (and hence ) from a dated radiostratigraphy using the above models, we must also constrain , a surface boundary condition of broad glaciological value but which only indirectly and slowly affects the basal thermal state. For a plausible range of values, both models are overall more sensitive to than to or h. Both 1-D ice-flow models are explicitly steady-state. Hence, following Fahnestock et al.  and MacGregor et al. , we restrict these models to only use reflections that are 9 ka old or less (the last three quarters of the Holocene epoch), when was stable across millennial time scales. Following MacGregor et al. , we require at least four dated reflections within each 1-km along-track bin to estimate and h as the best-fit parameters from a non-linear unconstrained minimization of χ2 statistic:
where N is the number of Holocene-dated reflections, and are the age and age uncertainty of the ith dated reflection, respectively, and is the modeled age of that same reflection, calculated using its depth and either Equation 1 or Equation 4. Confidence bounds (95%) for both model parameters are estimated using Δχ2 distributions. We grid along-track values using ordinary kriging with parameters similar to MacGregor et al. [2015a, 2016].
These two 1-D ice-flow models use the same data to produce distinct but related inferences of basal melting or motion. For example, a negative value of h (and hence > 1) is non-physical, indicating that is positive (implying basal melting) and possibly greater than [Fahnestock et al., 2001]. While these ice-flow models are relatively unsophisticated and non-unique, they are also relatively simple to evaluate and interpret. Dahl-Jensen et al.  introduced a more sophisticated 1-D ice-flow model that is essentially a combination of the above two models (“Dansgaard–Johnsen + melt”). While this model is more physically complete, we do not consider it here as we found separately that simultaneously solving for its four parameters (, , h and ub/us) results in less well-constrained solutions than the two models that we consider here.
Separate from the steady-state assumption, 1-D modeling of radiostratigraphy is not appropriate for the entirety of an ice sheet, because horizontal gradients in ice flow can alter locally observed depth–age relationships substantially [Waddington et al., 2007; Koutnik et al., 2016]. Hence, the strain history of the particles that form observed isochrones may differ from their apparent history at their present location. Following MacGregor et al. , we restrict our interpretation of radiostratigraphy-inferred values of , h and to the portion of the GrIS where we consider the local-layer approximation to be acceptable for reflections younger than 9 ka, i.e., the region where depth–age relationships may be represented reasonably by 1-D models that neglect horizontal gradients in ice flow. This region encompasses the majority of the GrIS (71% by area; Figure 5c,d).
Fast ice flow is widely accepted as an indicator of basal motion and hence of a thawed bed. Because creep deformation also contributes to ice flow, us is not equal to ub. A simple estimate of the maximum flow speed due to internal ice deformation only is therefore valuable for constraining where basal motion is likely contributing to the observed pattern of us and hence where the bed is thawed. We generate this estimate by calculating the surface speed due to ice deformation (udef) for an ice column that is entirely at the pressure-melting point (i.e., temperate), following Cuffey and Paterson  as
where Ā is the rate factor for temperate ice, τd = ρicegHα is the gravitational driving stress (Figure 1b), ρice is the mean density of the ice column (900 kg m−3), g is the acceleration due to gravity (9.81 m s−2) and α is the along-flow surface slope. For H, we use the gridded field and its uncertainty from Morlighem et al. . We assume that the uncertainty in udef is due to uncertainty in Ā and H only and propagate their uncertainties into lower- and upper-bound estimates of udef.
For temperate ice, a synthesis of reported values by Cuffey and Paterson  suggests Ā = 2.4 × 10−24 Pa−3 s−1, which is consistent with borehole-deformation measurements of the GrIS [Ryser et al., 2014a]. For polar ice experiencing simple shear, they also suggest an enhancement factor E ≥ 2 to account for fabric development, and for Greenland ice from the Last Glacial Period, E ≈ 3 [Paterson, 1991; Lüthi et al., 2002]. Such ice can greatly influence ice flow [e.g., Ryser et al., 2014a], although it generally constitutes a smaller portion of the ice column [MacGregor et al., 2015a]. Hence, we ignore E in our calculation of udef, but assign Ā a relative uncertainty of 25%.
Rignot and Mouginot  also calculated udef for the GrIS and determined that, when attempting to match us and udef in the northern GrIS interior, the best-fit value for Ā was 1.2 ± 0.2 × 10−24 Pa−3 s−1. This value is half of that which we assume and equivalent to a depth-averaged temperature of −3.8°C. However, our approach does not require the assumption that any particular region of the ice sheet is frozen at its base. Of course, the whole of the GrIS is not temperate (e.g., Figure 1), and the 3-D models also produce estimates of the englacial temperature. The key advantage of our particular approach, as compared to more sophisticated thermomechanical modeling, is simplicity of interpretation. We consider the large-scale flow field of the whole ice sheet only, rather than attempting to interpret the detail and significance of driving stress variability in any given region (e.g., Figure 1b, Sergienko et al. ).
To determine the direction of α, we first use a composite map of Greenland surface velocity (I. Joughin, pers. comm., 2015) developed using the same methodology described by Joughin et al. , i.e., it is derived from a combination of interferometric synthetic aperture radar (InSAR) analysis and speckle tracking of spaceborne imagery. This composite field is derived from 1995–2013 SAR data, and it is the same as that used by MacGregor et al. . In the slower flowing interior, ionospheric noise limits the quality of satellite-derived flow azimuths. Therefore, in this region (≤100 m a−1) we progressively weight the surface-velocity azimuth towards the direction of the surface-elevation gradient of the Greenland Ice Mapping Project (GIMP) digital elevation model [Howat et al., 2014]. This weighting depends exponentially on the magnitude and relative uncertainty of the local surface speed. We rotate the projected surface-velocity field onto this weighted flow azimuth, which improves representation of ice flow in the GrIS interior. This weighted flow azimuth is used to determine the magnitude of α.
Fast basal motion can effectively transmit large-scale (~1–20 H horizontally) bedrock topographic features into the ice sheet [Gudmundsson, 2003; Hindmarsh et al., 2006; Bell et al., 2007; De Rydt et al., 2013; Ryser et al., 2014b]. At the subaerial ice surface, this basal motion can manifest as transverse-to-flow surface undulations. Hence, the surface texture of an ice sheet is related indirectly to the slip ratio s = ub/udef. Where s >> 1 (equivalent to → 1), the surface texture is expected to be much rougher than where s ≤ 1 [De Rydt et al., 2013] and is an indirect indicator of a thawed bed.
The surface texture of the GrIS was explored by Scambos and Haran  using an Advanced High Resolution Radiometer (AVHRR) mosaic and in greater detail along a few transects by van der Veen et al. . Here we use the Moderate Resolution Imaging Spectroradiometer (MODIS) Mosaic of Greenland (MOG) to analyze the texture of the subaerial ice surface of the GrIS [Haran et al., 2013]. MODIS imagery is well suited to surface-texture analysis because of its relatively high spatial and radiometric resolution, and because of the high-pass filtering that can be applied to this imagery to enhance visibility of surface texture. MOG is posted at 100 m, with an effective resolution of ~200 m, so it offers a ~5-fold improvement in spatial resolution over the AVHRR mosaic studied by Scambos and Haran . We examine MOG at a small scale (~1:200,000) with a narrow, localized contrast stretch (~8-bit) and delineate the ice-sheet transition between a relatively smooth surface and prominent surface undulations. We assume that this transition in surface texture represents a boundary where s >> 1 downstream. This transition is traced twice, once as described above (“standard”) and a second time using a stricter visual criterion for the onset of surface undulations (“conservative”).
We synthesize the multiple methods of constraining the basal thermal state across the GrIS described above by simply assessing where each of the above independent methods produces a clear signal regarding this state. Table 2 lists the bounds used to discriminate between a frozen and thawed bed for each method. We initialize a 5-km gridded ice-sheet mask S to zero. For each method at each grid point, if a signal exists for a frozen (thawed) bed, then −1 (+1) is added to S. Because the radiostratigraphic constraints are not independent of each other, is the only radiostratigraphic constraint used to determine S.
If a given method does not yield an unambiguous signal regarding the basal thermal state, then S is not adjusted there based on that method. For example, where [(n + 1)/(n + 2)] < ϕ < 1, does not clearly distinguish between a frozen and a thawed bed. We do not weight any of the methods with respect to each other. Prior to but following the same procedure as for S, a separate mask is generated using the 3-D thermomechanical model outputs, each weighted equally. In this manner, each independent method contributes to an unbiased synthesis of the GrIS basal thermal state.
Based on confidence bounds or uncertainty estimates for each of the four methods described above and their discriminating characteristics (Table 2), two additional instances of S are generated: a cold-bias instance and a warm-bias instance (Scold and Swarm, respectively). We then generate a new mask (L) that synthesizes the agreement between the different methods and represents the likely thermal state of the bed of the GrIS. This mask is somewhat analogous to the analysis performed by Pattyn  using multiple instances of a thermomechanical model of the Antarctic Ice Sheet. L is also initialized to zero and then assigned −1 (+1), representing a frozen (thawed) bed where at least two of the three instances of S agree on the basal thermal state (sign of S), regardless of their degree of agreement in this state (magnitude of S). If only two instances of S agree, then the assignment is made only if the other instance does not suggest the opposite basal thermal state.
We assume that regions of the bed where L = 0 (uncertain) that are surrounded by likely thawed or frozen bed (L = −1 or +1, respectively) are more likely than not to possess the same basal thermal state as their surroundings. Following this reasoning, we reassign uncertain “holes” less than 10 grid cells in size (≤250 km2) with their surrounding basal thermal state. Similarly, in regions where L = 0, we reassign likely frozen or thawed “holes” of the same limited size to to L = 0.
To qualitatively evaluate the four methods’ inferences of the basal thermal state of the GrIS and our synthesis, we use the temperature–depth profiles measured in the six deep boreholes that have associated ice cores (Camp Century, DYE-3, GISP2, GRIP, NEEM and NorthGRIP). Additional full-thickness borehole-temperature profiles exist for the GrIS from its southwestern margin and elsewhere along the margin, and these are also included in our analysis. Figure 1 maps the locations of these boreholes in Greenland and Table 1 summarizes their key features, including their measured and pressure-melting-adjusted (“corrected”) basal temperatures (i.e., relative to the local pressure-melting point). Subglacial lakes are taken as indicators of a thawed bed, and the locations of the few known subglacial lakes beneath the GrIS are also shown in Figure 1. For completeness, we also show comparable inferences of basal temperatures for Greenland’s peripheral ice caps.
To evaluate radiostratigraphic inferences of basal motion, measured depth profiles of horizontal velocity are required. These profiles are derived by integrating borehole-inclinometry data collected at three boreholes: Camp Century [Gundestrup et al., 1993], DYE-3 [Dahl-Jensen and Gundestrup, 1987] and GISP2 [Bender et al., 2010]. Borehole-deformation studies by Lüthi et al.  and Ryser et al. [2014a], within and near Jakobshavn Isbræ, are too far (>50 km) from radiostratigraphic inferences of basal motion for such a comparison to be reliable.
Figure 3 shows the eight 3-D thermomechanical model estimates of GrIS basal temperatures and Figure 4 synthesizes these estimates. Figure 3 shows that there is a diversity of 3-D estimates of GrIS basal thermal state, but because they are weighted equally when synthesized, that diversity is muted in Figure 4. From Figure 3, there is no clear relationship between the predicted basal thermal state and whether the model assimilates modern observations or is initialized using various paleoclimatic forcings.
The region of largest agreement that the bed is thawed is in southwestern Greenland (ice-drainage systems – hereafter IDS – 6.2, 7.1 and 7.2, following Figure 1b), which also has the most extensive ablation zone, but this agreement also extends well into its accumulation zone. These models also predict that two southeastern drainage systems (IDS 4.1 and 4.2) are mostly thawed, as well as the northwestern portion of the GrIS that faces Baffin Bay (IDS 7.2 and 8.1) and the main trunks of Humboldt and Petermann Glaciers in northwestern Greenland (IDS 1.1). The large ice-drainage system (IDS 2.1) that includes the Northeast Greenland Ice Stream (NEGIS; hereafter the “NEGIS system”) is sometimes predicted to be thawed, but not as consistently as the aforementioned regions.
Only the northwestern extensions of the central ice divides are predicted consistently to be frozen to the bed, particularly between Camp Century and NEEM. Elsewhere, ice divides are generally predicted to be frozen, but not uniformly so, as in the vicinity of NorthGRIP and DYE-3. The northernmost reach of the GrIS is expected to be frozen (most of IDS 1.2, 1.3 and 1.4), as is the central half of the eastern GrIS coastline (IDS 3.1) and a substantial portion of the interior west of the central ice divide.
The NEGIS system and the GrIS south of ~65°N are the regions of greatest uncertainty in the 3-D modeled basal thermal state. Elsewhere, the location of the along-flow transition from a frozen to a thawed bed (if there is one) is rarely well constrained and often spans >100 km in the along-flow direction.
Figures 5 and and66 show the 1-D model estimates of , h and across the GrIS. These quantities are related to each other, as expected [Fahnestock et al., 2001]. Rather than restrict calculation of to regions where h < 0, as was done by Fahnestock et al. , here we simply estimate all quantities where sufficient traced Holocene radiostratigraphy exists and restrict our interpretation to the region where 1-D modeling of such radiostratigraphy is acceptable.
Evidence of significant apparent basal melting ( > 0) is restricted primarily to three portions of the GrIS. The largest of these three regions includes most of the NEGIS system, first identified by Fahnestock et al.  using nearly the same methodology. The other two smaller regions of significant apparent basal melt are in the northwest (NW; IDS 1.1) and southwest (SW; IDS 6.2) ice-drainage systems. They are of comparable size and have not been identified previously as possessing significant basal melting. Both regions are relatively near (<100 km) to frozen boreholes.
The NW region of basal melting is north of the ice divide along which both Camp Century and NEEM are located. It is composed of two sub-regions, one of which is approximately equidistant from those two boreholes and is well constrained by the spatial coverage of the radiostratigraphy. The other region is farther northwest and reaches the marginal limit of our 1-D modeling. Gridding suggests that this basal-melting signal nearly reaches the ice margin (Figure 5d), but we consider this inference to be speculative because of the sparse along-track coverage (Figure 6a). The SW region of basal melting is west of DYE-3, on the other side of the central ice divide. It is also remarkable because the basal temperature at DYE-3 was <−10°C (Table 1), <50 km southeast of this region’s eastern edge.
We find that < 0 in the immediate vicinity (≤ 5 km) of all interior boreholes. This pattern is qualitatively consistent with borehole-measured basal temperatures, except for NorthGRIP. We infer basal melting at a rate of less than 1 cm a−1 within 15 km of NorthGRIP, but not at the borehole itself. This result differs slightly from that of Dahl-Jensen et al. , who used a Dansgaard–Johnsen + melt model in conjunction with an along-divide radar-sounding transect and found that typically exceeded 0.5 cm a−1 along the flowline immediately upstream of NorthGRIP, which agreed with borehole observations. This discrepancy emphasizes the limitations of 1-D ice-flow modeling and the need to focus on larger scale patterns in these 1-D inferences of (and hence also h and ).
Although two new regions of significant apparent basal melting are found, < 0 for the large majority (87%) of the portion of the GrIS where we calculated this quantity, potentially suggesting widespread basal freeze-on (Figure 5a,d). While basal freeze-on has been inferred for several portions of the northern GrIS, including where we find < 0 [Bell et al., 2014; Wolovick et al., 2014], it is not anticipated to be as widespread as implied by these calculations for two reasons: 1. The model’s assumption of plug flow; 2 Because thermodynamic considerations and likely geothermal flux values (~55 mW m−2; Shapiro and Ritzwoller ) imply that widespread and sustained basal freeze-on at high rates (>5 cm a−1) is unlikely. The likely cause of the widespread inference of < 0 is that the magnitude of the mean vertical strain rate in the 9–0 ka portion of the ice column is greater than that of a Nye model where = 0. The simplest explanation for this pattern is that sufficient vertical shearing is occurring deeper within the ice column, as in a Dansgaard–Johnsen model, that the presumably low rate of basal melting or freeze-on beneath most (but not all) of the ice-sheet interior cannot be well constrained by modeling of isochrones younger than 9 ka.
The pattern of h (Figure 5b,e; Figure 6b–f) shows that vertical shear within the GrIS is remarkably heterogeneous. Regions where h < 0 and > 0 are well correlated, as expected. Dansgaard and Johnsen  originally developed their ice-flow model for Camp Century and estimated h = 400 m there. We find similar values in the vicinity of Camp Century, but also substantially larger values farther inland and nearby (<50 km away). Variability in h is expected, but has not been considered previously at this spatial scale. Following Waddington et al. , h ≈ 0.7 H is anticipated in the vicinity of an ice divide and h ≈ 0.25 H is anticipated elsewhere in the interior, where the ice sheet is in “flank” flow. We find that h is often larger at ice divides, but not uniformly so, particularly in the southern GrIS. At the confluence of multiple ice divides, >100 km east of NEEM, h does approach the expected range for divide flow (>2000 m or >0.5 H). The pattern of flank flow is variable and generally differs between adjacent IDSs. The basal shear layer thickness tends to be larger west of the central ice divides, except in regions of apparent basal melt.
The shape factor is related directly to h (Equation 3). As such, the GrIS-wide pattern and interpretation of closely resemble that of h and (Figures 5c,f and 6g–i). The advantage of considering over the other two 1-D modeled parameters is that interpretation of its values can constrain the location of both frozen and thawed beds, as opposed to only thawed beds. This analysis shows that the largest contiguous region of ambiguous basal thermal state is centered on the central ice divide in the southern GrIS, encompassing the uppermost reaches of IDS 4.1, 6.2 and 7.1. Elsewhere, the transition from an apparently frozen bed to a thawed one is often abrupt, occurring across a region of less than 50 km.
Figure 7a shows a histogram of values and conventional glaciological interpretations thereof. Where ≥ 1 indicates that basal motion is likely occurring and that the Dansgaard–Johnsen model’s assumption that basal melting is negligible is no longer valid. Conversely, if we assume that n = 1 (a potential lower limit value at ice divides; Pettit and Waddington, 2003]), then where < 0.67 ([(n + 1)/(n + 2)]) indicates that basal freeze-on is likely occurring. The gridded distribution is skewed left compared to the along-track distribution, likely because of the sparser radar coverage in regions of low along-track values (e.g., western interior). Both distributions peak above 0.67 but below 0.8 ([(n + 1)/(n + 2)] assuming n = 3). Assuming that a large but unspecified portion of the GrIS interior is indeed frozen and that its flow is adequately represented by the Dansgaard–Johnsen model, this observation implies that the effective value of n for the GrIS interior is somewhat less than 3. This result is consistent with the assumed rheology embedded in most ice-flow models, including Equation 6. Excluding observations near ice divides, where n is less clearly approximated by 3 [e.g., Pettit and Waddington, 2003; Gillet-Chaulet et al., 2011], does not affect these distributions significantly.
A comparison between borehole-measured values and those we infer from the radiostratigraphy (Figure 7b) shows that we tend to underestimate relative to borehole observations. This underestimation is most apparent at GISP2 (borehole: 0.85; radiostratigraphy: 0.77 ± 0.01), which is frozen at the bed. The most likely explanation for this underestimation is that ice flow was non-steady over the period represented by the entire ice column, as opposed to that which we infer from the 9–0 ka portion of the ice column only. A non-uniform rate factor in situ may also influence this difference. Given the paucity of borehole-measured values, we cannot perform an ice-sheet-wide correction to these values. Hence, while we use values to infer the basal thermal state from radiostratigraphy, we acknowledge that it may underestimate (overestimate) the portion of the bed that is thawed (frozen).
The filtered surface-velocity us, the temperate-column estimate of udef and the ratio of these two quantities (us/udef) are shown in Figure 8. Lower- and upper-bound estimates of us/udef = 1 are determined using the reported uncertainty for us and our modeling uncertainty for udef.
Except for NEGIS and small regions along the central ice divides, where flow azimuths are not as well constrained, this analysis infers the presence of basal motion within ~50–250 km of most of the GrIS margin. This pattern is clearer away from ice divides that reach the margin, and us/udef tends to be greater along the western margin than the eastern margin. This analysis infers basal motion at NorthGRIP, where the bed is thawed, but not at the other interior boreholes.
The along-flow uncertainty in the region where us/udef > 1 generally exceeds 200 km, especially in southwestern (IDS 6.2) and northwestern Greenland (IDS 1.1 and 8.1), near radiostratigraphically-identified regions of basal motion (Figure 5c,f). This large uncertainty is due primarily to uncertainty in ice thickness in the ice-sheet interior and to a lesser extent by the assumed relative uncertainty in the rate factor for temperate ice.
Figure 9 shows the contrast-stretched MOG and our traced outlines of the onset of surface undulations, based on both a standard and a conservative assessment of MOG surface texture. Several zoom-ins of MOG are also shown (same regions as Figure 6). Although NEGIS is a prominent feature of the surface of the GrIS and its shear margins are well defined [Fahnestock et al., 1993], in terms of the onset of surface undulations, it is otherwise undistinguished compared to elsewhere on the GrIS. The central ice divide is smooth, but the deep interior of the southern GrIS is relatively rough and the onset of surface undulations reaches the ice divide below 65°N, near DYE-3. In general, the standard and conservative estimates of the surface-undulation onset agree well, except in two relatively small regions southwest of the Summit region (GRIP/GISP2) and north of NEGIS.
Figure 10 shows the agreement between the four estimates of the GrIS basal thermal state, along with cold- and warm-bias versions of this agreement. Only thermomechanical models can estimate where the bed is both frozen and thawed for the entire ice sheet. Analysis of can constrain the location of frozen and thawed beds, but not for the entire ice sheet. Our surface-velocity and surface-imagery analyses constrain only where basal motion is likely. Hence, any synthesis of these four methods is biased toward inferring a thawed bed.
Figure 10a illustrates the variety of estimates of the frozen/thawed transition generated by the four methods. Although the large-scale structure of all these estimates broadly resembles a scalloped frozen core, the details of this structure and its total extent vary significantly. Nevertheless, in some portions of the GrIS, multiple methods agree well at small scales (<100 km), e.g., the frozen/thawed transition in west-central Greenland (IDS 7.2) as estimated independently by both the 3-D models and MOG surface undulations.
Table 3 summarizes the portion of the GrIS bed that is identified as frozen or thawed by each method. The SeaRISE synthesis suggests that a greater portion of the GrIS bed is frozen (56–80%) than that inferred from radiostratigraphy (19–41%). Only 10% of the GrIS is identified as thawed from radiostratigraphy, and this range is narrow (9–12%). The surface-velocity analysis has the greatest range of inferred thawed bed (18–76%), whereas the surface-texture analysis produces the largest standard estimate of the thawed extent of the GrIS bed (69%).
For both the standard and cold-bias agreement masks (Figures 10b and 10c, respectively), the extent of the scalloped frozen core is similar. This similarity occurs because only two of the methods can constrain indirectly where the bed is frozen. The primary difference between these two masks is the extent of the uncertain region (no or low agreement). The warm-bias mask (Figure 10d), which has a smaller scalloped frozen core, is more distinct from the standard mask and has generally greater agreement regarding the location of a thawed bed. For all agreement masks, most interior boreholes are represented correctly as frozen. However, there is no agreement regarding NorthGRIP’s basal thermal state for any mask, and in the standard and warm-bias masks it straddles the boundary between large regions of agreement that is bed is frozen (to the west of NorthGRIP) and thawed (to the east). Further, for all agreement masks, DYE-3 is identified as thawed, rather than frozen.
Figure 11 shows the likely basal thermal state of the GrIS, based on the above agreement and its confidence bounds. We find that 43% of the bed is likely thawed, 24% is likely frozen and the thermal state of the remainder (34%) is uncertain (Table 3). The hole-filling operation adjusts the size of these regions by less than 0.5% each. The distribution of each region amalgamates many of the characteristics described above: 1. The largest contiguous regions of likely thawed bed are in southwestern Greenland (primarily IDS 6.1, 6.2 and 7.1), along the northwestern coast (IDS 8.1) and within the NEGIS system (IDS 2.1), especially NEGIS itself. 2. The portion of the GrIS bed that is likely frozen lies between 68°N–80°N and is generally west of the central ice divides. 3. Smaller discontinuous marginal regions associated with major outlet glaciers are likely thawed (e.g., Helheim Glacier in IDS 4.1; Humboldt and Petermann Glaciers in IDS 1.1). 4. The thermal state of the bed is uncertain within large contiguous swaths of the GrIS interior that can be >100 km wide along-flow. 5. The majority of the ablation zone is identified as likely thawed (65%) and none of it is likely frozen.
Figure 12 shows the relationships between the likely basal thermal state and several fundamental ice-sheet properties. While the uncertain region complicates interpretation of these relationships and no property provides an unambiguous binary distinction between a frozen and a thawed bed, some noteworthy patterns emerge. Both increasing surface elevation and ice thickness result in a greater likelihood of a frozen bed, while increasing surface slope and speed imply that a thawed bed is more likely. Surface elevation and ice thickness co-vary, so the consistency of their relationships is unsurprising, as is the relationship between likely basal thermal state and surface speed. Above a surface slope of 50 m km−1 or speed of 100 m a−1, a frozen bed is unlikely. While driving stress is related to both surface slope and ice thickness, either thermal state is possible across the range of inferred values. For example, between 30–70 kPa, a frozen bed is more likely, but for values outside of this range a thawed bed is more likely. A similar pattern occurs for modeled surface mass balance, but its distribution is more complex. Between 15–30 cm a−1, a frozen bed is more than twice as likely as a thawed bed, but below 10 cm a−1, a thawed bed is most likely.
For ice-sheet wide applications, no single method considered here determines the thermal state of the bed of the ice sheet in a manner that is clearly more reliable than all others. 3-D thermomechanical models include more complex physical representations and more boundary conditions, but not all physical processes or boundary conditions are yet sufficiently well understood. Further, there is a commensurate increase in the challenge of interpreting results from such models. Our synthesis (Figures 10 and and11)11) demonstrates that significant uncertainty remains in the present basal thermal state of the GrIS, even with respect to ice-sheet properties that are often considered diagnostic (e.g., surface velocity, Figure 12). This uncertainty must be considered when evaluating processes that rely substantially on that state [e.g., Colgan et al., 2015]. In certain regions (e.g., southwestern Greenland, NEGIS, the onset regions of major outlet glaciers), we can be reasonably confident that the bed is thawed contiguously at the scales considered here (≥5 km), supporting the conclusions of earlier studies related to those regions [e.g., Fahnestock et al., 2001; Poinar et al., 2015; Tedstone et al., 2015; Rogozhina et al., 2016]. Confidence regarding a frozen bed is restricted mostly to the vicinity of boreholes and a contiguous region in northern and central Greenland.
To date, the majority of boreholes drilled through the GrIS are located where our synthesis suggests that the thermal state of the bed is reasonably well constrained (Figure 11). That argument is somewhat circular, because some of the 3-D thermomechanical models were tuned to match contemporary borehole observations. Only NorthGRIP and DYE-3 are located within regions that remain poorly constrained. The standard hypothesis explaining these discrepancies is that an incorrect geothermal flux in these regions leads to an incorrectly predicted basal thermal state, a problem that can be addressed by adjusting the geothermal flux [e.g., Greve, 2005; Seroussi et al., 2013; Rogozhina et al., 2016]. While certainly plausible, this hypothesis can be invoked only for the 3-D thermomechanical models, which directly apply the geothermal flux as a boundary condition, and not for the other three methods included in our synthesis. Yet, those other methods also produced varied responses in the vicinity of these two boreholes. This pattern suggests that the underlying cause of the NorthGRIP/DYE-3 discrepancy stems from more than just poorly known geothermal flux, and that some other property or process must also be invoked (e.g., spatially varying basal friction or rheology). Within ~200 km of NorthGRIP, this extended hypothesis is qualitatively consistent with the results of Christianson et al. . They reported observations of dilatant till within NEGIS but not outside of it, emphasizing the likelihood of spatially variable basal properties within the NEGIS drainage system.
Our synthesis evinces a pattern in northeastern Greenland that may have wider implications for our understanding of GrIS thermodynamics. The transition between a likely frozen bed and an uncertain basal thermal state follows the western and northern boundary of the NEGIS system (its interior ice divide) remarkably well (IDS 2.1, Figure 11). Both 3-D models and radiostratigraphic analysis agree regarding the large-scale structure of this transition, if not its exact position (Figure 10). There is also a large gradient in the Holocene-averaged accumulation rate across this ice divide [MacGregor et al., 2016], suggesting a strong coupling there between local accumulation rate and basal thermal state. Such a coupling could have arisen in at least two ways. First, ice in this region flows slowly (<10 m a−1, Figure 8a), so if this region has been quasi-stable for tens of millennia, temperature–depth profiles may be approximated using a 1-D steady-state model [Cuffey and Paterson, 2010]. For a given geothermal flux and surface temperature, such models predict that a lower accumulation rate results in a greater likelihood of a thawed bed. Second, a heretofore-unrecognized boundary in Greenland’s subglacial geology could exist at or near this ice divide [Dawes, 2009; Rogozhina et al., 2016]. This boundary could influence the local basal thermal state through a horizontal gradient in the geothermal flux, potentially leading to basal motion within the NEGIS system and ice drawdown, thus influencing the present position of the ice divide as determined from surface topography [Zwally et al., 2012].
The large region of uncertain basal thermal state also indirectly informs our understanding of ice-sheet thermodynamics. Undoubtedly, part of the extent of the uncertain region is due to artifacts or limitations in our methodology and synthesis. Nevertheless, our synthesis suggests that the basal thermal state typically transitions from frozen to thawed over a relatively wide region (often >100 km along-flow). The along-flow width of this uncertain region represents the distance over which ice must flow before its modeled state, englacial and surface properties consistently indicate a thawed bed. Such a transition likely exhibits large spatiotemporal variability that is only partly tied to local basal properties, in particular the geothermal flux [van der Veen et al., 2007; Brinkerhoff et al., 2011; Meierbachtol et al., 2015], and can also be influenced by model resolution in some regions [Aschwanden et al., 2016]. Heat advection (whether englacial or subglacial) and dynamic thinning or thickening can modify the local basal temperature gradient significantly at centennial time scales. In turn, the frozen/thawed boundary migrates. Such processes may be rate-limited and possess negative feedbacks [e.g., Huybrechts, 1996; Phillips et al., 2013; Wolovick et al., 2014]. Regardless of the specific factors and history that have led to the present GrIS basal thermal state, our synthesis emphasizes that the frozen/thawed transition is unlikely to be a singular, stationary, contiguous line around the ice sheet.
Evaluations of the role of the subglacial hydrologic system in the connection between surface melting and seasonal acceleration often hinge on the assumption that the bed is thawed, as is consistently observed by boreholes in the ablation zone and supported by our synthesis (Figure 11). However, virtually all such field studies for the GrIS have so far taken place in southwestern Greenland (IDS 6.1, 7.1 and 7.2; Table 1), within the largest contiguous thawed region of the bed. Many ground-based, remote-sensing and modeling studies also focus on this region [e.g., Price et al., 2008; Andrews et al., 2014; Poinar et al., 2015; Tedstone et al., 2015] or other extensively thawed regions, such as the NEGIS system [Karlsson and Dahl-Jensen, 2015]. While recent perturbations to local subglacial hydrology may be significant along the margin there [van de Wal et al., 2008], a previously unrecognized bias may be introduced by applying process-level insights from this particular region to the whole of the GrIS [e.g., Shannon et al., 2013; Tedstone et al., 2015]. Elsewhere along the margin of the GrIS, the bed is less clearly thawed, and rarely over such a large portion of an IDS. If the upstream portion of a given IDS is not thawed, then the subglacial water contribution of that upstream region is likely much smaller than in the southwestern GrIS. This contribution likely influences the onset and evolution of efficient/channelized subglacial drainage during the summer melt season. Hence, the seasonal evolution of subglacial hydrology elsewhere along the GrIS margin could differ substantially from that directly observed and modeled so far. Identification of presently frozen regions that are experiencing new surface melting or surface-to-bed connections will help bound the impact of this process on the evolution of the basal thermal state [e.g., Poinar et al., 2015; Willis et al., 2015].
MacGregor et al. [2015b] inferred that most of the southern GrIS, including the SW region noted in this study, contains ice at mid-depths that is significantly colder than at the surface, a pattern that could not be fully explained by horizontal advection of colder inland ice. This pattern may suggest the presence of widespread basal melting that is sufficient to draw down (vertically advect) colder ice, but this possibility was discounted in favor of changing surface boundary conditions in the past (temperature and accumulation rate). Here, through our synthesis, we identify a portion of this region as likely thawed, but it is not yet clear that basal melting there is fast enough to invalidate the original hypothesis regarding the cause of relatively cold ice there.
Radar bed reflectivity is commonly used to infer ice-sheet basal conditions, but we did not independently evaluate or include this metric in our synthesis due to the challenge of calibrating bed reflectivity across an entire ice sheet. Oswald and Gogineni  undertook the most extensive investigation of GrIS bed reflectivity so far, but their analysis was focused on the northern half of the GrIS only. Several regions and at least one borehole that they identified as possessing a wet bed are identified as either likely frozen or uncertain by our synthesis and existing borehole measurements (IDS 1.1, 1.2 and 8.2; Camp Century), but there appears to be better agreement regarding a thawed bed within the NEGIS system (IDS 2.1) [Rogozhina et al., 2016]. Further analysis of existing extensive radar data is clearly necessary to reconcile these differences, perhaps by applying existing methods to newer data [e.g., Oswald and Gogineni, 2008, 2012; Christianson et al., 2014] or newer methods to existing data [e.g., Schroeder et al., 2013, 2014, 2016; MacGregor et al., 2015b].
Using a suite of 3-D thermomechanical model runs, Pattyn  predicted that ~55% of the grounded Antarctic Ice Sheet is thawed at its bed, a value somewhat larger than the 43% we report here for the GrIS. Given the size of the uncertain region we identify (34%), these total values cannot be considered to be significantly different, but it is worthwhile to consider the difference between our synthesis of the GrIS’s likely basal thermal state (Figure 11) and that reported by Pattyn  for Antarctica. The GrIS pattern generally transitions along-flow from frozen to uncertain to thawed, whereas for Antarctica this pattern is often reversed, most notably in East Antarctica. This pattern is presumably due to the generally larger ice thicknesses and lower accumulation rates in East Antarctica as compared to Greenland (Figure 12). Pattyn  only considered 3-D models in his assessment, while we synthesized multiple independent methods. However, our synthesis of GrIS 3-D models (Figure 4) more closely resembles our final synthesis than the other methods (Figure 11), supporting the notion that evaluation of either multiple distinct models or multiple instances of a single model are prudent paths when seeking to constrain the basal thermal state of an ice sheet.
We evaluated the present basal thermal state of the GrIS using four independent methods: 3-D thermomechanical modeling and basal motion inferred from radiostratigraphy, surface velocity and surface texture, respectively. These methods vary in their degree of sophistication, but all are motivated by their ability to potentially discriminate between frozen and thawed beds. The resulting synthesis identifies distinct regions where the bed is likely frozen (24% by area) or thawed (43%), and where this basal thermal state remains uncertain (34%).
This first synthesis of the location of the ice sheet’s scalloped frozen core and thawed margins, along with the substantial region of uncertain basal thermal state, are a fundamental constraint on the ice sheet’s large-scale thermodynamics. We find that the three largest regions of the ice sheet that are likely contiguously thawed are the southwest margin, the northwest margin and NEGIS (both the trunk and coastal mouth). A frozen bed is likely west of the central ice divides and in the summit region.
Future work to constrain the basal thermal state of the GrIS should focus on the region that we identify as most uncertain. This region is mostly away from present ice divides, where the ice sheet is in flank flow. It is also where additional borehole measurements to the bed and remote observations would most improve our understanding of the basal thermal state of the GrIS. Such measurements concern not only the basal temperature itself but also the geothermal flux, as inferred from the vertical basal temperature gradients within the ice and underlying rock [e.g., Cuffey et al., 1995; Fisher et al., 2015; Rogozhina et al., 2016]. Improved knowledge of the relatively stable geothermal flux is at least as valuable for predictive ice-sheet modeling as knowledge of the present basal temperature and its vertical gradient, although the latter two quantities are more easily compared with observations.
NSF (ARC 1107753 and 1108058; ANT 0424589) and NASA (NNX12AB71G, NNX13AM16G, NNX13AK27G and NNX13AD53A) supported this work. We thank the organizations (Program for Arctic Regional Climate Assessment, Center for Remote Sensing of Ice Sheets, Operation IceBridge, SeaRISE) and innumerable individuals that both supported and performed the development, collection and processing of the radar data and numerical models used in this study. G. D. Clow was supported by the U.S. Geological Survey Climate and Land Use Change Program. S. F. Price was supported by the U.S. Dept. of Energy Office of Science’s Biological and Environmental Research Program. H. Seroussi was supported by NASA Cryospheric Sciences and Modeling Analysis and Prediction Programs, under a contract with Caltech’s Jet Propulsion Laboratory. We thank A. N. Mabrey for analyzing MOG surface texture, I. Joughin for providing the updated composite surface-velocity field, H. Thomsen for borehole-temperature data and L. C. Andrews for valuable discussions. We thank the editor, associate editor, M. Lüthi and two anonymous reviews for constructive reviews that substantially improved this manuscript. A mask of the likely basal thermal state of the GrIS (Figure 11) will be archived at the National Snow and Ice Data Center (NSIDC).