GBaTSv2: A revised synthesis of the likely basal thermal state of the Greenland Ice Sheet

. The basal thermal state (frozen or thawed) of the Greenland Ice Sheet is under-constrained due to few direct measurements, yet knowledge of this state is becoming increasingly important to interpret modern changes in ice flow. The first synthesis of this state relied on inferences from widespread airborne and satellite observations and numerical models, for which most of the underlying datasets have since been updated. Further, new and independent constraints on the basal thermal 15 state have been developed from analysis of basal and englacial reflections observed by airborne radar sounding. Here we synthesize constraints on the Greenland Ice Sheet’s basal thermal state from boreholes, thermomechanical ice-flow models that participated in the Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6), BedMachine v4 bed topography, Making Earth Science Data Records for Use in Research Environments (MEaSUREs) multi-year surface velocity mosaic v1, and multiple inferences of a thawed bed from airborne radar sounding. Most constraints can only identify where the bed is 20 likely thawed rather than where it is frozen. This revised synthesis of the Greenland likely Basal Thermal State version 2 (GBaTSv2) indicates that 33% of the ice sheet’s bed is likely thawed, 40% is likely frozen, and the remainder (28%) is too uncertain to specify. The spatial pattern


30
The basal interface of an ice sheet is a fundamental control upon its flow and response to external forcings. As such, the icesheet bed is a perennial focus of much glaciological fieldwork and modeling studies, especially its lithology, hydrology and morphology, along with spatiotemporal variability in those properties (e.g., Cuffey and Paterson, 2010). However, the relevance of most basal properties to modulating ice flow is often predicated on the basal temperature being at or very near the pressure-melting point, i.e., a "thawed" basal thermal state. In other words, the bed is only as significant to ice flow as its 35 temperature permits. If the bed is frozen and does not permit significant basal motion or subglacial water flow, then neither its roughness or rheology are likely to significantly influence ice flow at sub-centennial time scales. Resolving an ice sheet's basal thermal state is thus a prerequisite to interpretation of large-scale investigations of most other basal properties and processes. radiostratigraphy modeling, and surface-velocity and surface-texture analyses. The value of this synthesis lay not in its (in)certitude, but in its reduction of the substantial challenge of constraining basal temperature across an entire ice sheet to a simpler ternary determination of the likely basal thermal state of the GrIS. GBaTSv1 has served as a baseline for more sophisticated and localized interpretations of basal properties (e.g., Jordan et al., 2018;Chu et al., 2018;Oswald et al., 2018Bowling et al., 2019, context for other observations of the GrIS (e.g., Bons et al., 2018;Leysinger-Vieli et al., 2018;45 MacGregor et al., 2020; Maier et al., 2021;Karlsson et al., 2021), and as a conceptual framework for investigations of former ice sheets (e.g., Menzies et al., 2020).
Since the generation of GBaTSv1, most of the key datasets that underlie its synthesis have been updated, and some of its Project 6 (ISMIP6). An improved synthesis of GrIS thickness has been generated (BedMachine v4; Morlighem et al., 2017, 55 2021) relative to that used previously (BedMachine v1; Morlighem et al., 2014). A new, complete long-term surface-velocity field for the GrIS is available from the NASA Making Earth Science Data Records for Use in Research Environments (MEaSUREs) program (Joughin et al., 2016(Joughin et al., , 2017. Multiple new studies of airborne radar-sounding data have since been conducted to identify either basal water or deep englacial structures potentially related to a thawed bed (Panton and Karlsson, 2015;Oswald et al., 2018;Jordan et al., 2018;Leysinger-Vieli et al., 2018;Bowling et al., 2019). Finally, recent investigations 60 of basal roughness beneath the GrIS and the transmission of that roughness to the surface warrant a reevaluation of whether surface texture is a reliable indicator of non-negligible basal motion and hence a thawed bed Ignéczi et al., 2018;Cooper et al., 2019aCooper et al., , 2019b. Here we generate a new synthesis of the likely basal thermal state of the GrIS (GBaTSv2) using these new and updated datasets and refined methods. We then consider its differences relative to GBaTSv1 and its implications for interpretation for the present and future flow of the GrIS. 65 2 Data and methods

Direct observations of basal thermal state
As for M16, we consider "direct" observations of the GrIS basal thermal state (and that of Greenland's peripheral ice masses) to include both observations and inferences from deep boreholes, along with unambiguous evidence for subglacial lakes.
Except for NEEM (discussed below), we use the same borehole and subglacial lake observations included in M16 (their Table   70 1). We further include basal thermal state information from 14 additional boreholes (EastGRIP, discussed below, and 13 other marginal boreholes) and two additional active subglacial lakes reported by Bowling et al. (2019)  ibid. 1 As in M16, "corrected" means adjusted for pressure melting using the local ice thickness. "T" means that the basal temperature was not measured directly but that a thawed bed can be confidently inferred.
Since GBaTSv1, two key borehole observations of the GrIS interior have arisen from its two most recent deep boreholes: 80 NEEM and EastGRIP (Fig. 1). The ice thickness at NEEM is ~2538 m, indicating a pressure-melting point of -2.2ºC (assuming a decrease of 8.7 × 10 -4 K m -1 ; Cuffey and Paterson 2010, p. 406). Drilling was completed in 2012 and repeat logging of borehole temperature after the 2011 profile reported by MacGregor et al. (2015a) confirms a basal temperature of ~ -3.5ºC, inferred from the deepest englacial thermistor. However, subsequent logging directly at the base measured a higher temperature of -2.4ºC, presumed to be due to the presence of subglacial water (Colgan et al., 2022). Combined with the recovery of several 85 meters of refrozen, debris-rich ice from the bottom of the NEEM core (D. Dahl-Jensen, pers. comm., 2021), these observations indicate that the base of the NEEM ice core is thawed, rather than frozen as previously estimated by M16. This change in identified basal thermal state at NEEM also implies that the temperature threshold for assuming the bed is thawed in 3-D thermomechanical models should be lower than previously assumed by M16 (Sect. 2.7). 90 Figure 1: Reference map for GrIS study area with driving stress overlain (Sect. 2.5) and locations discussed in the text labeled. Ice-drainage basins are outlined and labeled following Mouginot et al. (2019). The reported basal thermal state of boreholes follows M16 and Table 1, except for NEEM and EastGRIP (Sect. 2.1). The four known subglacial lakes included in M16 are shown, along with two additional active subglacial lakes identified by Bowling et al. (2019).
While the basal thermal state at EastGRIP has not yet been directly measured by borehole thermometry, ice-core drilling 95 there is underway (80% of ice thickness as of September 2021). Preliminary interpretation of the core's depth-age relation indicates that the bed there is thawed, an approach that has correctly predicted the basal thermal state in the past (Dahl-Jensen et al., 2003). Recent phase-sensitive radar measurements also indicate that basal melting is occurring there (Zeising and Humbert, 2021), so we assume that EastGRIP is indeed thawed for this study (Table 1).

3-D thermomechanical modeling of basal temperature
Since GBaTSv1, it remains the case that only 3-D thermomechanical numerical models can estimate basal temperatures beneath the entire ice sheet. To do so requires explicitly solving coupled mass-, momentum-and energy-conservation equations using imperfectly known initial conditions, boundary conditions and constitutive relations. This challenge is met by multiple families of ice-sheet models, of which the most recent and suitable ensemble is ISMIP6 (Goelzer et al., 2020;Nowicki et al., 105 2020). For the GrIS, ISMIP6 constitutes a 21-member ensemble of nine different model families. For the purposes of generating GBaTSv2, such an ensemble is strongly preferred over multiple instances of a single model, as it permits evaluation of a wider range of models with varying ice-flow parameterizations and numerical schemes, whose outputs were homogenized prior to the ensemble analysis. Several of the models used in the SeaRISE ensemble are no longer developed actively (Nowicki et al., 2013), further motivating a transition to the ISMIP6 ensemble. While the choice of ensemble is new, similar trade-offs 110 exist as for the previous ensemble, i.e., variability in initialization and data-assimilation strategies and prescribed boundary conditions (e.g., Table A1 of Goelzer et al., 2020). As for GBaTSv1, we explicitly accept and welcome this diversity of model implementations, and here simply evaluate their agreement with one another. Table 2 lists the ten model instances (hereafter simply "models") from the ISMIP6 ensemble that we consider for GBaTSv2 and our rationale for their selection. Most groups participating in ISMIP6 submitted multiple models with different 115 spatial resolutions or stress-balance approximations; however, basal temperature outputs were not available for all model instances. We selected a single model from each participating group that we assessed to be the most physically complete (e.g., higher-order stress balance) or had the finest spatial resolution. Following M16, for each ISMIP6 model we only consider the modeled basal temperature on grounded ice ("litempbotgr" in ISMIP6 nomenclature) at the end of their 86-yr control runs ("ctrl_proj" in ISMIP6 nomenclature), which are intended to simulate the unforced state of the GrIS at the end of 2100 CE. 120 We assume that this temperature accommodates further thermodynamic relaxation following spin-up, without additional external forcing. Modeled basal temperature ( ! ) is corrected upward for pressure melting ( ′ ! ) using the contemporaneous modeled ice thickness ( ) and assuming a melting-point decrease of 8.7 × 10 -4 K m -1 (Fig. 2).  Figure 2: Modeled basal temperature ( ′ ! ) across the GrIS at the end of ten IMSIP6 control-run experiments, corrected for pressure-melting using each instance's ice-thickness field. The highest temperature range shown (darkest red) represents the range of basal temperatures we consider thawed for our "standard" temperature threshold. Symbology follows Fig. 1, except that the color for borehole symbols instead 135 follows the color scale for their reported or apparent values of ′ ! .
As in M16, the agreement in basal temperature between the ten models selected from the ISMIP6 ensemble is then combined for subsequent inclusion in a multi-method synthesis (Fig. 3). This pattern is qualitatively similar to that of M16 (their Fig. 4) but shows generally greater model agreement overall and more tightly defined thawed regions for some northern, eastern and southeastern outlet glaciers. Part of the key difference between this study and M16 lies in the selection of the 140 temperature thresholds for identifying a thawed bed. Our reinterpretation of the NEEM bed as thawed, despite a basal temperature >1 K below the pressure-melting point, suggests that M16's temperature thresholds were too conservative, i.e., they erred on the side of a frozen bed identification. We thus select -1ºC below the pressure-melting point as the standard temperature threshold for identifying a thawed bed (M16 used -0.05ºC), and increase the range considered for the cold-(-0.5ºC) and warm-bias (-1.5ºC) thresholds. This adjustment acknowledges greater uncertainty in basal thermal state from 145 directly measured borehole temperatures (e.g., Sect. 2.1), which implies greater ambiguity in interpretation of modeled basal temperatures. Figure 3: Agreement in modeled basal thermal state between the ten selected ISMIP6 control-run experiments (Fig. 2), assuming that the bed is thawed where ′ bed ≥ -1ºC and frozen where ′ bed <1ºC.

Basal melting from radiostratigraphy
M16 used one-dimensional (1-D) steady-state modeling of radar-observed Holocene (9-0 ka) depth-age relations to constrain the multi-millennial scale pattern of ice flow across a broad swath of the GrIS interior (69% by area), which can indirectly constrain its basal thermal state. The two primary models used to interpret these depth-age relations, "Dansgaard-Johnsen" (Dansgaard and Johnsen, 1969) and "Nye + melt" (Fahnestock et al., 2001), are both two-parameter representations of vertical 155 strain with differing underlying assumptions about local ice flow. However, for the purposes of constraining basal thermal state, they are fundamentally related. As shown by Fahnestock et al. (2001) and M16, a best-fit Dansgaard-Johnsen model that infers a negative basal shear-layer thickness (ℎ < 0) is qualitatively comparable to a best-fit Nye + melt model that infers a positive basal melt rate (̇ > 0), whereas vice versa implies non-negligible basal freeze-on (̇ < 0). M16 recast the basal shearlayer thickness of the Dansgaard-Johnsen model as a geometric shape factor for the horizontal ice flow of the bulk column.

160
This interpretation offered the potential to constrain not only where the bed is thawed ( > 1) but also where the bed is frozen, because the natural lower limit for should be ( + 1)/( + 2) ≈ 0.8 (Cuffey and Paterson, 2010, p. 310), where is the flow-law exponent and assumed to be 3 (Sect. 2.5). However, a surprisingly large fraction of the interpretable area (57%) displayed values below this limit, calling into question the assumptions underlying that interpretation of .
For GBaTSv2, we retreat from the possible interpretation of a frozen basal thermal state from radiostratigraphy and instead 165 focus only where these data clearly indicate basal melting and hence a thawed bed. This simplifies interpretation of Holocene radiostratigraphy to using the Nye + melt model only and provides a straightforward significance cutoff for interpreting a thawed bed, i.e., where ̇ > 0 cm yr -1 , which we conservatively increase to where ̇ ≥ 1 cm yr -1 (regions with red coloring in Fig. 4). Conversely, apparent ̇ values inferred from radiostratigraphy indicate large regions where ̇ < 0 (Fig. 4). However, as explained by M16, those values should be interpreted primarily as due to a limitation in interpretation of the Nye + melt 170 model in regions where there is non-negligible basal shear, rather than an indicator of widespread, rapid basal freeze-on. This caution is further supported by the independent modeling study of Dow et al. (2018), which indicates that the mean basal freeze-on rate across the GrIS is < 0.02 cm yr -1 .
While the focus of the interpretation has changed, the underlying dataset has not. Here we use the same GrIS dated radiostratigraphy dataset from MacGregor et al. (2015b) considered in M16, because no revision to the GrIS radiostratigraphy dataset yet exists. Uncertainty in ̇ is reflected by the range between its lower-and upper-bound values (̇m in and ̇m ax ), which are determined from the 95% confidence bounds for this model parameter in the Nye + melt model and are the same as for M16.

Basal water from airborne radar sounding
Since M16, multiple studies have mapped the apparent presence of basal water across the GrIS from analysis of airborne radarsounding data, including investigations of bed reflectivity (Jordan et al., 2018;Oswald et al., 2018;Bowling et al., 2019), the morphology of the ice-bed reflection (Bowling et al., 2019), and indirectly via the identification of disrupted basal ice (Panton and Karlsson, 2015;Leysinger-Vieli et al., 2018). Airborne survey coverage is often sparse in the GrIS interior, where large gaps persist that can be tens of kilometers wide; further, at finer scales (< ~50-100 km) there can be notable differences in the 190 inferred location of basal water between individual studies. However, at the scale of the whole ice sheet, the ensemble of the above analyses shows reasonable agreement and can therefore be credibly synthesized to interpret where airborne radar sounding has found evidence of basal water and hence a thawed bed. We merge four of these basal water datasets into a single mask of the likelihood of the presence of basal water ( !' ; Fig. 5). Where nuanced results were reported, indicating various degrees of confidence in the individual datasets by the study authors, we attempt to preserve that nuance when merging them. 195 We attempted to acquire the gridded basal water estimate of Oswald et al. (2018) (their Fig. 13) but were unsuccessful, so it is not included in our synthesis.  Panton and Karlsson (2015) or Leysinger-Vieli et al. (2018). (d) Merged inferences of presence of basal water from analysis of bed reflections or deep radiostratigraphy, respectively, in NASA airborne radar-sounding data ( !% ). A value of 1-4 for !% indicates low confidence, 5-9 indicates medium confidence, and >10 is assigned high confidence. Symbology follows Fig. 1, except that open symbols are used so that underlying inferences of basal water can be better shown. 205 We use the basal water identifications of Jordan et al. (2018) (their Fig. 6), which include along-track binary identifications of basal water from basal radar reflectivity analysis of Operation IceBridge (OIB) and pre-OIB NASA airborne radar-sounding surveys. We binned these identifications into a 5-km grid by the total number of identifications within the nearest grid cell (Fig. 5a). For larger subglacial water bodies, Bowling et al. (2019) synthesized evidence for subglacial lakes beneath the GrIS using multiple well-established criteria to analyze ice-bed reflections in OIB and pre-OIB NASA radar-sounding data. They 210 assigned four possible confidence levels to their identifications: "low", "medium", "high" and "very high" (their Fig. 3). To render these confidence levels compatible with the other datasets, we re-assigned these confidence levels to values of 1, 5, 9 and 10, respectively. We then add those values to the nearest 5-km grid cell (Fig. 5b). In this manner, the contributions to the synthesis of basal water estimates are roughly equalized.
Finally, we include two separate maps of disrupted basal ice by Panton and Karlsson (2015) and Leysinger-Vieli et al.
(2018) further analyze the structure of their identified basal plumes and conclude that they are most likely initiated by basal 220 freeze-on. While the significance of basal freeze-on is controversial (e.g., Dahl-Jensen et al., 2013;Bons et al., 2016;Dow et al., 2018), it remains possible that the genesis of these features could require locally sourced basal water and hence a thawed basal thermal state. In northern Greenland, the patterns of disrupted radiostratigraphy from Panton and Karlsson (2015) and Leysinger-Vieli et al. (2018) are qualitatively similar, so here we simply merge the maps from both studies. We assume that the putative basal water source that initiated these features still exists, and we neglect any horizontal displacement of location 225 of that source relative to their identified location (typically the apex). Following the nomenclature of Leysinger-Vieli et al.
(2018), we bin UDR/plume identifications to the nearest 5-km grid cell and add to them a value of either 1 (for small plumes) or 5 (large). For the UDR identifications of Panton and Karlsson (2015), we ignore regions where ice thickness is less than 1 km, due to a lower signal-to-noise ratio there. We assume all UDRs represent small plumes, except where the ratio of their height above the bed to the ice thickness exceeds ⅓ H, which is the threshold selected by Leysinger-Vieli et al. (2018) for 230 identification of a large plume.
Given sparse survey coverage of the GrIS interior, we assume that evidence of local basal water implies that a broader region of the adjacent bed possesses similar evidence for basal water but is as-of-yet unsurveyed. For each summed bin, we assign all eight adjacent bins the same value, effectively assuming that value for any individual 5-km grid cell is valid within a 15-km-square region centered on that grid cell (Fig. 5d). Similar strategies have been employed previously (e.g., Oswald et al., 2018), although ours is somewhat more conservative in that the regional extrapolation of the basal water signal has a fixed and finite range.

Minimum basal slip ratio
For GBaTSv2, we follow the method introduced in M16 (with minor modifications) to model the ice column's maximum 240 possible deformation speed ( ( ) ) under the end-member assumption that the whole of the ice column is temperate and hence as soft as possible, without directly invoking additional rheological processes (e.g., crystal-orientation fabric or damage).
Where the observed surface speed (|4⃑ * |) is greater than that hypothetical "speed limit", i.e., where the minimum basal slip ratio ( min = |4⃑ * |/ ( ) ) exceeds unity, this implies that non-negligible basal motion is occurring there and that the bed is likely thawed.

245
Similar to M16, we calculate ( ) by assuming that the shallow-ice approximation is appropriate for large-scale estimates of the deformation speed as (Cuffey and Paterson, 2010, p. 310): where the bulk density of the ice column ice is 900 kg m -3 , the rate of acceleration due to gravity is 9.81 m s -2 , the rate factor for temperate ice ( ̅ ) is 2.4 × 10 -24 Pa -3 s -1 and is the depth-averaged enhancement factor for the whole column, is 250 ice thickness and is surface slope in the ice-flow direction.
The three main modifications to M16 are the use of updated datasets, an adjusted filtering scheme, and a revision of the value of and uncertainty therein. First, we now use BedMachine v4 for (Morlighem et al., 2017(Morlighem et al., , 2021  areas of slower ice flow, which includes most of the ice sheet and is also where the basal thermal state is most poorly constrained. Third, M16 assumed a value of unity for and that the relative uncertainty in the product of ̅ was 25%, but Cuffey and Paterson (2010, p. 74) indicate that this value of is too low for polar ice undergoing simple shear. Instead, here we treat = 1 as a lower bound, assume a new default value of 2 and an upper bound of 4. Uncertainty in the value of is considered separately below.   Fig. 6 shows the filtered observed surface speed, modeled temperate-column deformation speed and the ratio of these two fields ( 670 ). Where 670 > 1, the ice column is inferred to exceed its speed limit due to deformation alone and that basal motion is occurring, implying a locally thawed bed. Assuming that regions where surface speed is poorly known limit the usefulness of this method, regions where the uncertainty in the surface speed ( C * ) is ≥ 25% of |4⃑ * | are not included in this method's assessment of the basal thermal state; these regions are located primarily along central and southern ice divides (Fig.   275 6a). Assumed and reported uncertainties in ( ) and |4⃑ * |, respectively, are used to calculate lower-and upper-bound values of 670 to then assess uncertainty in the basal thermal state agreement from this method. This analysis produces a substantially smaller region than M16 where min > 1, principally because ( ) is roughly twice as large as was previously assumed, due to the change in assumed value of .
Since M16  and it cannot be simply disentangled from its associated value and then corrected using an Arrhenius relation to a presumed temperate value, which is necessary for our method of calculating min . Further, the range of values that we now consider (1-4) produces substantially larger increases in ( ) with = 3 than using the Bons et al. (2018) rheological parameters, so we consider this range more suitable for calculating min .

Discontinued methods
As part of GBaTSv1, M16 mapped the onset of surface undulations across the GrIS from surface imagery, as they are suggestive -but not definitively indicative -of the onset of substantial basal motion and hence a thawed bed. Since M16, multiple additional studies have further explored both the nature of basal roughness beneath the GrIS (Cooper et al., 2019a), how that roughness is transmitted to the surface via basal motion Ignéczi et al., 2018), and made independent 295 observations of surface texture (e.g., Cooper et al., 2019b). When considered together with de Rydt et al. (2013), which formed part of the rationale for including this method in M16, we conclude that the onset of surface undulations can no longer be considered a reliable indicator of a thawed bed and we discontinue its use for GBaTSv2. Our rationale is explained further below.
Surface undulations due to ice flow over bedrock obstacles are expected to be more prominent where either basal 300 roughness is more pronounced or the ratio of basal motion to the deformation speed is greater ( , Sect. 2.5). Cooper et al. (2019a) found that basal roughness beneath the GrIS observed by airborne radar sounding at along-track scales of 200 m is typically greater within ~200 km of the ice margin than farther inland. Along-flow roughness is more likely to be efficiently transmitted to the surface than across-flow roughness (e.g., Ng et al., 2018), and the pattern of greater marginal roughness is less pronounced along-flow. Basal roughness at a 200-m horizontal scale is unlikely to generate significant surface undulations 305 where the ice sheet is generally several times thicker than that . However, in northwestern Greenland, rougher marginal areas also have a higher degree of self-affinity, suggesting they are also rougher at larger horizontal scales (Jordan et al., 2017). Overall, these studies imply that a priori we should expect more surface undulations closer to the ice margin due to increasing basal roughness there, independent of any change in basal thermal state. Separately, Ng et al. (2018) and  refined modeling of bed-to-surface transmission and further emphasize the primary role of topography in generating 310 modeled surface undulations that credibly reproduce observations, rather than those generated purely by a non-zero slip ratio.
The value of outlining the onset of surface undulations for GBaTSv1 was predicated on the dominant role of the latter mechanism only. Finally, Cooper et al. (2019b) found evidence of the surface expression of englacial features (e.g., disrupted basal units) and subglacial channels oriented along-flow in northwestern Greenland. These expressions are clearly not surface undulations that might be diagnostic of a thawed basal thermal state yet can be easily confused for them. Similarly, Kjaer et al.

315
(2018) and MacGregor et al. (2019) demonstrated that the presence of subglacial impact craters can be discerned partly from their surface expressions, and these structures are not yet conclusively associated with a particular basal thermal state. In summary, we conclude that it is no longer clear whether outlining surface undulations can be considered a reliable method for demarcating a low or negligible basal slip ratio ( → 0), for which no basal thermal state assignment can be made, from a high or non-negligible ratio ( ≫ 0) that clearly indicates a thawed bed. 320

Synthesizing basal thermal state estimates
We follow a similar methodology to M16 for generating GBaTSv2, with several minor adjustments. The thresholds for a positive identification of a particular basal thermal state are summarized in Table 3, including both the "standard" values and cold-and warm-bias values that consider uncertainty in each method; these are later used to assess the likelihood of a particular 325 basal thermal state. Same threshold, but estimated using |%⃑ | − # and ( ) (1 + ̅ 5 ) a Note these changes from GBaTSv1, which used -0.05ºC, 0ºC and 0.5ºC as the standard, cold-and warm-bias thresholds, respectively. b The standard (5), cold-(10) and warm-bias (1) thresholds are equivalent to the number of basal water identifications within each 5-km grid 330 cell synthesized by !% .
We synthesize the four methods of constraining the likely basal thermal state of the GrIS by first assessing where they each produce a clear signal regarding this state (Table 3; Fig. 7a). We initialize a 5-km gridded ice-sheet mask to zero. For each method, if that method indicates a thawed bed, then +1 is added to (Fig. 7b). Conversely, if the method indicates a 335 frozen bed, then -1 is added to . For the ISMIP6 ensemble, the agreement is considered significant only where at least 7/10 models agree that the bed is either frozen or thawed (Fig. 3), a more conservative assessment from GBaTSv1, for which only a plurality (more than half) of the eight SeaRISE models had to agree to reach the same assessment. All 3-D models are weighted equally, as are each of the methods. This process of generating is repeated using the cold-and warm-bias thresholds to generate cold and warm , respectively (Fig. 7c,d). For GBaTSv2, only one method can distinguish a frozen bed (3-D Figure 7: (a) Outlines of a thawed GrIS bed for the four methods considered in this study (Fig. 3, 4, 5d, 6c). For the ISMIP6 agreement, their outline denotes where at least 7/10 models agree that the bed is thawed. (b) Agreement between the four methods regarding the basal 345 thermal state using standard thresholds ( ; Table 3). (c, d) Cold-and warm-bias agreement ( cold and 9:;< , respectively) determined using each method's confidence bounds or ad hoc uncertainty estimates. Because only one method constrains where the bed is frozen (thermomechanical models), but all four constrain where it is thawed, the range of possible values is -1 (frozen) to +4 (all thawed).
Based purely on , cold and warm , we generate the likely basal thermal state mask ( ), which synthesizes their agreement 350 and is the primary GBaTSv2 product. is initialized to zero (uncertain basal thermal state) and then assigned +1 for a likely thawed bed where both and warm agree that the bed is thawed and cold does not indicate that the bed is frozen. Similarly, L is assigned -1 for a likely frozen bed where both and cold agree that the bed is frozen and warm does not contraindicate them. In other words, given our present understanding of the uncertainty of each method, we do not assign a likely basal thermal state if any of the three instances of contradicts the other two. We only consider the sign of and ignore its 355 magnitude. Where contains small "holes" (≤ 10 grid cells, equivalent to ≤ 250 km 2 ) in predominantly thawed regions, i.e., small regions with a different basal thermal state (uncertain or frozen), these are filled in as in M16 and assumed to be likely thawed. This process is then repeated for small holes in frozen regions, except those are assumed to be likely frozen. However, we do not repeat this process for holes in uncertain regions, as was done for GBaTSv1, because we infer that there is insufficient evidence to justify a particular assignment of basal thermal state there. These hole-filling procedures result in less than a 1% 360 difference in the total area assigned to each basal thermal state. Fig. 8 shows version 2 of the likely basal thermal state of the GrIS (GBaTSv2), based on the four methods and their synthesis described in Sect. 2. At the scale of the whole ice sheet, this synthesis is qualitatively similar to GBaTSv1, but there are notable regional differences highlighted below and summarized by ice-drainage basin in Table 4. The most prominent differences are 365 along the southern portion (≤ 68ºN) of the central ice divide (more contiguous regions of likely frozen bed in GBaTSv2), west of the central ice divide that lies between Summit and NorthGRIP (less confidence in a frozen bed in GBaTSv2), and within the drainage basin that includes the NEGIS (NE; more contiguous regions of likely thawed bed northwest of NEGIS in GBaTSv2). Similarities between the two versions include large contiguous regions of likely thawed bed along the southwestern and northwestern coasts (up to Melville Bay), and within the NEGIS ice-drainage basin. The "scalloped frozen core" described 370 by M16 is now potentially dissected between its southern and northern reaches, primarily due to the reduced agreement between 3-D thermomechanical models on the extent of the frozen-bedded region.  As compared to GBaTSv1 (M16), GBaTSv2 reports an uncertain basal thermal state at both NorthGRIP and DYE-3.

Results
However, it evinces greater uncertainty in the vicinity of NorthGRIP (especially farther west), but greater confidence that regions near DYE-3 are frozen. We interpret both changes as minor improvements in GBaTSv2 over GBaTSv1. However, we note two areas of concern in terms of GBaTSv2 misidentification of the basal thermal state, as compared to direct observations 385 ( Table 1). The first area is NEEM, which is not surprising given the change in its assessed basal thermal state. The second area is the Prudhoe Lobe of the GrIS in far northwestern Greenland where Palmer et al. (2013) identified two subglacial lakes from radar sounding. While GBaTSv2 at the lake locations is uncertain, most of the rest of this lobe is likely frozen.
The assigned likely basal thermal state of 46% of the GrIS changed between GBaTSv1 and GBaTSv2; for < 6% of the GrIS, the assigned state changed from likely frozen to likely thawed or vice versa. GBaTSv2 identifies more of the bed to be 390 likely frozen (+16%) and less to be likely thawed (-11%) than GBaTSv1. At first glance, this is surprising, because only 3-D thermomechanical models are used to identify a frozen bed in GBaTSv2. However, the loss of the discontinued method (onset of surface undulations) decreases the likelihood of a thawed bed identification, and the new method employed (basal water from airborne radar sounding) is inherently sparser in its more robust identification of a likely thawed bed (Fig. 4).

395
A comparison of Fig. 3 in this study to Fig. 4 from M16 suggests that changes in bed topography influence thermomechanical model agreement on basal thermal state, as the pattern of agreement appears more focused in some regions, particularly along major outlet glaciers. Fig. 9a,b shows this difference in agreement in basal thermal state between the SeaRISE and ISMIP6 thermomechanical models. While we observe a possible relation between change in ice thickness and basal thermal state in the vicinity of several outlet glaciers, particularly in southern Greenland, the pattern is more nuanced 400 across most of the ice-sheet interior. There is a noticeable divergence between northern and southern Greenland at around ~73ºN (Fig. 9b). This difference is not attributable to a new geothermal flux field, because most models from both ensembles use the older geothermal flux field derived from seismic data of Shapiro and Ritzwoller (2004), rather than a more recent field derived from aeromagnetic data (Martos et al., 2018) or machine learning (Colgan et al., 2022). Most ISMIP6 models used the BedMachine v3 bed topography, which on average results in thicker ice than the various bed topographies used by SeaRISE.

405
A local increase in reference ice thickness between SeaRISE and ISMIP6 (Fig. 9a) would presumably tend to increase agreement where the bed is thawed, as the pressure-melting point at the bed will decrease. However, these changes are poorly correlated (linear correlation coefficient = 0.08; Fig. 9c). Figure 9: (a) Change in ice thickness between the most commonly used syntheses for SeaRISE (Bamber et al. (2001) with modifications) and ISMIP6 (Morlighem et al., 2017) on the 5-km grid used in this study. (b) Change in agreement in modeled basal thermal state from SeaRISE (Fig. 4 of M16) to ISMIP6 (Fig. 3), where a positive (negative) difference indicates greater agreement that the bed is thawed (frozen). Values greater than ±100% indicate that model agreement in the basal thermal state changed significantly. (c) Histogram of change in ice thickness vs. change in model agreement of a likely thawed bed.

415
A more likely explanation for the change in modeled basal thermal state is a mean change in the spin-up surface mass balance (SMB) field between the two ensembles, i.e., higher SMB in southern Greenland and lower SMB in northern Greenland for the ISMIP6 ensemble, along with possible changes in modeled surface paleotemperatures. Higher snowfall rates in the dry snow zone over multiple millennia lead to increased downward vertical advection and an overall colder ice column, 420 which increases the likelihood of a frozen bed, and vice versa for lower snowfall rates. Unfortunately, the spin-up SMB fields used in SeaRISE and ISMIP6 models are more varied, so a simple comparison as in Fig. 9 for ice thickness cannot be generated easily to verify this hypothesis. Therefore, the root cause of the change in modeled basal thermal state remains not yet well understood.
While GBaTSv2 continues to be reported on a 5-km grid, it is increasingly clear that the basal thermal state can vary at 425 scales finer than that (e.g., Chu et al., 2018). Further, englacial thermal structure can be quite variable at finer scales than 5 km (e.g., Lüthi et al., 2002, Harrington et al., 2015Maier et al., 2019). Colgan et al. (2021) recently highlighted the role of bed topography in influencing geothermal flux at kilometer scales, a likely primary control on basal thermal state. However, this influence may be less important where there is negligible ice advection and basal temperature gradients are dominated by heat diffusion (Willcocks et al., 2021). The sum of these studies suggests that finer-resolution geophysical methods and models are 430 required to further specify the nature of Greenland's basal thermal state. This need could potentially be partly addressed by more intensive borehole investigations of regions where the basal thermal state is in question, especially in the deep interior of the GrIS and perhaps along existing flight lines where interpretations of airborne radar-sounding data disagree. Following the conclusion of OIB (MacGregor et al., 2021a), an opportunity now exists for an updated and more complete synthesis of basal water identifications from the data that mission collected, following existing methods (e.g., Jordan et al., 2018;Chu et 435 al., 2018).
We next consider the impact of the GBaTSv2 revision upon several existing studies that consider either the ice sheet's basal thermal state or used GBaTSv1 explicitly. Poinar et al. (2015) concluded that surface-melt-induced acceleration of ice flow is unlikely to propagate inland significantly within the SW basin of the GrIS, partly because they modeled that the bed is mostly thawed farther inland there. Although GBaTSv2 indicates decreased confidence in a thawed bed beneath the upper 440 reaches of SW basin of the GrIS, as compared to GBaTSv1, the extent of this change does not appear to impact the conclusions of Poinar et al. (2015). Maier et al. (2021) found that inferences of Greenland's regional basal rheology were insensitive to which portion of each ice-drainage basin was identified as likely thawed or uncertain in GBaTSv1. For most basins, the combined region identified as likely thawed or uncertain in GBaTSv2 is similar in extent to GBaTSv1's likely thawed region, so we expect that Maier et al. (2021)'s conclusions regarding Greenland's regional basal rheology are not substantially affected 445 by this revision. A similar assessment applies to Karlsson et al. (2021)'s estimate of GrIS basal meltwater production. They used GBaTSv1 to constrain the geothermal and frictional contributions to basal melt. Because frictional heating is concentrated within typically thawed outlet-glacier systems whose likely basal thermal state is mostly unchanged from GBaTSv1 to GBaTSv2, the frictional contribution to basal melt is unlikely to change significantly except for some eastern outlet glaciers south of ~75ºN, where it may decrease. For the geothermal contribution to basal melt, GBaTSv2 suggests that southern 450 drainage basins may produce less basal melt, while northern drainage basins may produce more. The net effect of the GBaTSv2 revision is likely that the total estimated GrIS basal melt is less than that reported by Karlsson et al. (2021), but it is unlikely that this change exceeds the ~20% relative uncertainty in their estimate.
It is unlikely that the basal thermal state of the GrIS will significantly affect its evolution over this century, which is the period considered by ISMIP6 (Goelzer et al., 2020). However, if anthropogenic climate forcing persists beyond this century 455 and continues for a substantial portion of this millennium, then GrIS retreat will likely be substantial and its present basal thermal state will have a progressively greater influence upon the nature of this retreat, because submarine melting will very likely outpace thermal diffusion at the bed (Aschwanden et al., 2019). Because observed GrIS mass loss presently tracks at the upper end of the range projected by the ISMIP6 ensemble and the socioeconomic impacts of rising sea levels are vast (Aschwanden et al., 2022), it remains essential to produce ice-sheet-wide assessments of basal boundary conditions, such as

Data availability
The core result of this study, version 2 of the likely basal thermal state of the Greenland Ice Sheet (GBaTSv2; Fig. 8b), is available at https://doi.org/10.5281/zenodo.5714527 (MacGregor et al., 2021b) and will also be later made available through 465 the National Snow and Ice Data Center, where further dataset-specific documentation will be provided. This dataset also includes the syntheses of other freely available datasets shown in Figs. 3, 4, 5d, 6c and 7b/c/d.

Conclusions
We have developed and presented the second version of the likely basal thermal state of the Greenland Ice Sheet (GBaTSv2). This second estimate is broadly similar to the first, although there is substantial regional variability therein and a greater 470 tendency toward a likely frozen basal thermal state. The large-scale similarity is likely due to the applied methods being mostly replicated from the first version, despite underlying updates to associated datasets and the discontinuation of one method. This new synthesis suggests that the bed of the GrIS is roughly equal parts thawed (33%), frozen (40%) or too uncertain to specify (28%). Although the use of an improved bed topography beneath the GrIS within 3-D thermomechanical models does not appear to be related to greater agreement in basal temperature within those models, we do observe more spatially focused 475 patterns of likely thawed bed within outlet-glacier systems that have been better mapped since GBaTSv1. The effect of these revisions upon existing studies that used GBaTSv1 is likely to be modest, but the influence of the basal thermal state upon ice flow is likely to increase if anthropogenic climate forcings persists beyond this century. Absent future investigations to directly measure basal temperature in new boreholes, to more extensively identify basal water from remote sensing and to map likely pathways for that basal water, the suite of methods we employed may be approaching a natural limit in its ability to resolve 480 the basal thermal state. Future syntheses should consider new, finely resolved yet ice-sheet-wide observations, which will most likely come from further campaigns or advances in airborne or satellite remote sensing.