Summit of the East Antarctic Ice Sheet underlain by thick ice-crystal fabric layers linked to glacial–interglacial environmental change Open Access

Ice sheets contain several types of single ice crystals, such as a perfect single pole, an elongated single pole and a vertical girdle crystal (Fig. 1). Ice crystals experience an anisotropic response to stress, whereby their principal axes, termed C-axes, rotate during deformation towards the direction of compression and away from the tensional direction (Fig. 2; Alley 1988; Azuma et al. 1999; Durand et al. 2009; Montagnat et al. 2014). Collectively, ice-crystal alignment leads to macro-crystallographic orientations, or ice-sheet fabric, which is the product of both ice rheology and stress. As the strain rate in anisotropic ice under simple shear (Azuma 1994; Durand et al. 2007; Kennedy et al. 2013) can be up to an order of magnitude greater than in isotropic ice (Azuma 1994), a feedback is possible whereby ice fabrics lead to preferential deformation and, in turn, further anisotropy (Alley 1988; Azuma 1994). The spatial variation of ice rheology is fundamental to ice-sheet flow modelling, yet it is normally averaged through depth as an isotropic medium, as in Glen’s flow ‘law’ (Glen 1954). Differences in the velocity of cold-based ice between observations (Thorsteinsson et al. 1999) and modelling (Gillet-Chaulet et al. 2006; Martin & Gudmundsson 2012), in which respective ice-surface elevations can vary by hundreds of metres, and ice velocities by tens of metres per year, may be attributable to the rheological impacts of ice fabric development. Crystal alignment in ice cores is analysed routinely, revealing marked and step-wise changes in the fabric of ice deposited across glacial–interglacial periods (Durand et al. 2007; Montagnat et al. 2014). Additionally, sonic logging (Gusmeroli et al. 2012) shows how variations in ice-sheet fabric is correlated highly with oxygen isotope ratios, implying a link with atmospheric parameters (Paterson 1991; Kennedy et al. 2013). The spatial distribution of ice-sheet fabric is largely unknown, however, as ice cores are heavily restricted in number and distribution across Antarctica and Greenland. As a consequence, the crystallographic structure of ice in large ice sheets, including its development through time and space, and its potential rheological consequences, remains uncertain.

Dome A is the highest ice divide on the East Antarctic Ice Sheet (Jouzel & Masson-Delmotte 2010), with surface elevations reaching over 4100 m above sea-level. During the 2007–09 period of the International Polar Year the Dome A summit region was surveyed by an ice-penetrating polarimetric radar (Sun et al. 2009; Fig. 3). Initial inspection of the non-polarized radar data revealed the first detailed images of the Gamburtsev Subglacial Mountains, which were found to contain well-preserved upland glacial geomorphological features and a clearly defined water drainage pattern originating from fluvial activity in preglacial times (Sun et al. 2009). Numerical ice-flow modelling suggests these mountains were the nucleation site of the East Antarctic ice sheet, and are likely therefore to have been continuously glaciated for at least 14 Ma (DeConto & Pollard 2003).

Because radio-waves reflect off ice fabrics as a function of orientation, polarimetric radars, which measure radio-wave reflections simultaneously in orthogonal axes, are able to detect gross levels of ice-crystal orientations. Here, we analyse polarimetric data from the Dome A radar survey, revealing insights into the deep-ice crystal structure at and around the ice divide, and its relation with subglacial topography. By linking the data stratigraphically to the Vostok ice core, using airborne radar data, we also provide an assessment of the temporal evolution of ice fabric.

In the case of a random ice-crystal fabric, the response of the electric field within the medium is independent of its orientation and ε₁ = ε₂ = ε₃. If ε₁ = ε₂ = ε_┴ and ε₃ = ε_//, an idealized perfect single-pole behaviour occurs (Fig. 1b). For a single-pole ice-fabric distribution, as the ice crystals are aligned in planes, distinct electrical properties will exist along and across those planes. If we denote ε_┴ and ε_// the permittivity perpendicular and parallel to the normal vector of the crystal planes, respectively, the permittivity tensor can be expressed as ε₁ = ε₂ = ε_┴ and ε₃ = ε_// (Fig. 1b).

In the case of a vertical girdle fabric, the crystal planes are randomly oriented in a vertical plane. In this situation, the two permittivities are (1) the average between ε_┴ and ε_// in the vertical plane where fabric is random, and (2) equal to ε_// in the normal direction, as follows: ε₂ = ε₃ = (ε_// + ε_┴)/2 and ε₁ = ε_// (Fig. 1d).

If the ice crystals have elongated single poles then ε₁ < ε₂ < ε₃ (Fig. 1c). The ellipse ratio in this fabric corresponds to the anisotropy of permittivity in the horizontal direction. The larger this ratio is, the higher the level of anisotropic permittivity.

