Temporally stable surface mass balance asymmetry across an ice rise derived from radar internal reflection horizons through inverse modeling

Abstract

Ice rises are locally grounded parts of Antarctic ice shelves that play an important role in regulating ice flow from the continent towards the ocean. Because they protrude out of the otherwise horizontal ice shelves, ice rises induce an orographic uplift of the atmospheric flow, resulting in an asymmetric distribution of the surface mass balance (SMB). Here, we combine younger and older internal reflection horizons (IRHs) from radar to quantify this distribution in time and space across Derwael Ice Rise (DIR), Dronning Maud Land, Antarctica. We employ two methods depending on the age of the IRHs, i.e. the shallow layer approximation for the younger IRHs near the surface and an optimization technique based on an ice flow model for the older IRHs. We identify an SMB ratio of 2.5 between the flanks and the ice divide with the SMB ranging between 300 and 750 kg m⁻² a⁻¹. The SMB maximum is located on the upwind side, ~4 km offset to today's topographic divide. The large-scale asymmetry is consistently observed in time until 1966. The SMB from older IRHs is less-well constrained, but the asymmetry has likely persisted for >ka, indicating that DIR has been a stable features over long time spans.

1. INTRODUCTION

Ice rises are numerous around the coastal ice-shelf belt of the Antarctic ice sheet, and they form where the fast-flowing ice shelves locally reground on topographic highs in seafloor topography. Ice rises protrude from the horizontal ice shelves with parabolically-shaped surfaces that include ice divides from which ice flows slowly towards the ice-rise boundaries (Martin and Sanderson, 1980; Matsuoka and others, 2015). Because ice rises decelerate ice flowing from the continent towards the ocean, they play an important role in defining the dynamics of the sheet-shelf system. Favier and Pattyn (2015), for example, showed that ice rises significantly delay grounding-line retreat during deglaciation and hence impact the timing of sea level changes. In order to better understand the deglaciation history and inevitably also ice-sheet dynamics in the future, there is high interest in collecting observational data that testify to ice-rise evolution (Matsuoka and others, 2015). One possibility for a better understanding is by analysis of the Raymond effect (Raymond, 1983), which exploits a characteristic flow pattern near the ice divide (e.g. Drews and others, 2015; references therein). We take a complementary approach by investigating the spatial and temporal variability of surface mass balance (SMB) on Derwael Ice Rise (DIR) in Dronning Maud Land, Antarctia.

The SMB in Dronning Maud Land is dominated by episodic but strong precipitation events with synoptic forcing contributing more than half of the total accumulation (Schlosser and others, 2010). During these events, moisture is transported in a south-west (SW) direction, deviating from the mean wind direction. In a dedicated modeling study for Dronning Maud Land, Lenaerts and others (2014) investigate the time period from 2001 to 2012. They find that strong precipitation events bring moist air from the Southern Ocean, which rises at the north-east (NE) flanks of ice rises. The orographic uplift causes higher SMB values on the upwind side. On the downwind side, on the other hand, the air column is depleted in moisture causing a precipitation shadow that may extend over adjacent ice shelves. With changing surface slope at the downwind side, near-surface winds accelerate and erode snow, which is subsequently deposited farther away, where winds slow down again. Snow redeposition by winds also occurs over smaller spatial scales whenever surface slope changes (King and others, 2004). Precipitation and drifting snow are the main atmospheric processes forming the large-scale SMB asymmetry across ice rises, which we investigate using radar.

