1. Introduction

Efforts are underway to illuminate the ice sheets with a polari-metric airborne system as a precursor to a potential space- borne mission (Reference DallDall, 2010). It is important to understand the polarization dependence of the backscattered signal and to successfully link it to physical properties within the ice. This delivers input for ice-sheet modeling, that enhances both the correct interpretation of paleo-ice-core records and the predictability of the ice sheets. We analyze polarimetric radio-echo sounding (RES) data from the Antarctic plateau, collected near the German summer station Kohnen in east Dronning Maud Land (DML), Antarctica. The data were acquired with an airplane on the ground, measuring the effect of varying polarization as the airplane followed a circular profile. The backscattered power is strongly related to the antenna orientation (i.e. incident polarization), and backscatter patterns change with depth. We use physical properties gained from a nearby deep ice core to discuss mechanisms for anisotropic scattering. Particularly, we focus on ellipsoidal air bubbles and anisotropic crystal-orientation fabric (COF). The ice core was drilled within the European Project for Ice Coring in Antarctica (EPICA) and is referred to as EDML. The drill site is located near a forking ice-divide in a flank-flow setting. An overview of the scientific program at Kohnen station and major results from the EDML ice core are given by Reference Oerter, Dmcker, Kipfstuhl and WilhelmsOerter and others (2009).

The polarization of the emitted radar signal is parallel to the airplane heading. In the upper 900 m the radar backscatter is strongest when the polarization is roughly perpendicular to the nearby ice divide. The orientation of maximum backscatter changes in direction by 90° at ˜900m depth, coinciding with the bubble/air-hydrate clathrate transition (under increasing pressure the air bubbles turn into a second solid phase of crystalline air hydrates, also called clathrates) and the transition from Holocene ice to ice from the last glacial period. The polarization dependence exhibits primarily a 180° symmetry. The aim of this study is to pinpoint a mechanism causing the anisotropy observed in the RES data, and to elucidate the capability of polarimetric RES surveys for the profiling of stress/strain regimes along the vertical.

RES measurements over ice sheets record reflections from dielectric non-uniformities. Density variations are dominant at shallow depths and alter the real and imaginary part of the complex dielectric constant (ε = ε’ — iε’’, where i² = -1). Changes in electrical conductivity influence the imaginary part and can be separated from density changes at intermediate depths and below, where density variations are negligible. Both effects often produce laterally coherent internal reflection horizons in the radar data. The dependence of the backscattered signal on the polarization of the incoming wave was noted in the early days of radio-glaciology. Reference DoakeDoake (1981) provides an overview. The backscattered signal appeared elliptically polarized, and showed a varying amplitude, depending on the respective azimuth orientation of receiving and transmitting antennas. Using solutions of the Maxwell equations in layered media with a tensorial form for the dielectric properties (Reference HargreavesHargreaves, 1977, Reference Hargreaves1978; Reference DoakeDoake, 1981), this could be explained in terms of birefringence and anisotropic reflection coefficients. Reference HargreavesHargreaves (1978) proposed as a mechanism for the birefringence the dielectric anisotropy of the single ice crystal in combination with preferred COF. The dielectric anisotropy has been confirmed in subsequent studies (Reference Fujita and MaeFujita and others, 1993, Reference Fujita, Matsuoka, Ishida, Matsuoka and Mae2000; Reference Matsuoka, Fujita, Morishima and MaeMatsuoka and others, 1997), which found a difference in the real part of the dielectric constant of for the components parallel (||) and perpendicular (⊥) to the crystal’s c-axis. The relative permittivity ranges between 3.1 and 3.2, with small dependences on temperature and frequency.

Numerous multi-frequency and multi-polarization measurements have been performed since then to infer various COF types from RES measurements (Reference Fujita and MaeFujita and Mae, 1993; Reference Liu, Bentley and LordLiu and others, 1994; Reference FujitaFujita and others, 1999, Reference Fujita, Matsuoka, Maeno and Furukawa2003, Reference Fujita, Maeno and Matsuoka2006; Reference Siegert and KwokSiegert and Kwok, 2000; Reference Doake, Corr and JenkinsDoake and others, 2002; Reference MatsuokaMatsuoka and others, 2003, Reference Matsuoka, Uratsuka, Fujita and Nishio2004; Reference Eisen, Hamann, Kipfstuhl, Steinhage and WilhelmsEisen and others, 2007; Reference Wang, Tian, Cui and ZhangWang and others, 2008; Reference DallDall, 2010; Reference Kravchenko, Besson, Ramos and RemmersKravchenko and others, 2011). Variations in density, electrical conductivity and changing COF are now accepted to be the major mechanisms producing RES horizons. Nevertheless, predominantly due to a lack of ice-core data, uncertainties about the nature of the anisotropic reflection mechanisms remain. In particular, effects from volume scattering have so far been largely neglected. Independent of the specific mechanism for the polarization dependence, the theory of birefringence and anisotropic reflections is well understood.