Polarimetric radars have orthogonal transmitting and receiving antennas, receiving four series of polarized echo (noted as HH, HV, VH and VV). Relative echo powers from these configurations vary through anisotropic media. The classic Stokes vector parameters (referred to below as S₀, S₁, S₂ and S₃) can be used to evaluate and quantify the polarization of such radio-waves as follows.

Here, $E→$ represent the electrical field vector. Its superscript ‘i’ means that it is the incident wave and ‘r’ indicates the reflection wave. $R¯$ is a 2 × 2 reflection coefficient matrix describing the relationship between the incident and reflecting waves.

For inducing reflections measured in the V antenna we obtain similar expressions.

S₀ refers to the total echo-power from the orthogonal antennae, and resembles the reflected power using in the conventional sole-channel radar system. S₁, S₂ and S₃ provide the difference in reflected power between the orthogonal electric fields in the standard Cartesian arrangement, a Cartesian system rotated by 45° and in a circular coordinate arrangement, respectively, thus providing measures of the anisotropy in ice permittivity.

Several authors have confirmed the relationship between the echo difference in orthogonal directions and the permittivity tensor of ice (Hargreaves 1977; Fujita et al. 2006; Matsuoka et al. 2012), which is linked to the second-order tensor of ice fabric (Matsuoka et al. 2009). Hence, in the absence of plausible alternative explanations, it can be assumed that the Stokes vector parameter description of the polarimetric radar reflections S₁, S₂ and S₃ expresses ice-sheet anisotropy caused by ice crystal fabric.

The polarimetric radar used a frequency of 179 MHz and recorded radio waves in orthogonal orientations over a square grid of fourteen 30 km-long transects with along-track sample acquisition every c. 17 m (Fig. 3). The radar employed dual receivers in orthogonal orientations (H and V) with a time-sharing sole transmitter (H or V). It comprised four assembly antenna pairs: HH, HV, VH and VV. An in-phase/quadrature (I/Q) demodulator can measure the amplitude, phase and instantaneous frequency of radar signals, so it is widely used in polarimetric radar systems. Normally each antenna pair records both I and Q components for detecting phase. Therefore there are eight data records in the polarimetric radar files necessary to measure information of both the amplitude and phase of reflected waves. Further details of the radar are provided in Sun et al. (2009).

The I and Q components of the four antenna pairs were imported from raw data files, from which the Stokes vectors of both H and V transmitters were calculated. The four parameters in the Stokes vector each have an explicit physical meaning. S₀ represents the average power of the EM waves. S₁ provides the echo difference between H and V components. Its value can be positive, zero or negative depending on whether linear polarization of H or V is dominating. S₂ provides the echo difference of orthogonal components between 45 and −45°. S₃ represents the echo difference of orthogonal components within a circular coordinate system. After normalization, S₁, S₂ and S₃ all range between −1 and +1.

S₀ has similar results to traditional (non-polarimetric) radar reflections. It can display clear isochrones and subglacial topography. In contrast S₁, S₂ and S₃ display obvious bands within the normal isochrone arrangements (e.g. red (positive) and blue (negative) layers in Fig. 4). The units in S₁ and S₂ are more distinct than in S₃ and show little difference between the H and V transmitters. By adding S₀, S₁ and the average power from H and V together, the positive layers link to reflections recorded in the V component, and the negative layers link to the H component. If the coordinates are fixed (i.e. as in a straight survey line), variation of the Stokes vectors in the H and V components within different layers of ice must be attributed to differences in the ice crystal fabric (i.e. different layers have different fabric and orientation tensors).

The Stokes vectors were calculated in each of the transects (Fig. 4), revealing four obvious englacial layered units in which discrete variations in S₁ were noticed, named T₁, T₂, T₃ and T₄, respectively (T₀ designates the bed). Beneath T₄, further units of alternating S₁ occur, but they are less obvious and not as easy to trace owing to the lack of returned power from deeper englacial reflectors. Broadly, positive S₁ units (red–yellow zones in Fig. 4) correspond to high levels of anisotropy in the horizontal plane (as the radar is insensitive to the anisotropy in the direction of transmission), and negative S₁ units (blue zones) are derived from either isotropic ice or low levels of anisotropy in the horizontal plane. As expected, the deeper S₁ ice fabric units thin gradually with depth owing to vertical compression (Fig. 4). Unit T₁ is the thickest of those measured, from the surface to ice depths of c. 800–1000 m. Based on crystallographic information from deep ice cores, this upper unit will have a nearly random (isotropic) fabric (Azuma 1994; Azuma et al. 1999; Fujita et al. 2006; Durand et al. 2009; Montagnat et al. 2014). Unit T₂ (positive S₁) is up to c. 400 m thick and contains an abundance of highly anisotropic ice related to a vertical girdle fabric distribution. T₃ (negative S₁) is a product of low levels of horizontal anisotropy, as the fabric clusters towards a broad single-pole distribution (Hargreaves 1977; Eisen et al. 2007; Matsuoka et al. 2009). T₄ (highly positive S₁) is similar in crystal arrangement to T₂. While both the thickness and depth of the units fluctuate horizontally substantially, we estimate that c. 20% of the ice column beneath Dome A, and up to c. 50% of the ice column at depths >1 km around this ice dome, are characterized by high levels of anisotropy (primarily from T₂, T₄ and T₆; Fig. 4).

