Geologic Provinces Beneath the Greenland Ice Sheet Constrained by Geophysical Data Synthesis

Present understanding of Greenland's subglacial geology is derived mostly from interpolation of geologic mapping of its ice‐free margins and unconstrained by geophysical data. Here we refine the extent of its geologic provinces by synthesizing geophysical constraints on subglacial geology from seismic, gravity, magnetic and topographic data. North of 72°N, no province clearly extends across the whole island, leaving three distinct subglacial regions yet to be reconciled with margin geology. Geophysically coherent anomalies and apparent province boundaries are adjacent to the onset of faster ice flow at both Petermann Glacier and the Northeast Greenland Ice Stream. Separately, based on their subaerial expression, dozens of unusually long, straight and sub‐parallel subglacial valleys cross Greenland's interior and are not yet resolved by current syntheses of its subglacial topography.


Introduction
The Greenland Ice Sheet (GrIS; not including peripheral ice masses) covers 1.1% of Earth's land surface and nearly 79% of Greenland itself.The ice sheet's large-scale flow pattern, its discharge into the surrounding oceans and the island's response to changing ice loading are all influenced by solid-Earth boundary conditions that are obscured by the GrIS and thus difficult to constrain.The solid-Earth boundary conditions that have received the most attention are its subglacial topography, geothermal heat flux, crustal geology and mantle viscosity (e.g., Adhikari et al., 2021;Alley et al., 2019;Braun et al., 2007;Colgan et al., 2022;Martos et al., 2018;Morlighem et al., 2017;Rezvanbehbahani et al., 2017;Rogozhina et al., 2016).These properties are potentially related, for example, a geologically young province with a rough topographic surface may have a higher geothermal heat flux and overlie a less viscous mantle, but these relations remain difficult to constrain precisely despite decades of investigation.
Of these properties, Greenland's subglacial geology is perhaps the most poorly constrained, despite its potential influence on ice flow (e.g., Winsborrow et al., 2010).Dawes (2009) (hereafter D09) produced the seminal map of Greenland's subglacial geology that is still widely used for both interpretation of island-wide studies and global maps of geologic provinces (e.g., Alley et al., 2019;Darbyshire et al., 2018;Gowan et al., 2019;Hasterok et al., 2022;Jones et al., 2021;Maier et al., 2022;Martos et al., 2018;Mordret, 2018;Pourpoint et al., 2018;White et al., 2016).However, the boundaries between major provinces that D09 identified were typically unconstrained across several hundred kilometers between similar synchronous formations known only from Greenland's ice-free periphery.Contemporaneous knowledge of magnetic anomalies was considered only for southern Greenland and seismologic inferences received limited consideration.Subglacial topography was not considered, despite being directly identified as a potential constraint on subglacial geology, nor were gravity data.D09 acknowledged these limitations by describing their resulting map as "conjectural" and mapping of subglacial geology as "the next major challenge for Greenland mapping."Since D09, the extent of marginal geologic mapping, seismic imaging, aerogeophysical surveys, satellite mapping of gravity and magnetic anomalies, and surface properties across Greenland have improved substantially.An opportunity exists to improve our understanding of Greenland's subglacial geology by synthesizing this new wealth of pertinent data and associated analyses.Here we map major boundaries in these data sets and reconcile these geophysical constraints with marginal geologic exposures and samples from beneath and within deep ice cores.We consider the implications of this new map of Greenland's geologic provinces for interpretation of Greenland's geologic history.