In birefringent media, waves are refracted depending on the incoming polarization. The solution of the corresponding wave equation results in two characteristic waves, the ordinary wave and the extraordinary wave, which are perpendicularly polarized and generally differ in wave speed. The ordinary wave is refracted according to Snell’s law, which does not hold for the extraordinary wave. In uniaxial media (i.e. two out of three principal components of the dielectric tensor are equal) an ordinary wave and an extraordinary wave superimpose. In biaxial media (i.e. all three principal components of the dielectric tensor differ) two extraordinary waves are generated. Typical polarimetric experiments on ice sheets emit an electromagnetic pulse and record amplitude variations of the reflected signal as a function of the azimuth antenna orientation. We discuss only the co-polarized case, which means that the emitting and receiving antennas remain parallel. As described by Reference HargreavesHargreaves (1977), the azimuth dependency of the backscattered power is a function of both anisotropic scattering and birefringence, which have different symmetries. In birefringence, the amplitude variation is caused by the superposition of the two characteristic waves with a phase shift. The phase shift depends on the depth-integrated difference in dielectric components perpendicular to the propagation direction. The latter causes a 90° periodicity, since the magnitude of the difference remains the same for a 90° rotation of the coordinate system. Anisotropic scattering, however, has a periodicity of 180°, since the corresponding rotation of the antennas at the surface results in an indistinguishable setup. Apart from some special cases (e.g. a fully isotropic COF or a perfect single-maximum COF aligned with the direction of propagation), polar ice is generally birefringent. Anisotropic scattering requires anisotropic gradients in the dielectric properties, for example due to a vertically varying COF, aligned and non-spherical inclusions, or a preferred roughness direction at the interfaces of internal layers. Birefringence and anisotropic scattering are related, and when combined lead to a complex amplitude variation with changing antenna orientation. This can be simulated for the different cases with matrix-based models (Reference DoakeDoake, 1981; Reference Fujita, Maeno and MatsuokaFujita and others, 2006; Reference Matsuoka, Wilen, Hurley and RaymondMatsuoka and others, 2009). However, sometimes either a 90 ^° or a 180 ^° periodicity in backscattered power is dominant; this allows workers to focus solely either on birefringence or on anisotropic scattering, respectively.

1.1. Site characteristics and previous studies

We consider a polarimetric dataset that was acquired along a circle (which changes the antenna’s azimuth orientation continuously) with an RES system operating in a co-polarized mode. In the following, we refer to this dataset as the circular RES profile. The location of the measurements is in Dronning Maud Land, ˜500 m downstream of the deep-drilling site for the EDML ice core (Fig. 1). The topography is smooth, with gradients of a few meters per kilometer. The area contains three main ice divides, which merge in a triple junction. The drill site is close to one arm of the forking divides, and the entire setting is in a flank-flow regime. Ice thickness at the drill site is 2774 m and local surface velocities average <1 ma^-1 (Reference Wesche, Eisen, Oerter, Schulte and SteinhageWesche and others, 2007).

Fig. 1. Map of surface elevation, ice flow (Reference Wesche, Eisen, Oerter, Schulte and SteinhageWesche and others, 2007; black arrows), and bedrock topography (Reference Steinhage, Nixdorf, Meyer and MillerSteinhage and others, 2001). The EDML drill site is at the origin of the light-gray arrows, indicating the principal directions of compression and dilatation.

The bedrock topography was mapped with a dense grid of RES profiles (Reference Steinhage, Nixdorf, Meyer and MillerSteinhage and others, 1999, Reference Steinhage, Nixdorf, Meyer and Miller2001) during pre-site surveys. After the ice core was drilled, Reference Eisen, Wilhelms, Steinhage and SchwanderEisen and others (2006) used data from dielectric profiling and a forward model to link several internal reflectors in the RES data to prominent changes in electrical conductivity. Two of those reflectors (referred to as h1 and h2) are discussed below. Using the same circular RES profile as presented here, a subsequent study (Reference Eisen, Hamann, Kipfstuhl, Steinhage and WilhelmsEisen and others, 2007) showed that another localized reflector at 2035 m depth is caused by a drastic change in COF which is observed in ice-core data. Linescan images from the ice core show the stratigraphy on a sub-centimeter resolution. The stratigraphy appears increasingly disturbed with greater depths. This suggests that the disappearance of internal reflection horizons below ˜2100 m depth may be caused by stratigraphic disturbances on a larger scale (Reference Waddington, Bolzan and AlleyWaddington and others, 2001; Reference DrewsDrews and others, 2009).

Reference Wesche, Eisen, Oerter, Schulte and SteinhageWesche and others (2007) derived an elevation model as well as surface velocities with corresponding stress/strain rates, based on a pentagon-shaped stake network. The flow regime is characterized by along-flow compression and lateral dilatation. The maximal principal strain-rate component has a magnitude of -1.85 × 10^-4a^-1 and points towards 24 ^° N. The minimal principal component points towards 114° N with a magnitude of 2.32 × 10^-4 a^-1. Reference UeltzhöfferUeltzhöffer and others (2010) investigated size, number and shape of the air bubbles in vertical cuts from the ice core. The transition of air bubbles into clathrates starts at ˜700m and progresses rapidly below 800 m. Typical effective bubble radii cluster around 0.1 mm, decreasing with depth, similar to the volume fraction, which ranges from ˜0.004 (200 m depth) to ˜0.001 (900 m depth). In vertical cuts, it was found that the bubbles are elongated, with eccentricities varying from 0.3 to 0.8. Due to the lack of continuous azimuth control of the ice core, no correlation with ice flow is possible. Reference BendelBendel (2009) did some further analysis on vertical cuts and found that number density and aspect ratio vary on centimeter scales, in particular for glacial ice, but also at shallower depths.

Crystal-orientation fabric studies on horizontal and vertical thin sections exploit the birefringent nature of the single crystals by using crossed polarizers to determine the c-axis orientation. The orientation distribution is often visualized using Schmidt diagrams, which project the piercing point of the c-axes through a hemisphere onto a plane. By fitting an ellipsoid that has the same moment of inertia as the distribution of piercing points (with imaginary masses attached), the COF types can be characterized with three eigenvalues, λ₁, λ₂, λ3, and their corresponding eigenvectors. The eigenvalues define the ellipsoid’s shape, the eigenvectors its orientation. As explained below, it can be assumed that the first two eigenvectors are in the horizontal plane, with the third eigenvector pointing vertically. Two specific COF types are the single-maximum and the vertical- girdle distribution. In a perfect single-maximum distribution, all c-axes point vertically (λ₁ = λ₂ = 0, λ₃ = 1). In a vertical-girdle COF (λ₁ < λ₂ < λ₃), the c-axes cluster in a deformed cone around the vertical axis. More details of the fabric measurements and the results are presented below.

