Substantial spatial variation in glacial erosion rates in the Dronning Maud Land Mountains, East Antarctica

The coast-parallel Dronning Maud Land (DML) mountains represent a key nucleation site for the protracted glaciation of Antarctica. Their evolution is therefore of special interest for understanding the formation and development of the Antarctic ice sheet. Extensive glacial erosion has clearly altered the landscape over the past 34 Myr. Yet, the total erosion still remains to be properly constrained. Here, we investigate the power of low-temperature thermochronology in quantifying glacial erosion in-situ. Our data document the differential erosion along the DML escarpment, with up to c. 1.5 and 2.4 km of erosion in western and central DML, respectively. Substantial erosion at the escarpment foothills, and limited erosion at high elevations and close to drainage divides, is consistent with an escarpment retreat model. Such differential erosion suggests major alterations of the landscape during 34 Myr of glaciation and should be implemented in future ice sheet models. The topography of the Dronning Maud Land Mountains, Antarctica, has become more pronounced and rugged since preglacial times due to higher glacial erosion at low elevations and lower erosion at high elevations, according to low-temperature thermochronology

A ntarctica remains the last fully glaciated continent. Its extensive ice sheet is both an important driving force of climate and weather patterns, as well as vulnerable to climate changes. The development of the ice sheet and its interactions with the climate have been the focus of intense study during the last few decades. One important factor during the initial formation of the ice sheet is the paleotopography at the Eocene-Oligocene boundary 1,2 . However, reconstructing this topography poses a challenge. On the one hand, the landscape should be expected to have been strongly modified by 34 Myr of glaciation, on the other hand, it is today largely hidden by the ice sheet and thus inaccessible to direct investigations. Some clues to the erosion of the landscape during the long-lasting glaciation come from the offshore record, which document increased sedimentation rates since the Eocene-Oligocene 3-6 . Early models used to reconstruct the initiation of the Antarctic ice shield often use the present-day topography, isostatically compensated to account for the absence of the ice mass in pre-glacial times, as a starting point 1 . Recent and more sophisticated reconstructions of the pre-glacial landscape take shelf sedimentation into account and also consider the effect of e.g., thermal subsidence and horizontal plate movements 7,8 . Converting offshore sedimentation into onshore erosion remains a challenge, though. For DML (Fig. 1, inset), offshore deposits correspond to an onshore average eroded thickness of c. 400-500 m 7,8 . A uniform erosion of 400-500 m over all of DML is highly unlikely though, as the landscape shows clear signs of strongly focused erosion, including the deep incision of Paleozoic paleosurfaces that is mainly ascribed to wet-based glacial erosion of a pre-existing fluvial or tectonic morphology prior to or during the Oligocene 9 . This morphology has later been preserved by cold-based ice sheets 9 . Attempts have been made to determine the precise provenance of the glacially-derived sediments off the coast of DML, though mostly focusing on the Weddell Sea 10,11 . Only limited attention has been given to the Lazarev Sea, although this area is probably the most important sink area for the DML mountains, as the mountain range acts as a barrier for the terrigenous sediments from the East Antarctic interior. Thus, the majority of the sediments along the DML margin in the Lazarev Sea probably is derived from the mountain range itself 6 .
So far, the onshore glacial erosion in DML has not been constrained in-situ in time and space at a regional scale. Quantitative approaches for in situ measurement of erosion include exposure age dating and low-temperature thermochronology. However, exposure age dating is mainly sensitive to the past ten million years 12 and is therefore not ideally suited for estimating the total erosion during long-lasting glaciations such as the 34 Myr long Antarctic glaciation. Low-temperature thermochronology, including apatite fission-track (AFT) and apatite (U-Th)/He dating (AHe), is applicable within a large age range (c. 400-2 Ma). These methods are, however, not particularly sensitive to very low temperatures (<40°C), so only substantial glacial erosion can be  recorded, such as during long-lasting glaciations. AHe data have recently been applied for constraining the exhumation within the Transantarctic Mountains ( Fig. 1) since the Eocene, recording up to c. 8.8 km of exhumation in the deepest parts of the Beardmore Glacier 13 . For the DML mountains, previous low-temperature thermochronology studies have speculated about strongly increased erosion and exhumation associated with the developing ice sheet [14][15][16] , but the Cenozoic exhumation has so far not been systematically constrained here. After a review of all available thermochronological data from large parts of DML, we have chosen a set of samples that are well suited for thermochronological modeling to investigate temporal and spatial changes in cooling rates 15,16 . The samples come from various geomorphic features, such as mountain tops, escarpment front, drainage divides, ice streams, and coastal areas. This allows us to illuminate the long-lasting glacial erosion history of different parts of the DML mountains since the Oligocene. Furthermore, this highlights the power of low-temperature thermochronology to document differential erosion along the length of the DML mountains, including regions of fjord-like incisions versus ice divides, as well as mountain tops versus foothills. This approach leads to a more holistic view of the differential glacial erosion of the DML mountains, eventually leading to a much better understanding of the paleotopography of the DML mountains at the beginning of the Oligocene glaciation. Our data document a negative correlation between elevation and erosion rate, indicating strong erosion in the foothills of the escarpment and a more-or-less stationary highland. Up to c. 1.5 km and c. 2.4 km of erosion is recorded in western and central DML, respectively, corresponding to time-averaged erosion rates of 45 and 75 m Myr −1 . The dataset clearly documents differential erosion along the DML escarpment, which is mainly attributed to variations in elevation, geomorphic setting, and the drainage area. Detailed thermal history modeling of low-temperature thermochronological data has proven to be useful for quantifying erosion in situ in areas with long-lasting glaciations.
Thermochronology, tectonic, geomorphic and glacial setting of the DML mountains. The DML mountains in East Antarctica form a c. 1500 km long, mostly coast-parallel, high-standing passive continental margin escarpment (Fig. 1). The escarpment is located c. 200 km inland of the present-day coastline, forming the boundary between the Antarctic polar plateau to the south and the coastal areas to the north. With mountain tops reaching c. 3300 masl, and deeply incised fjords down to c. 2000 mbsl (i.e., Jutulstraumen, Fig. 1), a total relief of more than 5000 m is present in DML. Bedrock exposures are dominated by Archean to early Neoproterozoic ( Supplementary Fig. 1), often granitic, apatite-bearing and zircon-bearing rock types that are very well suited for low-temperature thermochronological studies [17][18][19][20][21] . In the west, the Precambrian basement is unconformably overlain by remnants of late Carboniferous-early Permian sedimentary rocks of the Beacon Supergroup, now exposed at an elevation of c. 2000-2350 masl [22][23][24] . This unconformity is a part of a distinct paleosurface that can be traced over large parts of western DML (wDML) 9,22,25 , and possibly also to central DML (cDML), where mountain tops form a shallowly southward-dipping enveloping surface 16 . Furthermore, the same paleosurface can be identified as the Kukri erosional surface within the Transantarctic Mountains, separating the Devonian-Triassic section of the Beacon Supergroup from the Cambrian-Ordovician basement [26][27][28] . During the rifting of Gondwana, the Beacon sediments of wDML were in turn intruded and buried by up to 2 km of c. 183 Ma basaltic lavas, related to the Karoo mantle plume 15,25,29 . Substantial pre-Jurassic erosion of the Beacon sediments is recorded in southwest DML, where only a few meters of sediments separate the basement from the basalts 22 . AFT and AHe data suggest that the thickness of the basalts diminishes towards the east 15 , with a possible eastern boundary located between c. 4°E and 7°E 16 . While only remnants of the basalts are preserved in wDML, they are correlated to the c. 5 km thick Sabie River Basalt Formation of the Lebombo Monocline on the conjugate margin in SE Africa 30,31 . The escarpment in DML also formed during Jurassic rifting, as a result of dextral pull-apart rifting of East and West Gondwana in the late Jurassic 32,33 , before it successively retreated further inland over time.
Already prior to the Eocene-Oligocene ice sheet initiation 1 , the crest of the DML escarpment was a major drainage divide, with large parts of DML draining towards the inland and only a comparatively small coastal zone draining northwards 33 . As the DML mountains were oriented parallel to the coastline, they likely attracted large amounts of precipitation, which made them an important nucleation site for the initial Oligocene ice sheet 1 . The initial Antarctic ice sheets were wet-based and highly erosive, resulting in high offshore sedimentation rates 3,5,34 . After a short mid-Miocene climate optimum, individual ice caps grew further, amalgamated and expanded to the continental margin, resulting in the onset of full glacial conditions at c. 15 Ma 35 . By that time, Antarctica was isolated, and super-dry conditions persisted 34,36 . The ice sheet had developed a great thickness by then, with coldbased glaciers in upland areas and wet-based glaciers where they were thick enough to experience pressure melting at their base 34 . Since 13.6 Ma, the ice sheet had a similar extent as today, but continued to retreat and advance on Milankovitch cycles. At present, the DML mountains are crossed by several ice divides with very low erosion, as well as deep, highly erosive, fjord-like incisions, such as Jutulstraumen (Fig. 1).
Low-temperature thermochronological data are available from eight studies in western and central DML [14][15][16]25,[37][38][39][40] . The data that have been published so far include 203 apatite fission-track (AFT) analyses, 71 apatite (U-Th)/He (AHe) analyses, 14 zircon fission-track analyses, 11 zircon (U-Th)/ He analyses and 22 titanite fission-track analyses. Some AFT samples were re-analyzed recently to conform to modern analytical methods and calibrated experimental conditions 15,16 . The apatites of wDML record Carboniferous-Cretaceous fission-track ages and early Ordovician-Eocene (U-Th)/He ages, although the datasets are dominated by Cretaceous-Triassic AFT ages and Late Jurassic-Cretaceous AHe ages, respectively 14,15,25,40 . Neoproterozoic-Triassic zircon (U-Th)/ He ages are also recorded in the same area 15 . The published datasets have been interpreted to document burial of wDML by up to c. 2 km of continental flood basalts during the Early Jurassic, followed by Late Jurassic-Cretaceous rift-related cooling 15,25 . Thermal history modeling also documents accelerated cooling during the early-middle Miocene, possibly related to increased glacial erosion 15 or differential exhumation and flexural isostatic rebound as a result of the load of the developing ice sheet 14 . In cDML, apatites record late Carboniferous-middle Miocene (U-Th)/He ages 16,39 and middle Carboniferous-late Cretaceous fission-track ages. The sample suite is, however, dominated by early-Middle Jurassic AFT ages 16,38,41 . Late Devonian-Middle Triassic and late Ordovician-early Permian fission-track ages are also recorded for zircons and titanites, respectively 38,41 . The data have originally been interpreted to represent mainly slow post-Pan-African cooling with phases of accelerated rift-related cooling during the early Jurassic, late Jurassic and middle Cretaceous 38,39,41 . More recently, a scenario with late Paleozoic paleosurface formation and subsequent re-burial by 3-6 km of Beacon sediments has been put forward 16 .

