Uplift and tilting of the Shackleton Range in East Antarctica driven by glacial erosion and normal faulting

Unravelling the long‐term evolution of the subglacial landscape of Antarctica is vital for understanding past ice sheet dynamics and stability, particularly in marine‐based sectors of the ice sheet. Here we model the evolution of the bedrock topography beneath the Recovery catchment, a sector of the East Antarctic Ice Sheet characterized by fast‐flowing ice streams that occupy overdeepened subglacial troughs. We use 3‐D flexural models to quantify the effect of erosional unloading and mechanical unloading associated with motion on border faults in driving isostatic bedrock uplift of the Shackleton Range and Theron Mountains, which are flanked by the Recovery, Slessor, and Bailey ice streams. Inverse spectral (free‐air admittance) and forward modeling of topography and gravity anomaly data allow us to constrain the effective elastic thickness of the lithosphere (Te) in the Shackleton Range region to ~20 km. Our models indicate that glacial erosion, and the associated isostatic rebound, has driven 40–50% of total peak uplift in the Shackleton Range and Theron Mountains. A further 40–50% can be attributed to motion on normal fault systems of inferred Jurassic and Cretaceous age. Our results indicate that the flexural effects of glacial erosion play a key role in mountain uplift along the East Antarctic margin, augmenting previous findings in the Transantarctic Mountains. The results suggest that at 34 Ma, the mountains were lower and the bounding valley floors were close to sea level, which implies that the early ice sheet in this region may have been relatively stable.