2. Data Acquisition and Observation

The radar data were acquired in 2003 with an RES system designed for the airborne sounding of ice sheets (Reference NixdorfNixdorf and others, 1999). The short backfire antennas emit bursts with a pulse length of 60 or 600 ns at a center frequency of 150MHz. The electric field vector oscillates parallel to the airplane’s heading. For this study, the short pulse was used to increase the vertical resolution, at the cost of total penetration depth. The antennas are fixed on the wings. To investigate the response with varying incoming polarization, the airplane (a Dornier 228) taxied on the ground in a circle. In this co-polarized mode, the emitted and recorded polarization planes were always parallel and rotated continuously around the vertical while the airplane moved along the circle. Travel- time-to-depth conversion was done with the interpolation of tie points from ice-core dielectric profiling data with RES horizons (Reference Eisen, Wilhelms, Steinhage and SchwanderEisen and others, 2006). The angle heading is defined as the heading angle of the aircraft in respect to true north. It depicts the angle of the electric field vector.

The backscattered power along the circular profile is displayed for the depth interval of 200-1400 m in Figure 2a. The diameter of the circle is ˜100m and the profile consists of 170 traces with an average trace spacing of 1.8 m. Each trace is a tenfold stack of shots and no further processing has been applied. The dominant feature is the overlying sinusoidal variation with changing azimuth of the two antennas. The signal has a clear 180° periodicity and splits in two parts along the vertical: from ˜200 to 900m depth, the maxima develop at a heading angle of ˜206° N/26° N (note that the heading is tangential to the circle). Below 900 m, the maxima shift by 90°, and occur at a heading of ˜116° N/296° N. We define the two depth intervals as anisotropic reflections zones (ARZ) 1 and 2, respectively. The different directions are shown in Figure 2b. The strain ellipsoid marks the directions of compression and dilatation given by Reference Wesche, Eisen, Oerter, Schulte and SteinhageWesche and others (2007). The directions of minimal and maximal backscatter correspond well with the principal axes of the strain ellipsoid.

Fig. 2. (a) The backscattered power of the circular profile varies sinusoidally. The two horizontal lines marked h1 and h2 are internal reflection horizons linked to prominent peaks in dielectric profiling data by Reference Eisen, Wilhelms, Steinhage and SchwanderEisen and others (2006). (b) Top view of (a) marking the shot point coordinates of the circular RES profile (black dots). The dark-gray arrows depict the main flow direction of ice together with main wind direction; strong wind events turn to more northerly directions. The light-gray arrows mark the extrema in backscatter on the circle for the different zones (numbers are trace numbers of (a)). The electric field vector is tangential to the circle. The strain ellipsoid (dark ellipse) indicates the directions for the principal components of compression and dilatation, following Reference Wesche, Eisen, Oerter, Schulte and SteinhageWesche and others (2007). These directions correspond well with the backscatter extrema observed in the RES profile.

3. Mechanisms for Polarization-Dependent Backscatter

As described above, the polarization dependence of the backscattered radar signal can principally be induced by birefringence and by anisotropic scattering. Birefringence is a direct consequence of the ice crystal’s dielectric anisotropy. Anisotropic scattering requires anisotropic changes in the dielectric properties. The two effects can be distinguished on the basis of symmetry.

To better visualize the symmetry, vertically averaged intervals of 200 m thickness are shown in Figure 3. Due to absorption and spherical spreading, the averages within the individual intervals are biased towards shallower values, but this does not affect the symmetry. The 180 ^° periodicity in backscattered power appears overlaid by a 90 ^° pattern resulting in two small local maxima between the main maxima mentioned above. From these two symmetries we deduce that both birefringence and anisotropic scattering influence the azimuth dependency of the backscattered power. However, anisotropic scattering is dominant and causes the larger amplitudes of the 180° cycle. The interweaved smaller maxima have a 90° offset, so both effects appear to have similar principal axes, that do not change noticeably with increasing depth. Due to the dominance of the 180 ^° periodicity, we focus on potential anisotropic scattering mechanisms.

Fig. 3. Vertical averages of backscattered power (dBm, i.e. referenced to 1mW) in 200m intervals from 200 to 1400m depth plotted against heading (bottom axis) and trace number (top axis). The main maxima with the 180° periodicity are partly interweaved with small maxima with a 90° offset.

The appearance of two anisotropic reflection zones has also been observed in other studies. Reference Fujita and MaeFujita and Mae (1993) compared data from a 179 MHz ground-based radar with fabric measurements from the ice core at Mizuho station, Antarctica. They observed high amplitudes in the upper part of the ice column with the antennas parallel to the flowline, and lower amplitudes for the perpendicular arrangement. Below 800 m depth, a similar pattern appears, but in the reversed order. The Mizuho ice core does not reach the bedrock and ends at ˜700m. This depth is above both the bubble/clathrate transition (Reference Shimada and HondohShimada and Hondoh, 2004) and the transition from the Holocene into the last glacial period (Reference Nakawo, Ohmae, Nishio and KamedaNakawo and others, 1989) at this site. Using different frequencies, pulse lengths and antenna orientations in the same study area, Reference Fujita, Matsuoka, Maeno and FurukawaFujita and others (2003) confirmed the previous findings of changing anisotropy with depth. They proposed as mechanisms for the upper anisotropy, a varying cluster strength for the girdle-type fabric, and for the lower anisotropy layered strata composed of girdle and single-pole transitions. Using similar techniques Reference MatsuokaMatsuoka and others (2003) also observed the two different zones along a ground-based traverse from Dome Fuji, Antarctica, towards the Shirase Glacier drainage basin. For the lower anisotropy they supported the idea of transitions from girdle to single-pole COF as the primary mechanism. An airborne survey (Reference Matsuoka, Uratsuka, Fujita and NishioMatsuoka and others, 2004) identified high scattering in the deeper layers of convergent and compressional ice- flow regimes when the polarization was perpendicular to the flow direction. Using a matrix model, Reference Fujita, Maeno and MatsuokaFujita and others (2006) discriminated between birefringence and anisotropic scattering at Mizuho and Dome Fuji. In both cases aligned COF and changes thereof were assumed to be the main mechanism for the observed polarization dependence. Based on the previous studies, the scenario for the two anisotropic reflection zones near Mizuho station is:

• One of the dielectric tensor’s principal directions is close to vertical. A varying cluster strength of a girdle-type COF translates into a variation of the first eigenvalue in the horizontal plane. This increases the backscatter for an incoming polarization parallel to the corresponding principal axes. There are no significant changes in the second eigenvalue in the horizontal plane.

• At larger depths, girdle to single-pole transitions change the vertical, and the second horizontal principal component. The latter increases backscatter perpendicular to the direction observed in the upper interval.

If only one eigenvalue varies in the horizontal, a rotation of the principal axis with increasing depth may result in a similar pattern. The problem with the hypotheses is that only a single study (Reference Eisen, Hamann, Kipfstuhl, Steinhage and WilhelmsEisen and others, 2007) directly linked (an exceptional large) short-scale variation in COF from ice- core data to a reflector in RES data in the lower part of the ice sheet. For the middle to upper depths, vertically higher resolved COF data in conjunction with polarimetric RES measurements are needed to confirm this idea. Therefore, before discussing the effect of variable COF at the EDML drill site, we also consider other possibilities, such as anisotropic inclusions or layer boundaries with a directional roughness.

It can be speculated that patterns that arise on the surface from preferred wind directions are advected to larger depths while preserving their directionality, thus causing anisotropic interfaces between different layers. Reference BirnbaumBirnbaum and others (2010) investigated the formation of snow dunes in the EDML area. They found that strong wind events create aligned wind-shaped features, which can also be seen in the physical properties of shallow firn cores (down to 5 m depth). The main (katabatic) wind direction is depicted in Figure 2b, but strong wind events have a more northerly direction (Reference Reijmer and Van den BroekeReijmer and Van den Broeke, 2003) and may align with the backscatter extrema in the radar data. However, deeper ice seems to have no memory of structures caused by deposition at the surface. First analyses of the density and microstructure of firn (Horhold and others, 2012) indicate that impurities significantly influence, or maybe even control, density and to a large extent the pore-space structure in deeper firn. Characteristic surface features, such as sastrugi or snow dunes seem not to be imprinted in deep ice; in deeper firn at densities between ˜630 and 750 kg m^-3, density and porosity are rather correlated with the concentration of impurities. It seems reasonable that a preferred roughness will decrease with increasing depths, but the sinusoidal pattern observed in Figure 2 starts at ˜200-300m and sharpens its directional dependence with increasing depth. However, at shallower depths, the quality of the radar data is insufficient and wind-induced anisotropy may also be present in the upper 0-300 m. Polarimetric data from a radar system designed to illuminate the shallow layering of ice, and an improved understanding of densification processes are needed to further investigate this topic.

The effect of enclosed (and possibly distorted) air bubbles on the radar backscatter has been discussed in several studies. Using a spherical shape, a frequency of 35 MHz, the dipole approximation and a volume-type radar equation, Reference Robin, Evans and BaileyRobin and others (1969) suggested that air bubbles may have a small but noticeable effect on the backscattered power. Discussing birefringence, Reference HargreavesHargreaves (1978) estimated the effect from elongated air bubbles to be too small to explain experimental results. Reference Ackley and KeliherAckley and Keliher (1979) used data from the Cape Folger ice core, near the coast of East Antarctica, and estimated the effect from the elongated air bubbles to be larger than that from COF variations. They observed layered strata where the bubble elongation corresponded to density fluctuations, suggesting a link to deformation. Reference Alley and FitzpatrickAlley and Fitzpatrick (1999) analytically described the bubble elongation as an interplay between strain-induced deformation and diffusive restoration. Very elongated bubbles with the long axis being several times the short axis are seen most frequently in ice cores near the coast within high-shear regimes. The observed eccentricities in the EDML ice core are much smaller. Reference Fujita, Matsuoka, Ishida, Matsuoka and MaeFujita and others (2000) used new COF data and placed the effect from air bubbles on the backscatter lower than the dominating density, electrical conductivity and COF mechanisms.

Three indications motivate us to investigate the effect of air bubbles on the backscatter more closely: (1) Reference UeltzhöfferUeltzhöffer and others (2010) showed in vertical cuts that air bubbles at the EDML drill site are vertically deformed. This means that bubbles can maintain an ellipsoidal shape within the ice matrix, and there is reason to believe that they might also be deformed in the horizontal plane; (2) the coincidence of the clathrate transition zone with the transition between the anisotropic reflection zones (ARZ-1 and ARZ-2) and (3) the strain ellipsoid in Figure 2b indicates the orientation of deformed bubbles (if surface strain is considered as the exclusive mechanism for bubble deformation). The direction of the long bubble-axis coincides with the direction of maximum backscatter. In Section 3.1 we apply an analytical volume-scattering model to estimate the effect of discrete ellipsoidal scatterers in the ice matrix. In Section 3.2, COF data from the EDML ice core and the corresponding effect on scattering are discussed, in order to differentiate between the different mechanisms.

3.1. Volume-scattering model for a half-space with discrete ellipsoids

3.1.1. Model set-up