Province Delineation
We evaluated 19 geologic and geophysical data sets for this study to delineate major subglacial geologic provinces beneath the GrIS (Table S1 in Supporting Information S1).Most are evaluated as is with only minor standard corrections applied.A primary challenge in reconciling these data sets lies in their differing depth sensitivities, resolutions, interpolations and interpretations.The spatial resolution, precision and accuracy of each data set are discussed in detail in the cited references; deeper-sensing or satellite-derived data sets typically have coarser resolution.For simplicity, we only consider airborne data collected by NASA's Operation IceBridge and predecessor NASA surveys (MacGregor et al., 2021).Except for nearby subaerially exposed rocks, the Greenland Ice Sheet Project 2 (GISP2) bedrock sample and to lesser extent ice-core basal ice lithologies, all data sets are indirect indicators of subglacial geology, but none are uniquely diagnostic of a particular subglacial geologic province.However, combinations thereof have been used regularly to diagnose subglacial geology in Greenland and Antarctica (e.g., Brozena, 1995;Fahnestock et al., 2001;Jordan et al., 2023;MacGregor et al., 2019;Paxman, Austermann, & Tinto, 2021;Paxman, Tinto, & Austermann, 2021;Paxman et al., 2019;Rogozhina et al., 2016;Tinto et al., 2019).
We generally selected the most recent, comprehensive and accessible data set for each relevant method.The numerous existing inferences of Greenland crustal properties from seismic tomography are particularly challenging to reconcile, so here we consider only one recent study that applied state-of-the-art methods to a recent seismic data catalog to explore Greenland's crustal velocity structure (Darbyshire et al., 2018).While some geophysical inversions can coarsely resolve three-dimensional crustal structure, here we focus on inferences relevant to the delineation of apparent geologic provinces at their outcrop, that is, as they might be exposed subaerially if the GrIS and any surficial sediment cover were absent.
We first combine the D09 map of Greenland's subglacial geology and a more recent map of the island's subaerially exposed geology to delineate subaerial province extents (Pedersen et al., 2013;Watt, 2019).At a regional scale (∼1:500,000), we then manually trace apparent boundaries at >10-km spacing for each geophysical field separately.Where possible, this delineation is guided by each field's relation to mapped subaerial province boundaries in ice-free regions.The subglacial boundaries are based on where the magnitude, azimuthal trend or regionalscale pattern of each field show significant discontinuities, indicating a possible change in subglacial geology (e.g., Paxman et al., 2019).Based on the apparent significance of each discontinuity, geophysical boundaries are assigned either low or high confidence.For subglacial topography, we trace one boundary set that considers the ice-free rebounded subglacial topography, and a second set based on the ice sheet's surface expression of presentday subglacial topography.The latter set is informed by the pattern of mapped surface lineations (Section 2.2).
We then manually delineate new subglacial geologic boundaries based on our geophysical boundaries, their spatial coherence across both their ice-free exposures and apparent subglacial extent, our confidence therein, and consensus among this study's authors.Where multiple margin-originating geophysical boundaries are aligned within ∼50 km of each other, we identify the geologic boundary as high-confidence and assign it to either the highest-confidence geophysical boundary or the midpoint between several of them.We identify most apparent geologic boundaries that do not reach the ice-sheet margin as low-confidence.Where available data cannot be reconciled with known provinces, we identify the region as unresolved.Similarly, we do not assign subglacial geology beneath Greenland's smaller peripheral ice masses, because these are often poorly resolved by the geophysical data sets considered here and because several of them reside near known transitions in subaerial geology (e.g., Flade Isblink in far northeastern Greenland).

Flow-Aware Hillshade
We introduce one new method: the calculation and interpretation of a "flow-aware" hillshade of a digital elevation model (DEM) of the ice surface.This method is motivated by the observation and supporting theory that an ice mass flowing over variable bed topography transfers a signature of that subglacial topography to its subaerial surface (Gudmundsson, 2003;Ng et al., 2018), a phenomenon that is not presently considered to infer bed topography in the ice-sheet interior (Morlighem et al., 2022).To calculate this hillshade, we first filter the measured flow direction across the GrIS interior (Joughin et al., 2016(Joughin et al., , 2017) ) toward the direction of the surfaceelevation gradient instead, as measured by the Greenland Ice Mapping Project (GrIMP; Howat et al., 2022) and following MacGregor et al. (2022).We illuminate the GrIMP DEM using a standard hillshade algorithm (Greene et al., 2017), modified so that the illumination azimuth may vary for each pixel, which we select as 90°counterclockwise to the filtered flow direction.The advantage of the flow-aware hillshade is that it highlights surfaceslope features while cognizant of the frame of reference most relevant to an ice mass (its flow direction).Its appearance is similar to the "residual surface elevation" method introduced by Cooper et al. ( 2019) but the flowaware hillshade (hereafter simply hillshade) is easier to implement ice-sheet-wide.
The resulting hillshade highlights sub-kilometer-scale variability in GrIS subaerial surface slope (Figure 1).At this scale, variations in surface slope are related primarily to variability in subglacial topography, although this bed-to-surface transfer is complex and depends on the rate of basal sliding and rheology of the ice column (Ng et al., 2018).Contiguous surface slope features are often interpreted as subglacial valley networks in both Greenland and Antarctica (e.g., Ekholm et al., 1998;Livingstone et al., 2017;Ross et al., 2014;Scambos & Haran, 2002).Changes in surface slope may also be related to short length-scale variability in subglacial sliding (Cooper et al., 2019;Gudmundsson, 2003), but these are also detectable by surface-velocity observations and appear to be uncommon-shear margins notwithstanding.As an indirect but independent constraint on icethickness-scale subglacial topography, the hillshade is especially valuable for interpreting gaps between radarsounding surveys that are widespread across the GrIS interior (Morlighem et al., 2022).
The GrIMP hillshade highlights innumerable subtle but coherent surface lineations that are associated with ice flow across subglacial valleys (Figures 1b and 1d) and mountains (Figure 1e), although their interpretation is complicated by the geometry of the radar-lineation intersections and their relation to the ice-flow direction.Some of these lineations are also visible in Moderate Resolution Imaging Spectroradiometer (MODIS) mosaics (Haran et al., 2013), but many more are evident in this GrIMP hillshade.We manually traced prominent surface lineations across the GrIS, focusing on those that were longer than 10 km, wider than 1 km, detectable in the interior, or intersected radar-sounding surveys.This tracing is non-exhaustive but representative of the dominant lineation patterns.