Radar is a standard tool to map the internal stratigraphy of ice. Internal reflection horizons (IRHs) in the radar data originate from changes in density, conductivity and crystal orientation fabric (Fujita and others, 1999), and it is widely accepted that most of the IRHs are isochrones representing former surfaces. The depth and shape of IRHs depends on the SMB, the firn-densification rate and ice flow, leading to the formation of intriguing IRH patterns, which are not easily interpreted (e.g. Arcone and others, 2005; Ng and King, 2011). However, if all the time-varying mechanisms forming the IRH geometry can be disentangled, IRHs are an excellent archive to study the atmospheric and ice-dynamic history. For IRHs shallower than a few percent of the ice thickness, deformation by ice flow is small, and the IRH depth is primarily determined by densification and SMB. This approximation (also known as the shallow layer approximation; Waddington and others (2007)) allows mapping the SMB with higher-frequency radars by tracking the depth of laterally continuous IRHs (e.g. Eisen and others, 2008). IRHs at intermediate depths and below, on the other hand, have experienced strain thinning, deformation and are also affected by the SMB from farther upstream. In that case, the SMB can be reconstructed using a forward model with an SMB parameterization as input to simulate the ice flow and an inverse method that finds the optimal SMB parameters to fit the observed IRH geometry. In a similar manner it is also possible to infer ice flow velocities, although this problem is more ill-posed (Parrenin and Hindmarsh, 2007; Eisen, 2008; Ng and King, 2011). We follow a variation of the approach outlined by Waddington and others (2007), which has been used previously to derive the SMB from Martian polar layered deposits (Koutnik and others, 2009) and from a deep IRH extending over large parts of the Greenland ice sheet (Nielsen and others, 2015).

Our study area is DIR, situated in the Roi Baudouin Ice Shelf, one of the fringing ice shelves of eastern Dronning Maud Land, East Antarctica. The ice rise is quasi-circular and characterized by a distinct ridge in its center, oriented in a SW–NE direction and culminating 350 m above the surrounding ice shelf (Figs 1, 2). Between 2010 and 2012 we carried out several radar surveys across DIR to investigate the Raymond effect causing the formation of distinct arches in the IRHs underneath the ridge. The amplitudes of these arches testify to the ice-divide residence time, which is an important parameter for deglaciation scenarios. However, explaining the arch amplitudes is an underconstrained problem (Matsuoka and others, 2015). For DIR, fits were obtained using different ice rheologies (Martín and others, 2009) and thinning scenarios. The magnitude and spatial variability of SMB, on the other hand, was assumed to be time invariant. Under this set of conditions, it was hypothesized that the ice divide of DIR has existed as a local feature for >5 ka (Drews and others, 2015). However, IRH arches in the flanks and a systematic misfit between modeled and observed IRHs in the south-east (SE) high SMB regime remained unexplained. One goal of the study performed here is to extend the domain spatially, and to use the IRHs in the ice-rise flanks to investigate the assumption of a temporally stable SMB (independent of the IRHs arches beneath the ridge) using an inverse method. If the SMB pattern changes with time, this may trigger a divide migration that can falsely be attributed to, for example, dynamic changes in the surrounding ice shelves. The SMB history is also of interest for ice-core sciences, because ice rises are potential drill sites for the envisaged network of ice cores within the International Partnerships in Ice Core Sciences aimed at retrieving climate records over the past 2 ka and 40 ka (Brook and others, 2006).

Fig. 1. Map of DIR. The thick black line is the 2 MHz profile. The dashed and solid purple lines are the 400 MHz profiles of 2012 and 2013, respectively. The white line is the link to the borehole (blue star). The core site is located at the divide. Radarsat is used as a background image (Jezek and RAMP-Product-Team, 2002). The dotted lines are contours of 300, 350 and 400 m a.s.l. The inset in the bottom left corner depicts the location of DIR in Antarctica. NW and SE (mentioned in the following figures) refer to north-western and south-eastern ends of the profile, respectively.

Fig. 2. (a) Geometry of DIR. Bed and surface are in black and older IRHs detected with the 2 MHz radar are marked in red. The grey zone is the detection range of the 400 MHz radar. (b) Depth of the younger IRHs located in the grey zone of panel (a). The small data gap situated around +10 km is the link between data from 2012 and 2013.