Results and discussion
Thermal history models. The present study utilizes thermochronological data from the two latest studies from western and central DML (Supplementary Notes 1) 15,16 . AFT ages for the 55 samples included in this study span from c. 300 to c. 80 Ma. Mean track lengths range from 10.6 to 13.5 µm (Supplementary Notes 2 and Supplementary Data 1). Etch pit diameters (Dpar) were measured as a proxy for apatite annealing kinetics. The measured Dpars are between 1.2 and 1.8 µm, however, most Dpars are <1.6 µm, typical for fast-annealing, near-end member fluorapatites 42 . Included AHe ages range from c. 260 to 35 Ma (Supplementary Notes 3). The lowest AFT and AHe ages are commonly found at low elevations, specifically at the foothill of the DML escarpment, as well as on the eastern side of Jutulstraumen, while the oldest ages are found on mountain tops along the crest of the escarpment, as well as close to ice divides. For both regions, AFT and AHe ages generally increase with elevation. In wDML, the ages mostly increase to the east, possibly documenting the diminishing thermal influence from the continental flood basalts 15 . In cDML, the ages generally increase towards the hinterland, typical for retreating continental escarpments 16,43 .
The thermal histories for the 55 samples from wDML and cDML show large variations with respect to number of paths tested and found, style of post-Jurassic cooling, range of recorded temperatures at 34 Ma and timing of surface exposure. In addition to the 1000 acceptable paths found, up to 352 good paths (goodness-of-fit ≥0.5) were found. All thermal history models document pre-Carboniferous cooling related to paleosurface formation, followed by burial by late Carboniferous-early Permian Beacon sediments. The models from wDML additionally record early Jurassic near-surface conditions, and subsequent burial by the Jurassic continental flood basalts. The pre-glacial thermal evolution of the samples have been thoroughly addressed previously 15,16 , and are therefore not further discussed here.
At 34 Ma, temperatures of the individual thermal history paths range from c. -25°C to c. 67°C (Supplementary Data 1). The unweighted average temperatures span from c. -4°C to c. 60°C in cDML, and from c. -1°C to c. 57°C in wDML. The precision of the temperature estimates is given by the standard deviation of the temperature distribution. Here, the standard deviation ranges from c. 2°C to c. 18°C (1σ), and correlates negatively to the average temperature (Supplementary Notes 1, Supplementary  Fig. 3, and Supplementary Data 2). The post-Eocene thermal histories, and thus also the temperature distributions and average temperatures, are constrained mainly by 1) the timing of cooling through the partial annealing/retention zone and 2) the apatite fission track length distribution. The negative correlation between standard deviation and average temperature is therefore as expected, as the sensitivity of the measurements are higher close to the partial annealing/retention zones.
By comparing the thermal histories with the mean annual surface temperature curve ( Supplementary Fig. 8), useful information on the mechanism inducing the cooling can be obtained (Fig. 2a, b). More specifically, whether the cooling is linked mainly to simple climatic cooling or rather to erosional exhumation. Seven samples from cDML and four samples from wDML record average temperatures ≤5°C at 34 Ma, indicating (near-)surface conditions at the onset of the glaciation (Supplementary Data 2). The observed cooling in the thermal histories is therefore mainly climatic cooling. In cDML, these samples mainly come from mountain tops, but one sample is derived from the escarpment front of Conradfjella (sample JJ1746; Fig. 1). In wDML, no such trends are observed. On the contrary, the samples that are dominated by erosional cooling are characterized by high temperatures at 34 Ma, rapid post-Eocene cooling and little to no overlap between the individual thermal histories and the surface temperature curve (e.g., sample J7.2.94/4). Fifteen samples from cDML and nine samples from wDML show little to no overlap between individual thermal histories and the surface temperature curve, clearly indicating a component of erosional cooling. In central and northwestern DML, these samples are found within all geomorphic settings, but they are all from lower elevations (<1800 masl). The four samples from southwestern DML all record considerable erosional cooling. These four samples are from the same vertical profile in Heimefrontfjella (Fig. 1) and are the only samples from an area that drains towards the Weddell Sea. The resolution of the thermochronological methods used is known to be limited at very low temperatures, and thus we do not claim to be able to precisely determine the relative contribution of climatic versus erosional cooling. However, we suggest that this method can be used as a first approximation on whether the cooling history of a sample is dominated by climatic cooling (e.g., sample Z7.2.1; Fig. 2b) or erosional cooling (e.g., sample J7.2.94/4; Fig. 2a).
Regional variability of glacial erosion and erosion rates. Large spatial variations in erosion and erosion rates are observed both vertically and laterally across DML. The highest Oligocene paleotemperatures, and therefore also the highest amount of erosion and corresponding time-averaged erosion rates, are generally observed at low elevations in the foothills of the DML escarpment (Fig. 2c, d) Fig. 2 Outcome of thermal history models. a, b Thermal history models exemplified by samples from wDML. External modeling constraints are described in detail in Supplementary Notes 1. a Large difference between the mean annual surface temperature curve 66 and the modeled temperatures indicate substantial erosional cooling. b The similarity between the mean annual surface temperature 66 at the onset of the glaciation suggests that the sample experienced (near-)surface conditions already at the onset of the glaciation, and that the cooling is mainly attributed to climatic cooling. foothill erosion and a more-or-less stationary highland. Such an erosion pattern is indicated by the negative correlation between erosion rate and elevation that we observe in both western and central DML (Fig. 3). Instead, this is a classic pattern for a retreating passive-margin escarpment 45 . This is also supported by the observed age pattern in cDML, with older ages at the coast, younger ages at the escarpment front, and increasing ages towards the inland (Supplementary Figs. 4 and 5) 43 . The glacially modified escarpment, now positioned c. 200 km inland, probably resulted from continuous escarpment retreat after the Jurassic breakup and until the Miocene, when much of the mountains were covered by cold-based glaciers, preserving the topography and preventing the escarpment from further retreat. By assuming a typical escarpment retreat rate of c. 1 km Myr −1 45 , the escarpment may have been positioned c. 20 km further north at the onset of the glaciation than it is today.
Since the Miocene, the erosion was restricted to areas where the ice sheet was thick enough to experience pressure melting at its base, such as within the pre-existing fluvial drainage networks 9,34 or along ancient rift systems (i.e., Jutulstraumen 46 ). This is exemplified by the two samples from the vicinity of Jutulstraumen, which record some of the highest time-averaged erosion rates of wDML (samples H-10 and Z7.27.6, c. 33 m Myr −1 ). The present-day ice flow velocity of Jutulstraumen exceeds 750 m y −1[ 47 , making this the most extensive glacial outlet in our study area, and also the most important conduit for sediments from the continental interior and to the Lazarev Sea 6 . As the overall topography was preserved by the cold-based glaciers, such focused erosion resulted in over-deepening of the pre-existing incisions, making the relief more accentuated. It should be pointed out that the six ice stream proximal samples in this study are derived from a large area, making it difficult to fully disentangle the erosional effect of the ice stream from the other potential erosional mechanisms (e.g., escarpment retreat).
Besides samples from high elevations, very low erosion rates (<10 m Myr −1 ) are also measured for samples close to ice divides (Figs. 2 and 3). Latest from the middle Miocene, these regions were characterized by non-erosive cold-based glaciers. Erosion rates increase away from the drainage divides, which is evident when comparing erosion rates of samples from similar elevations, but at different distances and directions from the drainage divide, e.g., in southwestern DML. The erosion rates observed in Heimefrontfjella are markedly higher than the ones observed in the adjacent Kirwanveggen and elsewhere in DML ( Fig. 1 and Supplementary Data 2). A reason for this may be that Heimefrontfjella in western DML has a different glacial drainage system than the remaining DML. The glacial drainage of Heimefrontfjella is directed towards the Weddell Sea through the Endurance Glacier, while the remaining part of DML drains into the Lazarev Sea and the western Riiser-Larsen Sea (Fig. 1).
The samples from Schirmacheroasen differ considerably from the rest of DML (Fig. 3). These samples show large variations in erosion rates, yet they are all from similar elevations (50-150 masl). This is related to strong local variations in AFT ages, however, the reason for the age variations is currently not sufficiently understood.
Our data suggest up to c. 1.5 km of erosion in wDML, and up to c. 2.4 km of erosion in cDML in the areas with the highest erosion rates. Other samples record only minor erosion in the last 34 Myr. Erosion is thus spatially highly variable, and the irregular distribution of our samples makes it difficult to calculate spatiallyaveraged amount of erosion for all of DML. However, our data appear to be broadly compatible with previously suggested, spatially-averaged eroded thicknesses of c. 0.4-0.5 km for the entire DML 7,8 . Spatially-averaged onshore erosion estimates based on offshore sedimentation show the large-scale trends of the entire continent, whereas our study provides a higher degree of detail in specific locations, documenting the variability related to elevation and geomorphic settings by in situ measurements and modeling.
Implication for ice sheet models. Ice sheet models rely to a large extent on the sub-ice topography. Although the reconstructions of the paleotopography have become increasingly advanced lately 7,8 , they are limited in terms of resolution and accurate quantification of onshore erosion.
The current data from the DML mountains document major glacial modifications to the landscape. The 34 Myr of glacial influence has resulted in a more accentuated relief, but a lower mean elevation. Prior to the glacial initiation, we would therefore expect a higher potential for orographic precipitation in the coast-parallel DML mountains. Additionally, the crest of the escarpment retreated southwards over time. By implementing this into ice sheet models, the locus of initial ice sheet nucleation is likely to be affected and would shift closer to the present-day margin compared to models based on the presentday topography.
Not only onshore, but also offshore, for the marine ice sheet, a more detailed reconstruction of the topography can have important implications. Marine ice sheets (i.e., ice sheets resting on grounds below sea level 48 ) are sensitive to variations in ocean temperature and sea level 49 , and for such ice sheets, the grounding line probably represents the single-most important factor for the developing ice sheet as it mainly controls the ice sheet volume 50 . The position of the grounding line is in turn susceptible to changes in bedrock geometry, and a high resolution of the bedrock elevation models is therefore necessary for accurately tracking the grounding line. For the most recent reconstructions, a horizontal resolution of 5 km is applied 8 . Such resolution is suitable for reconstructing the continental interior, but for coastal regions, a km-scale resolution is recommended 50 . By including in-situ methods for estimating the onshore erosion, the evolution of the bedrock geometry can be pinned at discrete locations, improving the reconstructions of the migrating grounding line.
This study shows that low-temperature thermochronological methods can be applied for constraining variations in erosion over large areas. These methods are, however, not ideal for increasing the general resolution over large areas, due to the scarce distribution of outcrops on glaciated margins, as well as the time-consuming sampling and analyses. Therefore, we suggest that low-temperature thermochronological analyses should be focused in key areas along the coast (e.g., escarpment front and glacial outlets), and used as discrete pinpoints for the pre-glacial to glacial topographic reconstructions, which are often based on the present-day (sub-ice) topography 7,8 . We believe that a better representation of the pre-glacial topography can be obtained by combining these approaches.

Conclusions
This study has documented how low-temperature thermochronological data are useful for quantifying glacial erosion in regions with a long glacial history and major glacial erosion, such as in Antarctica. DML is a classic example of a glaciated passive continental margin, showing substantial variability in the amount of glacial erosion. Along the escarpment, limited glacial erosion is recorded at high elevations and in regions close to the ice divides. Our data indicate up to c. 1.5 km of glacial erosion in wDML, while up to c. 2.4 km of erosion is recorded in cDML. High erosion rates are recorded at low elevation, more specifically at the escarpment foothills and along glacial outlets. On the contrary, limited glacial erosion is recorded at high elevations and in regions close to the ice divides. Based on this erosion pattern, the present-day topography of the DML escarpment is mainly attributed to escarpment retreat, starting after the Jurassic Gondwana break-up, and almost ceasing during the middle Miocene transition from wet-based to cold-based glaciers. Such processes will over time result in a more pronounced morphology, leading to a higher relief, but lower mean elevation. The landscape of DML was therefore quite different at the onset of the glaciation than it is today: The escarpment was located 20-25 km closer to the coast, the mean elevation was higher, but the landscape was less rugged. Such changes in topography will most likely have an impact on the development of the East Antarctic ice sheet, especially with respect to surface temperature, the potential for orographic precipitation and the locus for ice sheet nucleation. Future studies modeling the Antarctic ice sheet should therefore take these findings into consideration when reconstructing the pre-glacial landscape.

Methods
Thermochronological data from two studies from western and central DML form the basis of this study 15,16 . No permissions were required for the sampling. In these two studies, the external detector method was applied for the AFT analyses, using standard etching conditions (5 M HNO 3 , 20 s at 20 ± 0.5°C). Detailed analytical techniques can be found in the original publications 15,16 . However, for the present study, we additionally calibrated the track length measurements 51 , using a length reduction factor of 0.998, to improve the resolution of the t-T models.
For the thermal history modeling, HeFTy 1.9.1 52 was used, focusing on the Cenozoic exhumation history. Thirty-four samples were modeled using both AFT and AHe data, while models for 21 samples are based on AFT data alone. For the samples where AHe data were available, three different helium diffusion models were tested [53][54][55] . The RDAAM model 55 appeared to work best for our data and was therefore preferred (Supplementary Fig. 6). For testing the robustness of our data and our thermal history models, all samples were modeled unconstrained (only start points and end points provided) and constrained, based on geological observations and previously published, independent thermochronological data (see details in Supplementary Notes 1 and Supplementary Fig. 7). The starting point for each model was selected based on previously published thermochronological data from the specific areas (TFT, K-Ar, 40 Ar-39 Ar), recording temperatures of c. 265 and 350°C between c. 450 and c. 975 Ma 40,56,57 . The samples were forced to nearsurface conditions during the Late Paleozoic when the paleosurfaces formed 9,16,22,25 , after that, reheating due to burial by up to 3 km of Beacon sediments was allowed 22,27 . The wDML samples were once more forced to near-surface conditions during the early Jurassic 22 , and then, reheating up to 300°C was allowed during the Jurassic to account for burial by the c. 183 Ma Karoo flood basalts. A surface temperature of -25°C was the end point for all samples.
The model runs were terminated after 1000 acceptable paths (goodness-of-fit ≥ 0.05) were found. Between c. 34,000 and 15 million paths were tested before the ending conditions were met. Although the samples have been modeled previously 15,16 , re-modeling has been done in order to 1) treat samples and data from different studies and distributed over a large region in the same and thus comparable way, 2) test the effect of the selection of different helium diffusion models for AHe data, 3) eliminate potential bias from the use of external constraints, and 4) base our calculations on a significantly high number of individual thermal histories (n ≥ 1000) found, rather than on thermal histories tested (i.e., n = 100,000 15,16 ).
The temperature of each individual, modeled t -T path at 34 Ma was then extracted and an unweighted mean temperature at 34 Ma was calculated for each sample. These generally reproduce well for the unconstrained and constrained model runs (Supplementary Notes 1 and Supplementary Fig. 7), although some samples show a slight shift in modeled temperature distribution (e.g., sample JJ1890, Supplementary Fig. 7). We consider the models with additional external constraints to be more realistic and continue to use temperature estimates only from these models henceforth. The thermal history models are documented in Supplementary Fig. 8. The unweighted mean temperatures have been converted to thickness of eroded overburden, using a geothermal gradient of 25°C km −1 for the basement rocks and Beacon sediments [58][59][60] . As up to 2 km of Jurassic flood basalts were present in wDML 15 , we assume that any overburden exceeding the vertical distance between the sample and the Late Paleozoic paleosurface in wDML, were basalts. As basalts are good insulators, an elevated geothermal gradient was applied for the basalts (50°C km −1 ). The selected geothermal gradient will affect our erosion estimates. By applying a geothermal gradient of e.g. 35°C km −1 for the basalts instead (minimum value for the Deccan Traps 61 ), this will only result in a difference of up to c. 400 m in estimated erosion in wDML. Erosion estimates based on the two alternative values for the geothermal gradient for the basalt thickness are reported in Supplementary Data 2. Our selection of geothermal gradient has been tested by adding the calculated overburden thickness to the sample elevations in a vertical profile in both wDML and cDML. The surface elevations of the preglacial landscape then reproduce within 200 m for the four samples in each vertical profile (Supplementary Notes 1).

Data availability
This publication is supported by multiple previously published datasets 15,16 . The data required for modeling can be obtained at https://doi.org/10.5281/zenodo.5634265.