Results
Figure 2 shows the traced boundaries for each geophysical field.For most province boundaries south of ∼70°N, there is general agreement on their location (boundaries within ∼100 km of each other), except for boundaries in the depth to the Mohorovičić discontinuity, which are generally meridionally oriented and rarely coincide with other boundaries.Farther north, some boundaries roughly coincide across multiple fields, such as in far northwestern and eastern central Greenland, but more common is a complex set of intersecting boundaries.This complexity favors our synthesis approach, rather than a focus on one geophysical field as in single-method studies.

10.1029/2023GL107357
Figure 3 shows our synthesis interpretation of apparent subglacial geologic provinces, Table S2 in Supporting Information S1 summarizes their geophysical expressions and Figure S2 in Supporting Information S1 compares each field's boundaries with our synthesis.We delineate 12 partly subglacial provinces; those whose extent differ substantially from D09 are discussed below (Figure S3 in Supporting Information S1).Five additional fully subglacial regions are discussed, four of which are delineated in this study.
For the Inglefield Orogen, we do not find evidence for a dissected eastward extension from Inglefield Land (∼78.5°N,68°W) across to northeastern Greenland, but rather multiple high-confidence geophysical boundaries at ∼50°W.Following D09 and Nutman et al. (2019), we distinguish the Victoria Fjord Terrane exposed at the head of Victoria Fjord (∼81.6°N,45.5°W) from the younger Rae Craton.For this terrane, we find some geophysical evidence for narrower eastern and western subglacial extensions and a more southerly extension from Victoria Fjord.For the main exposure of the North Atlantic Igneous Province along the Blosseville Coast (∼69°N, 25°W), multiple geophysical boundaries indicate a limited extension inland, as is the case for the Nuussuaq Basin in western Greenland.The geophysical expression of the predominant Rae Craton is varied, but our mapped boundaries suggest that its inland extent from the northeastern coast is generally ∼100 km west of the Petermann Megacanyon (Figure 3b; Bamber et al., 2013), although it may extend farther inland into northern central Greenland.The inland boundary of the Rinkian Orogen is approximately constrained along the northwestern coast, while the better constrained extents of the North Atlantic Craton, the Nagssugtoqidian, Ketilidian and Caledonian orogens are similar to D09.Conversely, we interpret the Independence Fjord Basin to extend farther inland to the west, and it appears to underlie the mid-section of the Northeast Greenland Ice Stream (NEGIS; Figures 3b and 3d).
We outline three larger regions whose subglacial geophysical expressions cannot yet be reconciled with the inland extension of margin-exposed provinces.In the central interior of northern Greenland, we outline a large region with relatively indistinct and lower-magnitude anomalies for most fields, which is bounded by known provinces with more distinct and widespread geophysical expressions (e.g., gravity and magnetic highs to its northwest and northeast).GISP2 is located near the apparent southeastern boundary of this region (Figure 3a), and its granitic bedrock appears consistent with the Rinkian Orogen (Weis et al., 1997).However, given limited geophysical similarity with the outcrops of that province or the regionally predominant Rae Craton, we opt not to ascribe this region of unknown geologic age and uncertain boundaries to any known province, and instead denote it generically as "Central Region A." Directly southeast of Central Region A, we outline a separate interior region with distinct gravity and magnetic anomalies (in some cases co-located), along with seismically slow crust, and rough subglacial and surface topography.This combination of properties suggests a volcanic origin and unknown age but likely younger than Central Region A (e.g., Fahnestock et al., 2001).This province is located just south of the onset region of the NEGIS, and ice flows across its northern boundary toward the southeastern flank of the NEGIS (Figure 3d).Both Rogozhina et al. (2016) and Jones et al. (2021) have previously inferred higher geothermal heat flux in the vicinity of this region as a potential contributor to the anomalous flow of NEGIS.This approximate region was also previously described as the "Central Greenland Magmatic Province" by Brozena (1995), based on its aerogeophysical properties, but we designate it agnostically and generically as "Central Region B."  (Forsberg & Jensen, 2015).(d) Earth Magnetic Anomaly Grid (Meyer et al., 2017).For panels C and D, Operation IceBridge airborne gravity and magnetic data are overlain (Cochran & Bell, 2016;Cochran et al., 2014;K. Tinto, 2020).(e) Isostatically rebounded BedMachine v5 topography assuming ice-free conditions (Morlighem et al., 2017(Morlighem et al., , 2022;;Paxman et al., 2022aPaxman et al., , 2022b)).(f) GrIMP hillshade with surface-exposed provinces and major faults (Harrison et al., 2011;Watt, 2019).Colors for provinces follow Figure 3a.(g) All traced geophysical boundaries.Background map is MODIS Mosaic of Greenland (Haran et al., 2013).
A third region borders both Central Regions A and B in eastern Greenland.Both its observed and apparent subglacial topography are notably rougher than those regions, and it shares few other geophysical characteristics with those regions or the Caledonian Orogen to its east (Figure 2).The simplest geologic explanation for this region may be that it is part of the Rae Craton (e.g., Thrane, 2002), but it is sufficiently distinct geophysically that we designate this third ambiguous region as "Central Region C." Tinto et al. (2015) identified a prominent ellipsoidal gravity, magnetic and smooth topographic high east of Petermann Glacier, northwestern Greenland (Figure 3c), which they interpreted as an "igneous body."It is clearly bounded on its western flank by the onset region of faster ice flow of Petermann Glacier, and along its northern and eastern flanks by the Petermann Megacanyon (Figure 3c).This ∼4,000 km 2 region is not large enough to constitute a separate province, so we nominally interpret it as part of the Inglefield Orogen and denote it generically as the "Petermann Anomaly."Based on basalt and porphyry erratics found in northwestern Greenland, Dawes et al. (2000) hypothesized the presence of an unexposed volcanic subglacial province beneath the northern GrIS.This putative province was outlined by D09 and-coincidentally-the Petermann  S2 in Supporting Information S1; Cohen et al., 2013Cohen et al., updated 2023 at https://stratigraphy.org; Watt, 2019).Bedrock and basal sediment lithology and ages from three interior boreholes and ice cores also shown (Blard et al., 2023;Christ et al., 2023;Weis et al., 1997).(b) Major glaciological and subglacial features discussed with surface speed overlain (Joughin et al., 2016(Joughin et al., , 2017)) Anomaly is at the center of that outline and nearly entirely within it (Figure S3 in Supporting Information S1).If this feature is indeed an igneous body and if ice flow over it was diverted farther west during a more extensive Pleistocene glaciation, then it is conceivably the source of those volcanic erratics and likely somewhat younger than the Inglefield Orogen.Separately, we do not find clear geophysical evidence to bound another possible unexposed province identified by D09 in northwestern Greenland that could explain the presence of redbed erratics found along the coast of Baffin Bay (Figure S3 in Supporting Information S1).For completeness, we still include this putative yet unlocated subglacial province in Table S2 in Supporting Information S1.
We traced 1062 surface lineations within the GrIMP hillshade (Figures 1a and 3b).Where they intersect radarsounding surveys, these lineations are mostly related to sparsely surveyed subglacial valleys.Many of these lineations connect farther downstream with the better-surveyed deeply incised troughs beneath most outlet glaciers around the GrIS periphery (Morlighem et al., 2017); along the northwestern coast facing Baffin Bay, several of these valleys appear to be >300 km long.Networks of these lineations are typically subparallel to each other across hundreds of kilometers, an unusual configuration that implicates past crustal tectonics or subglacial incision as a control on their morphology (Rose et al., 2014).Near the margin, the azimuths of these subglacial valley networks are often similar to those of previously mapped major faults in adjacent ice-free forelands (Figure 2f), further suggesting a relation to subglacial geology.

