2.1. Radar data acquisition and processing

We collected ice-penetrating radar profiles across DIR using two systems: a (higher-frequency, HF) 400 MHz radar (GSSI:SIR 3000) and a (lower-frequency, LF) 2 MHz radar system with resistivity loaded dipole antennas (Matsuoka and others, Reference Matsuoka, Pattyn, Callens and Conway2012). The former detects IRHs up to 50 m depth, while the latter is designed to detect deep IRHs as well as the bed. From hereon IRHs from the 400 MHz radar will be referred to as ‘younger’, and IRHs from the 2 MHz radar will be referred to as ‘older’.

Figure 1 displays the 32 km long HF profile, which is aligned perpendicular to the ridge. The acquisition of this profile was split into two field campaigns in austral summers 2012 and 2013. The profile collected in 2013 extended the profile from 2012 by ~13 km towards the grounding line, both profiles merge at km 10 in Figure 2b. The sampling rate was 10 Hz and the antenna was towed at ~3 m s⁻¹, resulting in an average trace spacing of 0.3 m; traces were recorded for 600 ns and split into 2048 samples. The profile was geolocated with a dual phase Trimble GNSS receiver, which was attached to the radar and which recorded positions at 1 s intervals relative to a non-moving base station at the dome of DIR. The base station coordinates were fixed using precise point positioning from the Canadian Geodetic Survey; the kinematic data were post-processed differentially using the GAMIT/GLOBK software (Herring and others, Reference Herring, King and McCluskey2013). Radar processing included dewow filtering, bandpass filtering and a depth-variable gain function. Using a semi-automated line-tracking routine in OpendTect, we identified five younger IRHs (Fig. 2).

The same profile was also observed with the LF system in the austral summer of 2010 (Fig. 1). We used the same processing techniques as for the HF data and picked the bed together with five older IRHs (Fig. 2a). The positioning of this survey was obtained using a Leica L1 receiver and was processed using the kinematic precise point positioning from the Canadian Geodetic Survey. We derived the depths of older IRHs from the two-way travel-time using a uniform propagation speed (168 m µs⁻¹) and by applying a normal-moveout (NMO) correction (Yilmaz, Reference Yilmaz1987) to account for the unknown ray-path trajectory between transmitter and receiver (which were separated by ~100 m in our case); because all older IRHs are below the firn/ice transition we added a bulk correction of 8.8 m (based on the density from ice-core data discussed in Section 2.2) to correct for the faster radar-wave propagation in firn.

Errors in the depth determination of the older IRHs are not easily quantified. The firn correction assumes a laterally homogeneous firn layer, which may not be the case for ice rises because the SMB and surface winds change over comparatively short horizontal distances. Such spatial differences in the depth/density relationship may bias the traveltime-to-depth conversion by several meters. Errors from the assumed radio-wave velocity in pure ice increase with depth but remain in the same order of magnitude (~5 m at 500 m depth). Errors in IRH picking are more random, and less likely to result in a constant bias. The largest and most ill-constrained error source presumably lies in the NMO correction, which presupposes horizontal IRHs and oblique (but straight) raypaths. In Antarctica, dips of IRHs are typically small and even if the density varies strongly with depth, effects of ray bending are equally small (Barrett and others, Reference Barrett, Murray and Clark2007; Drews and others, Reference Drews2016). However, the positions of the sleds carrying both receiver and transmitters were not well defined: depending on the surface slope, the towing rope was not always fully tensioned and sleds also drifted into an off-axis position compared with the heading of the snowmobile. Contrary to the radio-wave velocity errors, errors stemming from the NMO correction are more significant for shallower IRHs than deeper IRHs. Given the uncertainties in our error estimation, and given that some error sources increase in depth whereas other decrease with depth, we avoid assigning a depth-dependent error of the older IRHs and estimate the combined uncertainty to be ~±10 m.

2.2. SMB reconstruction using the shallow layer approximation

