Preserved landscapes underneath the Antarctic Ice Sheet reveal the geomorphological history of Jutulstraumen Basin

The landscape of Antarctica, hidden beneath kilometre‐thick ice in most places, has been shaped by the interactions between tectonic and erosional processes. The flow dynamics of the thick ice cover deepened pre‐formed topographic depressions by glacial erosion, but also preserved the subglacial landscapes in regions with moderate to slow ice flow. Mapping the spatial variability of these structures provides the basis for reconstruction of the evolution of subglacial morphology. This study focuses on the Jutulstraumen Glacier drainage system in Dronning Maud Land, East Antarctica. The Jutulstraumen Glacier reaches the ocean via the Jutulstraumen Graben, which is the only significant passage for draining the East Antarctic Ice Sheet through the western part of the Dronning Maud Land mountain chain. We acquired new bed topography data during an airborne radar campaign in the region upstream of the Jutulstraumen Graben to characterise the source area of the glacier. The new data show a deep relief to be generally under‐represented in available bed topography compilations. Our analysis of the bed topography, valley characteristics and bed roughness leads to the conclusion that much more of the alpine landscape that would have formed prior to the Antarctic Ice Sheet is preserved than previously anticipated. We identify an active and deeply eroded U‐shaped valley network next to largely preserved passive fluvial and glacial modified landscapes. Based on the landscape classification, we reconstruct the temporal sequence by which ice flow modified the topography since the beginning of the glaciation of Antarctica.


| INTRODUCTION
The geomorphology of the subglacial landscape beneath the Antarctic Ice Sheet (AIS) is an essential part to our understanding of its future stability (Sugden & Jamieson, 2018). The bed topography has been modified by extensive glacial and subglacial erosion on continental, regional and local scales (Jamieson et al., 2010). However, large ice sheets do not only erode the bed, but also preserve the subglacial environment in certain regions at the same time, depending on the flow velocity of ice, and on whether or not the ice sheet is frozen to the ground or sliding (Jamieson et al., 2014). In order to reconstruct the landscape evolution of the AIS, we thus need to understand present and former ice-stream activity and ice-sheet configurations (Siegert et al., 2005) since the glaciation of Antarctica. The reconstruction of the pre-glacial landscape and the determination of regions where the traces of former ice-flow activity are preserved helps to discriminate where glaciers have built up in the past (Sugden & Jamieson, 2018). To do this, we can refer to an increasing amount of radar data to map the ice-sheet base and examine the large-scale topography (km scale; e.g., Cui et al., 2020;Fretwell et al., 2013;Morlighem et al., 2020) and individual landforms (King et al., 2009).
Depending on the nature and scale of the structures investigated, this needs to be done at high spatial resolution in some locations (tens to hundreds of metres; e.g., Cui et al., 2020;Holschuh et al., 2020;King et al., 2009).
The pre-glacial topography steers ice flow and determines the distribution of selective erosion (Jamieson et al., 2014). Hence, to study the behaviour of the AIS, we depend on an understanding of the relationship between pre-glacial landscapes and the processes by which ice flow can modify them (Jamieson et al., 2014). It is necessary and important to understand this feedback mechanism if we are to improve our knowledge on the past and future ice-sheet behaviour (Jamieson et al., 2010;Oerlemans, 2010). In particular, constraints on ice-sheet and landscape development are key boundary conditions for numerical simulations of geomorphological development of the landscape beneath the AIS (Jamieson et al., 2010).
In this study, we use existing and new high-resolution airborne radar data to improve our understanding of the bed topography and landscape evolution in the source region of the Jutulstraumen Glacier (JG) in western Dronning Maud Land (DML; see Figure 1). JG is the largest and most important outlet glacier in DML (Callens et al., 2014;Høydal, 1996) and feeds the second largest ice shelf (Fimbul Ice Shelf) in DML (Rignot et al., 2013). We present various geomorphological structures at the glacier bed, and interpret them in terms of the regional erosional history. Our results give insight into the effects of past and present ice dynamics on the glacier bed and into landscape development in glaciated areas in general.