Discussion
A substantive change between our mapping of subglacial geologic provinces and D09 is that we do not find clear geophysical support for a dissected subglacial connection between the North Atlantic Igneous Province and the Nuussuaq Basin (Figure S3 in Supporting Information S1).If Central Region B is indeed magmatic, as hypothesized by Brozena (1995), it extends too far north to connect those peripheral provinces in a simple way.If there is a connection between those provinces and this unresolved region-perhaps associated with the Iceland hotspot for which numerous tracks have been hypothesized (cf., Martos et al., 2018)-then this region may support existing interpretations of a thinner lithosphere in this region due to thermal erosion by the hotspot (Nielsen et al., 2002;Steinberger et al., 2019).The hotspot track suggested by Martos et al. (2018), which was based primarily on magnetic data, trends farther north than most other tracks and may have been spatially coincident with Central Region B at ∼60 Ma.Together, their hotspot track and the distinct geophysical properties of Central Region B support the existing hypothesis that the transit of the Iceland hotspot beneath north-central Greenland produced anomalous lithospheric heat that influences present crustal properties (Alley et al., 2019).
A comprehensive evaluation of the relation between radar-observed subglacial topography and its surface expression is beyond the scope of this study.However, the remarkable similarity between simple running-mean anomalies of the radar-measured subglacial topography and surface elevation indicates that a non-negligible relation between kilometer-scale subglacial and subaerial topography is common across the GrIS (Figure 1).This relation is qualitatively consistent with that predicted from linearized small-perturbation theory for bed-tosurface transfer (Gudmundsson, 2003;Ng et al., 2018), that is, subglacial topographic perturbations that are sufficiently wide along-flow relative to local ice thickness commonly induce detectable surface expressions.That linearized theory is potentially exploitable to resolve subglacial topography at finer resolution (Ockenden et al., 2022(Ockenden et al., , 2023)), but more physically complete modeling has demonstrated that three-dimensional effects and non-linear ice rheology further complicate its interpretation (e.g., Sergienko, 2012).Regardless, it is clear that numerous subglacial topographic features across the GrIS interior-primarily subglacial valleys-are poorly resolved by existing surveys and interpolations.
Adjacent to the highlands of eastern Greenland, the ice surface of Central Region C and northeastern Central Region B suggests an underlying alpine-style subglacial geomorphology reminiscent of the Gamburtsev Subglacial Mountains in Antarctica (e.g., Figure 1e; Lea et al., 2023;Rose et al., 2013).However, the rough bed topography of this region is not yet well captured in current syntheses as it has only been sparsely surveyed so far (Morlighem et al., 2022).Traced surface lineations tend to be shorter there than elsewhere across the GrIS and cover a broad azimuthal range (Figure 3b).We suggest that the existing hypothesis that this region hosted part of the nascent GrIS during the Pliocene and models indicating that it remains ice-covered until the final phase of deglaciation are consistent with its apparently well-preserved subglacial morphology (Aschwanden et al., 2019;Bierman et al., 2016;Paxman et al., 2024).
Across the rest of Greenland, 6% (n = 64) of the traced surface lineations are both longer than 100 km and have a sinuosity less than 1.1 (Figure 3b), and all appear to overlie subglacial valleys (Figure 1).The presence of these remarkably long and straight valleys suggests low variability in bedrock resistance to channel incision across the low topographic relief of Greenland's interior (Lazarus & Constantine, 2013), a pattern that may be consistent with the predominance of Archean and Proterozoic provinces across the island.Further, there appear to be two roughly orthogonal preferential directions for many of the features/valleys, with the southwest-northeast azimuth being predominant (30-60°; Figure 3b).This pattern suggests that Greenland's tectonic history generated preferential pathways of lithological weakness for subsequent fluvial or subglacial erosion to exploit (e.g., Dühnforth et al., 2010).The converse scenario-whereby high-competence bedrock redirects valley incision-appears to be valid for at least part of Greenland, that is, the Petermann Megacanyon steers around the relatively smooth and presumably competent Petermann Anomaly (Figure 3c; Paxman, Tinto, & Austermann, 2021).
Whether this subset of long, linear valleys formed by fluvial or subglacial processes remains unclear.Where they are intersected by airborne radar sounding data, the cross-sectional morphology of these valleys shows a range of forms resembling both glacial and fluvial valleys elsewhere in the Northern Hemisphere, although definitively Ushaped cross-profiles are uncommon in the Greenlandic interior (Paxman, 2023).A fluvial origin is most likely given low long-term glacial erosion rates in the vicinity of GISP2 in central Greenland (Bierman et al., 2014), although subglacial meltwater may have also contributed to valley incision where the bed is thawed (e.g., Kirkham et al., 2022).

