The timing of fjord formation and early glaciations in North and Northeast Greenland

Understanding the behavior and long-term stability of the Greenland Ice Sheet (GrIS) is key for predicting its future course. However, our understanding of the timing and extent of early glaciations and their influence on long-term landscape formation in Greenland remains fragmented. Provisional signs of glaciation in Northeast Greenland date back to the Eocene-Oligocene transition (Eldrett et al., 2007; Tripati et al., 2008; Bernard et al., 2016). Glaciation resumed from the mid-Miocene (14–11 Ma; Helland and Holmes, 1997; Thiede et al., 1998; Winkler et al., 2002; Berger and Jokat, 2008) and intensified markedly from the late Miocene (Larsen et al., 1994; Butt et al., 2001; St. John and Krissek, 2002; Pérez et al., 2018). For the Pleistocene, some studies suggest a persistent GrIS for the past several million years (Bierman et al., 2014, 2016), whereas others suggest that Greenland was deglaciated almost completely for extended periods during the second half of the Pleistocene (Schaefer et al., 2016).

Largely ice-free conditions in Greenland have also been suggested for the late Pliocene–early Pleistocene, based on the marine Kap København Formation (KKF; Funder et al., 1984, 1985, 2001; Fig. 1B) with fossils that point to the presence of boreal forest-tundra and summer temperatures as much as 6–8 °C higher than at present (Bennike and Böcher, 1990; Penney, 1993; Bennike et al., 2010). The KKF sediments were deposited in a shallow marine setting (water depth ∼50 m), but are now found ∼40–230 m above sea level (a.s.l.) (Funder et al., 1984), suggesting that they have been uplifted since deposition at ca. 2.5 Ma. Above 165 m a.s.l., the KKF may have experienced post-depositional glaciotectonic dislocation (Funder et al., 1984), although most sediment successions are undisturbed (Bennike, 1990). Farther south at Scoresby Sund, the similar Lodin Elv Formation (LEF) extends from 24 to 62 m a.s.l. (Feyling-Hanssen et al., 1983), whereas three additional occurrences on Kap Rigsdagen (Funder and Hjort, 1980), Île de France (Landvik, 1994; Bennike et al., 2002), and Store Koldewey (Bennike et al., 2010) show intermediate elevations of 100 m, 35–80 m, and 110–130 m a.s.l., respectively (Fig. 1B).

In this study, we use the present elevations of the late Pliocene–early Pleistocene marine deposits to constrain the timing of isostatic uplift in response to fjord formation in North and Northeast Greenland. Previous efforts to estimate erosion-induced uplift in Greenland (Medvedev et al., 2008, 2013, 2018) have not constrained the timing of erosion or the inception and spatial evolution of the GrIS.

We use the concept of geophysical relief (Small and Anderson, 1998; Champagnac et al., 2007) as a proxy for erosion. This metric defines erosion from elevation differences between present topography and a smooth surface connecting high points in the landscape within a given radius. Previous applications of geophysical relief have estimated fjord erosion in Scandinavia using a radius of ∼2 km (Steer et al., 2012). However, because the largest fjords in North and Northeast Greenland are much wider (as wide as 50 km) than those in Scandinavia (up to 10 km wide), we use a radius of 10 km for fjords and a radius of 5 km for onshore regions, reflecting the smaller wavelength of onshore valleys compared to the much wider fjords. We note that wide bay-like regions such as the outer parts of Scoresby Sund would only be filled to sea level in their central parts, and can be regarded as intermediate regions between fjord and shelf. We have tested our approach on a synthetic landscape from a glacial landscape evolution model (Egholm et al., 2017), and were able to match the known total erosion within ∼6% (Fig. DR1 in the GSA Data Repository¹). Unlike previous studies exploring erosion-induced isostatic uplift (e.g., Medvedev et al., 2008; Steer et al., 2012), we also estimate shelf erosion. Here, we use a radius of 40 km, reflecting the large-scale troughs carved by ice streams into older sediments.

We calculate geophysical relief based on the digital elevation model BedMachine v3, which defines ice thickness and subglacial bed topography for Greenland, and the adjacent shelf bathymetry, with a grid spacing of 150 m (Morlighem et al., 2017). For a few regions (including Independence Fjord), synthetic fjord geometries define the bathymetry (Williams et al., 2017). Specifically, we use bathymetry and subglacial bed topography corrected isostatically for the loading of the present ice sheet (Fig. 1A). This correction results in ∼900 m of uplift below the summit of the GrIS, but has a negligible influence (<10 m) on the late Pliocene–early Pleistocene marine sediments (Fig. 1B). For calculating flexural isostasy, we assume a thin elastic plate model, solved using a spectral two-dimensional cosine transform (Makhoul, 1980), with zero-deflection gradients at the boundaries. Guided by studies of effective elastic thicknesses (EETs) in Scandinavia (Pérez-Gussinyé et al., 2004), we use an EET of 20 km, but explore also 40 km (Figs. 1 and 2; Figs. DR2–DR3). We use a Young’s modulus of 70 GPa, Poisson ratio of 0.25, and densities of 950 kg/m³, 1020 kg/m³, 2500 kg/m³, 2670 kg/m³, and 3300 kg/m³ for ice, water, sediment, eroded material, and mantle, respectively.