We present our approach of deriving the SMB from younger and older IRHs in Section 2, show the results in Section 3, discuss the robustness of the inversion in Section 4 and close by drawing conclusions if the spatial gradient reproduced by regional SMB modeling (Lenaerts and others, 2014) is clearly observable in the radar data, and whether this pattern is a robust feature in space and time.

2. METHODOLOGY

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, 2012). 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, 2013). 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, 1987) 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, 2007; Drews and others, 2016). 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, 2013). 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 (1965). 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, 2010), 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, 1960, 1962), 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, 2010):

(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, 2003). 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 (1999) and Pattyn (2010), 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 (1999):

(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, 1983). Such arches are clearly visible in the LF data between −2 and +3 km (Fig. 2a) and have been investigated by Drews and others (2015) 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, 2009), 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 (2015) 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, 1963) 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, 1983).

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 (1994, 1996) 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).

3. RESULTS

3.1. Younger IRHs and recent SMB

Figure 2b illustrates the younger IRHs picked from the 400 MHz radar data, all of which exhibit a similar pattern on the broader scale: The younger IRHs are generally deeper in the SE flank than on the north-west (NW) flank. The largest depths are observed 4 km away from the ice divide towards the SE. Superimposed on this broad-scale pattern, the IRH depths vary on kilometer scales and below. Two regions stand out: (1) near the divide where a local minimum is followed by a local maximum, and (2) km ~20 on the SE side, where a similar depth oscillation occurs.

The depth of the younger IRHs was used to infer the corresponding SMB, which naturally shows the inverted spatial patterns as described above (Fig. 3). Errors in SMB depends on the error on the depth/age scale and the RMS of the density/depth model of ±18.2 kg m⁻³, which results in an additional uncertainty on the depth estimation that varies with depth. We estimated the combined effect using standard error propagation resulting in SMB ranges between ±20% for the layers averaging over only a few years. This error does not include that the depth/density function can also vary laterally.

Fig. 3. Spatial distribution of the SMB across the DIR inferred from younger and deeper IRHs. The color code for the shallow IRHs is the same as in Figure 2. The shaded areas denote the SMB uncertainties. The solid black curve is the result of the optimization on all older IRHs. The grey shaded area represents the range of SMB derived while taking into account the depth uncertainty of the older IRHs.

The assigned ages of the younger IRHs are 17, 21, 27, 35 and 46 a before 2012. The SMB estimates of the corresponding bulk layers bounded by the younger IRHs thus represent a temporal average between 1995 and 2012 (averaging 17 a), 1991–1995 (averaging 4 a), 1985–1991 (averaging 6 a), 1977–1985 (averaging 8 a), and 1966–1977 (averaging 11 a), respectively. Layers averaging over shorter time intervals have larger relative errors due to uncertainties in the depth/age scale. The SMB for most layers varies from ~300 kg m⁻² a⁻¹ on the profile end-points in the flanks, to ~800 kg m⁻² a⁻¹ 4 km offset from the divide. This corresponds to variation by a factor of 2.5 over a horizontal scale of 10–20 km. The layer for the time interval 1991–1995 shows a suspiciously elevated SMB magnitude, which we discuss in Section 4.1. All other layers show similar signatures over time.

3.2. Older IRHs and the SMB history

The time-averaged SMB distribution that best fits the older IRHs is shown as a black curve together with the SMB estimates from the younger IRHs in Figure 3. The inferred SMB magnitude is roughly consistent with the results from the younger IRHs. However, the magnitude strongly depends on the parametrization of the depth/age scale, and is thus not reliable. The broad-scale pattern (with low values on both ice-rise flanks and increasing values towards the ice divide), on the other hand, is reproduced. On the SE side, the increase is not linear and deviates systematically from the SMB distribution inferred from the younger IRHs. On the NW end, the SMB inferred from the older IRHs also deviates significantly from the younger estimates and shows an anomalous increase over a few kilometers. Figure 4 shows the corresponding fit between the observed and modeled IRHs. The fit is better for the three shallowest, older IRHs than for the two deepest, older IRHs.