Conclusions
We improved knowledge of Greenland's subglacial geology by synthesizing multiple island-wide geologic, seismic, gravity, magnetic and topographic data sets that have been developed since Dawes (2009).This geophysically driven synthesis clarifies the apparent extents of Greenland's geologic provinces, but it still indicates that a greater variety of provinces exist in northern Greenland than southern Greenland.Three regions in central and northern Greenland possess sufficiently distinct geophysical expressions that they cannot yet be assigned to known provinces.While this new map is an inevitably incomplete representation of Greenland's subglacial geology, it provides a modernized framework for combined interpretation of the subglacial and solid-Earth properties of Greenland.Additionally, our observation of an extensive subglacial valley network beneath the GrIS, with many long and sub-parallel valleys, suggests that tectonic inheritance may have had a larger influence on Greenland's subglacial topography than previously recognized.This new synthesis of Greenland's subglacial geology has value for several potential applications.Foremost is the modern context it provides for new applications of existing methods (e.g., interpretation of new seismic, gravity or magnetic data sets), evaluation of aerogeophysical data sets beyond those we considered, or tracing the provenance of offshore sediments (e.g., White et al., 2016).Another application could be exploring relations between these apparent province boundaries and other subglacial properties of geophysical, glaciological, geomorphological or geochemical interest, such as geothermal heat flow (Jones et al., 2021), bedrock erodibility (Campforts et al., 2020), basal friction (Maier et al., 2022), drainage history (Jess et al., 2020;Keisling et al., 2020), or isotope geochemistry (Briner et al., 2022;Colville et al., 2011).While our results were based on conventional manual delineation of province boundaries and interpretation, machine learning techniques could also be applied to reveal such structures with reduced influence from expert biases (e.g., Colgan et al., 2022;Li et al., 2022;Rezvanbehbahani et al., 2017).