We are able to link radar isochrones from Kunlun Station to the Vostok ice core, allowing us to estimate a depth–age relationship at Dome A (Figs 5 & 6). The depth of the ice fabric units at Kunlun is then used to assign estimated ages to the fabric units (Table 2). We also used the stratigraphic connection to Vostok to translate the ice core dust and ΔT records to Dome A (Fig. 7).

Isochrone layers were picked along two airborne radar transects, collected by the Alfred Wegener Institute. The radar used operates at 150 MHz, alternating between 600 and 60 ns bursts, with internal stacking of shots prior to storage and a pulse repetition frequency of 20 Hz. The effective distance between 60 ns shots is c. 8 m, and for interpretation the shots are stacked two-fold (so that c. 8 m becomes c. 16 m). The two-way travel times of the layers were converted into ice depth using the Kovacs formula (Kovacs et al. 1995) and ice density information (Barkov 1970; Ekaykin et al. 2010), with an accuracy of c. 5%. The radar data were differentiated and amplitude balanced to aid the extraction of the layers. The age–depth relationship at Vostok Station was taken from Parrenin et al. (2004). The layers were traced using a commercial software package from Landmark, which displays picked horizons at cross-over points and thus allows an evaluation the quality of the data to ensure the correct isochrone is identified along the transect. Importantly, layers were traced along both transects from Dome A and compared at their intersection at Vostok Station. Agreement between the layers at Vostok Station provides confidence that the correct layers were followed to Dome A, and that the depth–age relationship there is reliable. The error for picking both the ice-surface reflection and an isochrone is about five sample points, which corresponds to 66.65 ns or c. 5.6 m in ice thickness.

From the stratigraphic correlation between Vostok and Dome A, we estimate that the base of T₂ is c. 105 kyr, which roughly corresponds with the end of the Eemian interglacial period. T₃ is aged between c. 105 and c. 136 kyr, and there is a good correspondence between the peak of S₁ within T₃ and the peak of the Eemian. T₄ is between c. 136 and c. 190 kyr, similar to MIS 6. Two faint units below T₄, named T₅ and T₆, are possibly consistent with previous glacial/interglacial periods.

Broadly, positive S₁ units (T₂, T₄) correspond well to glacial periods and negative S₁ units (T₃, T₅) are associated with interglacials (Fig. 7). This suggests that gross crystal fabric across Dome A is influenced by past environmental change at the glacial–interglacial transition.

In a vertical column of ice near the ice divide, the orientation of ice crystals will gradually tend to a predominant (vertical) orientation under compression. In doing so, the ice fabric develops from ‘random’ to ‘single pole’. Horizontally, perfect single pole and ‘vertical girdle’ crystal alignments represent the two end-member extreme fabrics (Figs 1 & 2a) that could be expected in idealized ice-divide and ice-sheet margin regions, respectively. The resulting development of ice fabric in layers T₁–T₅ is illustrated in Figure 8. An alternative, theoretical, explanation is that under shear the crystal axis experiences rotation in line with the stress field (Fig. 2b). Both explanations, summarized in Table 3, would yield similar polarimetric radar measurements.

It is known that soft layers (glacial ice) are more easily deformed under stress than hard layers (interglacial ice) and isotropic ice (e.g. Paterson 1991; Thorsteinsson et al. 1999; Durand et al. 2007; Matsuoka et al. 2012; Kennedy et al. 2013). The deformation leads to the glacial ice layers more easily experiencing crystal realignment (Fig. 2a), which is what is recorded by the Stokes vectors measurements (Fujita et al. 2006; Wang et al. 2008; Matsuoka et al. 2012). Hence, the polarimetric radar measures the thickness and distribution of layers with discrete ice fabrics.

