1. Introduction
Ice cores provide input for both climate and ice-sheet models. Impurities, stable isotopes and enclosed gases give insight into atmospheric conditions in the past. Measurements of physical properties, such as crystal orientation fabric (COF), improve the parameterization of the flow law of ice. Drill locations at domes or divides often simplify the interpretation of ice-core records, as horizontal ice flux is small, the catchment area is well defined and flow disturbances are minimized (Reference Waddington, Bolzan and AlleyWaddington and others, 2001). However, commonly used approximations (e.g. the shallow-ice approximation) do not apply in these areas since all the components of the stress tensor are significant. This hampers, for example, reliable derivation of the age–depth scale with the corresponding annual-layer thickness. Our aim is to provide these parameters in a comprehensive approach (using satellite, airborne and ground-based techniques) in which surface and bedrock topography, accumulation rates and the internal stratigraphy are combined with a numerical model specifically tailored for ice flow at divides. This allows us to draw conclusions about the ice-flow history and thus to characterize current and past glaciological properties of a potential ice-core site in the context of the International Partnerships in Ice Core Sciences (IPICS).
The IPICS 2k/40k initiative (Reference Brook, Wolff, Dahl-Jensen, Fischer and SteigBrook and others, 2006) aims to increase the spatial and temporal coverage of climate records with a series of ice cores to intermediate depths. Its focus is on the beginning of anthropogenic impacts in the more recent past as well as, among others, on the transition from the Last Glacial Maximum to the Holocene. Drill sites within the IPICS 2k/40k arrays should provide a vertical resolution that allows us to resolve annual and seasonal cycles. The high-accumulation ice ridge Halvfarryggen, Antarctica, is a candidate for these arrays. This Y-shaped ridge faces the northern Weddell Sea and extends the grounded ice seawards, separating the Ekström and Jelbart Ice Shelves to the west and east, respectively (Fig. 1). The German overwintering station Neumayer III is located on the Ekström Ice Shelf with an air-line distance of ∼80 km to the dome of Halvfarryggen. The ice thickness of the marine-based ridge is ∼900 m. The internal stratigraphy beneath the Halvfarryggen ice divides is arched upwards. This effect was first described by Reference RaymondRaymond (1983) as a consequence of the nonlinearity in Glen’s flow law in which ice behaves more stiffly at lower deviatoric stresses. Therefore, the downward motion of internal layers beneath the divide is slower than the downward motion farther away from the divide. The upward arching typically ceases at a distance of a few ice thicknesses from the divide position and has been observed in many radargrams across ice divides in Antarctica (e.g. Berkner Island (Reference SandhägerSandhäger, 1995; Reference Steinhage and BlindowSteinhage and Blindow, 1996), Siple Dome (Reference Nereson and RaymondNereson and Raymond, 1996), Fletcher Promotory (Reference Vaughan, Corr, Doake and WaddingtonVaughan and others, 1999) and Roosevelt Island (Reference Conway, Hall, Denton, Gades and WaddingtonConway and others, 1999)).
Fig. 1. (a) Overview of the survey area (adapted from Reference NeckelNeckel, 2010). The black box marks Halvfarryggen which is zoomed in on the right. The inset delineates the location (red box) within Antarctica. (b) MODIS image of Halvfarryggen (Reference Scambos, Bohlander and RaupScambos and others, 2002 updated 2011). The Y-shaped ice divides appear clearly and are partly accompanied by double-ridge features (green arrows).
We describe the upward arching of layers (defined as a deviation from a smooth curve fitting the flanking isochrones 5 km from the divide) with the term isochrone arch, which assumes that the internal layers as seen via radar are isochrones. Other publications use the term Raymond bump, referring to Reference RaymondRaymond (1983). To distinguish between airborne and ground-based radar profiles we use the abbreviations RES (radio-echo sounding) and GPR (ground-penetrating radar), respectively.
At an ice divide, an ice-flow model that can account for the longitudinal stresses present is necessary. The disadvantage of having to use a more computationally advanced model can be turned into an advantage by recognizing that the more advanced model can also provide additional information regarding flow patterns and flow history. The latter can be exploited to infer the paleo-ice-sheet dynamics in the surrounding region (e.g. the Roosevelt Island Climate Evolution Project (Reference ConwayConway and others, 2012)). In this context, a number of dedicated modeling studies have focused on near-divide flow (Reference HindmarshHindmarsh, 1996; Reference HvidbergHvidberg, 1996; Reference Nereson, Hindmarsh and RaymondNereson and others, 1998a,Reference Nereson, Raymond, Waddington and Jacobelb, Reference Nereson, Raymond, Jacobel and Waddington2000; Reference Nereson and RaymondNereson and Raymond, 2001; Reference Nereson and WaddingtonNereson and Waddington, 2002; Reference Pettit and WaddingtonPettit and Waddington, 2003; Reference Pettit, Jacobson and WaddingtonPettit and others, 2003, Reference Pettit2011; Reference Jacobson and WaddingtonJacobson and Waddington, 2005; Reference Martín, Hindmarsh and NavarroMartín and others, 2009a, Reference Martín and Gudmundsson2012). However, the correct reproduction of the amplitude vs depth distribution of the isochrone arch remains a challenge. For example, the synclinal shape of isochrones at larger depths, which is visible in some radargrams and also referred to as a double bump, has not been replicated with an isotropic rheology. The two-dimensional (2-D), anisotropic case has been considered by Reference Pettit, Thorsteinsson, Jacobson and WaddingtonPettit and others (2007) and Reference Martín, Gudmundsson, Pritchard and GagliardiniMartín and others (2009b), both showing that aligned COF generally increases the arch amplitude relative to the isotropic case. Reference Martín, Gudmundsson, Pritchard and GagliardiniMartín and others (2009b) use an initially random fabric, which develops with ongoing ice dynamics. This model reproduces the synclinal folding and other features such as the flanking synclines of the isochrone arch or the subtle surface concavities observed in some satellite images next to the divides (Reference Goodwin and VaughanGoodwin and Vaughan, 1995). Reference Martín, Gudmundsson, Pritchard and GagliardiniMartín and others (2009b) hypothesize that these features are a consequence of anisotropic ice flow and appear at ice divides that are frozen to the bed and have been stable over multiples of their characteristic timescales (i.e. ice thickness divided by accumulation rate (ice eq.), ∼900 years for Halvfarryggen). Reference Gillet-Chaulet and HindmarshGillet-Chaulet and Hindmarsh (2011) presented a three-dimensional (3-D), isotropic model and compared the results to radar surveys around triple junctions of Antarctic ice ridges (Fletcher Promontory and Berkner Island (Reference Hindmarsh, King, Mulvaney, Corr, Hiess and Gillet-ChauletHindmarsh and others, 2011)) which are in a comparable setting to the Halvfarryggen ice ridge. They predicted that the arch amplitude should be largest beneath the dome and decline with increasing surface slope beneath the individual divide arms. The observations partly match the modeled results, but the increase in arch amplitude with increasing proximity to the dome could not be confirmed.
We draw on the results of these previous studies to interpret the internal structure of Halvfarryggen. In the following we present a multi-method data acquisition together with a comprehensive data synthesis for a numerical model, that will constrain the glaciological conditions of Halvfarryggen in the light of ice coring.
2. Methods and Data
We investigated the internal layering using GPR and RES profiles. Primary results of the observations are the spatial variation in accumulation, the bedrock topography, and the 3-D structure of the isochrone arch around the triple junction. Existing digital elevation models (DEMs) for Halvfarryggen were synthesized in order to obtain a single DEM describing the surface of the dome. The current positions of ice divides were derived from satellite imagery, and the internal stratigraphy is compared with an ice-flow model. This section comprises a discussion of the radar profiles, followed by the derivation of surface topography, divide location and accumulation rates, and a description of the applied numerical model.
2.1. Ground-based radar
Radar systems emit an electromagnetic pulse and record the signal’s reflections caused by dielectric contrasts encountered along the ray. Changes of dielectric properties in ice are primarily caused by variations in density, electrical conductivity and COF. Internal changes in density and conductivity are linked to depositional events on the former surface and as such have an isochronous character. Spatially coherent changes in the dielectric properties appear as internal reflection horizons in the radargrams, which display amplitude vs two-way travel time, or corresponding depth, of the individually recorded traces side by side.
A field campaign in 2007 mapped the upper 100 m of Halvfarryggen (Fig. 2) while a follow-up campaign in 2010 used a low-frequency radar to image the geometry of the isochrone arch in the vicinity of the triple junction. The data from the 2007 field campaign presented here were acquired with a commercial shielded 100 MHz, bistatic RAMAC GPR from Malå Geoscience. Shots were triggered approximately every meter in constant offset mode. Simultaneous measurements with a Trimble GPS receiver were used to geolocate the profiles. The raw data were high-pass filtered (lower cutoff frequency 5 MHz), bandpass filtered (lower/upper cutoff frequency 65 MHz/150 MHz) and amplified with an automatic gain control window of 75 ns. More details of the acquisition system, general post-processing and isochrone tracking are provided by Reference Eisen, Nixdorf, Wilhelms and MillerEisen and others (2004). The profiles were arranged in a star-like pattern with six 13–20 km long legs, centered ∼1 km to the northeast of the dome position (Fig. 2).
Fig. 2. Overview of the airborne RES profiles and ground-based GPR (higher-frequency GPR color-coded with accumulation). Elevation contours are based on the merged DEM. The background image is MODIS data (Reference Scambos, Bohlander and RaupScambos and others, 2002 updated 2011). Sections of profiles shown in other figures are marked with white lines.
The 2010 survey used a low-frequency radar with resistivity loaded dipoles (Reference Funk, Gudmundsson and HermannFunk and others, 1994) and a transmitter as described by Reference Narod and ClarkeNarod and Clarke (1994). The antenna length can be adapted for nominal frequencies of 3.75, 15 or 30 MHz. A Tektronix THS 730 A storage oscilloscope recorded the data as an eightfold stack. The radar was towed with a snowmobile at ∼10 km h−1. A Trimble GPS was mounted on the snowmobile. Shots were taken at a constant time interval and interpolated via the GPS to 20 m postings during the post-processing. This corresponds to a horizontal stacking of about two shots per trace. The 15 MHz profiles were arranged pentagon-like, centered at the triple junction as determined (Section 2.4) from the MODIS (Moderate Resolution Imaging Spectroradiometer) Mosaic of Antarctica (MOA) (Reference Haran, Bohlander, Scambos and FahnestockHaran and others, 2005). The 30 MHz profiles were star-like and centered on the apex of the isochrone arch seen in RES profile 063102a (Figs 2 and 3). Signal processing included static correction, bandpass filtering and background removal (subtracting an ensemble average of traces from the individual traces to reduce background noise) where necessary. The background removal and the values of the cut-off frequencies varied between the individual profiles in order to maximize contrast of internal layering in the upper 100–300 m.
Fig. 3. Radargrams of airborne RES profiles. The horizontal blue line at the top of profile 063102a indicates the overlap with a near-parallel ground-based radar profile in Figure 6. Isochrone arches are clearly visible in all three radargrams. Double bumps first appear in profiles 063102a and 0631202b at ∼510 and 370 m below the surface, respectively.
2.2. Airborne radar
Airborne RES profiles were collected in the area of interest over a number of field seasons during the years 1997–2011. The airborne radar system emits a pulse at a center frequency of 150 MHz with an alternating pulse duration of 60 and 600 ns, resulting in a theoretical vertical resolution of 5 and 50 m, respectively. The short pulse aims at imaging the internal structure of ice sheets, whereas the long pulse is used to measure ice thickness. A radar altimeter records variations in flight altitude. Other components and the technical specifications of the system are described by Reference NixdorfNixdorf and others (1999) and Reference Steinhage, Nixdorf, Meyer and MillerSteinhage and others (2001). Standard signal processing of the airborne data involved ten- or twofold stacking, differentiation, a Hanning tapered low-pass filtering (tapering window between 100 and 150 MHz) and automatic gain control with a 100 ns time window. Reference Steinhage, Nixdorf, Meyer and MillerSteinhage and others (1999) provide further details of the individual processing steps and the software used. Topographic corrections were done with a local DEM (Section 2.3) after variations of flight altitude were compensated using the radar altimeter. Twofold stacked data have a horizontal shot spacing of ∼20 m, depending on the flight speed. The horizontal accuracy of the geolocation of the shots with the airplane’s internal navigation system depends on the spatial baseline to the available reference station, but is usually considered to be accurate within 5–100 m (Reference NixdorfNixdorf and others, 1999). The travel-time to depth conversion was done assuming a constant speed of electromagnetic waves in ice of 169 m μs−1 (Reference Bogorodsky, Bentley and GudmandsenBogorodsky and others, 1985). To correct for higher propagation velocity in the top firn layer, a constant offset of 13 m was subtracted, based on averaged firn-core measurements in Dronning Maud Land (Reference Steinhage, Nixdorf, Meyer and MillerSteinhage and others, 2001).
2.3. Surface topography
Antarctic-wide DEMs are often based on spaceborne altimetry data (e.g. Reference Bamber, Gomez-Dans and GriggsBamber and others, 2009), which perform well in areas with low surface slope but deteriorate at the steeper margins of the Antarctic ice sheet (Reference Griggs and BamberGriggs and Bamber, 2009). In our area of interest the surface slope is relatively large and the Antarctic-wide DEMs deviate up to hundreds of meters (Reference Wesche, Riedel and SteinhageWesche and others, 2009). To describe the topography near the summit of Halvfarryggen we merged two local DEMs, CW (Reference Wesche, Riedel and SteinhageWesche and others, 2009) and RD (Reference Drews, Rack, Wesche and HelmDrews and others, 2009a), together with a local excerpt (BB) from the coastal DEM mosaic derived by Reference BindschadlerBindschadler and others (2011). All three DEMs cover different parts of Halvfarryggen and show smaller standard deviations to ground-truth data than the available Antarctic-wide DEMs. Characteristics of the individual DEMs and the specifics of the merging process are presented in the Appendix.
The grounded part of the resulting DEM was compared to all available ground control points (kinematic GPS along the radar profiles in Fig. 2, ∼120 km of airborne laser altimetry, and Ice, Cloud and land Elevation Satellite (ICESat) data (GLA 12, Release 24; Reference ZwallyZwally and others, 2003)). The comparison yielded a mean and standard deviation of −0.4 ± 11 m. In transitional areas of the initial DEMs, mismatches of 5–25 m are visible. At the northwestern edge of Halvfarryggen they increase to 50 m. However, in the vicinity of the dome no large steps occur and the comparison with ground-truth data is within 10 m.
2.4. Divide location
Ice divides are defined by the change in slope at the top of an ice rise. Because of the overall low slopes near divides, defining the exact position of the ice divide is not trivial. In DEMs, ice divides run perpendicular to the contours and vertical errors transfer into a horizontal misplacement of the divide location. This may be on the order of kilometers, since the slope near the divides is only around 0.002–0.006 (Fig. 4, red lines).
Fig. 4. Different picks (×) for the ice-divide location, based on the DEM and different MODIS scenes. The red continuous line marks the GPS transect (which is only near-perpendicular to the divide). See Figure 2 for location of this profile (071111).
In satellite imagery, on the other hand, ice divides become evident through a change in backscatter (which depends on the surface slope). We manually picked the divide from individually contrast-enhanced satellite images, by following the larger scale, in the direction of the divide coherent border between shaded and illuminated areas. We excluded shorter-scale, irregular contrast variations across the divide. These are evident in the entire scene and probably linked to subtle variations in bedrock geometry. Due to these shorter-scale variations and the double-ridge features (Fig. 1b, green arrows), the inferred divide location does not always coincide with the strongest gradient in gray values.
The derived line has a systematic offset, compared to those based on other scenes. This offset depends on the illumination angle. To estimate this effect, we picked the divides from different scenes with partly opposing look angles (MODIS(1) clipped from the MOA (Reference Haran, Bohlander, Scambos and FahnestockHaran and others, 2005), MODIS(2) from the MODIS image archive (Reference Scambos, Bohlander and RaupScambos, 2002 updated 2011)), and a clipped image from the RADARSAT mosaic (Reference JezekJezek and others, 2002). The results were compared with the ‘true’ divide locations derived from GPS transects (note that the GPS transects are not always perfectly perpendicular to the assumed divide position). An example is shown for the two MODIS images in Figure 4. The RADARSAT image does not have enough contrast in that area for a successful evaluation.
For our dataset, the satellite imagery provides a more accurate location of the ice divide than using the DEM, which does not have the required vertical accuracy. As a result, points based on the MODIS(1) image were selected because they were closest to the GPS-derived divide locations. Based on a comparison with the other satellite-inferred divide positions and a comparison with the GPS data, we estimate the horizontal misplacement to be around ±750 m. The uncertainty increases near the triple junction, where the divides generally appear less distinct.
2.5. Accumulation
In order to derive the spatial variability in accumulation, we followed a standard procedure of accumulation mapping (e.g. Reference EisenEisen and others, 2008). The relative depth of two continuous internal layers was converted into total accumulation using density and age–depth data from the 85 m deep firn core B38 in the center of the profiles (Fig. 2). Details of the firn-core analysis are given by Reference Fernandoy, Meyer, Oerter, Wilhelms, Graf and SchwanderFernandoy and others (2010). The age–depth relationship was found from the firn core via layer counting. The relation between two-way travel time and depth was established using the empirical formula given by Reference Kovacs, Gow and MoreyKovacs and others (1995). It converts density to relative permittivity, which in turn determines the velocity profile in firn.
Error estimates for this approach are composed of errors in determining the depth of the reflectors, the cumulative mass above the reflectors, and their respective age. Reference Fernandoy, Meyer, Oerter, Wilhelms, Graf and SchwanderFernandoy and others (2010) estimated the dating uncertainty for the B38 core to be around ±1 year. The errors for the assigned depth and cumulative mass depend mainly on the density profile. For the B38 core this was measured via the attenuation of γ-rays (Reference WilhelmsWilhelms, 2000) with a relative error of <1%. The spatial variation in density remains unknown, but based on Sorge’s law (Reference BaderBader, 1954) we assume that on short time and spatial scales neither the age–depth profiles nor the density–depth profiles change considerably. Other studies, summarized by Reference EisenEisen and others (2008), give error estimates within 5% for this approach, mostly caused by uncertainties in dating.
To derive accumulation we chose two horizons which could be traced over the longest distance in the radar data. At the core site B38, these have depths of 51.6 and 63.0 m with corresponding ages of 27 and 34 years before (January) 2007, respectively. Several tracked horizons within the radar profile are shown in Figure 5. Accumulation estimates were calculated from the surface to the first layer, surface to the second layer and between the first and second layers, each of them averaging different time periods. The results mostly differ by a constant value. For the B38 firn core, Reference Fernandoy, Meyer, Oerter, Wilhelms, Graf and SchwanderFernandoy and others (2010) reported an accumulation rate of 1257 ± 347 kg m−2 a−1) averaged over a time period of 47 years. They observed a minimal trend towards higher accumulation rates in recent years.
Fig. 5. Internal reflection horizons picked from the 100 MHz GPR data parallel to RES line 063102a (Fig. 2). The isochrone arch develops about 30–60 m below the surface. The apex below this depth has an offset with respect to the maximum surface topography of ∼500 m to the east.
2.6. Modeling internal stratigraphy and the age–depth relationship
To determine the dome’s stratigraphy and the age distribution, we used a transient thermomechanical full-Stokes model that considers anisotropic rheology (Reference Martín and GudmundssonMartín and others, 2012). The model is an improved version of that presented by Reference Martín, Hindmarsh and NavarroMartín and others (2009a). We describe the model briefly here and we refer the reader to those papers for details. The 2-D ice-flow model assumes plane strain across the ridge, a frozen ice/bedrock interface, and an outflow boundary condition with zeroth-order anisotropic shallow-ice approximation that conserves the total ice mass. The temperature model assumes a uniform geothermal flux at the base, no horizontal temperature gradient at the margins and a uniform surface temperature. The primary model input is the surface and bedrock topography and a spatially non-uniform accumulation rate as described above. The averaged surface temperature (−17.9°C) is taken from the AWS11 automatic weather station near the summit (the position overlaps with the B38 triangle in Fig. 2). The model domain was chosen along the profile 063102a (Fig. 3) which intersects the southern divide close to the dome and exhibits double bumps in the lower third of the ice column. The accumulation rate along this profile is based on the underlying, quasi-parallel GPR line (Fig. 2) and was smoothed with a linear regression. The rheology employed is an anisotropic extension of Glen’s flow law (Reference Martín, Hindmarsh and NavarroMartín and others, 2009a). The dependence of rate factor with temperature follows Reference Dahl-JensenDahl-Jensen (1989) for a rheological index n = 3. To extend this temperature dependence to other values of n, we multiplied the rate factor by a constant. For a given n, we chose the constant that produces a better fit between the steady-state modeled surface and the surface topography derived in Section 2.3. The model was initialized, similarly to that of Reference Martín and GudmundssonMartín and Gudmundsson (2012), with a flat surface over the domain and an ice fabric which changes linearly from isotropic at the surface to a vertically aligned single-maximum at the base. The model was run to ten times the characteristic time of the divide, after which the age–depth scale does not change significantly (apart from the stagnant bottom layer).
The model results must be interpreted considering the applied simplifications, which we now discuss: The 2-D approach neglects the geometry for the triple junction and cannot account for the slanted alignment of the RES path across the divide (angle of intersect ∼45°). We expect that the domain is, to a certain extent, imparted by a component of along-ridge flow and, as a result, mass is being redistributed across the divide plane. To account for that, we assumed that the along-ridge flow is smaller than the across-ridge flow (and hence that the along-ridge stratigraphy is homogeneous in a length scale smaller than the ice thickness). Accordingly, we projected bedrock topography and surface accumulation for the model input in a plane perpendicular to the ridge (Reference Martín, Hindmarsh and NavarroMartín and others, 2009a). The same is done for the radar stratigraphy for comparison.
The model assumes that the ice fabric is induced by deformation. It does not include polygonization or migration recrystallization. The dynamical effect of these processes is unclear and here we assumed that their overall contribution is comparatively small. The explicit effect of temperature on the fabric development is small compared to the implicit temperature dependency linked to the ice deformation. The largest uncertainty in calculating the temperature profile is the value of the geothermal heat flux, G. We address its influence on the age–depth estimate with a sensitivity analysis. No direct measurements for G exist in our area, and based on studies by Reference Shapiro and RitzwollerShapiro and Ritzwoller (2004) and Reference Maule, Purucker, Olsen and MosegaardMaule and others (2005) we consider values for G ranging between 30 and 70 mW m−2.
The model uses a constant accumulation through time. Changes in accumulation on a timescale comparable to or smaller than the characteristic time of the divide (∼900 years for Halvfarryggen) may give rise to over- or undersized anticlines in the stratigraphy of the divide area (Reference Martín, Hindmarsh and NavarroMartín and others, 2006). Since the rheology is calculated so that the model reproduces the observed isochrones, this could lead to errors in the estimation of the rheological parameters. The value of the rheological index near ice divides may differ from the modeling standard n = 3 (e.g. Reference Pettit and WaddingtonPettit and Waddington, 2003 and addressed more recently by field measurements (Reference Gillet-Chaulet, Hindmarsh, Corr, King and JenkinsGillet-Chaulet and others, 2011)). We estimate the effect of a varying n on the age–depth relationship with a sensitivity analysis for n = 2,3,4.
3. Results
3.1. Wind- and topography-induced accumulation pattern
The accumulation pattern derived from the 100 MHz data is depicted along the GPR profiles in Figure 2 by the color-coded GPR lines. The western flank of Halvfarryggen has lower accumulation rates than the eastern side, with the ice divides marking the boundary between the two regimes. A map interpolated from point measurements (Reference RotschkyRotschky and others, 2007) confirms this general trend. However, for that study only one data point located on the northwestern side of Halvfarryggen was available, resulting in significantly lower absolute accumulation estimates than observed here. The analysis of the B38 and another firn core (Reference Fernandoy, Meyer, Oerter, Wilhelms, Graf and SchwanderFernandoy and others, 2010) located further south revealed higher accumulation rates at Halvfarryggen, with a decreasing trend towards the north, and we confirm the gradient from east to west, which was already hypothesized by Reference Fernandoy, Meyer, Oerter, Wilhelms, Graf and SchwanderFernandoy and others (2010).
The automatic weather station AWS11 (run by the Institute for Marine and Atmospheric Research Utrecht) is located near the drill location of the firn core B38 (Fig. 2) and has recorded data continuously since its deployment in 2007. The preferred wind direction originates from 110° (geographic) north, i.e. east-southeast, indicating a strong influence of westerly low-pressure systems. Therefore, it seems reasonable to assume that the recent accumulation pattern is dominated by orographic snowfall, which continuously deposits more snow on the eastern flank than the western flank, resulting in differences of up to 50%. A similar gradient across the divide of Siple Dome has been observed by Reference Nereson, Raymond, Jacobel and WaddingtonNereson and others (2000) who also conclude that more snow is deposited on the windward side due to orographic uplift.
3.2. Geometric divide characteristics
The highest points in the GPS profiles and the divide position, as picked from MODIS data, are shown in Figure 6. While the southern and northwestern arm agree well in both datasets, the northeastern divide based on MODIS appears shifted ∼1.5 km farther north compared to the peaks in the GPS profiles. The dome derived from the GPS data appears 1 km southeast of the MODIS-based triple junction. We can only speculate on the origin of these differences, but the ice divides appear less clear in the MODIS imagery as they approach the triple junction. It is therefore likely that the corresponding picks are misplaced, despite the apparent correspondence with the underlying DEM.
Fig. 6. The spatial characteristics of the isochrone arch near the triple junction as seen in the low-frequency GPR profiles. The color code depicts the change in depth of an exemplary internal reflection horizon. The apices (red dots) are visible on the southern and northeastern branch of the ice divides. The maxima in surface elevation of the GPS transects are marked with bold black crosses; the divide is marked with black dots (the faint gray lines indicating the horizontal uncertainty).
Double-ridge features, as mentioned in the introduction, appear in the satellite imagery as faint lines, parallel to the divides. Reference Goodwin and VaughanGoodwin and Vaughan (1995) analyzed examples from Landsat imagery of Fletcher Promontory and concluded that these double ridges are caused by subtle concavities in the surface topography, which they measured with GPS profiles running perpendicular to the ice divide. In the case of Halvfarryggen, parallel lines to the divide are visible in full-resolution satellite images, but their existence and specific location are subjective and by no means distinct.
In MODIS(1), a linear surface feature is visible ∼12 km south of the triple junction on the eastern side of the divide. In the RADARSAT mosaic, both sides of the southern divide in this area are accompanied by dark lineations, and the northwestern divide has a parallel lineation on its western side, ∼13 km north of the triple junction. In MODIS(2) both sides of the northeastern divide are accompanied by parallel lineations, starting ∼14 km north of the triple junction. All double ridges disappear in the vicinity of the triple junction, where the divide itself is also less clear in all images. In that area, the merged DEM and the GPS data, which were acquired together with the radar profiles (Figs 2 and 3), do not show concavities in the surface elevation (e.g. the profile in Fig. 4 crossing the northwestern divide, or the profile in Fig. 6 crossing the southern divide).
3.3. Horizontal and depth-varying architecture of the isochrone arch
A prominent internal reflection horizon with an average depth of 260 m could be linked across different RES profiles spanning an area of about 25 km × 30 km. The spatial interpolation (Fig. 10, further below) illustrates that the isochrone arches beneath the southern and the northeastern divide are clearly visible throughout this area. Below the northwestern divide, the isochrone arch is only weakly pronounced and vanishes a few kilometers away from the triple junction.
The spatially denser low-frequency GPR data in the dome’s vicinity (Fig. 6) show a similar picture, as no upward arching is observed below the northwestern divide. Near the dome, an offset between the apex of the isochrone arch and the peak in surface topography (as determined from the GPS data) becomes apparent. This is also observed in the 100 MHz profile that runs near-parallel to RES line 063102a. In that profile, the initial development of the isochrone arch can be observed in shallow ice (Fig. 5). The upward arching of layers starts at around 30–60 m below the surface. The apex near the dome shows the previously mentioned offset by ∼500 m to the east.
The shape of the isochrone arches at larger depths varies between the different RES lines. Double-peaked arches are visible in multiple RES lines crossing the southern and the northeastern divide. Flanking synclines can be observed in some profiles on the western side of the ice divide at ∼25–75% of the total ice thickness. The internal layers on the eastern side decline more strongly and generally do not show a flanking syncline. It becomes increasingly difficult at larger depths to connect the different RES profiles and we cannot interpolate the synclinal shape of the isochrone arch spatially in the lower third of the ice column.
3.4. Model results
The steady-state geometry is compared to the RES layering in line 063102a in Figure 7. The modeled surface topography places the ice-flux divide position ∼1000 m to the east of the highest point in the DEM profile. The topographic divide does not coincide with the flux divide due to the asymmetry in accumulation. The observed elevation differences are within the uncertainty of the satellite-based elevations. The apices of the modeled isochrones are tilted towards the east in a line perpendicular to the bedrock slope. The tilt in the RES layering is less pronounced and in the opposite direction. Further away from the divides the strong descent in internal layering, particularly at the eastern flank, is not reproduced by the model, while the depth below which the double-peaked Raymond arches occur matches the observation.
Fig. 7. Model output (G = 50 mW m−2, n = 3) along RES line 063102a close to the dome (see Fig. 3 for location). The modeled isochrones (black lines) are shown in comparison to the measured RES layers (red lines). The temperature field is displayed in the background and predicts −7.3°C at the ice/bed interface underneath the divide.
The model (for n = 3 and G = 50 mW m−2) predicts ∼11.4 ka BP ice at 90% of the ice thickness beneath the divide (Fig. 8a). The corresponding annual-layer thickness is ∼7 cm (Fig. 8b). In 30–70% of the ice thickness, the ice age is 400–4000 ka BP, with a corresponding layer thickness of 1–0.3 m. To reinforce our age–depth estimate and to give an idea of maximum bounds depending on particular rheological conditions, we show in Figure 8 the age estimate under the divide using a Dansgaard–Johnsen model (Reference Dansgaard and JohnsenDansgaard and Johnsen, 1969) applied at a distance of five times the ice thickness from the divide. The age–depth scale is traced to the divide along isochrones in the radar data. The analytical approximation in the flanks depends on a parameter representing the changes in vertical strain rate from constant at the top to linear at the bottom. We show variations of age values for the solutions with its minimum and maximum value (Reference Martín and GudmundssonMartín and Gudmundsson, 2012).
Fig. 8. (a) Comparison of age–depth estimates along RES line 063102a. Pink and blue lines are based on the analytical approximation by Reference Dansgaard and JohnsenDansgaard and Johnsen (1969) and depict the age–depth at the western and eastern flanks (8 km away from the divide), respectively. The width of the colored area is linked to the uncertainty in choosing the critical depth for the onset of vertical shearing (extreme cases correspond to 0 and 75% of the ice thickness). The extrema of both age–depth relationships are traced to the divide along internal RES layers (green crosses) so that the horizontal distance between the crosses indicates the uncertainty. The black line is the output of the full-Stokes model beneath the divide. (b) Corresponding estimates for the vertical layer thickness.
The sensitivity analysis (Fig. 9a–c) shows that the age–depth scale does not depend strongly on the specific value of the geothermal heat flux as long as it is <70 mW m−2. Also, for the depth interval where control points are available, the different values for the rheological index are of minor importance for the age–depth scale.
Fig. 9. Sensitivity analysis for (a) varying rheological indices (for G = 50 mW m−2), (b) varying geothermal heat flux (for n = 3) and (c) the temperature profile within the ice for the results in (b).
4. Discussion
Choosing a drill site for an ice core depends on the scientific goals, which vary from retrieving an undisturbed climate record to measuring the physical properties of ice (e.g. to better understand the interplay of COF development and ice deformation). Regional climate patterns and the trade-off between annual-layer resolution and temporal coverage are important for the climate record, while studies concerning ice deformation depend more on the ice-flow regime. We propose that the Halvfarryggen site offers possibilities for both applications.
The regional climate patterns and the corresponding moisture transports in the Ekström Ice Shelf area are influenced by cyclonic activities in the circumpolar trough and katabatic winds with a south–north orientation. This has been investigated in studies using the long-term meteorological observations at the Neumayer stations (Reference Schlosser, Reijmer, Oerter and GrafSchlosser and others, 2004, Reference Schlosser, Oerter, Masson-Delmotte and Reijmer2008; Reference König-Langlo and LooseKönig-Langlo and Loose, 2007). The radar data evaluated here identify the summit region of Halvfarryggen as a high-accumulation site with a spatial variability due to a preferred wind direction from the east-southeast and the symmetrical alignment of Halvfarryggen towards the north (Fig. 2). The accumulation gradient is associated with a variability in age–depth scales and corresponding annual-layer thicknesses (Fig. 8a and b): the western flank has older ice while the eastern flank has a larger annual-layer thickness. Due to the isochrone arches, the depth interval of ∼40–75% of the ice thickness at the divide combines both aspects as ice beneath the divide appears older than at both flanks while the annual-layer thickness is larger (the depth interval corresponds to a time frame of about 700–4000 years BP). A full Holocene record is to be expected down to ∼90% of the ice column at the divide. Ice from the last termination is potentially conserved at larger depths. In spite of the simplifications discussed in Section 2.6, the age–depth estimates of the numerical model compare well with the analytical solutions in the flanks which have been traced to the divide along (isochronal) radar layers. Hence, in cases where few data on ice thickness are available, the general notion of dating internal layers in the flank-flow regime via Dansgaard–Johnsen-type approximations and transferring the dates to the actual divide position is justified by our comparison. However, although our results are consistent, the error created by such approaches can be considerable. Moreover, the extrapolation along isochrones works only for clearly identifiable layers. Especially deep layers in the flanks often cannot be tracked unambiguously to the divide. The simplified dating is thus often only possible for the upper 50% of the ice thickness. The results do not strongly depend on the specific value of the geothermal heat flux. Similarly, the value of the rheological index, in the context of an anisotropic rheology, does not significantly alter the age–depth scale (Fig. 9a–c).
The modeled ice flow qualitatively reproduces the internal layer architecture near the dome (Fig. 7). The results suggest that the ice stratigraphy is the consequence of stable conditions for at least three to five times (about 2700–4500 years for Halvfarryggen) the characteristic time of the divide. This time is needed to produce a sufficiently aligned vertical girdle-type COF in the top half of the ice thickness and a single-maximum COF pattern towards the base. The latter eventually gives rise to the formation of double-peaked isochrone arches due to a subtle tilt in the COF alignment with respect to the vertical (Reference Martín and GudmundssonMartín and others, 2012). The opposing tilt in apex positions of modeled and observed stratigraphy can either be interpreted as unresolved 3-D effects or as a migration of the southern divide towards the west. The latter will be on a timescale much shorter than the characteristic time (∼900 years for Halvfarryggen), because the Raymond stack is not (yet) aligned with the (migrated) divide position. The localized offset between the highest point in surface topography and the apex of the isochrone arch near the dome (Figs 5 and 6) may support this hypothesis, as it also indicates a westerly migration of the dome (within the last 10–20 years). However, such a migration is also likely to affect the position of the divides away from the dome, below which the offset is not observed. Raymond stacks below the northeastern divide are both vertically aligned and tilted (Fig. 3), and the double-ridge features, which are visible in satellite imagery (Fig. 2b) and interpreted by Reference Martín, Gudmundsson, Pritchard and GagliardiniMartín and others (2009b) as an indicator for long-term stability, do not support a migration of the divides further away from the domes on a timescale larger than the characteristic time. Based on the available dataset, we therefore have no conclusive evidence to explain this observation, and 3-D, anisotropic modeling is needed to clarify the flow dynamics around triple junctions. The need for such a model is also illustrated by the spatial representation of the internal reflection horizon in the RES data (Fig. 10). It shows that the southern divide is quasi-2-D, meaning that the corresponding isochrone arch is little affected by ice flow in directions other than perpendicular to the divide. To a lesser extent this is also the case for the northeastern branch. The muting of the isochrone arch below the northwestern branch, however, highlights the importance of along-ridge flow, which is more pronounced along this divide due to the steeper surface slope.
Fig. 10. (a) Fence diagram of RES data near the summit. The green surface with shaded relief is interpolated from a picked horizon with an average depth of 260 m below the surface. The isochrone arch is visible as a shaded ridge. The blue dots mark the location of the ice divides. (b) RES profiles (blue and red) used for the interpolation. The red lines mark the RES profiles shown on the left. In the background is the merged DEM with 50 m contours.
The high accumulation rates make Halvfarryggen primarily a candidate for an ice-core drill site as defined within the IPICS 2k array. The east–west positioning of the drill location is important, since the accumulation gradient results in a gradient of the age–depth distribution, which is further amplified by the isochrone arches beneath the divides. We propose a drill site along the southern divide where the along-flow component is small. The scientific output in this area is threefold:
-
1. Due to the predicted anisotropy in COF, ice-core measurements have the potential to deliver evidence for the role of COF in ice deformation (especially of isochrone arches).
-
2. The corresponding climate record is expected to be highly resolved in the time frame needed to quantify anthropogenic influences. At the divide, the annual-layer thickness for ∼2000 year old ice ranges around 0.5 m.
-
3. At larger depths ice from the last termination is potentially conserved, also fulfilling requirements of the IPICS 40k array.
A caveat lies in the RES stratigraphy, which appears not to be continuous in the lowest 5–10% of the ice column throughout the area. Whether this is an issue of system sensitivity or potentially linked to stratigraphic disturbances (Reference Waddington, Bolzan and AlleyWaddington and others, 2001; Reference DrewsDrews, 2009b) has yet to be determined. The latter seems unlikely near the dome area, where there is little horizontal ice flow. The numerical model predicts warmer ice at the base beneath the divides where the isochrone arches develop. In these areas we also observe less backscattered power (Fig. 3), indicating that the loss of RES layering in the lowest part of the ice column is more likely linked to increased signal absorption. Future airborne RES profiles should investigate this further. Ground-based radar measurements should be directed to spatially better resolve the accumulation gradient, which at the same time will overcome the limitations of the current DEM with the help of additional GPS profiles. This will facilitate the precise identification of a drill site according to the yet to be determined scientific priorities.
5. Summary and Conclusions
We have derived a DEM and located the ice divides with a combination of different remote-sensing techniques. The internal layering was investigated by means of airborne RES and ground-based GPR surveys. Analysis of the top layers reveals a strong spatial variation in accumulation which appears to be driven by a preferred wind direction and changing surface slope. Starting 30–60 m below the surface, the layers beneath the dome arch upwards due to the Raymond effect. At the dome, we observe a small offset of the apices compared to the highest point in the surface topography. Raymond stacks around the triple junction are not symmetric, vary in amplitude and appear partly tilted. The modeled age–depth scale at the divide is robust with respect to the specific value of the geothermal heat flux and is also only little influenced by the magnitude of the rheological index. It compares well with analytical solutions in the flanks which have been traced to the divide area along isochrones in the radar data. However, errors in transferring the age from the flanks to the divide can be considerable (Fig. 8), and the technique generally fails at larger depths due to discontinuities in the radar stratigraphy. Based on the appearance of the double-peaked isochrone arches and the comparison with the anisotropic ice-flow model, we conclude that the divide has been stable for at least 2700–4500 years. The tilt of the Raymond stack potentially indicates a more recent, westerly migration of the southern divide, but the role of a temporally changing accumulation as well as the role of along-ridge flow is unclear and hampers a solid interpretation. The inclusion of an anisotropic rheology predicts (at the divide) a girdle-type COF in the upper half of the ice column, which evolves into a single-maximum distribution at larger depths. Physical properties derived from ice-core measurements in this area could confirm or belie this prediction and thus contribute to the further development of an anisotropic flow law. The derived age–depth scale favors Halvfarryggen as a drill site for studies which require a highly resolved climate record for up to 4 ka BP. At larger depths, ice from the transition into the last glacial period may be preserved.
Acknowledgements
Preparation of this work was supported by the Emmy Noether program of the Deutsche Forschungsgemeinschaft (DFG) grant EI 672/5-1 to O. Eisen and a scholarship of the Evangelisches Studienwerk e.V. Villigst to R. Drews. We thank V. Helm for the airborne laser scanner data, and C.H. Reijmer and M.R. van den Broeke for providing the data from AWS11. We thank M. Funk for providing details of the low-frequency radar antenna design. Initial ideas for comparison of data and models at Halvfarryggen were based on discussions between R. Hindmarsh and C. Martín during a stay at the British Antarctic Survey, funded by the Scientific Committee on Antarctic Research (SCAR) fellowship scheme 2006–07 and an Emmy Noether scholarship of DFG grant EI 672/1-1 to O. Eisen. Comments from E. Pettit, J. MacGregor and an anonymous reviewer greatly helped to improve the manuscript.
Appendix: Derivation of Surface Topography
The CW model was interpolated from data of airborne radar altimetry, complemented with spaceborne laser altimetry (ICESat), and ground-based GPS surveys. In particular, it incorporates the ground-based GPS data from 2007 gridded to 1 km × 1 km. The RD model was derived by differential synthetic aperture radar (SAR) interferometry. It has the benefit of a higher horizontal resolution (50 m × 50 m grid), at the cost of partly unknown atmospheric contributions and a strong dependency on the coherence of the processed image pairs. The BB model is based on photoclinometry with data from the Enhanced Thematic Mapper Plus (ETM+) instrument on board the Landsat-7 satellite. It is gridded to 30 m × 30 m. The primary pitfall of this method is albedo variation caused by factors other than changing surface slope.
In order to work with a single DEM, the individual DEMs were sampled to a common 100 m × 100 m grid and compared with independent ground-truth data. For RD and BB the GPS data from the surveys in 2007 and 2010 were used. In both GPS surveys, a local reference station was positioned in the center of the survey area to enable a differential post-processing with baselines smaller than 30 km. Details of the differential processing of the GPS data are provided by Reference Wesche, Riedel and SteinhageWesche and others (2009). The estimated vertical accuracy from crossover points within the individual surveys is in the sub-meter range. The 2010 GPS elevations have a 1.3 ± 0.2 m offset compared to the 2007 GPS elevations. This may be caused by a real change in surface elevation or by an imprecise static solution for the local reference station, which depends on longer baselines to fixed reference stations. Additionally, profiles from an airborne laser scanner, which are not incorporated in the CW model, are used for the DEM validation. The vertical accuracy of the laser scanner data is discussed by Reference HelmHelm and others (2007) and lies within the sub-meter range.
The ground control points from the GPS surveys and the airborne laser scanner enable the identification of weak areas (deviation >25 m) in the individual DEMs. For example, around the dome, the BB and CW model perform well, whereas the RD model deviates significantly (>40 m) from the laser scanner and GPS data. However, in the larger vicinity (>10 km) of the dome, the performance of the RD model improves. Based on the initial comparison with the control points, the individual models were masked in weak areas and then merged to a single DEM. To minimize mismatch at the margins, a 2-D cosine tapering function was applied at the boundaries and the RD model was additionally fitted to the BB model near the dome with a first-order polynomial function.