| Regional setting
The present thick ice-sheet cover in DML inhibits a detailed understanding of its geological history since Precambrian times. Only the combination of geological sampling of the sparse outcrops/nuntaks (Jacobs et al., 1998(Jacobs et al., , 2017 combined with aerogeophysical data (Mieth & Jokat, 2014) made it possible to retrieve some first-order information of the regional geology ( Figure 2).
These investigations indicate numerous geological affinities between DML and South Africa. Both fragments re-collided between 650 and 500 million years ago (Jacobs & Thomas, 2004) along the East Africa-Antarctic Orogen to form a new southern supercontinent named Gondwana. Massive volcanism started around 180 million years ago both in East Africa and in DML (Elliot, 1992), which finally led to the fragmentation of Gondwana (Jokat et al., 2003;Leinweber & Jokat, 2012;Mueller & Jokat, 2019).
The Jutulstraumen Graben represents a deep topographic depression which cuts through the otherwise continuous DML mountains extending for 1500 km from 15 W to 30 E (Jacobs, 1991;Jacobs et al., 1992). This mountain chain includes the Borgmassivet and Sverdrupfjella, which form the flanks of the Jutulstraumen Graben, and the Kirwanveggen separating the western Jutul-Penck Graben from the main Jutulstraumen Graben (Riedel et al., 2012;Steinhage et al., 1999) in the east (Figure 1c). Since the graben system runs along several geological units of different ages it is still debated whether the graben formed along a geological weakness of a failed rift formed prior to the Gondwana break-up (Ferraccioli et al., 2005) or has a different origin. JG and its tributary, the Penck Graben Ice Stream, mainly follow the tectonic lineament ( Figure 1b,c), which provides a northward-oriented draining system for the East Antarctic Ice Sheet.
South of the exposed mountain range, DML hosts a complex subglacial landscape. Näslund (2001) proposes a step-by-step genesis, starting with the development of relief by subaerial weathering and erosion by mountain glaciers and river systems resulting in an alpine landscape. Afterwards, erosion of the mountain range by large ice streams is followed by the growth of the East Antarctic Ice Sheet in mid-Miocene times (Holbourn et al., 2005;Shevenell et al., 2004).
Further east, between the ice streams that flow around the Sør Rondane Mountains, the plateau inland of the escarpment has been interpreted to preserve features of a pre-existing fluvial landscape (Eagles et al., 2018). Näslund (2001) presented helicopter-based RES data collected along a 260 km profile across the Jutulstraumen Glacier and from a small area around Jutulsessen Nunatak with about 70 km of data. Näslund (2001) used these prevalent data to reconstruct the long-term landscape development in western and central DML. Around the same time, Steinhage et al. (1999) and Steinhage (2001) gathered more than 50 000 km of airborne RES data in DML as part of the pre-site survey for the European Project for Ice Coring in Antarctica (EPICA). Additional aerogeophysical datasets were acquired in the early 2000s in the Jutulstraumen region (Ferraccioli et al., 2005) and in western DML (Mieth & Jokat, 2014;Riedel et al., 2012). These data have supported detailed interpretations of the geological history of central and western DML, compilations of continental-scale bed topography and subglacial-lake distributions (Huybrechts et al., 2000;Fretwell et al., 2013;Goeller et al., 2016;Morlighem et al., 2020), and targeted studies of bathymetry beneath DML's ice shelves (Eisermann et al., 2020;Smith et al., 2020). To date, however, they have not been used for a regional-scale interpretation of subglacial landscape-forming processes that would complement and update the early interpretations of Näslund (2001). All JURAS-2018 radar profiles were recorded in a frequency range of 180-210 MHz and at a flight altitude equivalent to 365 m above ground. For radar data processing, we used the CReSIS Toolbox (CReSIS, 2020). The processing comprises algorithms for pulse compression and Synthetic Aperture Radar (SAR) processing via F-K migration and array processing. Our SAR-processed radargrams have a trace interval of 15 m and a range resolution of 4.3 m. Franke, Jansen, Binder, et al. (2020) give a detailed description of the radar acquisition and processing. The error in calculating the ice thickness is 24.7 m estimated from the RMS of the sum of the crossover mean at 62 intersections. This estimate takes also a 1% error of the dielectic permittivity into account (for further details see Franke, Jansen, Binder, et al., 2020).

| Ice thickness and radar data
The UWB radar data were acquired at the onset of streaming flow, as defined by ice surface velocities of 5-100 m a À1 , between the Troll (Norway) and Kohnen (Germany) stations ( Figure 3). In the south, the survey covers an area within which ice flow of JG accelerates. The southern survey lines are oriented perpendicular to the direction of ice flow. The northern survey lines are oriented at $60 to the southern survey lines. They cover parts of the eastern valleys leading to the Jutulstraumen Graben, which also corresponds to an area of accelerating surface velocities and parts of the higher elevated topography to the east of it.
Complementary ice thickness radar data acquired with the AWI EMR (Electromagnetic Reflection System) radar  and also depicted in Figure 3 were gathered over the last three decades throughout DML (see the section about previous radar surveys and Steinhage et al., 1999, Steinhage, 2001, and Riedel et al., 2012. The profiles of the JURAS-2018 survey were designed to map the ice-bed interface at high resolution while also filling substantial gaps in the existing RES data. Most of these data gaps occur over deep topographic depressions in which the bed could not be detected by older radar systems. F I G U R E 1 Overview of the geographical and glaciological setting of the JG drainage basin. (a) LIMA Landsat image mosaic (Bindschadler et al., 2008), ice-surface velocity (Mouginot et al., 2019) and BedMachine Antarctica (BMA) bed topography (Morlighem et al., 2020). Panels (b) and (c) zoom in on the ice-flow velocity and bed topography maps for the study region. White arrows schematically indicate the convergent ice flow upstream of the Jutulstraumen Graben. The black outline represents the JG drainage system boundary after Mouginot and Rignot (2017). It has to be noted that in our study we focus on the region south of the Maud Belt [Color figure can be viewed at wileyonlinelibrary.com]

| UWB radar ice thickness and bed topography
We determine ice thickness by converting the two-way travel time (TWT) of the radar wave between the ice surface and bed reflection to depth. To do this, we use the depth-dependent electromagnetic wave speed in firn and ice (Winter et al., 2017) with a constant electromagnetic wave speed of 1.689 Â 10 8 m s À1 (corresponding to a relative dielectric permittivity of 3.15). The electromagnetic wave speed is the same as used for synthetic aperture radar (SAR) processing (see Franke, Jansen, Binder, et al., 2020). In addition, we use a firn correction factor of 13 m (Steinhage, 2001). The bed topography is calculated by subtracting the ice thickness from an ice surface digital elevation model (REMA DEM; Howat et al., 2019). For the creation of an improved bed topography, we use the SAGA GIS multi-level B-spline interpolation module (Conrad et al., 2015) where line coverage is dense and ordinary kriging where line coverage is sparse. The grid size of our improved bed topography is 1 km.

| Valley morphology
Knowledge of the morphology, distribution and spatial orientation of subglacial landforms is important for investigating the ice-sheet dynamic history of glaciated and previously glaciated landscapes. For this study, we analyse the geometry of single valleys along the JURAS-2018 radar transects to determine the developmental stage for a glacial valley (Hirano & Aniya, 1988). We derive three geometrical parameters to categorise the valley type (see Next to the trivial parameter D, we use the parameters R and the V-index to obtain information about the geometry of the valley. The valley ratio factor R informs whether the valley is wide or narrow, whereas the V-index helps to discriminate between U-shaped and Vshaped valleys. Here, we bear in mind that the V-index is a relative measure; that is, negative values (V-index < 0) represent a tight V-   Wang et al. (2020)). The position of DML, East Antarctica is shown in part (a) relative to Africa and India. DML terrane boundaries in part (b) are interpreted from scattered outcrops and interpolated on the basis of regional aeromagnetic data (Jacobs et al., 2017;Mieth & Jokat, 2014).

| Basal roughness
The analysis of valley geometries is somewhat subjective because it relies on a starting interpretation to select the set of features to analyse as valleys. To complement these results, we also analyse basal roughness to quantify the bed topography in a systematic way. Topographic roughness is an indicator of subglacial conditions during present and past ice-flow activity and a control on current ice-sheet dynamics Rippin et al., 2014). Here, we consider a spectral method (Cooper et al., 2019;Li et al., 2010) in combination with the root-mean-square (RMS) height (Shepard et al., 2001) to quantify subglacial roughness.
Spectrally derived roughness is the relative vertical and horizontal variation of the ice-bed interface derived from the nadir bed F I G U R E 3 Flight tracks with new AWI UWB JURAS-2018 (red) airborne radar data in the Jutulstraumen drainage system (black outline Mouginot & Rignot, 2017). Previous airborne AWI EMR tracks (Steinhage, 2001;Steinhage et al., 1999) are highlighted in blue. Icesurface velocities (colour code), after Mouginot et al. (2019), are superimposed over the JG catchment. The EPSG projection code for this and all following maps is epsg:3031 [Color figure can be viewed at wileyonlinelibrary.com] reflection. We distinguish between two parameters: ξ, which reflects the vertical irregularity of the bed and provides information about the dominant vertical amplitudes; and η, which describes the dominant horizontal wavelength. For the calculation of these two parameters, we follow Franke et al. (2021). For this study, we use a moving window (MW) length of 512 samples, which corresponds to an approximate ground distance of 7.68 km. We use this value to balance between the detection of large valleys in the roughness analysis as a dominant feature (when the MW is too large) and the fact that valleys may not be detected as roughness when the MW is too small. The majority of radar lines are oriented parallel to each other and approximately perpendicular to ice flow. This alignment promotes a more robust comparison of roughness values than would be possible with a random survey layout.
In addition, we make use of the RMS height. The RMS height χ (or the standard deviation of heights about the mean) is defined by where n represents the number of sample points, z(x i ) the height of the bed at position x i , and z the mean height of the detrended profile.
Here, we perform a bed roughness analysis using the JURAS-2018 radar data. The combination of spectral vertical and horizontal roughness estimates enables us to infer the general appearance of the structure at the bed surface, while the RMS height provides information about the absolute elevation differences. The three methods require a uniform sample spacing, which is 15 m in the JURAS-2018 radar data.

| Isostatic adjustment
We follow the simple approach of Rose et al. (2015) to simulate the isostatic rebound in order to get an idea of the ice-free paleotopography. Here, we only simulate the bed elevation prior to the start of Antarctic glaciation (34 million years ago). The Airy-Heiskanen model (Airy, 1855) predicts that different topographic heights are fully compensated by variations in crustal thickness. We use this model to predict the elevation change the overburden ice will cause on the bed topography. For the calculation, we use densities of 915 kg m À3 for ice (ρ i ) 2750 kg m À3 for the crust (ρ c ) and 3330 kg m À3 for the mantle (ρ m ). For a single location or grid cell, we can calculate the uplift response h r using This approach does not account for the entire complexity of processes involved in isostatic rebound after overburden ice removal (Rose et al., 2015). Furthermore, it is possible that the isostatically rebounded topography is underestimated due to erosion and removal of sediment during the last 34 million years, and the original surface may have formed at a higher elevation (Sugden et al., 1995). Nonetheless, we obtain a rough idea about the paleo-topography and are able to use simple statistical methods to infer the characteristics of the landscape.

| Hypsometry
An isostatically corrected bed topography allows us to analyse the frequency distribution of elevations (hypsometry) to characterise the landscape morphology (Brozovi c et al., 1997). The aim of this method is to generate indications of whether or not the landscape has a fluvially or a glacially imprinted character. A widespread approach to hypsometry is to create a simple histogram of the frequencies in different elevation bins from a topographic DEM.
In this context, Brocklehurst and Whipple (2004) argue that hypsometry of large arbitrary regions potentially masks the detail in topographic variation in comparison to smaller regions. In addition, hypsometric analyses are sensitive to the resolution of the DEMs on which they are based. The method is most often applied for ice-free alpine regions where the topography is resolved with a grid cell resolution of 30-50 m (Brocklehurst & Whipple, 2004). By contrast, continental-scale bed topography DEMs of ice-covered areas in Antarctica have a resolution of 500-1000 m (e.g., Fretwell et al., 2013;Morlighem et al., 2020) and are interpolated over large distances (up to tens to hundreds of kilometres).
To overcome some of these drawbacks, we follow Creyts et al. (2014) and apply hypsometry on the radar-derived ice thickness data of the JURAS-2018 survey after applying an isostatic correction (Equation 2). For this study, we use the results of a hypsometric analysis to compare different landscapes and to (i) infer the degree of glaciation (Brocklehurst & Whipple, 2004), (ii) discriminate between fluvial and glacial landscapes (Sternai et al., 2011) and (iii) discuss modes of glacial landscape evolution (Jamieson et al., 2014;Li et al., 2010). To avoid a particular sub-region from dominating the signal in the hypsometry, we analyse the entire dataset as well as sub-regions that reflect a particular geomorphological setting.

| Water flow routing
We apply a simple water flux scheme to simulate the general water flow network, assuming a completely ice-free isostatically compensated topography. By doing so, we aim to approximate features of the paleo-fluvial system in our survey region. We make use of the following SAGA GIS algorithms (Conrad et al., 2015): (i) the fill sinks module, to fill surface depressions with a minimum angle of 0.1 and (ii) the catchment area algorithm, to calculate the flow accumulation, which represents the number of accumulated cells. In this study, we only perform the water routing for the drainage basin of JG.

| Radar data and bed topography
The additional ice thickness data of the JURAS-2018 survey allowed us to create an improved bed topography DEM of our survey region ( Figure 5). It was possible to fill significant data gaps where the older radar system did not show the bed reflection. However, some data gaps imaging the bedrock topography are still present in the new radar data. The distribution of the gaps is shown in Figure 5b. The majority is located over deep topographic depressions at the onset of the Jutulstraumen Graben. In the upstream region, data gaps occur in regions with comparable or thinner ice cover than in downstream regions with continuous bed reflections ( Figure 6). This indicates that higher englacial attenuation is causing the signal loss or that the bed reflection itself is weaker. Figure 6 shows in Figure 5b and radargrams in Figure 6). Despite these difficulties, it is still possible to estimate these valleys' general topographies on the basis of the shape of the adjacent bed and internal layering. We use this information in conjunction with the new ice thickness data to create an outline of the valley network associated with the Jutulstraumen Graben (see Figure 5c). A comparison between existing bed topography and new AWI UWB ice thickness data with a subsequent interpretation of the aforementioned data gaps enables the following observations: 1  Here, the ice-flow velocity is slightly lower than in Zone a, while bed topography is more elevated than in the west. In contrast, the horizontal roughness is similar between Zones a and b. Shorter wavelengths are more dominant in the eastern part of the northern radar set than in the western part.

| Hypsometry
Analysis of the elevations reveals that, in addition to their contrasting basal roughnesses, Zones a and bare hypsometrically distinct in Morlighem et al. (2020), where the authors compare the Bedmap2 and BedMachine Antarctica v01 bed topography).
U-shaped valleys are distinctive features that suggest selective linear erosion on a local to regional scale, and are thus an indication for the presence of warm-based ice at the time of erosion (Sugden & John, 1976). The deepest topographic depressions in the JG drainage area are located where present-day ice-flow velocity and acceleration are highest (Mouginot et al., 2019). This suggests that the region (black outline in Figure 7d) has experienced erosion over a long period. The wide but shallow U-shaped valleys in the upstream region (red outline in Figure 7d) are located close to the ice divide, where the ice column is thickest and current ice-flow velocities are low. This indicates that these valleys were probably not eroded under the present ice-sheet configuration. The locations, shapes and sizes of these valleys suggest they most likely formed in regions of tectonic or lithologic weakness, or by the exploitation of valleys in a pre-existing fluvial network (e.g., Baroni et al., 2005;Jordan et al., 2013;Rose et al., 2015).
As an ice sheet grows to an intermediate size over a pre-existing fluvial landscape, selective linear erosion will deepen and widen existing river valleys (Sugden, 1968(Sugden, , 1978. Ongoing ice-sheet growth may lead to the development of a cold base, concomitant with a reduction of sliding and increase in ice-bed coupling (Jamieson et al., 2010). As a consequence, ice flow will be dominated by internal deformation and ice-flow velocity will reduce (Creyts et al., 2014). At the base of the ice sheet, the consequence of all this may be a change from an erosive alpine environment to a passive setting, in which the bed undergoes little further modification. This kind of evolution may have occurred rapidly over the upstream part of the northern radar dataset (blue outline in Figure 7d), which reflects the lowest degree of glacial erosion in its average high elevation and the dominance of shallow, V-shaped valleys.
The relative orientation of the radar line to the valley pathways is the largest source of uncertainty in the valley geometry analysis.
Oblique crossings will make valleys appear artefactually too wide (a potential candidate could be the rightmost valley in Figure 7 (V 2 )).
Artefacts like this will affect the determinations of valley ratio (R) and  and Ross et al. (2014) at the West Antarctic Ice Sheet.

| Basal roughness
Most importantly, the analysis of spectral vertical roughness and RMS height supports a certain aspect of the valley geometry analysis, increasing the robustness of their interpretation. It has to be noted that the focus of our roughness analysis and the choice of our MW is that we are able to discriminate the different landscapes. The spatial clusters of deep valleys correlate with high vertical roughness and RMS height, whereas the locations of shallow valleys correlate with low vertical roughness and RMS height. Furthermore, we observe that in Zone b the combination of narrow, V-shaped valleys is consistent with the dominance of short horizontal wavelengths. We are aware that both roughness values and valley geometries depend on the orientation of the radar profiles. However, since both analyses are based on the same radar profiles, we can assume that they experience the same effects due to the orientation of the profiles.

| Hypsometry
The main differences in the hypsometric curves for Zone a and

| Landscape classification
Using the isostatically corrected bed topography (Figure 5d) as an approximation of the pre-glacial surface, we note that the deepest and widest U-shaped valleys extend down to present sea level ( Figure 5d). The Jutulstraumen Graben is of tectonic origin, dating back to the ice-free Jurassic, and so must have been subject to a long period of fluvial erosion after its formation (Näslund, 2001). This suggests that the bed in Zone a would have been subsequently modified by selective linear glacial erosion (Sugden & John, 1976).
The deep and wide U-shaped subglacial valleys in Zone a are likely candidates for pre-glacial fluvial activity and subsequent selective erosion after the onset of glaciation (Sugden & John, 1976). Hence the U-shaped valleys of Zone a are likely to follow the pre-existing fluvial system. It is likely that, regardless of the regional evolution in the extent and thickness of ice and the location of ice divides, some portion of the region's ice has always flowed through the Jutulstraumen Graben towards the northern continental margin. The initial formation of the graben itself and the subsequent long-term erosional activity could explain the deeply incised relief. The graben may have acted as an 'attractor' for higher ice dynamics, which then increased erosion rates.
The analysis of basal roughness and hypsometry for the area east of the main trunk of JG (Zone b in Figure 8) points towards a largely preserved fluvial landscape whose origin precedes extensive glaciation Rose et al., 2015). This is supported by the predominant appearance of shallow and narrow valleys at the glacier bed and the less variable and generally more elevated topography (and the resulting hypsometry). The mix of V-and U-shaped valleys suggests that the bed has been locally modified by glacial erosion at some point in the past. However, the high elevation topography likely resulted in slower flow and cold-based conditions (Sugden & John, 1976) for a large part of its history, resulting in a limited degree of landscape modification. Erosion under cold-based ice is several orders of magnitude slower than under temperate ice (Stroeven et al., 2002;Thomson et al., 2010). Jacobs et al. (1995)

| Basal thermal regime
Geothermal heat flow (GHF) at the base of the Antarctic Ice Sheet is a parameter which is subject to substantial uncertainty because it is difficult to measure in situ. Consequently, geophysical models diverge greatly. Van Liefferinge et al. (2018) compare existing GHF datasets Fox Maule et al., 2005;Martos et al., 2017;Shapiro & Ritzwoller, 2004) to predict the probable distribution of basal ice that has been frozen over the last 1.5 million years. Their modelled area includes the ice divides at the upstream part of the JG drainage basin.
Based on the individual GHF datasets, their modelling reveals it to be likely that the ice base has reached the pressure melting point at least once over the last 1.5 million years. Moreover, it has been shown that measured GHF in eastern DML is locally much higher than predicted by geophysical methods (Talalay et al., 2020), consistent with evidence for temperate bed conditions from the recovery of refrozen subglacial water at the EPICA-DML drill site (Weikusat et al., 2017;Wilhelms et al., 2014

| Landscape erosion and preservation
The bed morphology of the JG drainage system shows patterns of three successive and/or separate stages of landscape modification ( Figure 9). We regionally extend our interpretation of the structures we have classified in our study and define the following interpretation for landscape erosion and preservation: (1) a preserved fluvial landscape shaped by fluvial erosion, minimal glacial erosion and preserved thereafter; (2) a preserved glacial landscape, which was more strongly modified by selective linear erosion at some time in the past compared to (1); and (3) an active glacial landscape, in which the erosion has been long-lived and is still ongoing.
The glacial landscape close to the valley system of and further upstream (black and red outlines in Figure 9) was probably formed by a regional ice-cap glacier system during a period when Antarctica was influenced by a warmer climate than present (Holmlund & Näslund, 1994;Young et al., 2011). A former ice divide may have separated the two regions during this period (yellow dashed line in Figure 9). At some point, as ice coverage increased (Holmlund & Näslund, 1994), erosion terminated in the U-shaped valleys in the south and west (red dashed arrows in Figure 9), but continued further north where ice surface velocities are still fast today (black dashed F I G U R E 9 Summary of our landscape interpretation for the JG drainage system. The active glacial system (black dashed outline) is located in the immediate surroundings of the JG main trunk where the ice flow accelerates to join it. The preserved fluvial system (blue dashed outline) is located to the east and the preserved inactive glacial system (red dashed outline) to the west and south of the JG main trunk. The respective outlines do not represent physically fixed boundaries but instead areas of dominant landscape characteristics. Ice-flow velocity is indicated with contour lines in a 5 m a À1 interval. Hence the ice flow velocity of the large region south and east of the 5 m a À1 contour line is less than 5 m a À1 [Color figure can be viewed at wileyonlinelibrary.com] arrows in Figure 9). The precondition of the initially large-scale graben likely created a positive feedback causing faster ice flow, more meltwater and, thus, enhanced erosion. In addition, at locations where the ice sheet grows beyond a critical size, we expect basal melting and the development of a subglacial water network (Creyts et al., 2014). In our study area, this might explain the increased erosion rates that locally formed exceptionally deep valleys. Erosion mechanisms, such as plucking and abrasion, are mainly controlled by basal sliding, which requires in turn a lubricated bed (Hallet, 1996). By contrast, refreezing of subglacial water in cold-based areas with a thinner ice thickness (Creyts et al., 2014) is probably also a factor in the preservation of the fluvial landscape in the east. Jamieson et al. (2010) use an ice-sheet and erosion model to investigate the subglacial landscape evolution of Antarctica since the Oligocene (34 million years ago). These authors use three different stages for their simulations (see Table 1 in Jamieson et al., 2010) where: (a) ice sheets grow and shrink at similar spatial scales to northern hemisphere Pleistocene ice sheets (34-14 million years before present); (b) cooling shifts the system towards development of a continental polar ice sheet (14-13.6 million years ago); and (c) the continental polar ice sheet stabilises and subsequently experiences only minor changes (13.6 million years ago to present). In order to establish a rough timing for the landscape evolution in the JG drainage system, we compare the erosion rates modelled by Jamieson et al. (2010) with our landscape classification. For the comparison we mainly rely on the results shown in fig. 6 in Jamieson et al. (2010).
The model predicts that high erosion rates at the Jutulstraumen Graben developed at the end of stage (a) and persisted there ever since. We note that, at the beginning of stage (a), high erosion rates were restricted to the high-elevation areas of the DML mountains.
The fluctuating regional ice caps may not have seen the development of a stable ice divide around the JG drainage basin like today's.
Instead, they may have drained via ice streams running both northwards through the Jutulstraumen Graben and towards the southwest.
We find that modelled erosion rates in the lower-elevation regions (e.g., in the Jutulstraumen Graben and further upstream) must have been higher than in the highlands. This is consistent with evidence for reduced erosion in Zone b. We note that the erosion rates in our survey region during stages (b) and (c) correspond to those of today's ice-dynamic setting. That is, constant high erosion rates are to be found only in the Jutulstraumen Graben.

| SUMMARY AND CONCLUSIONS
The landscape of Antarctica has been evolving for tens of millions to 100 million years. Since the colder climate in Antarctica around 14 million years ago fluvial erosion only played a minor role in the landscape modification (Sugden & Jamieson, 2018). The thick ice cover has protected the underlying landscape where the ice is stagnant and cold based. In our work, we demonstrated that previous fluvial and glacial activity can be preserved beneath thick ice masses in western DML, also highlighting the longevity of the local dynamic setting of ice flow.
Their geomorphological characteristics reflect the processes and icedynamic setting under which they were created (Jamieson et al., 2014;Sugden & Jamieson, 2018), which can serve as constraints for modelling past ice dynamics .
We extend the radar data inventory for central DML, and in particular for the immediate surroundings of the JG onset region to characterise the modification stages of the subglacial landscape. We use the existing bed topography data in combination with new highresolution UWB radar data to analyse the morphology of the bed in the JG drainage system. For our analysis, we rely on the distribution of bed elevations (hypsometry), parametrisations of the roughness of the bed, estimates of subglacial valley types, and analyses of their slopes and connectivity for characterisation of a pre-glacial fluvial system.
The glacial landscape in the JG drainage system can be divided into three different regimes: 1. The area of contemporary high ice-flow velocities marking JG, which represents an alpine landscape that has been extensively Graben is presently the only passage for ice to flow towards the icesheet margin. The ice-dynamic setting must have been different to the present one to create the U-shaped valleys, which reach up almost to the ice divide in the southwest of the survey area. The high probability that the ice base of our survey region is at pressure melting point suggests a higher potential for erosion than if frozen to the bed. However, erosion in this area is limited as long as ice-flow velocities are low. The landscape evolution established in this study is in good agreement with modelled erosion rates in Antarctica over the last 34 million years by Jamieson et al. (2010).