Neogene‐Quaternary Uplift and Landscape Evolution in Northern Greenland Recorded by Subglacial Valley Morphology

The landscape hidden beneath the Greenland Ice Sheet remains one of the most sparsely mapped regions on Earth, but offers a unique record of environmental conditions prior to and during widespread glaciation, and of the ice sheet's response to changing climates. In particular, subglacial valleys observed across Greenland may preserve geomorphological information pertaining to landscape and ice sheet evolution. Here we analyze the morphology of a subglacial valley network in northern Greenland using bed elevation measurements acquired during multi‐year airborne radio‐echo sounding surveys. Channel profile morphologies are consistent with a primarily fluvial origin of the network, with evidence for localized modification by ice and/or meltwater. Gravity and magnetic anomalies suggest that the spatial organisation of the valley network is influenced by regional‐scale geological structure, implying a long‐lived and well‐established hydrological system. We also document two knickzones in the valley longitudinal profile and terraces above the channel floor in the lower course of the network. These observations, combined with stream power modeling, indicate that northern Greenland experienced two episodes of relative base level fall during the Neogene (∼150 m at ca. 12–3.7 Ma and ∼380 m at ca. 8.2–2.8 Ma) that resulted in channel profile adjustment via incision and knickzone retreat. The timing of the inferred base level fall correlates with other onshore and offshore records of uplift, denudation, and/or relative sea level change, and we suggest that tectonic and/or mantle‐driven uplift played an important role in the genesis of the modern landscape of northern Greenland.