The analysis presented here is based on Reference Tsang, Kubasci and KongTsang and others (1981). Volume-scattering models of this kind have been developed primarily for satellite applications, for example to classify the backscatter from canopy. We adapt the model for our purpose by changing the boundary conditions and focusing on a nadir-looking case. The geometry is shown in Figure 4. It includes two half-spaces 0 and 1, with permittivity ε₀ for air and ε₁ for ice, separated at the boundary z = 0. Ellipsoids with permittivity ε_s are randomly positioned in the lower half-space. The approach is based on the radiative transfer theory, which deals with specific intensities, I, rather than the electromagnetic fields, E (see, e.g., Reference Tsang, Kong and DingTsang and others, 2000, for an introduction). The specific intensities are 4 × 1 vectors. The first two components are the vertically (index v) and the horizontally (index h) polarized intensities l_v and l_h. Although for nadir incidence the distinction of vertical and horizontal becomes ambiguous (Reference Doake, Corr and JenkinsDoake and others, 2002), we keep the notation to remain consistent with other models. The indices v and h thus represent polarization types corresponding to orthogonal orientations of the (copolarized) antennas in the horizontal plane. The third and fourth components depict the correlation between vertical and horizontal polarization and enable the description of any polarization state. Using the modified Stokes vectors, the specific intensities are given by

Fig. 4. Sketch for the model of discrete ellipsoidal scatterer, modified after Figure 1 of Reference Tsang, Kubasci and KongTsang and others (1981).

where I = [1, 0, 0, 0] corresponds to the fully vertically polarized and I = [0, 1, 0, 0] to the fully horizontally polarized case; * depicts the complex conjugate and Re and lm refer to the real and imaginary part, respectively. The changes of the incoming intensity with depth are caused by extinction and scattering. The extinction is described with the extinction matrix, as a sum of absorption (index a) and scattering (index s) away from the receiver. We consider absorption to be isotropic. The phase matrix, operates for scattering from direction θ’, ϕ′. into the direction θ, ϕ All matrices are 4 × 4. Shape, dielectric contrast and spatial orientation of the ellipsoids define the individual components of . They can be calculated in a Rayleigh approximation using scattering dyads (Reference Tsang, Kubasci and KongTsang and others, 1981). The Euler angles, α, β, γ, define the orientation of the randomly positioned bubbles. The rotation around the z-axis corresponds to the angle α. For a definition of all angles see Figure 2 of Reference Tsang, Kubasci and KongTsang and others (1981). The principal axes of the ellipsoids are defined as a, b and c; a and b are in the horizontal plane when α = β = γ = 0. Each element is averaged over a common orientation density distribution, ψ(α, β, γ, Δα, Δβ, Δγ), which quantifies the degree of alignment for the bubbles. For simplicity, ψ is assumed to be box-like in all three angles, i.e. within the corresponding Δ-values all angles are equally probable, while outside the Δ-values the probability is zero. The fully aligned and the fully random cases for α correspond to Δα = 0° and Δα = ±90°, respectively. For non-spherical shapes, is non-diagonal and the different components in the specific intensity, I, mix, leading to depolarization. The governing equation is

(1)

The integration is over the solid angle dΩ’ = sin(θ’)dθ’ d?’. The angles are displayed in Figure 4, with ϕ being the azimuthal angle in the horizontal plane.

We first treat the general case of oblique incidence, and then consider the backscattering for a nadir-looking sensor as a special case. The integro-differential equation (1) can be solved using perturbation theory. The different terms of the perturbation series correspond to the different instances of multiple scattering. For backscattering, the incoming direction (index i) equals the scattered direction (index s). In terms of the normalized wavevectors this means with (with μ. being the magnetic permeability), where f is the frequency. The backscattering coefficient, σ_γ,δ , quantifies the ratio of incident specific intensity (direction with polarization δ = v or h) to backscattered specific intensity (direction - with polarization γ = v or h):

(2)