Fig. 4. Comparison between modeled isochrones and the observed older IRHs (red). Black isochrones are the results of the reference run, which included all five older IRHs for the inversion. The light blue curves are inversion results where the IRH closest to the respective blue curve was omitted during the inversion.

4. DISCUSSION

4.1. SMB inferred from the younger reflection horizons

The broad-scale SMB pattern with higher SMB values on the upwind- and lower SMB values on the downwind side is consistent with observations on other ice rises (Nereson and others, 2000; King and others, 2004; Drews and others, 2013). We have not considered here that the depth/density profile derived from the ice core may vary laterally due to changing wind and SMB regimes. Because we have no depth/density data in the ice-rise flanks, we do not investigate this further, but note that a systematic density variation across DIR will alter the SMB magnitude and spatial distribution derived here. However, given the observational evidence from multiple ice rises and additional evidence from atmospheric modeling (Lenaerts and others, 2014), it seems unlikely that this effect plays a major role in shaping the spatial pattern observed here.

King and others (2004) noted a substantial SMB variability on sub-kilometer scales on Lyddan Ice Rise, which is also seen here. They linked it to changes of the local topography using a snow-transport model. Such an effect can, for example, also explain the strong oscillation km ~20 (Fig. 3), which corresponds to a characteristic break-in-slope at the surface (Fig. 2). The oscillation near the ice divide has been investigated in more detail by Drews and others (2015) who found that this small-scale SMB anomaly significantly increases the arch amplitudes of the older IRHs beneath the divide. The origin of this SMB oscillation is not fully clear but could be a combination of wind erosion at the crest with subsequent re-deposition at the downwind side. Regional atmospheric models such as RACMO (e.g. Lenaerts and others, 2014) do not currently have sufficient spatial resolution to resolve such effects, which appear to occur on sub-kilometer scales.

We derived SMB estimates for five time intervals reaching back to 1966. All layers show a similar SMB pattern in space and time, apart from the layer averaging over 1991–1995. To rank the significance of these elevated values, we use simulations from Lenaerts and others (2013), who investigated snow fall anomalies in Dronning Maud Land using the regional climate model RACMO starting from 1979. Based on their SMB results, we calculated averages for the time periods 1991–1995 and 1985–1991 for the IceSAT Basin 6 to which DIR belongs. We find that for the younger time interval, the SMB is indeed higher than for the older time interval. However the relative difference is only 6% and thus barely significant. Even though the regional precipitation of DIR may not be fully captured by the larger-scale atmospheric models, it seems unlikely therefore, that the 1991–1995 outlier fully reflects a temporal SMB anomaly. Because the temporal average of this layer is only over 4 a, uncertainties in the depth/age scale for the IRHs at the layer boundaries are a more likely reason for the observed anomaly. We conclude that the 1991–1995 outlier should be interpreted with care, and ongoing analysis on the age/depth scale (such as using major ions) will help to constrain error sources better in the future. Independent of the SMB magnitudes, all SMB estimates show the same spatial dependence, indicating the precipitation patterns have been temporally stable over the past 46 a. To go further back in time, we will next discuss the results from the older IRHs.

4.2. SMB inferred from older reflection horizons

4.2.1. Sensitivity of the inversion

To test the robustness of the inversion, we run three experiments and compare the results with the SMB profile inferred from the older IRHs presented above. The first experiment is to vary the roughness coefficient (Eqn (5)), κ, within one order of magnitude in both directions to show the impact of the regularization on the results (Fig. 5a). Unsurprisingly, a lower κ results in a higher frequency content. However, using a higher κ seems to improve the large discrepancy between the younger and the older SMB estimates observed on the NW end of the profiles. Table 1 shows the cost function (normalized to the reference run) for the different sensitivity tests. A lower κ improves the matching because a larger part of the uncertainty due to the forward model is transferred to the SMB distribution (resulting in an overfitting of the data). Nevertheless, the roughness term does not change the bulk gradient of the SMB. In this respect, the method can be considered robust.