The spatial uplift pattern we infer from the elevated marine deposits is opposite to the uplift expected from mantle convection and the Icelandic plume thermal anomaly (Steinberger et al., 2015). Previous work suggests that mantle dynamics has resulted in 200–800 m of uplift in easternmost Greenland since 5 Ma, with decreasing values toward the north (Steinberger et al., 2015). However, the low elevations of the LEF at Scoresby Sund (24–62 m a.s.l.) suggest that most of this uplift had happened by 2.5 Ma. So, mantle dynamics can in principle explain the current elevation of the LEF, whereas this is not the case for the KKF farther north, with its distal location to the Icelandic plume thermal anomaly and negative dynamic topography (Steinberger et al., 2015).

Estimates of eustatic sea level during deposition of the late Pliocene–early Pleistocene marine deposits are uncertain. Reported global values for the mid-Pliocene thermal optimum (ca. 3.3–3.0 Ma) range from +10 to +40 m, but these reconstructions are affected by local non-eustatic sea-level changes such as residual isostatic adjustments associated with recent glaciation (Raymo et al., 2011). Particularly for regions proximal to the GrIS, such effects could diminish or even reverse the reported mid-Pliocene sea-level values (Raymo et al., 2011). Consequently, we do not consider eustatic sea-level changes explicitly for the interpretation of the elevated marine deposits.

By considering geophysical relief as a proxy for erosion, we get by far the most erosion in fjords and large valleys and less erosion in the adjacent areas, with estimated fjord erosion reflecting local variations in relief (Fig. 1). The dramatic relief in the inner parts of the Scoresby Sund fjord system results in erosion estimates of >3500 m (Fig. 1, profile B-B′), whereas the lower relief in Independence Fjord toward the north results in erosion estimates of <1500 m (Fig. 1, profile A-A′). Overall, the estimated erosion varies between 0 and 3678 m for onshore regions and fjords, with an average of ∼350 m, and between 0 and 613 m for the shelf, with an average of 167 m.

The erosion estimates derived from geophysical relief indicate up to ∼940 m of flexural isostatic uplift, with the largest values around Scoresby Sund (Fig. 2). A net surface uplift is achieved when the estimated erosion is smaller than the isostatic uplift, and accounts for up to ∼935 m of surface uplift in the high regions surrounding Scoresby Sund.

In order to compare our estimates of erosion-induced uplift with the present elevation of the marine deposits, we evaluate the isostatic deflection in 30 × 30 km windows around each locality (Figs. 2B–2F, histograms). We find good agreement between the erosion-induced uplift and the present elevation of the KKF (Fig. 2B), whereas there is a mismatch for the other localities. The mismatch increases toward the south, reaching >450 m for the LEF (Figs. 2C–2F).

Our estimates of erosion-induced isostatic uplift are some 15% lower than previous estimates from the region (Medvedev et al., 2008, 2013, 2018). Compared to these previous findings, our approach leads to more conservative estimates of erosion, assuming only a modest component of erosion (<500 m) in regions outside fjords and large valleys, which is what we expect for glacial erosion (Sugden, 1974; Egholm et al., 2017; Strunk et al., 2017; Andersen et al., 2018).

We find that both complete fjord formation and additional modest erosion in adjacent areas are needed in order to explain the amount of uplift observed for the KKF since 2.5 Ma (Fig. 2B), even when considering a potential sea-level change of some 10–40 m. This suggests that most glacial erosion in North Greenland occurred within the last 2.5 m.y. Assuming that our estimated erosion has taken place at a constant rate since deposition of the KKF, fjord erosion rates amount to >0.5 mm/yr in many places, with rates locally approaching 0.8 mm/yr (Fig. 3). Estimated erosion rates are significantly lower at high elevations between fjords, where the geophysical relief method predicts limited erosion (Figs. 3B and 3C). If we assume that fjord incision is limited to glacial periods—suggested to constitute a maximum of 45% of the last 2.5 m.y. in western Greenland (Strunk et al., 2017)—erosion rates would have been >1 mm/yr. A younger age of ca. 2.0 Ma for part of the KKF (Bennike et al., 2010) would raise the rate further to ∼1.4 mm/yr. We stress that these rates are minimum estimates, as the KKF and related sequences could have extended to higher elevations previously, as indicated by their upper erosive surfaces (e.g., Feyling-Hanssen et al., 1983; Funder et al., 1984).

