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

The timing and extent of early glaciations in Greenland and their co-evolution with the underlying landscape remain elusive. Here we use the concept of geophysical relief to estimate fjord erosion and the subsequent flexural isostatic response to erosional unloading in Northeast and North Greenland between Scoresby Sund (70 ° N) and Independence Fjord (82 ° N). The timing of erosion and isostatic uplift is constrained by marine sediments of late Pliocene–early Pleistocene age that are now exposed on land between ~24 m and 230 m above sea level. By determining the timing of fjord formation, we can establish the early history of the Greenland Ice Sheet. We find that the Independence Fjord system must have formed by glacial erosion at average rates of ~0.5–1 mm/yr since ca. 2.5 Ma in order to explain the current elevation of the marine Kap København Formation by erosion-induced isostatic up-lift. In contrast, fjord formation in the outer parts of Scoresby Sund commenced before the Pleistocene, most likely in late Miocene, and continued throughout the Pleistocene by progressive inland fjord formation. Our results demonstrate that the inception of the Greenland Ice Sheet began in the central parts of Northeast Greenland before the Pleistocene and spread to North Greenland only at the onset of the Pleistocene.


INTRODUCTION
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 longterm 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(Bierman et al., , 2016, whereas others suggest that Greenland was deglaciated almost completely for extended periods during the second half of the Pleistocene (Schaefer et al., 2016).
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(Medvedev et al., , 2013(Medvedev et al., , 2018 have not constrained the timing of erosion or the inception and spatial evolution of the GrIS.

METHODS
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 , and were able to match the known total erosion within ~6% (Fig. DR1 in the GSA Data Repository 1 ). 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 de-fine 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 3 , 1020 kg/m 3 , 2500 kg/m 3 , 2670 kg/m 3 , and 3300 kg/m 3 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 sealevel changes explicitly for the interpretation of the elevated marine deposits.

RESULTS
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 erosioninduced 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).

Glacial Erosion and Isostatic Uplift
Our estimates of erosion-induced isostatic uplift are some 15% lower than previous estimates from the region (Medvedev et al., 2008(Medvedev et al., , 2013(Medvedev et al., , 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 deposi tion. 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).

Inception and Evolution of the Greenland Ice Sheet
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 Feoxide finger print ing 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 ero-sion in North Greenland near the KKF prior to 2.5 Ma implies that widespread and fjordforming 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).

CONCLUSION
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 topog raphy, 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.

ACKNOWLEDGMENTS
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.