Fig. 5. Sensitivity analysis with respect to (a) the regularization during the inversion (i.e. magnitude of κ) and (b) uncertainties in the age/depth relationship of the older IRHs, which is estimated with a Nye time scale. In both panels, the thin grey curve is the reference run using unperturbed parameters.

Table 1. Table summarizing the results of the sensitivity analysis with respect to changing the roughness coefficient, the age estimate of the older IRHs and the number of older IRHs involved in the inversion

Shown is the minimum value of the cost function divided by the number of older IRHs included in the inversion and normalized by the value of the ${J_{{\rm ref}}}(\dot a)$ of the reference run.

The second experiment evaluates the influence of the age/depth relationship. The older IRHs were dated by applying a Nye time scale to the mean IRH depth calibrated with a time-averaged SMB of 550 kg m⁻² a⁻¹. Because this method may induce significant errors, we investigate the impact of systematic shift of 10% on the age/depth scale. As expected, this produces an overall shift in the SMB amplitude (Fig. 5b). The analysis of the cost functions favours older age at constant depth. Again the SMB profiles are similar and the result of the reference run stays within the bounds of the two extreme scenarios.

With the last experiment, we test the dependency of the model on the number of older IRHs. The optimization is run five times, each time discarding one of the older IRHs used in the inversion. Each run gives an SMB distribution, which we use to reconstruct the removed older IRH. On Figure 4, each modeled isochrone is obtained by a forward run with the SMB distribution inferred from the optimization on the other older IRHs. The results are similar to the reference run and demonstrate the quality of the inversion. Furthermore, the cost function values for the best match show that each run gives approximately the same results illustrating the global coherence of the method to handle multiple IRHs.

A peculiar feature of the inverted SMB is the apparent increase on the domain ends, particularly on the NW end. This is likely linked to the fact that the SMB is generally less-well constrained by the older IRHs at the boundaries, which has also been discussed for a similar approach by Nielsen and others (2015): Conceptually one can imagine that the SMB value at the last grid point is not constrained at all by the deeper IRHs because the corresponding ice particle exits the domain before it reaches the depth of the older IRHs. Close to the ice divide, on the other hand, the trajectory of the corresponding ice particle passes through all of the older IRHs, which means that the inverted SMB at that location is well constrained by observations. The length of the interval at the boundaries where the SMB is ill-constrained is determined by the amount of horizontal advection: fast ice particles will exit the domain without reaching the depth of the older IRHs. For DIR, flow velocities inferred from GNSS point measurements (cf. Fig. 11c in Drews and others, 2015) are larger in the NW than the SE. To a certain extent this explains, why we see stronger deviations of the inverted SMB compared with the SMB based on younger IRHs in this area.

4.2.2. Efficiency of the inversion scheme

We used Legendre polynomials to parameterize the spatial dependency of the SMB. These polynomials form an orthogonal basis and lead to a better conditioning of the optimization problem, thus necessitating less iterations to converge to the optimal solution. This step deviates from other approaches (Waddington and others, 2007; Nielsen and others, 2015) who use a piece-wise linear function for the SMB parameterization. In that case, the number of model parameters equals the number of grid points. In our case, we have used 100 m grid cells in the horizontal over 25 km, which results in 250 coefficients for the inversion when using the piece-wise linear approach. Using Legendre polynomials has reduced the number of coefficients by a factor of 10. However, the relatively simple flowline geometry of DIR has allowed this, and it needs to be checked whether using Legendre polynomials is also feasible for more complex and 3-D geometries.

4.2.3. Limitations of the forward model

While the inversion works satisfactorily, the forward model may be hampered by different simplifying assumptions. First, the steady-state assumption is questionable. Drews and others (2015) propose that the evolution of the Raymond arch beneath today's divide required modest thinning, however, during that process the ice divide remained laterally stable. Since Lenaerts and others (2014) show that the asymetric SMB distribution is mainly a function of the divide position, we argue that the steady-state assumption, which was imposed here, does not significantly alter the inferred SMB distribution (but may affect its magnitude). Second, we have assumed that ice is frozen to the bed. Drews and others (2015) have calculated the thermal regime of DIR and inferred ice temperatures at the base to be ~5°C. We therefore believe that basal sliding can safely be neglected.