The large mismatch we find between our estimates of erosion-driven uplift and the present elevation of the LEF in Scoresby Sund suggests that at least some of the predicted erosion must have taken place prior to deposition of the LEF at ca. 2.5 Ma. Either the majority of the erosion and fjord formation had already occurred in the broader region by 2.5 Ma, or extensive fjord erosion had propagated a significant distance inland from the Lodin Elv locality at the time of deposition. In order to explore the implications of this latter view of fjord propagation, we split our isostatic uplift estimates into four components, based on the erosion from four separate domains, with varying proximity to the LEF (Fig. 4, regions 1–4). From this, it becomes evident that erosion of the outermost regions 3–4 would have contributed considerably to uplift of the LEF, and a significant portion of this erosion must have taken place prior to 2.5 Ma. In contrast, erosion in regions 1–2 would not have contributed to uplift, but may have induced minor subsidence at the LEF locality, and therefore we cannot resolve the timing and amplitude of this erosion. However, based on numerical modeling studies (e.g., Pedersen et al., 2014; Egholm et al., 2017), we find it plausible that glacial erosion and fjord formation have propagated inland since glacial erosion commenced and as the fjord system has evolved.

As noted previously, mantle dynamics can in principle explain the current elevation of the LEF (Steinberger et al., 2015). This would point to a larger component of pre-Pleistocene glacial erosion in the Scoresby Sund region, with more erosion from regions 3–4 taking place prior to 2.5 Ma (Fig. 4). Finally, we note that a higher sea level during deposition of the late Pliocene–early Pleistocene marine deposits would imply that less erosion has occurred since 2.5 Ma, whereas a lower sea level would imply the opposite. However, the current lack of an accurate late Pliocene–early Pleistocene sea-level reconstruction prevents us from assessing these effects further (Raymo et al., 2011).

Our study implies that extensive glacial erosion and fjord formation in North and Northeast Greenland did not commence synchronously. In North Greenland, extensive glacial erosion was limited to the Pleistocene, with the relief of the Independence fjord system developing after deposition of the KKF at ca. 2.5 Ma. In contrast, extensive glacial erosion commenced earlier farther south in Scoresby Sund, Northeast Greenland, with only limited glacial erosion taking place after 2.5 Ma near the LEF location. This implies that the main source of ice-rafted debris (IRD) material dating back before the Pleistocene was Northeast Greenland, or regions outside of our study area. This view is consistent with recent work using Fe-oxide fingerprinting pointing to Northeast and Southeast Greenland as the main sources for Eocene–Oligocene IRD (Tripati and Darby, 2018), and other studies suggesting extensive glacial erosion in particular in Northeast and Southeast Greenland prior to the Pleistocene (e.g., Larsen et al., 1994; Eldrett et al., 2007; Bernard, et al., 2016; Bierman et al., 2016). The lack of extensive, relief-forming glacial erosion in North Greenland near the KKF prior to 2.5 Ma implies that widespread and fjord-forming glaciations did not commence in this region before 2.5 Ma. Farther south in Scoresby Sund, our results imply that widespread and prolonged glaciations commenced earlier, with a main phase of glacial erosion proximal to the LEF taking place prior to 2.5 Ma and probably as early as the late Miocene (ca. 7 Ma), or even the Eocene-Oligocene transition, as indicated by the IRD records offshore Northeast Greenland (e.g., Larsen et al., 1994; St. John and Krissek, 2002; Tripati and Darby, 2018).

We have evaluated the timing of extensive glacial erosion and fjord formation in North and Northeast Greenland using elevated marine sediments of late Pliocene–early Pleistocene age. By constraining the timing of fjord formation, we constrain the timing of widespread and prolonged glaciation in these regions. We find that in North Greenland, the full relief of the Independence Fjord system must have developed since ca. 2.5 Ma in order to explain the current elevation of the marine KKF by the flexural isostatic response to erosional unloading. This implies glacial erosion rates in the fjords >0.5 mm/yr over a 2.5 m.y. time period, or rates >1 mm/yr if erosion is limited to glacial periods which constitute a maximum of 45% of the last 2.5 Ma. Consistent with patterns in present topography, our study suggests that extensive and prolonged glaciation did not extend to North Greenland before 2.5 Ma, whereas glaciation and glacial erosion commenced before 2.5 Ma farther south around Scoresby Sund.

This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement 745669. Pedersen and Larsen thank the Villum Foundation’s Young Investigator Programme for their support (grants 15467 and 023440, respectively). We thank Sergei Medvedev and two anonymous reviewers for detailed and constructive comments.