The S₁ units align well among the individual transverses and cross profiles. We picked the top and bottom surfaces of all the S₁ units in each line (where S₁ became 0) to make three-dimensional wall and surface diagrams (Fig. 9), which display the extent and thickness of each fabric horizon revealing two obvious features. The first is significant S₁ variations in both ice depth and unit thickness. The second is changes to ice fabric caused by flow over macro-scale bed topography (Robin & Millar 1982). The fabric units are dominant across the left side of a major subglacial valley. On the right side of the subglacial valley, a cluster of relatively small peaks exists. Unit T₂ pinches out immediately above these peaks (Fig. 9). In this area, an arcuate belt with no S₁ signal covers the cluster of peaks. We infer that, above these smaller peaks, a pure compressive stress state (uniaxial or pure shear) dominates, causing the fabric to be characterized by a single-pole arrangement (Azuma 1994; Matsuoka et al. 2012). This induces isotropy in the horizontal direction, causing the S₁ to be near zero (c. 25 km mark in Fig. 4d), leading to local disappearance of the fabric units. Away from the peaks, simple shear will increase gradually, leading to the redevelopment of fabric units that can be seen on both sides (Fig. 9d).

Our polarimetric radar data reveal widespread ice-sheet anisotropy at Dome A. However, polarimetric radar only captures relative contrasts in electric permittivity between two orthogonal directions in the horizontal plane rather than the actual values. Interpretation of our radar data has been guided by typical ice-fabric distributions sampled in ice cores sites (e.g. Durand et al. 2007 at Dome C; Drews et al. 2012 at EPICA DML; and Montagnat et al. 2014 at NEEM). Ice-core data only show subtle changes in fabric during climatic transitions, however, in contrast with the marked transitions we image in radar data. This suggests that widespread but locally subtle transitions in ice fabric are likely to have a large effect on the permittivity of the ice sheet. Hence subtle transitions in ice fabric could potentially have a larger effect on ice flow dynamics than currently considered.

Our observations at Dome A are likely to hold for other ice divides, including central Greenland, as noted above. At Dome F in East Antarctica, using a dual frequency radar, simple-shear-dominated ice flow over large subglacial obstacles has been shown to lead to preferential recrystallization along bands of relatively soft ice (Fujita et al. 1999), meaning that the rheological consequences of initial anisotropy will evolve downstream. At Dome C, Siegert & Fujita (2001) analysed 60 MHz radar, and discovered that the form of the depth-to-power relationship of the returned radar signal could be used as a proxy for the along-track development of ice-crystal fabric. A similar investigation revealed how ice fabric may evolve between Ridge B and Lake Vostok (Siegert & Kwok 2000).

One situation where the downstream modification to ice crystal fabric and rheology may be highly relevant is in the preservation of major blocks of basal ice, interpreted to be caused by widespread and thick basal accretion (Bell et al. 2011). Stress-induced anisotropy-may lead to preferential englacial flow along bands, leading to smooth englacial flow over rough topography and, potentially, a partial flow ‘decoupling’ of the bulk of the ice sheet from basal ice units frozen to the bed. This new and alternative explanation for the so-called ‘accretion ice’ identified atop the Gamburtsev Subglacial Mountains (and other ‘fuzzy’ lower-englacial features seen in both Greenland and Antarctica) is consistent with the lack of significant surface uplift and velocity anomalies that a high rate of basal freezing (of the order of the surface rate of ice accumulation) should lead to, and is consistent with meteoric radar layering above the basal ice that drapes smoothly over it (see Wrona et al. 2017). Such flow may preserve the basal ice units as ancient ice blocks, which means that the rate of ice accumulation by freezing will be lower, and the age of ice considerably older (possibly by millions of years), than first thought.

Polarimetric radar measurements, linked stratigraphically to the Vostok ice core, demonstrate that ice-sheet crystal alignment at the centre of large ice sheets is a function of both climate-induced impurities acting to ‘soften’ the ice during periods of glaciation and internal stresses within the ice sheet, acting to preferentially deform the soft ice, yielding ice-sheet anisotropy. These observations at and around Dome A are likely to hold for other ice divides, including central Greenland, where soft anisotropic ice has been measured in the lower portion of the NEEM ice core (NEEM Community Members 2013).

The next steps to this area of research are two-fold: first, to study ice crystal fabric developments using polarimetric radar around other ice divides, to understand the similarities to, and differences from, the situation at Dome A; and second, to investigate the degree to which the ice fabric uncovered at Dome A relates to modification in ice rheology. Ice-sheet modelling is capable of such analysis, through assimilation of three-dimensional internal layer structures into the ice flow parameter. This coupling between data and modelling is possible, but requires high-resolution radar data acquisition across study regions.

This paper is based on the work supported by the National Natural Science Foundation of China (40874060, 41176167, 41276200 and 41676181), the 973 projects of China (2012CB957702), the 863 projects of China (2011AA040202) and the UK Natural Environment Research Council (NE/D003733/1). This work was carried out while Bangbing Wang was a visiting research fellow at the British Antarctic Survey. We thank Richard Hindmarsh, Hugh Corr and Ed King for help and advice in developing the work for this paper. We also thank Neil Ross and Adam Booth for offering constructive and helpful reviews.