The third caution concerns the plain strain assumption. Our flowline model does not account for convergence or divergence of ice into or out of our modeling domain. We cannot easily verify this assumption with observational data, because ice on ice rises flows too slowly to be reliable detected by remote sensing methods such as interferometric synthetic aperture radar (e.g. Neckel and others, 2012). The surface topography of DIR, however, is relatively well constrained by ground-based GNSS surveys (Drews and others, 2015) and the profile was aligned perpendicular to the surface contours (i.e. along the main-flow direction; Fig. 1). We therefore expect little ice flow across our modeling domain based on the surface data. More worrisome is the bed topography, which exhibits a comparatively flat plateau near the ice divide but sloped bed in the SE flank. While the plateau near the divide is well constrained with other radar profiles, we do not have sufficient parallel radar profiles to determine if this bed slope conforms to the plain strain assumption or not (cf. Fig. 1 in Drews and others, 2015). If it does not, this may explain why the SMB inferred from the older IRHs shows a different spatial pattern between km 5 and km 12 in Figure 3 compared with the SMB inferred from the younger IRHs in this area. Along those lines, Drews and others (2015) also noted that in their forward model, the oldest two IRHs in the SE flank were systematically too deep, which is also the case here. Because the inverse method does not improve the fit compared with a static approach, this renders further support that unaccounted ice-dynamic effects alter the IRH geometry. Given that the unexplained side arch is in close vicinity, those effects may also be of transient nature.

We conclude that the differences in SMB estimates derived from the younger and the older IRHs, respectively, are most likely linked to an underconstrained inversion at the ice-rise boundaries, and to limitations of the forward model, which does not account for potential 3-D effects. Nevertheless, both approaches indicate that the broad scale SMB asymmetry across DIR has been consistent over time for both the younger and the older IRHs. The time intervals older than 1966 are not well constrained here, because we do not have a good handle on the ages of the older IRHs. However, in the study of Drews and others (2015) ice at similar depth is estimated to be older than what is the case here. We therefore tentatively conclude that the age of the oldest IRH provided here (~ka before 2012) is a lower bound.

5. CONCLUSION

Combining geophysical data and modeling, we reconstructed the SMB across DIR in space and time using the shallow layer approximation for younger IRHs, and an inverse method for older IRHs. For the latter, we applied the shallow ice approximation in a plain-strain geometry along a flowline measured with radar, and parameterized the SMB spatially using Legendre polynomials.

The SMB reconstructed from the younger IRHs covers five time intervals dating back to 1966. In all of these intervals the SMB is asymmetric, with higher values on the SE flank, and lower values on the NW flank. This is consistent with observations on other ice rises (Nereson and others, 2000; King and others, 2004; Drews and others, 2013) as well as with atmospheric modeling (Lenaerts and others, 2014). The asymmetric distribution is related to orographic uplift of air masses, which induces an increase of precipitation on the upwind side and a deficit on the downwind side. The observed SMB maximum is located on the upwind side, ~4 km offset to today's ice divide. We find no significant temporal variability.

We also found a time-invariant SMB distribution that matches the geometry of the older IRHs in the ice-rise flanks. This SMB pattern shows the same broad-scale asymmetry across DIR as observed for the younger IRHs. However, some systematic deviations occur. We attribute those to unaccounted 3-D effects in our forward model, and to the fact that the SMB is in principle ill-constrained by observations at the domain boundaries. Because the oldest IRH is likely older than ka we conclude that the SMB pattern across the DIR has not significantly changed over the last millennium. This is consistent with the large Raymond arch observed on DIR, which requires stable conditions for its development (Drews and others, 2015).