Data Availability Statement
Most data sets used in this study are openly archived (Table S1 in Supporting Information S1), and some through QGreenland v3.0.0 (Fisher et al., 2023).Version 2 of the Grid files of the Total isostatic response to the complete unloading of the GrIS is publicly archived (Paxman et al., 2022b).All data sets that were previously unavailable (seismic tomography from Darbyshire et al., 2018) and those we generated (GrIMP flow-aware hillshade and associated code, traced surface lineations, traced boundaries in geophysical fields and the new subglacial geologic province map), are publicly archived (MacGregor et al., 2024).

Figure 1 .
Figure 1.Flow-aware hillshade of the GrIS derived from the GrIMP DEM.(a) Map of the whole island with manually traced lineations (Figure S1 in Supporting Information S1 shows uninterpreted hillshade).(b-g) Zoom-ins of selected regions with bed-elevation anomaly Δz b (subtracted from 5-km running mean) shown where |Δz b | ≥ 50 m in all 1993-2019 NASA radar-sounding data (CReSIS, 2024; thin white lines).The scale is the same for all six panels.Selected 50-km-long segments of radar-sounding tracks (thick white lines) are shown in the rightmost column, along with the surface-elevation anomaly.Cyan arrows on maps show mean ice-flow direction along those segments; the mean ice thickness in meters (H) for each segment is shown in the lower right of each profile.