. Regional setting. (a) Ice surface velocity (Joughin et al., 2010). The 50 m/yr velocity contour is marked in black to demarcate regions of fast-flowing ice. Solid lines mark the ice drainage divides (Zwally et al., 2012). Major outlet glaciers referred to in the text are labeled. White lines mark the subglacial valley network mapped in this study. (b) Bed topography and bathymetry (Morlighem et al., 2017). The sea level (0 m) contour is marked in black. Subglacial channel/canyon networks referred to in the text are labeled. A group of features in the Greenland landscape that can potentially yield new insights into this debate are the subglacial "trench", "trough", "canyon", "valley", and "channel" systems that have been documented by a number of recent studies (Bamber et al., 2013;Cooper et al., 2016;Dam et al., 2020a;Livingstone et al., 2017). Channel networks have been mapped upstream of the Hiawatha, Humboldt, and Petermann glaciers and Ilulissat Icefjord, at Uparuaqqusuitsut, and in Gåseland-Milne Land (Figure 1b). They appear to exert a control on contemporary ice flow (Cooper et al., 2016) (Figure 1a) and also facilitate the connectivity of low-lying bed from the ice sheet margin deep into the Greenland interior (Figure 1b), which may result in ice margins being continuously exposed to marine-driven melt as they retreat (Morlighem et al., 2017).
However, the history of these channel networks remains unclear. Some studies have suggested that the morphology of the networks indicates that they are inherited paleo-fluvial drainage features dating from prior to widespread glaciation (Bamber et al., 2013;Cooper et al., 2016). Moreover, geological evidence around the Greenland margins has been used to infer that canyons present in the modern landscape ( Figure 1) may have first developed as sediment pathways as early as the Mesozoic (Dam et al., 2020a(Dam et al., , 2020b. Conversely, numerical modeling has been used to suggest that the canyon systems may have been excavated by repeated sudden glacial outburst flood events during the Pliocene-Pleistocene (Keisling et al., 2020a), and/or by the erosive action of subglacial meltwater (Chambers et al., 2020). These contrasting hypotheses for channel incision are not necessarily mutually exclusive; it is likely that individual channels have, to varying extents, been influenced by a combination of fluvial, meltwater, and glacial erosion over the course of their history (Cooper et al., 2016;Keisling et al., 2020b;Livingstone et al., 2017).
Elsewhere on Earth, channel profiles have been highly informative in aiding our understanding of the tectonic, geological, erosional, and climatological history of a landscape, and in constraining the timing and pattern of uplift and denudation (Miller et al., 2013;Roberts & White, 2010). Analysis of the geomorphological and geological characteristics of the subglacial valley networks beneath the GrIS may therefore offer a valuable insight into the formation and antiquity of the modern Greenland landscape.
For this study, we aimed to conduct the first systematic appraisal of valley topography and geomorphology for the channel network located in northern Greenland (Figure 1a), which includes the Petermann "mega-canyon" described using early compilations of RES data (Bamber et al., 2013). We focussed on this system because it is the most extensive subglacial valley network in Greenland and because much of its course is situated beneath the slow-moving, cold-based interior of the GrIS (Figure 1) (Joughin et al., 2010;MacGregor et al., 2016), meaning there is potential for preservation of inherited landscape features. Our analysis was applied to the extensive archive of RES, gravity, and magnetic data acquired by Operation IceBridge and historical airborne geophysical surveys across Greenland (Figure 1c).
The aims of this study are threefold: (a) to describe the geomorphology of the subglacial channel system in northern Greenland as imaged by RES data, (b) to use gravity and magnetic data to assess the relationship between the valley network and the regional geological structure, and (c) to develop a model of the temporal evolution of the channel system and constrain the history of valley incision. Together, these aims will allow us to better understand the process(es) involved in the formation and evolution of the valley network, and to constrain the history of uplift, denudation, and ice sheet behavior in northern Greenland.

Valley Identification and Network Extraction
To map valley morphology in northern Greenland, we utilized the freely available Operation IceBridge (OIB; and Center for Remote Sensing of Ice Sheets (CReSIS;1993) RES databases (Figure 1c). Although the most recent BedMachine DEM (version 4) has a nominal horizontal resolution of 150 m (Morlighem et al., 2017), in the ice sheet interior the grid commonly relies on interpolation between RES profiles with large and irregular spacing. Consequently, extracting topographic profiles directly from the BedMachine v.4 DEM can result in (a) aliasing, smoothing, or loss of relatively subtle, yet important, details present in the original RES data, and (b) introduction of spurious gridding artifacts into the extracted profile. OIB bed elevation data have been shown to be reliable through repeat-track surveying, and are typically able to resolve bed features down to a horizontal length-scale of ∼30 m with a vertical uncertainty of 10 m . Using the original 1. Elevation of the valley floor. 2. Mean average of the elevations of the valley rim on both sides of the channel. 3. Incision depth, defined as the difference between (2) and (1). 4. Valley width. For consistency, and since the "bankfull" datum for the channel is unclear, we defined the valley width as the horizontal distance between the valley rims. Since the radar lines will usually cross the channel at an oblique angle and therefore overestimate the true valley width, we did not determine the width for any location where the acute angle between the flight line and the channel was <60°. We also applied a 13% uncertainty to the derived widths in locations where this angle was >60° (given that sin(60°) ≅ 0.87).
Elevation data were extracted with their in situ values (i.e., situated beneath the GrIS), and were also adjusted to account for the isostatic rebound resulting from (a) the complete removal of the modern GrIS ( Figure S4 in Supporting Information S1), and (b) the incision of near-coastal fjords following glacial inception (Medvedev et al., 2013;Pedersen et al., 2019) ( Figure S5 in Supporting Information S1). For the isostatic rebound calculations, we used an elastic plate flexure model assuming a spatially variable effective elastic thickness of the lithosphere (Steffen et al., 2018), a non-viscous fluid mantle, and uniform ice, bedrock, and mantle densities of 920, 2,700, and 3,330 kg m −3 , respectively. The thickness of the modern ice load is well-constrained (Morlighem et al., 2017), and we estimated the thickness of glacially eroded fjord material following the approach of Pedersen et al. (2019) (Text S1 in Supporting Information S1). However, there is significant uncertainty in the magnitude and timing of fjord erosion relative to subglacial valley development, which we accommodate by calculating a maximum erosion scenario where the fjords are fully infilled, and a minimum erosion scenario where subglacial valley elevations have been unaffected by the isostatic response to fjord incision ( Figure S5 in Supporting Information S1). We adjusted the valley elevations by the mean of these two scenarios, where the fjords are 50% infilled, with error bars defining the minima and maxima. Following the isostatic correction for ice loading and erosional unloading, rebounded valley floor elevations (hereafter referred to as "ice-free" elevations) were used to construct a longitudinal profile, where elevation is plotted as a function of distance along the valley.

Potential Field Anomalies
Where relevant, we compared the mapped valley network to the regional magnetic and gravity anomalies to examine the influence of geological structure on channel routing and morphology. We derived these potential field anomalies from circum-Arctic magnetic and gravity compilation grids (Gaina et al., 2011), and also used magnetic and gravity anomaly line data from OIB flights where available. The Bouguer anomaly, which removes the gravity effects of the ice surface and bed topography, was computed by subtracting the Bouguer correction from the free-air anomaly, assuming standard air, ice, and rock densities of 0, 920, and 2,700 kg m −3 , respectively (Parker, 1972). We also performed 2D forward modeling along selected sections of the airborne gravity anomalies to assess the relationship between particular topographic features of interest and the subglacial geology.

Methods: Longitudinal Profile Modeling
In this section, we test whether the observed valley longitudinal profile in northern Greenland is consistent with initial formation by fluvial incision, and whether it can be used to constrain regional base level and incision history. A commonly-used framework for modeling the spatiotemporal evolution of fluvial longitudinal profiles is to assume that the rate of change of profile elevation (z) with distance (x) and time (t) over geological timescales is primarily controlled by the balance between uplift (U) and erosion (E) (Whipple & Tucker, 1999) We exploit the well-established empirical formulation for erosion in the fluvial domain (Whipple & Tucker, 1999) where A(x) is the upstream drainage area, and A(x) m is a proxy for the integrated discharge at any point along the channel (Whipple & Tucker, 2002). The first term on the right-hand side of Equation 2 is advective and represents channel incision. The second term represents diffusion of bedrock slopes, but we have followed previous studies and assumed that over the long spatial and temporal scales being considered here, diffusion is insignificant relative to erosional advection, and can be neglected (Roberts, Paul, et al., 2012;Stephenson et al., 2014). By combining Equations 1 and 2, we have: This is the 1D stream power incision equation. It contains three geomorphological parameters that affect channel profile evolution: v is the advective coefficient, and m and n are dimensionless exponents. The value of each of these parameters remains debated and appears to be dependent on the regional climatic, lithological, and tectonic setting (Duvall et al., 2004;Harel et al., 2016).
Our approach is to simplify the stream power model while maintaining sufficient complexity to capture the major processes operating on the spatial and temporal scales of interest. We first constrain the values of the erosional parameters (Table 1), assuming that they remain constant in space and time; quantification of spatiotemporal parameter variability is beyond the scope of this study. We also note that this modeling framework only applies to time-dependent fluvial incision, and not to catastrophic flood events or glacial erosion.

m/n Ratio
Estimation of the ratio between the exponents m and n involves integrating the stream power equation from base level to a given upstream observation point (Perron & Royden, 2013). If we assume that the profile is in steadystate (dz/dt = 0), and U and v are spatially uniform, Equation 3 can be re-written as and A 0 is a reference drainage area. is a transformed horizontal co-ordinate measured from the drainage basin outlet in an upstream direction. When is plotted against elevation, Equation 6 describes a straight line if the appropriate m/n ratio is chosen. We therefore calculated along the longitudinal profile, with the upstream drainage basin area calculated using a hydrological flow routing algorithm on the rebounded (ice-free) BedMachine v.4 DEM ( Figure S6 in Supporting Information S1) (Schwanghart and Scherler, 2014). We then used least-squares regression to determine the m/n ratio that provides the strongest linear fit for the -elevation plot (Perron & Royden, 2013).
Previous studies show that when observed and predicted river profiles are compared, a global misfit minimum occurs at or close to n = 1 (Paul et al., 2014;. However, values of n > 1 have been reported, for example, in regions of low mean annual temperature (Harel et al., 2016). Since this likely applies to the interior of Greenland, we tested the sensitivity of our model results to values of n greater than 1 (and the corresponding values of m for our best-fitting m/n ratio).

Advective Coefficient and Vertical Erosion Rate
If diffusion is neglected, Equation 2 can be simplified and rearranged to give an expression relating the parameters v, m, and n for a particular profile The advective coefficient, v, encapsulates the ability of a river to erode, and is modulated by a number of factors including precipitation, lithology, channel width, and sediment flux (Duvall et al., 2004;Stock & Montgomery, 1999). Equation 8 shows that v depends not only on the vertical erosion rate (E), but also on the exponents m and n (Croissant & Braun, 2014;. This formulation allows us to determine v for realistic combinations of m, n, and E. While m/n is constrained from integral analysis of the channel profile (Section 4.1.1), it is more challenging to determine E empirically in northern Greenland. However, a recent global compilation of basin-averaged denudation rates derived from cosmogenic Beryllium-10 ( 10 Be) nuclide concentrations shows that the median value of E is between 50 and 200 m/Myr for arid, temperate, cold, and polar climate zones (Harel et al., 2016), which likely encompasses the realistic range of possibilities for the climate of northern Greenland prior to glaciation (Dolan et al., 2015;Funder et al., 2001). In other valley and canyon systems in North America, reported longterm (Myr timescale) average incision rates are typically on the order of 50-100 m/Myr (Bender et al., 2018;Crow et al., 2014;Karlstrom et al., 2008;Miller et al., 2013). On this basis, we performed our calculations for E = 50-100 m/Myr, and used sensitivity analysis to assess the relative influences of E, m, and n on longitudinal profile morphology and evolution.

Knickzone Retreat and Base Level History
If we assume that uplift is a function of time only, the 1D stream power equation can be reduced to In this formulation of the stream power equation, longitudinal profile evolution can be viewed in terms of the development and migration of knickzones (regions of steep channel gradient). A fall in base level causes the river to incise down into the bed, forming a knickzone that subsequently propagates upstream. The time taken for a knickzone to retreat to a given position along the longitudinal profile after an episode of base level fall (τ) can be estimated by integrating Equation 9 assuming that the upstream drainage area can be approximated as A ∼ x 2 (Pritchard et al., 2009) where x k is the distance from the channel source to the knickzone. We assume that the total length of the channel (L) is ∼1,700 km, which reflects the distance between the highest-elevation point (closest the drainage divide) and the modern calving front, which aligns with the drainage basin outlet and likely approximates the location of the paleo-coastline prior to glacial fjord incision (Medvedev et al., 2013).
Our approach was to use forward modeling (via Equation 9) to determine the base level history (U(t)) that can reproduce the shape of the observed longitudinal profile. Given the absence of constraints on the initial shape of the longitudinal profile, we assumed a typical steady-state power law profile graded to sea level, although we acknowledge the uncertainty associated with the possible initial state of the landscape. For simplicity, the uplift term (U) in Equation 9 can be conceptualized as representing a change in base level, which may be triggered by uplift/subsidence of the land surface at or near the river outlet and/or a fall/rise in local sea level. We assume that episodes of base level change occur at a linear rate and are spatially uniform. Previous studies have used an inverse method to derive a continuous regional uplift history for a large family of river profiles (Czarnota et al., 2014;Roberts & White, 2010). However, given the limited observational data, lack of direct age constraints, and uncertainties associated with the stream power parameters, we are unable to confidently resolve the additional spatiotemporal complexity associated with the inverse method. The aim instead was to identify the simplest base level history that matches the observed longitudinal profile, and assess the sensitivity of the results to the model parameters (Table 1).

Valley Morphology
The subglacial valley network in northern Greenland comprises a trunk channel ("main stem") that is ∼1,500 kilometres long between the first and last mapped points, and a number of comparatively short (<300 km-long) tributaries that meet the trunk channel along its course ( As noted in Section 4.2, the lower end of the longitudinal profile was extended out to the modern calving front, which coincides with the modeled drainage basin outlet (Figure 2a), increasing the total length of the profile to 1,700 km ( Figure 2f). The total drainage basin area is 5 × 10 5 km 2 , which represents ∼25% of the area of the Greenland landmass. In total, we identified 85 unique locations where the trunk channel is intersected and imaged by a RES line (not counting repeat flight tracks). We acknowledge that the mapped channel network in the interior will not be complete since smaller upland channels in particular may remain unresolved due to gaps in the RES data coverage or overprinting by localized subglacial erosion and deposition.
Along the channel network, valley cross profiles exhibit predominantly V-shaped geometries ( Figure 2). The mean incision depth in the trunk channel is ∼300 m, and the mean width is ∼4.9 km. Maximum incision depths of ∼700 m and widths of ∼7 km are attained in the lower course of the valley network (Figures 3c and 3d). The channel top width increases approximately linearly with distance downstream, with least squares regression indicating a rate of valley widening of 3.8 ± 0.27 m/km, albeit with a degree of scatter (R 2 = 0.58; Figure 3c). The longitudinal profile also exhibits two regions of overdeepening, where the bed descends into an enclosed "minimum" (Figure 3e; pink bars). The first is a ∼150 km-long section at 720 profile-km, and the second is a ∼80 km-long section at ∼1,130 profile-km. The valley floor in both regions is overdeepened by 100-200 m relative to the adjacent profile (Figure 3e), resulting in a pronounced "high" in the depth of channel incision (Figure 3d). Both locations appear to also coincide with steep gradients in the magnetic anomaly (Figure 3b), and the lower overdeepening is marked by a step increase in the Bouguer anomaly from −100 to −75 mGal ( Figure 3b) and a cluster of anomalously narrow channel widths (Figure 3c).
The valley profile is also characterized by two knickzones where the bed slope increases appreciably and the profile "steps" between regions of relatively uniform elevation and low gradient (Figure 3e; gray bars). The upper knickzone is situated at 910 ± 20 profile-km, with a magnitude of ∼150 m, and the lower knickzone is situated at 440 ± 40 profile-km, with a magnitude of ∼380 m (Figure 3e). Although smaller-scale undulations also exist along the channel, they are not systematic and likely reflect small uncertainties associated with RES-derived bed elevations and/or minor localized modification by (e.g.,) erosion and sedimentation.  An exception is found at 370 profile-km, where a narrow section of the longitudinal profile is ∼200 m higher than the surrounding channel floor (Figure 3e). This vertical offset is too large to be attributed to uncertainties in the RES data, and the bed is well-imaged in two separate radar profiles. Such a large, abrupt jump in elevation is unexpected for a valley longitudinal profile. Inspection of the Bouguer gravity anomaly over this section of the channel reveals a 3-4 mGal "low" relative to the immediate north and south (Figure 4e). Gravity modeling indicates that the wavelength and amplitude of the "low" are consistent with a shallow mass deficit caused by a ∼200 m infill of low-density (2,200 kg m −3 ) material within the channel relative to the background bedrock density of 2,700 kg m −3 (Figure 4b). This would place the base of the low-density infill at an elevation comparable to the adjacent parts of the profile (Figure 4f). The implication is that this section of the channel is infilled by Bouguer gravity anomaly and ice-free bed topography profile through the channel dam. Red star marks a reference location where the Bouguer anomaly and bed elevation in profiles a and b equate, indicating a consistent surface rock density. The Bouguer anomaly "minimum" is 3-4 mGal lower than the regional trend (see dashed red lines in profiles a and c and regional trend in panel e). This is consistent with a ∼200 m infill of low-density (2, low-density sedimentary material, which appears to dam the valley. We therefore temporarily neglect these two points for the longitudinal profile modeling. Towards the northern (downstream) end of the trunk channel, RES profiles show that the valley side is marked by a number of rises/scarps and flats/benches ( Figure 5). Some of these terraces are paired, while others are discernible on just one side of the valley. Although there are a number of scarp-like features visible, the elevations of the most clearly resolved terraces are grouped into two clusters, indicating that there are two separate surfaces along the valley sides ( Figure 5). These two terraces (referred to as T1 and T2), are characterized by very shallow gradients, dipping down to the north at ∼0.05°.

Longitudinal Profile Evolution
The -elevation plot is best linearized by an m/n ratio of 0.32, with an R 2 value of 0.90 ( Figure 6). Because knickzone migration disturbs the equilibrium of the channel profile, we applied the straight-line fitting to the segment of the profile above the downstream knickzone, which is clearly defined in the -elevation plot and more suitable for estimating m/n (Figure 6a). Values of m/n between 0.30 and 0.40 yield similarly good linear fits in the -elevation plot (Figure 6b). Moreover, sensitivity analysis showed that values of n between 1.0 and 1.2 produced a good fit with the observed longitudinal profile morphology, but when n > 1.2 the model is unable to satisfactorily match the geometry of the knickzones ( Figure S7 in Supporting Information S1). We therefore assumed 1.0 < n < 1.2 and 0.3 < m/n < 0.4; and hence 0.3 < m < 0.5, which is in good agreement with the typical range of m values reported for bedrock rivers (Harel et al., 2016;Stock & Montgomery, 1999).
Using Equation 8, with a measured average channel slope of 7 × 10 −4 m m −1 , a catchment area of 5 × 10 5 km 2 , m/n = 0.35 ± 0.05, n = 1.1 ± 0.1, and an assumed erosion rate of 75 ± 25 m/Myr, the corresponding range of v is 0.73-44 m (1−2m) Myr −1 (Figure 6c). Although no prior constraints for v exist in Greenland, values determined for numerous fluvial valley and canyon networks in neighboring North America, Iceland, and Scandinavia fall in the range 10 0 -10 2 m (1−2m) Myr −1 Miller et al., 2013;Pedersen et al., 2018;Stucky de Quay et al., 2019), which indicates that the parameters derived in this study are realistic. Moreover, these v values are consistent with those associated with crystalline gneisses, granitoids, and well-indurated metasedimentary lithologies (Stock & Montgomery, 1999), which reflects the likely subglacial geology in the Greenlandic interior (Henriksen et al., 2009).
Based on the inferred geomorphological parameters, our forward modeling shows that the shape of the observed longitudinal profile can be closely matched by a scenario in which an initially steady-state fluvial profile experiences two separate phases of base level fall (e.g., via land uplift; Figure 7). The model yields a good match with the upstream curvature of the longitudinal profile (Figure 7), and also accounts for the presence of two separate knickzones, with the magnitudes of the drops in base level constrained by the height of the knickzones (∼150 m for phase 1, ∼380 m for phase 2). For the parameter space encapsulated by m/n = 0.3-0.4, n = 1.0-1.2, and  (Figures 7 and 8). The observed knickzone gradients are consistent with the two phases of base level fall having durations of ∼0.2 and ∼1.5 Myr, respectively (Figure 7). The total time taken from the onset of the base level fall to the knickzone reaching its current location is therefore 3.7-9.6 Myr (upper) and 2.8-5.5 Myr (lower; Figure 7g).  Figure 5). (f) Final model profile compared to the observed channel profile (blue inverted triangles) and terrace (colored crosses) elevations with vertical error bars. Brown line marks the elevation of the channel rim. (g) Inferred uplift history. Solid line marks the uplift history for the scenario shown in panels a-e, which uses the mean average model parameters listed in the inset. Gray shaded region shows the range of uplift histories associated with the assumed range of the parameters m/n, n, and E (see Figure 8). Blue shaded region illustrates the uncertainty in the timing of uplift given the potential halt in profile evolution following major ice sheet expansion at ca. 2.7 Ma. The geological timescale labeled along the top of the panel is for the potential scenario where profile evolution ceased at 2.7 Ma. Sensitivity testing shows that the influences of m/n and n on the inferred timing of knickzone formation are relatively minor compared to that of E (Figure 8). This is because changes in m and n are partially offset by changes in v (see Equations 8 and 10), whereas E is not modulated by any other parameters. The range of knickzone ages quoted above therefore primarily reflects the uncertainty in the value of E. We also note that there is an uncertainty of up to ±100 km associated with the length of the channel and the distance from the knickzones to the channel mouth (Figure 3e). However, the positive covariance of these two parameters results in an uncertainty ellipse (Figure 8a) and consequently a relatively minor effect on the inferred knickzone age compared to the uncertainties associated with the other geomorphological parameters (Figures 8b-8g). Although there are also likely to be other systematic uncertainties that are difficult to capture numerically, the inferred knickzone ages provide a first-order estimate of the timescales involved in the development of the observed longitudinal profile. The two-phase base level fall model can also potentially account for the formation of the terraces in the downstream section of the profile. The location and elevation of the T1 terrace at 400-500 m above sea level under ice-free conditions is broadly consistent with having formed due to channel downcutting following the first phase of base level fall (Figure 7). The T2 terrace may in turn represent the position of the abandoned channel at 200-300 m following incision after the second phase of base level fall (Figure 7). We also note that the channel rim between 450 and 650 profile-km is ∼200 m lower than along the remainder of the profile, and appears to closely align with the elevation of the initial (uplifted) graded profile and the T1 terrace ( Figure 7f). Finally, we highlight that the fluvial incision model cannot account for the overdeepenings at 720 and 1,130 profile-km (Figure 7f), which require an alternative mode of formation.

Fluvial Channel Incision
Our geomorphological analysis reveals a number of lines of evidence that support a fluvial origin for the subglacial valley network in northern Greenland. The V-shaped valley cross profiles (Figures 2b-2e) are typical of fluvial incision (Bamber et al., 2013), and the approximately linear increase of channel width with downstream distance (Figure 3c) is consistent with the scaling often observed for fluvial bedrock channels (Finnegan et al., 2005;Langston & Temme, 2019;Yanites & Tucker, 2010). The planform organisation of the channel network is also typical of fluvial systems, with dendritic branching of tributaries and junction angles of 40-80° (Figure 2a) (Grau Galofre & Jellinek, 2017). The longitudinal profile of the trunk channel exhibits the concave-up form that is a hallmark of fluvial channels, with steep bed slopes in the upland region and a gently-sloping lower course (with the exception of the knickzones), which appears to reach an ice-free base level close to present-day sea level, albeit it with some uncertainty associated with the magnitude of the flexural response to glacial erosion within the Petermann fjord (Figure 3e; Figure S5 in Supporting Information S1).
Our stream power incision modeling demonstrates that the two knickzones identified in the longitudinal profile are consistent with the adjustment of the profile (via incision) to two separate episodes of base level fall (Figure 7). This result lends further support to the interpretation that the channels were formed by an erosional regime that was responsive to changes in base level, which is consistent with a fluvial system. We also note that the lower course of the channel is incised into high terrain east of the Petermann Glacier when there is a lower-lying route to the coast available to the west (upstream of the Humboldt Glacier; Figure 1). Moreover, such a deeply incised (∼700 m; Figure 3d) gorge is atypical of the lower course of a river system. These observations suggest that the channel network was initially established at a time when the terrain was lower, and that subsequent incision and canyon formation were triggered by a fall in base level (e.g., uplift of the land surface). Other potential signatures of base level fall and channel downcutting are the terraces in the valley walls of the profile's lower course (Figure 5), which may represent the "abandoned" locations of the original valley floor prior to each episode of incision and knickzone retreat (Figure 7). An alternative scenario is that the terraces in the valley walls are instead "lithological steps" formed by differential erosion along (sub)-horizontal bedding planes. Identifying the presence of fluvial sediment deposits versus stratified bedrock atop the terraces would help to differentiate between these two scenarios, but this is challenging given the lack of outcrop exposure and the inability of the airborne gravity and magnetic data to resolve subtle density or magnetic susceptibility contrasts over these small length-scales.
However, the elevations of the terraces and their locations relative to the knickzones are consistent with what would be expected for the incision and knickzone retreat model (Figure 7). In addition, gravity anomaly modeling over the anomalously low section of the channel rim between 450 and 650 profile-km (Figure 3e) indicates that relatively thin (<300 m) low-density sedimentary drapes are present across this flat-lying region, which also exhibits low magnetic anomalies ( Figure S8 in Supporting Information S1, Figure 9). The rim may therefore represent the former location of a broad floodplain dating from prior to the first phase of base level fall and channel downcutting, in a region where low channel gradients were conducive to sediment deposition. The apparent contiguity of the channel rim in this area and the T1 terrace to the north (Figure 7f) suggests that both may share a common mode and timing of formation, indicating that the terraces formed as a consequence of base level fall and channel incision.
It is also possible that knickzones arise due to geological structure, for example, where a channel crosses an uplifting fault block or a contact between lithologies of differing erodibility (Wolpert & Forte, 2021). Although  (Gaina et al., 2011). Contour interval is 80 nT. (b) Bouguer gravity anomaly calculated by correcting the OIB airborne free-air anomaly data (Tinto et al., 2019) for the gravity effects of the ice surface and bed topography (assuming air, ice, and rock densities of 0, 920, and 2,700 kg m −3 respectively). Contour interval is 10 mGal. The mapped subglacial channel network is marked in white, with colored circles showing the incision depth along the trunk channel. The two knickzones are labeled, as is the lower overdeepening, and the linear geological boundary (red dashed line), landslide dam (blue circle), and intrusive igneous body (magenta dashed outline) (Cooper et al., 2019;Tinto et al., 2015) described in the main text. Heavy black line marks the 50 m/yr ice velocity contour (Joughin et al., 2010). the lithology of Greenland's interior is difficult to constrain due to thick ice cover, we found that the mapped knickzones do not correlate with any abrupt changes in the magnetic or Bouguer gravity anomalies (Figure 3b) that might be indicative of a major change in magnetic susceptibility and/or bedrock density (i.e., lithology). They also do not appear to be marked by any clear change in channel width or morphology (Figure 3c), and are unrelated to the locations of major interpolated geological terrane boundaries (Henriksen et al., 2009). Although smaller-scale (<10 km) geological structure cannot be resolved with the available spacing of geophysical data, observations from elsewhere on Earth show that long-wavelength (>10 km) knickzones appear to be controlled by regional uplift and rarely correlate with lithological boundaries (Cook et al., 2009;Wapenhans et al., 2021).
However, we note that geological structure does appear to exert an influence over the spatial organisation of the valley network. Between 77°N and 78°N, the trunk channel is abruptly deflected by ∼90° where it encounters a sharp, linear magnetic and Bouguer anomaly "high" (Figure 9). When this ∼200 km-long magnetic/gravity lineament terminates, the channel reverts back to its previous trend. A combined gravity and magnetic anomaly "high" is indicative of a dense, magnetic lithology (e.g., crystalline basement or a mafic igneous body) that may be particularly resistant to erosion and consequently causes the diversion of the channel despite the lack of an associated topographic expression. Similarly, the channel is deflected around a concentric ∼50 km-diameter magnetic and Bouguer anomaly "high" at around 80°N, that corresponds with notably flat and smooth bed topography, and has previously been interpreted as the signature of an igneous intrusion (Tinto et al., 2015). These observations imply that the channels have established an equilibrated relationship with the regional geological structure, which is consistent with the long-term development of a fluvial network.

Channel Modification by Glacial Erosion and Meltwater Action
Upstream of the Petermann Glacier, almost the entire trunk channel profile is characterized by low ice velocities (<20 m/yr; Figure 3a), a "likely frozen" basal thermal regime, and negligible basal melt rates (Joughin et al., 2010;MacGregor et al., 2016), implying that substantial channel modification beneath the contemporary GrIS is unlikely. However, it is possible that the valley network was modified by glacial processes in the past, for example, when the GrIS was more restricted and/or oscillating across the land surface (Christ et al., 2021;Keisling et al., 2020b;Schaefer et al., 2016).
We have demonstrated that the stream power model, which simulates fluvial incision, is able to replicate several features of the observed longitudinal profile, with the exception of the overdeepenings at 720 and 1,130 profile-km (Figure 3e). Overdeepenings or undulating profiles are a hallmark of erosion by glacial ice and/or where pressurized flow under ice can drive subglacial meltwater up-slope (MacGregor et al., 2000;Sugden & John, 1976). As well as high channel depths of up to 400 m (Figure 3d), the lower overdeepening (at 720 profile-km) is characterized by anomalously narrow valley widths (Figure 3c). Since glacial erosion itself cannot decrease the width of a valley, we propose that the narrow widths in this region are an inherited feature of the channel. A narrowed section of a valley is indicative of a more resistant bedrock lithology, which results in lower lateral erosion rates (Finnegan et al., 2005;Schanz & Montgomery, 2016). This region is also characterized by a distinct magnetic anomaly "high" and elevated valley rim terrain (Figures 9 and 10), which lend further support to a lithological control on valley width.
The overdeepening itself likely formed at a later stage, when ice and/or subglacial meltwater was flowing through the existing paleo-fluvial channel system. On encountering a narrowed conduit with a reduced cross-sectional area, convergent ice/meltwater would be required to erode vertically downwards in order to conserve the volume flux, resulting in the development of an overdeepening in the profile (MacGregor et al., 2000). This interpretation is consistent with the approximate halving of the channel width from ∼6 km to ∼3 km and approximate doubling in incision depth from ∼200 to ∼400 m through this region (Figures 3 and 10), with the channel cross-sectional area remaining broadly constant. The overdeepening is also located immediately downstream of a confluence between the trunk channel and a well-defined tributary (Figure 10f), which would facilitate the convergence of ice and/or meltwater.
The relative roles of erosion via ice versus meltwater in overdeepening development in this setting are unclear. We note that the channel cross profile is more resemblant of a V-shape associated with flowing water than a U-shape associated with glacial ice (Figure 10b). The appearance of the overdeepened channel resembles features documented in formerly glaciated regions that are interpreted to have been cut by subglacial meltwater, such as inner gorges and tunnel valleys (Beaud et al., 2016;Jansen et al., 2014), although the valley has a lower widthto-depth ratio than is typically observed in subglacial meltwater channels (Grau Galofre & Jellinek, 2017). Alternatively, the large depths of incision and the V-shaped geometry may be attributed to rapid erosion by meltwater sourced from glacial outburst flood events (Keisling et al., 2020a), although it is unclear if such events are able to overdeepen the valley profile by hundreds of meters.
Gravity anomalies suggest that the dam at 370 profile-km is comprised of low-density sedimentary material (Figure 4) that we speculate may have been deposited by bedrock landsliding, which is common in steep-sided canyons and may be partly responsible for the observed large valley widths and channel asymmetry in this part of the profile (see Figures 3c and 4a-4c). Landsliding in canyons typically results in river diversion and incision to  (Morlighem et al., 2017). Contour interval is 100 m. Blue dashed lines mark the edges of the region of elevated terrain where the channel is anomalously deep and narrow and is directly superimposed on a magnetic anomaly "high" (see panel b). Magenta circle marks the confluence between the main channel and a tributary upstream of the overdeepening. bypass the blockage (Crow et al., 2014), or the formation of a lake upstream of the dam (Fan et al., 2020). Since the channel in northern Greenland is not re-routed around the dam, either (a) a significant lake existed upstream of the dam in the past, or (b) the dam formed after (and/or as a consequence of) GrIS expansion, meaning there was no induced fluvial response. Constraining the timing of dam formation and channel blockage would therefore assist our understanding of fluvial and/or meltwater drainage pathways during past intervals in which the region was fully or partially ice-free.
To summarize, while the morphology of the channel network overall appears to be consistent with initial development via fluvial erosion (Section 6.1), there is evidence that parts of the network have been modified and overprinted by erosion driven by a combination of ice, subglacial meltwater, and/or outburst floods during the Quaternary (Keisling et al., 2020a;Livingstone et al., 2017). Although untangling the roles of glacial and meltwater erosion is beyond the scope of this study, closer examination of valley morphology in the future may provide further constraints on the routing of past glacial and meltwater networks in Greenland, and in turn the dynamics and extent of past ice sheets.

Timing and Cause of Base Level Fall
Our longitudinal profile modeling indicates that the inherited (pre-glacial) fluvial network experienced two episodes of base level fall, incision, and knickzone retreat; here we discuss the possible age and causes of these incision events.
A number of lines of terrestrial and marine evidence suggest that Greenland witnessed an expansion from early ephemeral mountain ice caps to a more contiguous ice sheet during the Pliocene, before briefly retreating during the mid-Pliocene warm period (3.29-2.97 Ma), and subsequently growing again, culminating in a fully-expanded GrIS by ca. 2.7 Ma (Bierman et al., 2016;Knutz et al., 2019;Solgaard et al., 2013;Tripati & Darby, 2018). Our modeling suggests that the retreat of the upper and lower knickzones to their present locations took 3.7-9.6 Myr and 2.8-5.5 Myr from the onset of the first and second episodes of base level, respectively. Since it is unclear whether fluvial incision had largely ceased by the time of glacial expansion at ca. 2.7 Ma, or continued during icefree intervals during the Quaternary, we extend the upper age limit for the two base level fall events by 2.7 Myr (Figure 7g). This gives inferred geological age ranges for the onset of the two base level fall events of 12-3.7 Ma and 8.2-2.8 Ma.
Base level fall during the Neogene is most likely to have been caused by (a) a fall in local sea level, and/or (b) uplift of the land surface. Global Cenozoic records indicate a maximum eustatic sea level fall of ∼50 m during the mid-to-late Miocene (ca. 15-8 Ma), which is primarily attributed to the expansion of ice in Antarctica (Miller et al., 2020). Although the timing of this sea level fall overlaps with the inferred age of the upper knickzone, the amplitude is much lower than the base level change recorded by the longitudinal profile (Figure 7). This implies that while eustatic sea level fall may have contributed in part, uplift of the land surface is a more plausible cause of the majority of the recorded ∼530 m of base level fall. We also note here that the ice-free longitudinal profile appears to reach a base level that approximates present-day sea level, albeit with a degree of uncertainty associated with the correction for the isostatic response to fjord incision (Figure 3e). We suggest that if the uncertainty in the pre-glacial elevation of the valley lower course could be reduced, this may provide an important constraint on Neogene paleo-sea level in northern Greenland.
Apatite fission track data, thermal history modeling, and geological observations including high-elevation planation surfaces and erosional unconformities have been used to infer two episodes of near-surface cooling in east and west Greenland at 11-10 Ma and 7-2 Ma (Bonow et al., 2006a(Bonow et al., , 2006bJapsen et al., 2014). These cooling events are assumed to have been induced by Neogene uplift, which triggered denudation and valley incision Japsen et al., 2014). The timing of our two inferred base level fall events overlaps with the ages of these two reported cooling events. The inferred magnitudes of uplift in east and west Greenland are ∼1 km at 11-10 Ma and 1-2 km at 7-2 Ma Japsen et al., 2014), compared to the ∼150 and ∼380 m amplitudes of the knickzones in northern Greenland. This indicates that if these events were contemporaneous, the amplitude of uplift was considerably smaller in northern Greenland than in east and west Greenland.
Additionally, the inferred age of the first phase of base level fall overlaps with that of a prominent seismic angular unconformity on the northeast and northwest Greenland shelves with an estimated age of middle to late Miocene (ca. 15-10 Ma) (Berger & Jokat, 2009;Døssing et al., 2016;Japsen et al., 2021). The Miocene-Pliocene sediment units above this erosional unconformity are characterized by strongly prograding clinoforms, which are indicative of relative sea level fall and/or increase in sediment input (Berger & Jokat, 2009). This observed stratigraphy is consistent with our inferences of base level fall and enhanced fluvial incision onshore at this time.
We note that there is an ongoing debate as to whether the high-elevation peneplains in east and west Greenland should be interpreted as markers of Neogene uplift events , or if they formed in situ as an emergent feature of glacial erosion (Egholm et al., 2017). It has also been suggested that the highlands are instead remnants of ancient Caledonian topography, and have experienced protracted downwearing since orogenesis and no major Cenozoic cooling (uplift) (Jess et al., 2018;Pedersen et al., 2012). The modern rejuvenation of these highlands may instead be the result of glacial erosion and associated flexural isostatic uplift (Medvedev et al., 2013). However, this is unlikely to have been a primary driver of knickzone development in northern Greenland, given that fjord formation likely only commenced in the early Pleistocene (ca. 2.5 Ma) (Pedersen et al., 2019) and the calculated erosion-driven uplift near the downstream end of the valley profile is < 150 m ( Figure S5 in Supporting Information S1), compared to the total uplift of ∼530 m implied by the knickzones.
The interpretation derived from the elevated planation surfaces and thermal history modeling is that the landscape of east/west Greenland was graded to sea level during the Paleogene-Neogene, and then experienced two phases of kilometer-scale uplift and valley incision into the peneplains Japsen et al., 2014). In our modeling we assumed an initially steady-state landscape in northern Greenland, that subsequently experienced two phases of base level fall, resulting in two migrating knickzones without peneplain formation. Although both approaches infer two phases of uplift/base level fall during the Neogene, a degree of uncertainty persists as to the nature of the paleotopography prior to base level fall and incision, as well as the cause(s) of the uplift and its spatial variability.
One possibility is that uplift was induced by tectonic activity. The timing of base level fall appears to correlate with changes in global plate motion, as well as local bathymetric changes such as the widening of, and onset of seafloor spreading in, the Fram Strait in the late Miocene (Engen et al., 2008), implying that uplift was triggered by plate tectonic forces, flank uplift, and/or changes in intraplate stress (Døssing et al., 2016;. If the fluvial incisional response induced in northern Greenland was contemporaneous with uplift in east and west Greenland, the implication is that the uplift was widespread. However, specifically how mechanisms such as changes in plate motion cause consistent and widespread vertical motion far from plate margins remains unclear. Alternatively, Iceland Plume activity has been invoked as a possible cause of Pliocene uplift in east Greenland Steinberger et al., 2015), and may also be wholly or partially responsible for the broader pattern of Neogene uplift across Greenland. Surface displacement may have resulted from pulses of anomalously warm, buoyant mantle plume material that are inferred to have occurred during the Cenozoic (Parnell-Turner et al., 2014). The comparatively smaller magnitudes of uplift inferred in northern Greenland relative to east and west Greenland may imply a spatial gradient in uplift radiating northwards and westwards from Iceland/east Greenland.
Another process that could cause uplift of the land surface in this setting is the isostatic response to ice sheet (un) loading. However, it is unclear whether there is a plausible ice load history that could have induced one or more widespread uplift events on the order of 100s of meters-and also enabled a fluvial incisional response to the uplift-over the timescales indicated by our longitudinal profile modeling. Although beyond the scope of this study, we suggest that in the future a combined ice sheet -glacial isostatic adjustment -landscape evolution model may be able to test this possibility.
Although the causes and distribution of Neogene uplift remain to be fully constrained, we have demonstrated the value of analysis of subglacial valley networks extracted from the original RES data rather than purely from continental DEMs. As a final point, we reiterate that subglacial valley networks have been described in other regions of Greenland (Figure 1b). Although these channels may have been subjected to more substantial modification by glacial erosion than in northern Greenland (Cooper et al., 2016;Livingstone et al., 2017), they may retain informative signatures from pre-and syn-glacial landscape development. Further analysis of these other families of channels may be an important future target in developing a fuller understanding of the distribution of fluvial and glacial incision, landscape antiquity, and patterns of base level change and uplift that may have been facilitators of glacial inception across Greenland.

Conclusions
In this study, we used airborne geophysical data to analyze subglacial valley morphology in northern Greenland, and developed a simple numerical model of longitudinal profile evolution, which we used to infer base level history. We draw the following conclusions: 1. Mapped valley morphologies are consistent with a primarily fluvial origin of the network, with evidence for localized valley modification resulting from the erosive action of glacial ice, subglacial meltwater, and/or outburst floods. 2. The spatial organisation of the channel network and locations of profile overdeepenings are influenced by the regional geological structure, implying a long-lived and well-established hydrological system and highlighting the value of valley network analysis in constraining subglacial geology. 3. Longitudinal profile modeling indicates that the channel was likely active over multi-million-year timescales, with the locations of knickzones consistent with two episodes of base level fall: ∼150 m at ca. 12-3.7 Ma and ∼380 m at ca. 8.2-2.8 Ma. This interpretation is consistent with the elevations of terraces in the walls of lower course of the channel, which we suggest may represent relict valley floors from a time prior to base level fall, incision, and knickzone retreat. 4. The timing of the inferred base level fall correlates with other onshore and offshore records of Neogene uplift, denudation, and/or relative sea level change. Although the nature of the early Neogene topography and the mechanism(s) responsible for uplift remain unclear, tectonic and/or mantle-driven uplift and the resulting fluvial incisional response appear to have played an important role in the genesis of the modern landscape of northern Greenland.

Data Availability Statement
Data related to this article can be found at https://doi.org/10.5281/zenodo.5648042. The data set includes the northern Greenland subglacial valley morphology database generated from analysis of the radio-echo sounding data, and an ESRI shapefile of the interpreted subglacial channel network. The geophysical datasets used in this study are archived and freely available at the National Snow and Ice Data Center (NSIDC). Operation IceBridge: https://doi.org/10.5067/GDQ0CUCVTE2Q; BedMachine: https://doi.org/10.5067/VLJ5YXKCNGXO; MODIS Mosaic of Greenland: https://doi.org/10.5067/9ZO79PHOTYE5. Operation IceBridge radar, gravity, and magnetic profiles displayed in this manuscript are labeled with the relevant flight numbers. . Drainage basin calculation and analysis were performed using TopoToolbox software (https://topotoolbox.wordpress.com/) (Schwanghart & Scherler, 2014). Figures were prepared using the Generic Mapping Tools (GMT) version 5 software package (Wessel et al., 2013). We would like to thank Peter Japsen, a second anonymous reviewer, and the editor Mikaël Attal for their helpful and constructive reviews, which greatly improved the final manuscript.