Introduction
Antarctica's bedrock topography is an important boundary condition that influences the dynamics of the overlying ice sheet [Gasson et al., 2015]. In particular, Antarctic ice sheet stability in regions proximal to the grounding line is heavily dependent on the local ice thickness and bedrock elevation and slope . Near-coastal regions of Antarctica where the ice sheet is marine-based (i.e., the bed is below present-day sea level) are particularly susceptible to rapid grounding line retreat via marine ice sheet instability [Schoof, 2007] and calving mechanisms such as hydrofracturing and ice cliff failure .
The Recovery catchment in East Antarctica is located on the eastern margin of the Weddell Sea ( Figure 1). It is one of the largest and yet least explored marine-based sectors of the East Antarctic Ice Sheet (EAIS) that may be susceptible to such instability mechanisms [Le Brocq et al., 2008]. The area of the catchment is 1.5 × 10 6 km 2 ; it drains~10% of the EAIS and contains~5 m of sea level equivalent, which is approximately equivalent to the entire West Antarctic Ice Sheet [Rignot et al., 2008]. The regional ice velocity field shows that ice flow in the catchment is focused through three major outlet glaciers-Recovery, Slessor, and Baileywhere velocities reach almost 1000 m/yr at the grounding line [Rignot et al., 2011]. These glaciers are the major arteries that drain the EAIS into the Filchner Ice Shelf (Figure 1).
Situated between these outlet glaciers are the Shackleton Range and Theron Mountains, over which ice velocities are less than 10 m/yr. The strong bimodality in the ice velocity field is reflected in, and caused by, the bedrock topography ( Figure 1). The summits of the Shackleton Range and Theron Mountains protrude above the EAIS as nunataks at up to 1.8 km above sea level, while the bed at the floor of the Recovery, Slessor, and Bailey troughs is as deep as 2.5 km below sea level; the ice thickness in these glaciers exceeds 3 km [Fretwell et al., 2013]. The troughs trend E-W and are 300-500 km long and 50-100 km wide. This fjord-like landscape renders the Recovery catchment particularly susceptible to ice sheet retreat in a warming world [DeConto and Pollard, 2016]. However, if the subglacial landscape has evolved significantly in the past 34 Ma, the response of this sector to climatic and oceanic change in the past may have been very different compared to that of its modern configuration.
The timing and mechanism(s) responsible for the uplift of the Shackleton Range and Theron Mountains and the subsidence of the Recovery, Slessor, and Bailey troughs remain outstanding questions. Apatite fission track and (U-Th)/He dating indicate multiple phases of denudation and burial of the Shackleton Range in the Mesozoic before final uplift and formation of the present landscape since EAIS inception at 34 Ma [Krohne et al., 2016]. Previous studies have suggested that the uplift of the mountain ranges and subsidence of the troughs are inherently coupled. Sugden et al. [2014] hypothesize that post-Eocene uplift of the Shackleton Range was driven by the regional isostatic response to glacial overdeepening and erosion within the Recovery and Slessor troughs. Furthermore, they speculate that the observed tilt of the Shackleton block, with the highest elevations along the southern escarpment (Figures 2 and 3) occurred because excavation of the larger Recovery trough caused more flank uplift than the smaller Slessor trough.
The Recovery, Slessor, and Bailey glaciers likely exploited preexisting fault systems that separate the metamorphic basement of the Shackleton Range [Tessensohn et al., 1999b;Will et al., 2009Will et al., , 2010 from the Palaeozoic Beacon sediments and Jurassic dolerite sills exposed in the Theron Mountains and the isolated  [Rignot et al., 2011]. EAIS = East Antarctic Ice Sheet. The inset denotes the major Antarctic ice divides [Rignot et al., 2008], the area bounded by the red colored polygon is the Recovery catchment, and the black box denotes the area shown in the main figure. (b) Bedrock topography. The red lines show the major onshore basement faults that bound the Shackleton Range and Theron Mountains [Marsh, 1985]. The submarine Thiel trough is bounded by the Filchner Rift (magenta lines, ticks point to the downthrown side) [Jordan et al., 2013. The black dashed box marks our main study area.
Whichaway Nunataks [Brook, 1972] (Figure 2). These fast-flowing glaciers are bounded by ice surface lineaments that reflect the trend of major subglacial fault systems ( Figure 1) [Marsh, 1985]. These lineaments trend parallel to E-W trending approximately 500 Ma thrust faults in the Shackleton Range [Tessensohn et al., 1999b], regional aeromagnetic lineaments interpreted as reflecting major basement faults , and inferred half-graben basins upstream of the Slessor Glacier Shepherd et al., 2006]. These observations support the hypothesis that the deep subglacial troughs are structurally controlled [Patton et al., 2016]. Jurassic extension and horst-and-graben formation have been recognized in the adjacent Weddell Sea Rift System [Jordan et al., 2016, and references therein] and also onshore, in particular in Dronning Maud Land where the Jurassic Jutulstraumen Rift has been imaged [Ferraccioli et al., 2005a[Ferraccioli et al., , 2005b. The Shackleton Range and Theron Mountains may therefore represent fault-bounded horst blocks that experienced tectonic uplift and tilting [Skidmore and Clarkson, 1972].
The relative roles of erosion-driven and tectonic uplift in driving Shackleton Range and Theron Mountains uplift have yet to be quantified, in contrast to the Gamburtsev Subglacial Mountains [Ferraccioli et al., 2011;Paxman et al., 2016] or the Transantarctic Mountains (TAM) [Stern et al., 2005], where the relative roles of these processes have been addressed with the aid of quantitative modeling. In this study, we use 3-D  [Marsh, 1985] are marked by the red dashed lines. The arrow marks the direction of grid north in the adopted polar stereographic projection. The inset shows the location of the study area within Antarctica. (b) Field photograph of a peneplanation surface exposed on Stephenson Bastion in the Shackleton Range (location marked by blue star in Figure 2a). (c) Profile X-Y across the Recovery catchment (VE = x50). Bedrock (black line) and ice surface (blue line) topography were assembled from three RES flight lines. Ice flow direction is out of the page. The green dashed lines highlight the tilting of the mountain blocks. The schematic red lines mark the position of the faults (the arrows show the inferred sense of dip-slip motion).

10.1002/2016JB013841
flexural isostatic models to quantify for the first time both the mechanical unloading associated with normal border faults and erosional unloading associated with Cenozoic glacial incision in the Shackleton Range and Theron Mountains region. Our results have significant implications for understanding the evolution of paleotopography in this part of East Antarctica and for assessing how these changes in topography may have influenced the early history of the EAIS.

Geophysical Data Sets
This study utilizes bedrock elevation and free-air gravity data acquired during a number of recent airborne geophysical surveys over the previously unexplored Recovery catchment.

Bedrock Topography
We collated onshore ice thickness data from a series of recent airborne radio-echo sounding (RES) surveys over Coats Land and the Recovery catchment, including ICEGRAV (2013) [Ferraccioli et al., 2014], Operation IceBridge (2009-2012 [Leuschen et al., 2010[Leuschen et al., , updated 2016, and a 2001/2002 survey of the upper reaches of the Bailey Ice Stream and Slessor Glacier [Rippin et al., 2003;Bamber et al., 2006;Shepherd et al., 2006]. We also include direct ice thickness measurements that were previously incorporated into Bedmap2 [Fretwell et al., 2013] (data coverage is shown in Figure S1 in the supporting information).
RES profiles reveal that the Theron Mountains exhibit a lightly dissected mesa-like topography, whereas the Shackleton Range is more heavily incised ( Figure 2). The mesas in the Theron Mountains resemble those observed westward of the TAM, which are interpreted as the result of the Ferrar dolerite sills capping Beacon Supergroup sedimentary rocks [Ferraccioli et al., 2001[Ferraccioli et al., , 2009Studinger et al., 2004]. Both lithologies are also present in the Theron Mountains [Brook, 1972], hinting at a common mode of formation. A number of plateau surfaces are also exposed at up to 1.8 km above sea level in the Shackleton Range (Figure 2) [Skidmore and Clarkson, 1972;Kerr and Hermichen, 1999]. The plateau surfaces in the Shackleton Range cut different geological units; they do not reflect a stratigraphic dip slope but instead are surfaces that experienced erosion and were subsequently uplifted. These plateaux have been interpreted as remnant fragments of the Devonian Kukri Peneplain, a flat, once-continuous undulating erosion surface which is extensively observed in the TAM [Stern and ten Brink, 1989;Fitzgerald, 1994;Tessensohn et al., 1999a]. along five equally spaced (~20 km spacing), parallel lines. The profiles were isostatically rebounded to remove the effect of present-day ice sheet loading using an elastic plate model with T e = 20 km, which corresponds to our regional T e estimate (section 3.1). Only lines 1, 3, and 5 are shown for clarity. We gridded the new flight line data together with the existing direct ice thickness measurements from Bedmap2 [Fretwell et al., 2013] using a 2 km grid mesh with a continuous curvature tensional spline algorithm [Wessel et al., 2013]. The grid was masked to remove interpolated values more than 10 km from the nearest data point. These grid nodes were replaced by ice thickness values from the Bedmap2 compilation; while there are no direct ice thickness estimates in these areas in Bedmap2, approximate ice thicknesses have been computed using satellite-derived gravity field models [Fretwell et al., 2013]. We then subtracted the ice thickness grid from the surface digital elevation model (DEM) [Fretwell et al., 2013] to produce a bedrock DEM (Figure 3). Offshore bathymetry data were taken from Bedmap2.
We took the spectral average [see, e.g., Bassett and Watts, 2015] of an ensemble of five profiles crossing the mountain ranges. The ensemble average enhances the "common" features of the topographic profiles, such as the tilted plateau surface and the deep U-shaped glacial troughs, while at the same time suppresses the effects of the more localized dissection of the plateau surface by cirques and rivers. It can be seen that the Shackleton Range is on average tilted by 1.2°to the north and the Theron Mountains are tilted by 0.8°to the south (Figure 3).

Free-Air Gravity Anomaly
Our new free-air gravity anomaly (FAA) grid for the Recovery catchment ( Figure 3) was generated from flight line data from the Operation IceBridge [Cochran andBell, 2010, updated 2016] and ICEGRAV 2011 and 2013 surveys [Ferraccioli et al., 2014]. Continuation to 500 m above the bedrock elevation, crossover analysis, and leveling of the lines was performed, and a satisfactory standard deviation of 1 mGal at crossovers between intersecting flight tracks was achieved. Gravity data were gridded at 2 km horizontal spacing using a continuous curvature tensional spline algorithm [Wessel et al., 2013]. The grid was masked to remove interpolated values more than 10 km from the nearest data point. Uncertainties in the FAA grid were estimated to be ±2 mGal. The FAA grid was used in conjunction with the bedrock topography grid to estimate the regional flexural rigidity of the lithosphere (section 3.1).

Effective Elastic Thickness Estimation
The isostatic response of the lithosphere to (un)loading may be computed by modeling the lithosphere as a flexed elastic plate overlying an inviscid (non-viscous) fluid [Watts, 2001]. The amplitude and wavelength of the isostatic response are determined by the effective elastic thickness (T e ), a proxy for the integrated strength of the lithosphere [Watts and Burov, 2003]. We employed two independent methods to determine the appropriate T e for the Recovery catchment to see whether the T e values converged.

The 2-D Forward Modeling
The bedrock topography along two RES flight lines exhibits intermediate-wavelength (100-500 km) warping characteristic of plate flexure in response to surface loading (Figures 4a and 4b). In these profiles, the Bailey trough is downwarped toward the elevated Theron Mountains mesa. We envisage that this topography is largely the product of regional erosion of material from within the Bailey trough and normal fault action. Unloading of the material within the trough is equivalent to the loading of a flat sheet by the mesa. The topography is therefore analogous to the loading of a seamount on the ocean floor, except the mesa displaces ice rather than water. We do not explicitly model the mechanism of loading; our 2-D forward model comprised a distributed load approximating the shape of the mountain range (with a topographic density of 2670 kg m À3 ) (Figures 4a and 4b), which was applied to a thin elastic plate with a uniform T e overlying an inviscid fluid mantle (with density 3330 kg m À3 ) (equation (1)). The density of the material displaced by the load and infilling the flexure was that of ice (915 kg m À3 ). We modeled the topography for a series of T e values between 10 and 50 km. The wavelength of flexure is consistent with T e values of 20-30 km. The best fitting T e values (24 and 25 km for the two models) were determined using the root-mean-square misfit (Figures 4a and 4b).

The 3-D Inverse (Spectral) Modeling
The gravitational admittance is a transfer function that describes the relationship between the FAA and the bedrock topography for a range of load sizes (wavelengths) on an elastic plate with a given T e . The observed admittance was computed by taking Fourier transforms of our newly compiled bedrock topography and FAA grids over a 900 km × 900 km window ( of T e values [Watts, 2001] and compared to the observed admittance (Figure 4c), providing an estimate of the average T e value across the region. This method recovers a best fitting T e value of 11 km. However, taking a Fourier transform of data sets with limited lateral extent causes spectral leakage into the result, downward biasing the recovered T e [Kirby, 2014]. A calibration that accounts for the consequent underestimation of T e [Kalnins and Watts, 2009] was used to correct the recovered T e to 18 km ( Figure S2).
Despite the uncertainties associated with the interpretation of the admittance for topography/gravity data sets of limited lateral extent, such as spectral leakage and the fact that the window-based spectral estimates reflect a wide range of spatial and temporal loads [Kirby, 2014], the corrected T e value of 18 km is broadly consistent with the 24 and 25 km results from our 2-D forward models. In our subsequent flexure (c) Comparison of the observed free-air gravitational admittance (red dots with standard error bars) with model curves for T e = 0, 5, 10, 20, and 40 km. The admittance recovers a best fitting T e of 11 km, which is calibrated to 18 km (see Figure S2 in the supporting information).

10.1002/2016JB013841
calculations, we used a 3-D elastic plate model with a uniform T e of 20 km, which is intermediate between our estimates, and tested the sensitivity of our model by running the calculations for T e values between 5 and 50 km ( Figure S3).

Calculation of Erosional Unloading 3.2.1. Spatial Distribution of Erosion
In order to determine the 3-D distribution of eroded material, we used a peak accordance method [Stern et al., 2005;Champagnac et al., 2007]. This approach involves the identification of peaks and flat-topped surfaces in the bedrock topography that are assumed to have not experienced any erosion and the interpolation of a smooth surface between them. The resulting "peak accordance surface" represents the restoration of the eroded material to the topography without accounting for the associated isostatic response. The difference between the accordance surface and the bedrock topography is the eroded material.
We identified over 80 flat-topped surfaces in the vicinity of the Shackleton Range and Theron Mountains in RES flight lines (Figure 2), Google Earth satellite imagery, and field photographs ( Figure 2). We assumed that these now high-elevation plateau surfaces originally formed a contiguous, low elevation, and flat landscape prior to incision into it and thus have not experienced erosion since the onset of continental glaciation at 34 Ma. This assumption is supported by very low (0.10-0.35 m/Myr) long-term cosmogenic nuclide-derived erosion rates on the plateau surfaces in the Shackleton Range [Fogwill et al., 2004;Sugden et al., 2014] and also by a lack of glacial modification of these surfaces [Kerr and Hermichen, 1999]. If the peaks have been lowered since 34 Ma, the amount of erosion will be an underestimate. In addition, we used a spatial filter to identify local highs in the bedrock topography DEM within a circular moving window with a fixed radius of 15 km [following Champagnac et al., 2007;Paxman et al., 2016]. Peaks where the present-day ice velocity exceeds 10 m/yr, and have therefore likely experienced significant erosion, were discarded. The remainder were assumed to have experienced negligible erosion since 34 Ma; a smooth surface was interpolated between the peaks and flat-topped surfaces to produce a peak accordance surface that was assumed to exist just prior to the onset of glaciation at 34 Ma ( Figure 5).
The accordance surface was constructed by (1) adjusting the DEM to account for the loading of the presentday ice sheet using our preferred elastic plate model with a T e of 20 km (see section 3.2.3), (2) sampling the adjusted DEM at the location of each peak, and (3) interpolating between peaks using a continuous curvature tensional spline [Wessel et al., 2013]. The eroded material was calculated by subtracting the ice-free bedrock topography from the peak accordance surface ( Figure 5).
The assumption that the difference between the peak accordance surface and the bedrock topography is entirely due to removal of material by glacial (post-34 Ma) erosion is probably reasonable within the Theron Mountains and Shackleton Range themselves. However, this may not be the case over the large troughs, where some of the difference may also be attributable to tectonic subsidence due, for example, to movement on the border faults. For this reason there is uncertainty in the amount of material eroded from the troughs. We envisage two end-member scenarios for the amount of material that has been eroded from the troughs: 1. Minimum erosion scenario-tectonic subsidence has contributed to trough depth. In this scenario the peak accordance surface is dipped over troughs ( Figure 5), representing a preexisting depression caused by mechanical subsidence on border faults (section 3.3). We dipped the surface such that when the contributions of erosional unloading and fault motion were summed (see section 3.3), the modeled trough depth matched the observed trough depth. It is therefore assumed that the tectonic subsidence was not infilled with sediment. 2. Maximum erosion scenario-subsidence of the trough floors below sea level is entirely attributable to glacial erosion. In this scenario, the peak accordance surface is stretched across the tops of the troughs ( Figure 5) and the difference between the accordance surface and the bedrock topography is all glacially eroded material. Under this maximum erosion scenario, it is assumed that any fault movement predated glaciation and the resulting subsidence of the hanging wall blocks was completely infilled with sediment Shepherd et al., 2006].
Assuming an average eroded rock density of 2300-2700 kg m À3 (reflecting sedimentary and basement rock end-members), the minimum and maximum estimated mass of eroded material were 1.0 × 10 18 kg and 1.5 × 10 18 kg, respectively.
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013841

Offshore Sediment Estimates
The estimated mass of eroded material was compared to the mass of sediment located offshore on the continental shelf. Isopach maps for the Weddell Sea shelf north of the calving front ( Figure 5) have been constructed by interpolating sediment package thicknesses measured from seismic reflection lines [Huang et al., 2014]. Sediments are divided into preglacial (145-34 Ma), transitional (34-14 Ma), and full-glacial (14-0 Ma) sequences based on correlation of seismic stratigraphic facies across lines and age constraints from sediment cores [Huang et al., 2014]. Because of the uncertainties associated with the volume, provenance, and postdepositional reworking of sediment, we determine a maximum and minimum total 34-0 Ma sediment volume under the following assumptions: 1. Material eroded from the Recovery catchment is now located on the southeastern Weddell Sea shelf (eastward of 50°W and south of 75°S), including the Crary Fan [Diekmann and Kuhn, 1999] ( Figure 5). However, the Support Force Glacier and (during glacial periods) the Foundation Ice Stream also drain into the southeastern Weddell Sea via the Filchner Ice Shelf (Figure 1), so some fraction of the sediment will have been derived from this catchment. We assume that between 50 and 100% of the detrital sediment entered the Weddell Sea via the Recovery, Slessor, and Bailey glaciers, since they drain a larger area than the Support Force and Foundation glaciers. 2. Five to fifteen percent of the total offshore sediment is pelagic (biogenic) rather than detrital and therefore was not derived from onshore [Wilson et al., 2012]. 3. The average bulk density of offshore sediment is between 2100 and 2300 kg m À3 [Wilson et al., 2012]. This accounts for uncertainties in the amount of pore space between the grains (i.e., the degree of mechanical compaction) and their lithology.
By computing the total volume of sediment and applying these assumptions, the mass of 34-0 Ma Recovery catchment-derived detrital sediment in the Weddell Sea basin was determined to be 0.66-1.6 × 10 18 kg. The mass of rock eroded from onshore (1.0-1.5 × 10 18 kg) is therefore within the range of uncertainty of the mass of offshore material. It might be expected that the mass of eroded material exceeds the mass of offshore sediment, since material has likely been lost from the shelf and reworked in the Weddell Gyre or by contourite currents, and some sediment may have instead been deposited within interior sedimentary basins [e.g., Shepherd et al., 2006] during the early stages of EAIS development. Therefore, even our maximum erosion scenario does not obviously overestimate the amount of post-34 Ma erosion from the region. Due to a lack of constraints, we do not incorporate post-34 Ma sediment deposition onshore or beneath the Filchner Ice Shelf or the associated isostatic response in our models. Although this leads to an unrealistic gradient in sediment thickness/erosion at the calving front ( Figure 5), sensitivity testing indicates that onshore flexural uplift is insensitive to the amount of offshore erosion/deposition ( Figure S3).

Flexural Isostasy
The flexural isostatic adjustment (w(x,y)) to erosional unloading and sediment loading ( Figure 5) was computed by solving the general equation for the (un)loading (h(x,y)) of an elastic plate overlying an nonviscous fluid [Turcotte and Schubert, 1982].
is the flexural rigidity as a function of spatial dimensions x and y. Density terms represent the load (ρ load ), the material infilling the flexure (ρ infill ), the material displaced by the (un)loading (ρ displace ), and the mantle (ρ mantle ). We assumed values of 9.81 ms À2 for the acceleration due to gravity (g), 100 GPa for Young's modulus (E), and 0.25 for the Poisson ratio (v). By solving equation (1), we calculated the flexural uplift/subsidence due to the removal of the modern-day ice sheet (assuming an ice density of 915 kg m À3 ), the removal of the eroded material (assuming an average eroded material density of 2500 kg m À3 ), and the loading of offshore sediments (assuming an average sediment density of 2200 kg m À3 ).
The unloads were our updated ice thickness grid (section 2.1) and the grid of eroded material (section 3.2.1); the offshore sediment load was the postglacial (34-0 Ma) sediment isopach of Huang et al. [2014]. The modeled flexure is most sensitive to T e , since this governs the amplitude and wavelength of the flexural response; sensitivity testing was carried out by computing the flexure for T e values between 5 and 50 km, a typical range of values for the continental lithosphere [Watts, 2001]. Since crustal-scale faults may introduce a discontinuity in the plate where the flexural rigidity is effectively zero [Watts, 2001], we also tested a "broken plate" scenario, where T e was decreased to zero at the plate break along one or more of the major faults bounding the Shackleton Range and Theron Mountains. We used a fast Fourier transform method [Watts, 2001] to solve equation (1) analytically for spatially uniform T e scenarios and a numerical centered finite difference technique [e.g., Stewart and Watts, 1997] for spatially variable T e scenarios.

Calculation of Mechanical Unloading
The bedrock topography of the Recovery catchment, with broad valleys bounded by faults and uplifted flanks, is typical of extensional terranes. Vening Meinesz [1950] proposed a model for the uplift of rift flanks as a consequence of failure of the lithosphere by normal faulting. In this case, slip on a normal fault causes unloading of the footwall block by removal of the hanging wall; the result is flexural isostatic rebound and uplift of the footwall. Concomitant replacement of footwall crustal rock by the mantle causes isostatic subsidence of the hanging wall block. Long-term normal fault displacement may therefore be modeled as the flexural isostatic adjustment to the rigid uplift/subsidence of the footwall/hanging wall blocks, assuming that the lithosphere retains a finite flexural rigidity during extension [Weissel and Karner, 1989]. The resulting topography resembles a half-graben, and the footwall is flexurally uplifted. Uplift on the shoulders of normal Journal of Geophysical Research: Solid Earth
To determine the contribution of mechanical unloading associated with the border faults to Shackleton Range and Theron Mountains uplift, we modeled the displacement across each fault as the flexural isostatic adjustment to the rigid uplift and subsidence of the footwall and hanging wall [following Weissel and Karner, 1989] (Figure 6). The amount of flexure depends on the elastic thickness, thickness of the faulted layer (the crust), material densities, and dip and heave of the faults. We used our preferred uniform T e scenario (20 km) and a crustal thickness of 35 km . The assumed densities of the crust, infill (air), and mantle were 2670, 1, and 3330 kg m À3 , respectively. For simplicity, we assumed that each fault is continuous, dips at 60°toward the downthrown side, and exhibits a constant vertical displacement (throw) along-strike. We tested the sensitivity of the results to the elastic and crustal thicknesses and the fault dip and found that only T e strongly influences the distribution of flexure ( Figure S4). The amount of extension (heave) was tuned so the modeled displacement matched the observed topography next to the fault(s).
We also incorporated the diffusion of the scarp due to mass-wasting processes, which is given by [Watts, 2001]: where is the flexural response function, h t (k) is the topography after time t, h 0 (k) is the initial topography, κ is the "subduing coefficient," and k is the wave number (the computation is carried out in the frequency domain). This equation assumes that erosion is a diffusive process that transports mass from the uplifted side of the fault to the subsided region (and the resulting flexural isostatic adjustment is computed). The result is to smooth the edge of the fault-generated topography so it is more exponential in form. Values of t and κ were chosen so the modeled scarp slope matched the observed slopes of the mountain ranges.
The (diffused) flexure was calculated in 2-D along a series of 1000 km long profiles (with 10 km horizontal spacing) trending orthogonal to the faults ( Figure 6). These profiles were then gridded to produce a 3-D map of flexure driven by mechanical unloading (Figure 6). The displacement on the four major border faults was superimposed in various combinations to produce the total 3-D fault-driven displacement. We estimated the throw on the faults bounding the Recovery and Bailey troughs by measuring the difference in elevation of the bedrock on either side of the troughs. Elevation differences of 600-700 m provide first-order estimates of the cumulative long-term throw on the faults, assuming that the flexure associated with erosional unloading is approximately symmetrical either side of the troughs ( Figure 5).

Erosion-Driven Uplift
We calculated the contribution of erosional unloading to Shackleton Range and Theron Mountains uplift for our minimum and maximum erosion scenarios. For our preferred T e scenario of 20 km, we find that erosion in the Recovery, Slessor, and Bailey troughs has driven on average between 600 m (minimum erosion scenario) and 800 m (maximum erosion scenario) of flexural uplift throughout the Shackleton Range and Theron Mountains ( Figure 5). This represents~40-50% of the total elevation of the mountain blocks. The greatest amount of flexural uplift (~1 km in the maximum erosion scenario) occurs along the southern flank of the Shackleton Range, which is bounded by the Recovery trough. The Recovery trough is deeper and wider than the Slessor trough, resulting in a larger-magnitude and longer-wavelength erosional unload [Sugden et al., 2014]. However, this differential erosional unloading only confers a maximum northward tilt of 0.2°on the upper surface of the Shackleton Range, compared to the observed tilt of 1.2°N (Figures 5  and 7).
The flexure onshore is very insensitive to the amount and distribution of sediment offshore ( Figure S4). Offshore sediment loading accounts for <3% of uplift/subsidence in the Shackleton Range, Theron Mountains, and bounding glacial troughs. We suggest that this is because the locus of sediment loading Journal of Geophysical Research: Solid Earth

10.1002/2016JB013841
(the southeastern Weddell Sea) is too distal for significant isostatic uplift/subsidence to be transmitted to the inland fjord system, even if a flexurally rigid (T e = 50 km) lithosphere is assumed.
We tested the sensitivity of the model to the assumed T e scenario (Table 1 and Figure S3). However, we found that while the pattern of erosion-driven flexure is sensitive to the assumed T e , no value between 5 and 50 km was able to produce a satisfactory agreement with the observed magnitude and wavelength of mountain uplift. Intermediate T e values of 20-30 km give the best agreement with the observed wavelength of tilting, but the modeled tilt is only~0.2°. We also tested a scenario where the elastic plate was broken along faults bounding the Shackleton Range and Theron Mountains to investigate whether this could reproduce the observed tilting of the mountain ranges ( Figure S3). However, the difference between continuous-and broken-plate flexure is relatively minor except for regions very close to the faults. Under this scenario, hanging wall subsidence was not infilled, meaning that the contribution of erosion was reduced. The match between observed and modeled topography is worse than the maximum erosion scenario, but the relative contributions of erosional and mechanical unloading remain similar.

10.1002/2016JB013841
We find that irrespective of the assumed erosion and T e scenario, erosion-driven flexure accounts for~40-50% of the total elevation (and only~0.2°of tilting) of the mountain blocks ( Figure 7). The misfit between the modeled and observed topography is small on the flanks of the mountain ranges bounding the Slessor trough, but increases toward the flanks of the Recovery and Bailey troughs, where flexure underestimates the topography by up to 800 m (Figure 7).

Fault-Driven Uplift
Erosional unloading due to the removal of material (rock) from the troughs cannot account for the total observed elevations of the northern Theron Mountains and southern Shackleton Range or the observed tilt of the mountain surfaces (Figure 7). We therefore invoked mechanical unloading due to the unloading of the footwall by normal faults bounding the Bailey and Recovery troughs. For the maximum erosion scenario, the depressions created due to subsidence of the hanging wall blocks were assumed to be filled to sea level. For the minimum erosion scenario, the subsidence was not filled. The throw on the faults was estimated as 600-700 m (section 3.3). Incorporating mechanical unloading on these two major fault systems significantly improved the match between the observed and modeled flexural uplift and tilting of the Shackleton Range and Theron Mountains (Figure 7). Modeled tilts agree very well with the 1.2°N and 0.8°S tilting of the Shackleton Range and Theron Mountains, respectively (Figure 7). We find that erosional unloading accounts for 40-50% of the uplift of the mountains and mechanical unloading accounts for a further 40-50%. There is a small residual misfit; some of the topographic signature is likely the result of nonflexural processes, such as brittle deformation on faults. The maximum erosion scenario, where tectonic subsidence is infilled, produces a better overall fit with the observed topography that the minimum erosion scenario (Figure 7 and Table 1). This suggests fault activity and subsequent sedimentation predated glaciation (section 5.2).
In order to calculate a 34 Ma paleotopography, we restored the eroded material to the ice-rebounded topography and computed and subtracted the associated isostatic response (Figure 8). This calculation was based on the assumption that fault activity mostly predated the onset of Antarctic glaciation at 34 Ma (section 5.2) and has therefore not contributed to mountain uplift or trough subsidence since glacial inception. Since we determined a maximum and a minimum erosion scenario, which differ in their respective assumptions of how deep the troughs were prior to glaciation, we present a minimum and a maximum paleotopography (Figure 8).
Our key finding is that the model scenario that produces a best fit between process-oriented model topography and the observed modern topography requires major contributions from both erosion-and mechanically driven flexure (as well as slope diffusion, which is needed to explain the regrading of the fault scarps). Both processes, operating together, are necessary to achieve a satisfactory agreement with the observed elevation and tilt of the Shackleton Range and Theron Mountains. None of the processes alone can satisfactorily explain these observations. The model results are summarized in Table 1.

Other Mechanisms for Shackleton Range and Theron Mountains Uplift
Glacial erosional unloading combined with mechanical unloading provides a simple and elegant model for the asymmetric uplift (tilting) of the Shackleton Range and Theron Mountains. Our results indicate that the  flexural effects of glacial erosion have driven 40-50% of mountain uplift in this region near the Antarctic margin, which is similar to previous findings in the TAM [Stern et al., 2005]. Ongoing glacial erosion of the Recovery, Slessor, and Bailey troughs and associated flexural isostatic uplift also provides a simple explanation for the emergence of the Shackleton Range from beneath the EAIS at 2.5 Ma [Sugden et al., 2014]. Are there other processes that could account for the observed asymmetric pattern of uplift of the Shackleton Range and Theron Mountains (Figure 2)?
Bedrock uplift in the Recovery catchment could be linked to rift flank uplift on the margin of the Jurassic Weddell Sea Rift System. Such a mechanism has been suggested for the early Cenozoic uplift of the TAM on the flank of the West Antarctic Rift System [Stern and ten Brink, 1989;ten Brink et al., 1997]. However, the TAM are very wide (~300 km) for a rift flank, which is in part attributed to a major inferred contrast in T e across the lithospheric boundary between East (T e = 85 km) and West Antarctica (T e = 5 km) [ten Brink et al., 1997]. In contrast, this study indicates that the Recovery catchment is characterized by T e values of 20 km. Moreover, the highest elevations of the Shackleton Range are >300 km from the Filchner Rift, which marks the easternmost extent of the Weddell Sea Rift System [Jordan et al., 2013 (Figure 1). In addition, the trends of the faults inferred flanking the Shackleton Range are approximately orthogonal to the Weddell Sea Rift System. Together, this suggests that although Jurassic rift flank uplift and passive margin Another option is that the observed asymmetry in the topography is the result of spatially variable erosion rather than spatially variable uplift. However, this is unlikely to be the case, since the tilt is observed in the mesa/plateau surfaces (Figure 2), which have experienced negligible incision and cut different geological units [Kerr and Hermichen, 1999;Sugden et al., 2014;Krohne et al., 2016]. The presence of surfaces that all tilt away from the region of unloading and have slopes that are not the same as geological dip slopes is strong evidence for flexural tilting [Watts et al., 2000]. Krohne et al. [2016] have proposed a model in which thick (up to 3.4 km) sedimentary basins formed in the region following the opening of the Weddell Sea. Post-Jurassic sediments are not observed in the regional outcrops, suggesting that these sequences have been eroded. Could erosion of this overburden, which is not considered in our models, have driven isostatic uplift of the Shackleton Range and Theron Mountains? Erosion of sediments would indeed drive isostatic uplift of the underlying bedrock. However, removal of a spatially uniform overburden is not capable of driving spatially variable (asymmetric) bedrock uplift as is observed. Moreover, if the top of the sedimentary basin was close to sea level [Krohne et al., 2016], sediment erosion could not uplift the top of the bedrock (i.e., the basin floor) to above sea level, because the unload (the sediment) is less dense than the material it displaces (the mantle).

Timing of Fault Activity
While the timing of glacial incision is well constrained to the last 34 Ma [Coxall et al., 2005;Thomson et al., 2013;Krohne et al., 2016], the timing of fault activity remains a source of uncertainty. The inferred faults that bound the Shackleton Range and Theron Mountains lie approximately parallel to major Pan-African age thrust faults and proposed crustal-scale transpressional shear zones recognized within the Shackleton Range itself and in western Dronning Maud Land [e.g., Jacobs et al., 2015]. Recent thermochronology studies indicate a period of significant exhumation in the Shackleton Range area at approximately 190-180 Ma, which is attributed to the onset of crustal extension in the Weddell Sea Rift System [Jordan et al., 2013 and widespread mafic magmatism associated with Ferrar Large Igneous Province [Krohne et al., 2016]. A renewed period of exhumation at approximately 120-100 Ma is attributed to a change in spreading direction in the oceanic crust north of the Weddell Sea Rift System, which may have triggered oblique transtension onshore [Krohne et al., 2016].
Because 120-100 Ma is the most recent episode of exhumation prior to glaciation at 34 Ma [Krohne et al., 2016], it is the most likely time at which the faults inferred to bound the Shackleton Range and Theron Mountains were last active. The time between the conclusion of fault activity and the onset of glaciation was likely relatively short in order to preserve the (asymmetric) topography associated with faulting ( Figure 2). Although the faults may have been moving since 34 Ma, as has been inferred in other regions of East Antarctica [Cianfarra and Salvini, 2016], there is no geological evidence for this in the Shackleton Range region. Furthermore, the presence of fluvial valley slopes flowing toward the Slessor trough close to sea level would appear to rule out significant post-34 Ma tectonic uplift [Sugden et al., 2014]. However, valley incision can have a strong positive feedback on the growth and life span of major range-bounding normal faults in extensional systems [Olive et al., 2014]. The offsets on the range-bounding faults in the Shackleton Range region may therefore, in part, be the result of glacial or preglacial (fluvial) erosion within the troughs. We speculate that the ongoing process of erosion-driven isostatic uplift is accommodated on these faults, as has been suggested for the Lambert Glacier region [Phillips and Läufer, 2009]. The resultant unloading of the footwall by the hanging wall would also contribute to the total flexural uplift, highlighting that the faults were likely the cause and effect of uplift.

Landscape Evolution
The landscape evolution of the Shackleton Range region since Gondwana breakup was likely dominated initially by rifting in the Weddell Sea (commencing at approximately 180 Ma ), uplift of the passive continental margin, and dissection of the landscape by continental river systems [Sugden et al., 2014;Krohne et al., 2016]. With the locus of uplift along the continental margin, it is likely that Jurassic-Cretaceous river systems initially flowed eastward. At some stage, the passive margin was breached at the Journal of Geophysical Research: Solid Earth 10.1002/2016JB013841 location of the present-day confluence of the Recovery, Slessor, and Bailey glaciers; this could have occurred prior to or after glaciation. The modern ice streams exploit this breach today-it is the point through which the entire catchment drains into the Filchner Ice Shelf (Figure 1).
Plate reorganization at approximately 120-100 Ma triggered activity on the faults inferred to flank the Shackleton Range and Theron Mountains and drove exhumation of the region [Krohne et al., 2016]. Our models indicate that the fault systems that drove the majority of mountain uplift were those bounding the Recovery and Bailey troughs. However, magnetic modeling indicates that further upstream the Slessor Glacier is underlain by a sediment-filled half-graben . Because the topography prior to faulting is unconstrained, our models cannot be used to estimate the total amount of uplift on the faults. Our models suggest that the amount of uplift driven by the Recovery and Bailey faults was~600-700 m greater than by the Slessor faults. As well as following the location of the preexisting fault systems (see below), the Figure 9. Cartoon showing the proposed landscape evolution of the Shackleton Range region. The geometry of the region has been simplified to two elongate, subparallel mountain ranges bounded by three subparallel troughs. (a) Fault activity in the Jurassic/Cretaceous uplifted the Shackleton Range (SR) and Theron Mountains (TM) blocks. The troughs were filled with sediment. Large river networks drained the continental interior, flowing westward toward the Jurassic age passive margin. (b) After 34 Ma, the region was covered by the early East Antarctic Ice Sheet (EAIS). Valley floors subsided due to the loading effect of the ice sheet. (c) By the Quaternary, the EAIS had grown to continental scale, and three large ice streams had excavated large overdeepened troughs. The location of these troughs was controlled by the preexisting fault structure and river networks. Ice flow is from east to west. Erosional unloading in the troughs drove isostatic bedrock uplift of the Shackleton Range and Theron Mountains, causing the peaks to emerge from beneath the EAIS as nunataks [Sugden et al., 2014]. As a result of ice sheet loading and glacial erosion, the floors of the glacial troughs now lie up to 2.5 km below sea level. Assuming that fault activity had ceased by the Late Cretaceous, significant (500-1000 m) topography must have existed in the Shackleton Range region prior to glaciation (Figures 8 and 9). During or shortly after faulting, the grabens bounded by the faults were likely infilled with sediment [Krohne et al., 2016] and rivers exploited the structurally controlled topography and cut the valley floors to base level [Sugden et al., 2014] (matching our "maximum erosion scenario"; Figure 5). These river networks (Figures 8 and 9) would have flowed westward if the passive margin had been breached by this time; near the head of the Recovery trough, rivers may have drained east into the Recovery Lakes [Bell et al., 2007] (Figure 8).

Journal of Geophysical Research: Solid Earth
From 34 Ma onward, the landscape has been shaped significantly by glaciation. The modern landscape of the Recovery catchment bears the hallmarks of selective linear erosion [Sugden and John, 1976] by the EAIS. For selective linear erosion to occur, an existing (lower amplitude) topographic feature must have existed prior to glaciation [Sugden and John, 1976;Wilson et al., 2012]. The Recovery, Bailey, and Slessor glaciers therefore likely exploited preexisting depressions controlled by the faults flanking by the mountains and occupied by rivers prior to glaciation (Figure 9). The focusing of ice through the preexisting troughs initiated a strong positive feedback, whereby the troughs were rapidly overdeepened by fast-flowing, warm-based erosive ice, while the peaks of the neighboring mountain blocks were protected (and isostatically uplifted) beneath slowmoving, cold-based nonerosive ice [Kessler et al., 2008]. We estimate that~2 km of rock has been eroded from the Recovery, Slessor, and Bailey troughs. This implies long-term average vertical erosion rates of~0.06 mm/yr, which is consistent with observed erosion rates beneath modern polar glaciers [Koppes et al., 2015].

Implications for Past Ice Sheet Dynamics
The evolution of the bedrock topography of the Recovery catchment has significant implications for the dynamics and stability of past Antarctic ice sheets. With more subdued topographic relief at 34 Ma (Figure 8), topographic steering of the ice sheet would have been less effective during the early stages of glaciation. Early ice sheets may therefore have simply overridden the mountains and troughs. As the bed within the troughs was progressively overdeepened, topographic steering will have become more effective, allowing ice and subglacial erosion to be focused through the troughs [Kessler et al., 2008]. Gradually, ice thicknesses and flow velocities will have increased in the troughs and decreased over the mountain blocks [Sugden et al., 2014].
By correcting for erosion and erosion-driven uplift, we have shown that the Shackleton Range and Theron Mountains were~700 m lower at the time of EAIS inception at the Eocene-Oligocene climate transition (34 Ma) than today ( Figure 8). Furthermore, the Bailey, Slessor, and Recovery trough floors were likely close to sea level at this time (Figure 8), compared to almost 2.5 km below sea level today. This paleotopography therefore provides a new input for models of early ice sheet initiation and evolution. Crucially, the implication of the reconstructed topography is that the early ice sheets were less responsive to climate and ocean forcing, because the bed was not significantly overdeepened below sea level.

Conclusions
In this study, we have presented new compilations of radar and gravity data over the previously largely unexplored Recovery catchment and used the data sets to quantify for the first time the roles of erosion-driven and tectonic uplift. The 2-D forward and 3-D inverse (spectral) modeling indicates that the Recovery catchment is characterized by T e values of around 20 km. Our 3-D flexural models show that erosion-driven uplift has driven a substantial amount (~700 m) of post-Eocene uplift of the Shackleton Range and Theron Mountains, augmenting the previous study of Sugden et al. [2014]. However, the model results show that erosion alone cannot account for the elevation or the tilt of the Shackleton Range and the Theron Mountains. We propose that the Recovery, Slessor, and Bailey glaciers are structurally controlled. The glacially overdeepened troughs were superimposed on preexisting fault-bounded half-grabens that may have been active during Jurassic rifting and Cretaceous intraplate faulting as proposed from independent recent thermochronology studies [Krohne et al., 2016]. Overall, our results indicate that the Shackleton Range and Theron Mountains were likely~700 m lower and the bounding valley floors were close to sea level at the Eocene-Oligocene climate transition at 34 Ma. This has important implications for developing more robust models of the dynamics and stability of the early EAIS.