(see Fig. 2 The first-order approximation of Eqn (1), marked with superscript (1), corresponds to single scattering. For the nadir-looking geometry (θ_oi = θ_i = 0; π + ϕ_i = ϕ_oi = 0) we obtain

The scattering coefficient is a function of several parameters (e.g. frequency, volume fraction, incidence angles, eccentricities, alignment of bubbles) through the phase matrix. The transmission between air (index 0) and ice (index 1) for the different polarizations is labeled with the power transmission coefficients, t_vij and t _hij, where i and j can take the values 0 or 1. Since the air/ice interface is assumed to be isotropic, they play a minor role. The variables k_v and k _h are the first two diagonal elements of the extinction matrix, . Extinction includes off-angle scattering and absorption in ice. The latter is considered with the complex part of the permittivity: ε’’₁ = σ/(2πf ε₀).

We also calculated a second-order solution (Reference DrewsDrews, 2009) which treats multiple scattering. However, apart from depolarization effects, this is of minor importance for the parameters used here. The analytical solution for higher-order terms is only feasible for spherical inclusions (or spheroids), and a fully numerical approach or other approximations are needed if higher-order terms are desired.

3.1.2. Model outcomes

The polarization dependence of the backscatter coefficient is simulated by keeping the antenna orientation fixed, and rotating the bubbles in the horizontal by changing the values for the Euler angle, α. The sinusoidal variation of the backscattering coefficients is shown in Figure 5a for the case of fully aligned ellipsoids (Δα = Δβ = Δγ = 0). The absolute value of the phase matrix elements is mainly dependent on the ratio of bubble volume to wavelength. The magnitude of the backscatter coefficient is also influenced by absorption. For the examples presented here, we assume an average conductivity of 15 μS m^-1, as determined by dielectric profiling measurements (Reference Wilhelms, Kipfstuhl, Miller, Heinloth and FirestoneWilhelms and others, 1998; Reference WilhelmsWilhelms, 2000) from the EDML ice core. The measurements are scaled to the center frequency of the RES system (Reference Eisen, Wilhelms, Steinhage and SchwanderEisen and others, 2006). The difference of vertical and horizontal backscatter coefficients depends on the bubble eccentricity and the degree of alignment determined by the orientation density function.

Fig. 5. Model outcomes for the backscatter coefficients (hh, vv with solid curves and vertical axis on the left, hv with dashed curve and vertical axis on the right). (a) α rotates in the horizontal plane while β and γ remain fixed. No orientation density distribution is applied.This corresponds to the fully aligned case (Δα = Δβ = Δγ = 0). (b) Increasing disorder is simulated for α = β = γ = 0 and increasing Δα. Parameters for (a) and (b) are ε1 = (3.2+0.0018i)ε₀, ε_s = ε₀, Θ_oi = φoi = 0, vol = 0.0001 (volume fraction), a = 0.15 mm, b = c = 0.1 mm.

The effect of order and disorder in bubble alignment is shown in Figure 5b. The Euler angles α = β = γ = 0 are kept constant with Δβ = Δγ = 0. The disorder in the horizontal plane increases with increasing Δα. For the fully random case (Δα = 90°) the polarization dependence ceases.

3.2. Changing crystal orientation fabric

Crystal orientation fabric data of the EDML ice core were measured in 2005 in a -20°C cold room using an automatic fabric analyzer system (Reference Wilson, Russell-Head and Sim HMWilson and others, 2003). The EDML core was drilled between 2001 and 2004 and was stored at -30°C after transportation to the storage facilities. Thin sections were prepared using a microtome from samples cut horizontally (˜0.5 × 50 × 50mm³) and vertically (˜0.5 × 50 × 100 mm³). The automatic fabric analyzer measures the crystal orientation of these samples within 30min. This is possible due to the full automation of a polarization microscope (Reference Wilen, DiPrinzio, Alley and AzumaWilen and others, 2003) with a CCD (charge- coupled device) camera acquiring images of the thin section between rotating crossed polarizers from nine different viewing angles. The derived eigenvalues characterizing the COF distribution should, in theory, be independent of the cutting direction. However, in particular for eigenvalues λ₂ and λ3, an offset between the horizontal and vertical cuts is evident (Fig. 6a, difference between filled and unfilled symbols). We believe this artifact is introduced by the automatic exclusion of grains with a c-axis orientation that is difficult to measure (e.g. the c-axes being parallel to the section plane). Since the COF is non-random but aligned, these cases occur in different frequencies when the cutting direction is changed. The estimated uncertainty for the magnitude of the individual eigenvalues is 0.1, but confidence in relative changes of eigenvalues determined with the same method can be considered higher.

Fig. 6. (a) Crystal orientation fabric data from the EDML ice core in terms of eigenvalues (λ1=triangle, λ2=square, λ3=circle; filled symbols for horizontal cuts, open symbols for vertical cuts). The depths of two prominent internal reflection horizons in the radar data are marked h1 and h2. The difference, Δλ_i , for horizontal cuts measured in 0.9m intervals is displayed in (b) and (c). The top axes mark the corresponding power reflection coefficients (PRCs), as calculated from a two-layer approximation given in Eqn (3). Changes parallel to the propagation direction (i.e. changes in λ₃) do not influence the backscatter. (d) Non-sea-salt (nss) Ca²⁺concentration as published by Reference FischerFischer and others (2007) serves as a proxy for impurities within the ice which decrease in concentration during the transition from the last glacial into the Holocene (700–900m depth).

The COF in EDML shows a gradual transition from isotropically oriented c-axes at shallower depths (down to ˜450 m) to a broad girdle fabric (˜450-1700 m depth). The eigenvalues in Figure 6a reflect this transition by evolving from three similar values into three distinct values (λ₁ < λ₂ < λ3). We attribute this to the evolving influence of uniaxial extension and lateral flow dilatation. The girdle fabric narrows from 1700 m down to ˜2000m depth, indicated by converging λ₂and λ₁ This is possibly caused by a transition from dominating horizontal uniaxial extension to bed-parallel simple shear.

Fabric data are measured at intervals of 50 m or less using vertical cuts. Starting at 450 m depth, the measurements are complemented with two horizontal cuts within 0.9m distance every 50 m. We use this dataset as a proxy for short- scale variations in COF. The difference, Δλ_i (i = 1,2,3), within the 0.9m interval for the three eigenvalues is plotted for increasing depth in Figure 6b and c. In the range of the first anisotropic reflection zone (ARZ-1), the variation of the smallest eigenvalue, λ₁, is noticeable, and decreases towards larger depths. This corresponds to the previously discussed variation in girdle strength for ARZ-1. The opposite behavior becomes apparent for the two other eigenvalues, λ₂ and λ₃, which eventually swing in antiphase parallel to ARZ- 2’s appearance in the radar data. This implies the tendency of girdle to single-pole transitions in ARZ-2. The relative changes of the eigenvalues considered here are close to the detection limit of the fabric measurements. However, at least at larger depths, the general tendency appears systematic and seems to fit with the effects seen in the radar data, as discussed in Section 3.3. We therefore carry on, but bear in mind that more accurate COF data, preferably at higher vertical resolution, are needed to solidify the later conclusions drawn from this dataset with respect to the interaction with radar waves.

We use a layer approximation to estimate the effect on the backscatter. A premise is that the changes in COF as observed in the ice core are laterally coherent on a scale that is much larger than the illuminated area of a single RES pulse. For larger depths, this has been shown by Reference Eisen, Hamann, Kipfstuhl, Steinhage and WilhelmsEisen and others (2007), who linked a single RES reflector to changes in COF and traced it laterally for a few kilometers. For the intermediate depths this cannot easily be done, because changes in COF are less distinct. However, the evolution of COF is related to strain rates, and there are also indications that a varying impurity content influences the crystal size and thus the formation of COF (Reference ThorsteinssonThorsteinsson, 1996). Both effects favor a laterally coherent stratification of COF. Prominent internal reflection horizons (h1 and h2) are visible at 1068 and 1312 m depth in Figure 2a. Reference Eisen, Wilhelms, Steinhage and SchwanderEisen and others (2006) linked both horizons to (isotropic) reflectors caused by changes in electrical conductivity and we would expect – if at all - brightness variation in the polarimetric data with the symmetry of birefringence. The internal horizon h2 appears to be independent of the incoming polarization, whereas h1 shows some brightness variations, with minima and maxima coinciding with the characteristics of ARZ-2. Insufficient data are available to fully explain this. It is possible that impurity content in ice impacts the COF type. Changes in electrical conductivity may be accompanied by laterally coherent changes in COF, leading to a polarization dependence of the corresponding reflection horizon. It seems that the layer approximation is certainly simplified but not unrealistic.

In order to characterize the dielectric properties of the individual layers, we need to know the principal directions associated with the eigenvalues. This is difficult, since no continuous azimuth control in the horizontal plane is available for the ice core. However, the third eigenvector seems to be in the vertical, which becomes apparent in some Schmidt diagrams by an accumulation of points in the center of the girdle. We therefore tentatively assume that λ₃ points along the z direction (see also the discussion of Reference Eisen, Hamann, Kipfstuhl, Steinhage and WilhelmsEisen and others, 2007). The directions of λ₁ and λ₂ remain undefined in the horizontal plane, but based on the vertical alignment of minima and maxima in Figure 3, we assume that they do not change with increasing depth.

Reference Fujita, Maeno and MatsuokaFujita and others (2006) described a method to incorporate the eigenvalues into the dielectric tensor. Power reflection coefficients (PRCs) determine the ratio of incoming and scattered fields. They can be calculated in a planar two- layer approximation for changes in the permittivity, ε’, as illustrated by Reference ParenParen (1981). Along the principal axes of λ _i (i = 1,2,3), the PRC is given by

(3)

The values of r resulting from changes in eigenvalues in the 0.9 m intervals (upper boundary depicted with z₁, lower boundary depicted with z₂) are displayed in the top axis of Figure 6b and c.

For nadir incidence, the polarization of the electromagnetic wave is in the vertical plane parallel to the airplane’s heading. Therefore, changes in the third eigenvector do not contribute to the scattering process in this approximation. The stronger variation of λ₁ in the 0.9 m intervals within ARZ-1 fades in ARZ-2 and simultaneously the variation of λ₂ increases. This directly translates into the polarization dependence of the PRCs. In ARZ-1, backscatter is largest when the polarization is parallel to the principal direction of λ₁; conversely, in ARZ-2 the maximum in backscatter is shifted by 90° and aligns with the principal direction of λ₂. The differences in the corresponding PRCs reach up to 10 dB.

3.3. Ellipsoidal air bubbles and depth-varying COF

We estimated the effect of ellipsoidal air bubbles with a volume-scattering model, and the effect of depth-varying COF with reflections from a stratified medium. For ARZ-1, it is not straightforward to distinguish the effect of elongated bubbles and varying cluster strength in COF in terms of their orientation, since they most likely operate in similar directions. The deformation of air bubbles in the horizontal plane, as expected from strain measurements on the surface, is indicated with the strain ellipsoid in Figure 2b. The direction of the long axis is almost parallel to the polarization plane with maximum backscatter in ARZ-1. The azimuth angle of the COF measurements is not constrained; however, Reference Eisen, Hamann, Kipfstuhl, Steinhage and WilhelmsEisen and others (2007) inferred that the girdle at ˜2025 m depth is in the vertical plane parallel to the nearby ice divide. A recent modeling study for COF evolution at this site by Reference Bargmann, Seddik and GreveBargmann and others (2011) confirms this result, in the sense that the girdle is aligned with the principal axes of the strain ellipsoid. One of the axes is near-parallel to the ice divide. It is therefore reasonable to assume that both effects of elongated air bubbles and a varying cluster strength in COF superimpose and are aligned with the strain regime on the surface.

The estimations of both scattering coefficients and PRCs have their shortcomings, due to the applied simplifications. The volume-scattering model neglects characteristics of the finite radar pulse, as only the center frequency and scattering contributions from incoming nadir incidence are taken into account. The system’s bandwidth for the short pulse is ˜17MHz, and an interplay with higher- frequency components already increases the scattering coefficient by a few decibels, since the frequency enters the equations for the phase-matrix components with a power of four. The beam widening with increasing travel time may also give rise to scattering components with nonnadir incidence. This generally increases the polarization dependence. The absolute value for the scattering coefficient also heavily depends on the absorption, which dominates the extinction, since the phase-matrix components are small. Potential effects of anisotropic absorption, ora non-diagonal extinction matrix have been neglected. The bubbles were assumed to be randomly distributed. However, it is known that bubbles cluster in so-called ‘cloudy bands’, which leads to a non-uniform bubble distribution within the ice (Reference SvenssonSvensson and others, 2005). Reference Ackley and KeliherAckley and Keliher (1979) took this idea further and estimated the effect of air bubbles in a way similar to our approach for the depth-varying COF. The bubble stratification was modeled with a scalar mixing formula for air and ice with a different depolarization factor for each layer due to varying bubble eccentricities. The corresponding anisotropy and backscattered power with the approach for stratified media is larger than that estimated with the volumescattering model using similar parameters. More data for air-bubble statistics in ice are needed - especially in the horizontal plane - in order to decide whether surface or volume scattering is a more suitable description. In any case, the absolute values of the scattering coefficients have a high uncertainty. This is also the case for the PRCs derived for the variations in COF: reflections from within the ice are subject to interferences, and the sample interval of the COF data is not short enough to identify potential reflection layers reliably. Changes in COF with depth in the range of the wavelength (˜1 m in ice) may be larger than those portrayed in Figure 6. The interfaces, however, may not be planar, but could be slanted and/or rough. This substantially decreases the magnitude of the PRCs. This is to say that the absolute values of scattering coefficients and PRCs should be handled with care. Moreover, quantitative comparison of backscattered power from volume and surface scattering requires consideration of transmit and receive beam patterns, since the scattering coefficient is normalized per unit volume, while the PRCs are normalized per unit area subtended. In view of these difficulties, we do not attempt to identify the dominant backscattering mechanism on the basis of absolute backscattered power. We rather focus on the relative polarization dependence. The example for volume scattering in Figure 5a yields a ˜1 dB difference in a somewhat favorable case of fully aligned air bubbles with typical size parameters for intermediate depths. In the upper half of ARZ-1, the PRCs due to changes in eigenvalue λ₁ appear several decibels stronger than those caused by the changes in eigenvalue λ₂ (Fig. 6b and c), whereas in the lower half of ARZ-1 they appear almost equal. For ARZ-2, the PRCs due to changes in λ₂ outweigh the PRCs due to changes in λ₁ by >10dB. Based on the relative polarization dependence, we estimate that in ARZ-2, and probably also in ARZ-1, changes in COF are the dominant mechanisms for the observed anisotropy.

We briefly discuss the shading of ARZ-2 by ARZ-1. The increased backscatter parallel to the direction of maximal compression diminishes the transmitted power with the corresponding polarization in that direction. Therefore, polarization-dependent variations in backscattered power from ARZ-2 are to be expected, even without a specific anisotropic backscatter mechanism in that zone. To estimate this effect, we consider a lossless, stratified medium containing n equidistant layers with thickness Δz and anisotropic PRCs r_x and ry. The ratio, |E_t,x|²/|E_t,y|², of the transmitted power after two-way penetration for the two polarization types is

(4)

which ignores multiple reflections as they are small. To estimate an upper boundary we assume n = 900 layers with strongly anisotropic PRCs of r_x = -40dB and r_y = - 50dB with a layer thickness Δz = 1 m. This extreme case results in a ratio of -0.7dB, which is smaller than the ratio observed in ARZ-2.

4. Conclusion

We observe a variation in signal strength depending on the antenna orientation. The dominant symmetry has a 180° periodicity, which we assign to anisotropic scattering from within the ice. To a lesser extent we also observe a 90 ^° periodicity due to birefringence caused by the macroscopic alignment of the COF. The polarization-dependent backscatter changes its direction with depth at ˜900m, coinciding with the clathrate transition and the transition from Holocene ice to ice from the last glacial (700-900 m depth). Due to a lack of data, we neglected the effect of a potentially (wind-induced) directional layer roughness and only investigated backscatter mechanisms above 900 m (ARZ-1) and below 900 m (ARZ-2) in terms of elongated air bubbles and a depth-varying COF. For ARZ-1, we propose that both effects superimpose and operate in similar directions, predetermined by the principal axes of the stress regime at the surface. We assume that bubbles are compressed parallel to the axis of maximal compression and elongated along the axis of lateral extension, causing minimum and maximum in backscatter in these directions, respectively. Similarly we hypothesize that the girdle in the COF is aligned along the same axes; this idea is supported by recent fabric-evolution modeling results in this area (Reference Bargmann, Seddik and GreveBargmann and others, 2011). The eigenvalues in ARZ-1 suggest that the horizontal width of the girdle-type COF varies with depth, increasing the backscatter perpendicular to the direction of maximal compression. Most likely, the variations in COF dominate the observed anisotropy, since the induced polarization dependence through elongated bubbles is small. However, a potential stratification of air bubbles has been neglected and higher- resolved COF measurements, as well as measurements with respect to orientation and eccentricities of the air bubbles in the horizontal plane, are needed to quantify the effects reliably. The volume scattering for anisotropic inclusions becomes more important at shallower depths when higher frequencies are used and the volume fraction of air is larger.

In ARZ-2, we propose that changes in eigenvalue λ₂ and a stable eigenvalue, λ₁ , increase the backscatter in the direction normal to that observed in ARZ-1. This corresponds to the girdle to single-pole transitions suggested by Reference Fujita, Matsuoka, Maeno and FurukawaFujita and others (2003) and Reference MatsuokaMatsuoka and others (2003). In contrast to ARZ-1, the variations in COF appear more clearly in the ice- core data.

So far, it is unclear how the eigenvalues of COF change vertically in meter and sub-meter intervals and how the evolution is influenced by the impurity content of the ice. Given the available evidence, it is possible that the transition in anisotropic backscatter caused by changing variational patterns of COF indicates the climate transition at the end of the last glacial (the relationship with the clathrate transition at similar depths appears to be merely coincidental). To visualize the relationship of increasing dust concentration with the changing COF variations, Figure 6d displays the non-sea-salt Ca²⁺ concentration (Reference FischerFischer and others, 2007) as an impurity proxy alongside the COF measurements. For ARZ-1, the concentrations are low compared to all other depths. The exact mechanism on a grain-size scale still needs to be investigated, and it seems that the flank-flow setting is important, as Reference Fujita, Maeno and MatsuokaFujita and others (2006) identified anisotropic scattering near Mizuho station (converging ice-flow region), but not near Dome Fuji (located near a dome summit). More radar data are needed to determine whether the appearance of the ARZs is a universal feature in flank-flow regimes.