The younger IRHs were linked to a nearby firn core drilled in 2012 (Hubbard and others, Reference Hubbard2013). This 100 m long core provides vertical profiles of temperature, density and age. The density was measured for 48 discrete samples at various depths along the core. The age/depth scale is based on the layer counting of δ ¹⁸O stable isotopes. Because the HF profile is located ~400 m from the firn core, an extra HF profile was collected to link both sites within 6–7 m of the firn-core site (white line in Fig. 1).

To link the core with the younger IRHs, the two-way travel time was converted to depth using the density-depended radar-wave propagation speed parameterized by the mixing formula of Looyenga (Reference Looyenga1965). The continuous density/depth scale needed for the velocity-depth profile stems from a semi-empirical model of firn compaction (Eqn (1); 1, Arthern and others, Reference Arthern, Vaughan, Rankin, Mulvaney and Thomas2010), which was tuned to fit the discretely measured samples of the core. Assuming (1) steady-state conditions (∂ρ/∂t = 0), (2) a surface subsidence rate of 1.375 m (snow) a⁻¹ (based on the upper firn-core layers) and (3) that the density does not vary laterally (Bader, Reference Bader1960, Reference Bader1962), we parameterize the density:

(1) $$\displaystyle{{{\rm d}\rho} \over {{\rm d}t}} = \left\{ {\matrix{ {{c_0}({\rho _{\rm i}} - \rho ),} \hfill & {\rho \le 550\;{\rm kg}\;{{\rm m}^{ - 3}}} \hfill \cr {{c_1}({\rho _{\rm i}} - \rho ),} \hfill & {\rho \gt 550\;{\rm kg}\;{{\rm m}^{ - 3}},} \hfill \cr}} \right.$$

where ρ _i = 917 kg m⁻³ is the density of solid ice and ρ is the local density; c ₀ and c ₁ are rate parameters whose values are defined using the Nabarro–Herring relation (Arthern and others, Reference Arthern, Vaughan, Rankin, Mulvaney and Thomas2010):

(2) $${\left\{ {\matrix{ {{c_0} = \alpha \dot ag\exp \left( { - \displaystyle{{{E_{\rm c}}} \over {RT(z)}} + \displaystyle{{{E_{\rm g}}} \over {R{T_{{\rm av}}}}}} \right)} \cr {{c_1} = \beta \dot ag\exp \left( { - \displaystyle{{{E_{\rm c}}} \over {RT(z)}} + \displaystyle{{{E_{\rm g}}} \over {R{T_{{\rm av}}}}}} \right)} \cr}} \right.}.$$

Here, $\dot a$ is the SMB (550 kg m⁻² a⁻¹), g = 9.81 m s⁻² is the gravitational acceleration, R = 8.314 J mol⁻¹ K⁻¹ is the gas constant, E _c = 60 kJ mol⁻¹ is the activation energy for self-diffusion of water molecules through the ice lattice, E _g = 42.2 kJ mol⁻¹ is the activation energy for grain growth and T _av = −14°C is the mean annual temperature based on the temperature at 10 m depth in the core. T(z) is the temperature at depth z based on temperature profile of the core. The density model is tuned to minimize the difference (in a least-squares sense) with the 48 density samples of the core through adjustment of the parameter α and β. These parameters have been estimated to be 0.02 and 0.015 m s² kg⁻¹, respectively. The model was integrated vertically using a surface density of 400 kg m⁻³ and a pure ice density of 917 kg m⁻³. To estimate the SMB corresponding to a layer bounded by two younger IRHs, we divide the total amount of ice in that layer by the age difference of the IRHs at the boundaries (years were assigned based on summer peaks in δ ¹⁸O). For the youngest IRH, we consider the surface as a boundary of age 0.

The absolute values of $\dot a$ derived here are inflicted with uncertainties, such as the error on the firn-core dating and errors in the density/depth model impacting both the depth determination of the IRHs and the determination of the cumulative mass above the IRHs. The uncertainty on the layer age is estimated to be ±1 a for the shallow depths considered here. This estimate is based on a preliminary comparison of the ice-core isotope profile with corresponding annual peaks of major ions. The error on the density is defined as the RMS difference between the density/depth model and the 48 density samples of the firn core. The combined error on the derived SMB values is then calculated using standard error propagation.

2.3. SMB reconstruction from older IRHs

Using a computationally efficient forward model, we simulate the older isochrone geometry by using a SMB with a parameterized spatial dependence. Optimal parameters for the SMB distribution are found using an inverse method. These steps are detailed in the following paragraphs.

2.3.1. The forward model

For the forward model we make the following assumptions: (1) the flank-flow regime is adequately captured by the shallow-ice approximation (SIA), (2) ice is frozen to the bed, (3) the radar profile is aligned with ice flow, justifying a 2-D geometry with plain strain and (4) DIR is in steady state (i.e. no thinning/thickenning). The validity of these simplifications is discussed in Section 4.2.3.

Using a horizontal resolution of 100 m, we initialize the model with the observed geometry and solve for the velocities without time dependence. Both the ice flow field and the age field are solved using an Eulerian approach (Rybak and Huybrechts, Reference Rybak and Huybrechts2003). We redefine the depth scale as ζ = (b + H − z)/H where z is the elevation above the bed, b is the bed elevation and H is the ice thickness. At the surface, ζ = 0, while the bottom of the ice mass becomes ζ = 1. Following Hindmarsh (Reference Hindmarsh1999) and Pattyn (Reference Pattyn2010), the horizontal flow field is given as a shape function based on the balance flux. In the scaled coordinate system, the horizontal velocity u becomes:

(3) $$u(x,\;\zeta ) = \displaystyle{{n + 2} \over {n + 1}}\overline u (\dot a(x))(1 - {\zeta ^{n + 1}}),$$

where x is the horizontal coordinate, $\overline u $ is the balance velocity, n = 3 is the Glen's Law exponent and $\dot a$ is the SMB. The vertical velocity w is derived from the horizontal flow field assuming incompressibility. We use an analytical solution following Hindmarsh (Reference Hindmarsh1999):

(4) $$\eqalign{w(x,\;\zeta ) = & - \left[ {\displaystyle{{{\zeta ^{n + 2}} - 1 + (n + 2)(1 - \zeta )} \over {n + 1}}} \right]\dot a(x) \cr & + (1 - \zeta )u(x,\zeta )\nabla H + u(x,\zeta )\nabla b.} $$

Each of the velocity components is a function of the SMB $\dot a(x)$ . SIA-based models are fast and appropriate for flank dynamics. However, two ice-dynamic patterns in the IRHs derived from the LF data require special attention: (1) IRHs arch beneath the ice divide. This is due to the non-linear ice rheology, which stiffens ice under low deviatoric stresses encountered at the bed beneath ice divides (also known as the Raymond effect; Raymond, Reference Raymond1983). Such arches are clearly visible in the LF data between −2 and +3 km (Fig. 2a) and have been investigated by Drews and others (Reference Drews2015) to infer ice-rise evolution using a full Stokes model. (2) The secondary arch appearing in the SE flank (between +3 and +5 km).

Simulating the Raymond effect underneath the ice divide is possible using thermomechanically-coupled full-Stokes models. However, matching the arch amplitudes not only requires an anisotropic rheology (Martín and others, Reference Martín, Gudmundsson, Pritchard and Gagliardini2009), but may also involve transient processes, like surface thinning/thickening. Given these complicating factors (and the computational costs that go along with it), the IRH arches beneath the dome are ignored and the analysis is performed on the flanks of the ice rise only. Because Drews and others (Reference Drews2015) speculate that the secondary arch in the SE flank is a remainder of a previous divide (or a divide migration), this area is equally masked out.

2.3.2. Optimization

In the inverse procedure, we determine the SMB distribution that minimizes the difference between the depths of older IRHs and modeled isochrones. This problem is under constrained as neither the SMB of the last millennium nor the ages of the older IRHs are known. To estimate the age of the older IRHs, we apply a Nye time scale in the ice-rise flanks (Haefeli, Reference Haefeli1963) with a time-averaged SMB (550 kg m⁻² a⁻¹, based on the SMB inferred from the younger IRHs) to the mean depth of the older IRH in the region of interest. The average depths of the older IRHs range from 180 to 341 m and the corresponding ages range from 432 to 1030 a before 2012. This dating procedure was chosen because the presence of the Raymond arches prevents the use of the Nye time scale under the divide (Raymond, Reference Raymond1983).

The SMB is initialized with a constant value $\dot a(x) = 550$  kg m⁻² a⁻¹ for all x. A spatial dependence of the initial guess is intentionally neglected, so that the SMB reconstructed from the older IRHs is independent of the SMB reconstructed from the younger IRHs. We then invoke an optimization procedure to determine the spatial distribution of $\dot a(x)$ so that the modeled isochrone ( $z_{\rm j}^{\rm m} $ ) matches the older IRH ( $z_{\rm j}^{\rm o} $ ); z _j is the depth of both the modeled isochrones and the older IRH j. m and o denote modeled and observed quantities, respectively. As we work towards a global solution (i.e. time-averaged distribution of $\dot a(x)$ ), all older IRHs are included to get an optimal solution. The mismatch in depth between the modeled isochrone and the observed IRHs is formulated as a least-squares problem for which we seek $\dot a(x)$ that minimizes the objective function,

(5) $$J(\dot a) = \sum\limits_{\,j = 1}^{{n_{{\rm irh}}}} \displaystyle{1 \over {\sigma _{\rm j}^{\rm o}}} {\kern 1pt} \sum\limits_{i = 1}^{n_{\rm j}^{\rm o}} \omega (i){\kern 1pt} {\rm \parallel} z_{\rm j}^{\rm o} (i) - z_{\rm j}^{\rm m} (\dot a,i){{\rm \parallel} ^2} + {\kappa ^2}\displaystyle{{{{\rm d}^2}(\dot a)} \over {{\rm d}{x^2}}},$$

where n _irh is the number of older IRHs involved, $n_{\rm j}^{\rm o} $ is the number of data points of older IRH j, ω(i) is a mask set to 1 everywhere except under the divide (−2 to +5 km, where Raymond effect is most pronounced), where ω(i) = 0. Because the internal variability of the older IRHs increases with depth, we normalized the mismatch by the variance of the observed IRH ( $\sigma _{\rm j}^{\rm o} $ ) to ensure that each older IRH gets the same weight. The last term is a roughness term (see below).

Given an SMB distribution $\dot a(x)$ , the algorithm calculates the flow and age fields (according to Eqns (3) and (4)) providing the modeled isochrones. The minimization problem is solved with a vector-valued approach with an error vector composed out of the two terms of Eqn (5). The first vector components are the terms of the summation and the last components are the regularization terms. The error vector is used to compute a preconditioned conjugate gradient (computed numerically using small variations in $\dot a(x)$ along the flowline). The subspace trust-region method based on the interior-reflective Newton method (trust-region-reflective algorithm) described by Coleman and Li (Reference Coleman and Li1994, Reference Coleman and Li1996) then determines the modified $\dot a(x)$ -profile for the next iteration. This step is implemented in the MATLAB function lsqnonlin. The iterations stop when the variation of $J(\dot a)$ is below an arbitrarily small threshold (in our case when updates between floating point numbers become indistinguishable).

The spatial pattern of $\dot a(x)$ is parameterized using a linear combination of the first 25 Legendre polynomials, assuming that the SMB is continuous and differentiable in space. The benefits of this parameterization are discussed in Section 4.2.2. The roughness term (given by the second order derivative in Eqn (5)) tunes the smoothness of the SMB distribution and avoids wiggling and overfitting. We chose the roughness factor, κ, so that the roughness term has a comparable variance with the first term of the cost function (Eqn (5)). In the reference run, we chose κ = 3 × 10⁴. This value is discussed in a sensitivity analysis, which also takes into account uncertainties of depth and age for the older IRHs (Section 4.2.1).