Rapid and accurate polarimetric radar measurements of ice crystal fabric orientation at the Western Antarctic Ice Sheet (WAIS) Divide deep ice core site

. The Crystal Orientation Fabric (COF) of ice sheets records the past history of ice sheet deformation and inﬂuences present-day ice ﬂow dynamics. Though not widely implemented, coherent ice-penetrating radar is able to detect anisotropic COF patterns by exploiting the birefringence of ice crystals at radar frequencies. Most previous radar studies quantify COF at a coarse azimuthal resolution limited by the number of observations made with a pair of antennas along an acquisition plane that rotates around an azimuth centre. In this study, we instead conduct a suite of quad-polarimetric measurements consisting of 5 four orthogonal antenna orientation combinations at the Western Antarctic Ice Sheet (WAIS) Divide Deep Ice Core site. From these measurements, we are able to quantify COF at this site to a depth of 1500 m at azimuthal and depth resolutions of up to 1 ◦ and 15 m . Our estimates of fabric asymmetry closely match corresponding fabric estimates directly measured from the WAIS Divide Deep Ice Core. While ice core studies are often unable to determine the absolute fabric orientation due to core rotation during extraction, we are able to unambiguously identify and conclude that the fabric orientation is depth-invariant to at least 10 1500 m , equivalent to 7400 years BP (years before 1950), and coincides exactly with the modern surface strain direction at WAIS Divide. Our results support the claim that the deformation regime at WAIS Divide has not changed substantially through the majority of the Holocene. Rapid polarimetric determination of bulk COF compares well with much more laborious sample-based COF measurements from thin ice sections. Because it is the former that ultimately inﬂuences ice ﬂow, these polarimetric radar methods provide an opportunity for accurate and widespread mapping of bulk COF and its incorporation into ice ﬂow 15 models.


Introduction
There is a growing need to understand the dynamics of ice sheets and how they will respond to future climate change (IPCC, 2013). The flow of ice sheets is governed by the balance between the gravitational driving stress, basal resistance to sliding, and the internal deformation of ice (Cuffey and Paterson, 2010). Past flow history influences the ice crystal orientation fabric (COF), which, in turn, influences the present-day anisotropic ice viscosity and flow field. Because ice crystals effectively reorient themselves to minimise resistance when subjected to stress, the COF of ice is reflective of long-term strain at time scales proportional to the depth-age relationship (e.g., Alley, 1988), and has consistently been observed to align with ice deformation (Matsuoka et al., 2012). The COF of ice is also thought to be influenced by perturbations in climate on a yearly timescale (Kennedy et al., 2013). Additionally, abrupt vertical changes in COF are often indicative of paleoclimatic transitions (e.g., 25 Durand Montagnat et al., 2014;Paterson, 1991). Therefore, an examination of the present-day COF can reveal past changes in the stress-strain configurations associated with historical ice flow (Brisbourne et al., 2019).
Ice core analyses represent the traditional method to quantify COF within ice sheets, and remain the only direct means of ground-truth observation. However, sites suitable for ice coring are often restricted to slow-moving (< 50 m a −1 ) sections of ice sheets, and therefore only reveal a subset of the dynamics between ice flow and COF. These slow-flowing areas most likely 30 do not encapsulate the dynamics and physical processes responsible for ice sheet stability and sea level rise. Furthermore, ice core analyses are often unable to resolve the absolute direction of fabric asymmetry due to the rotation of the core in the barrel during or following extraction (Fitzpatrick et al., 2014). In contrast, ice-penetrating radar offers an alternative method to calculate anisotropic COF patterns through exploiting the birefringence of polar ice without the practical limitations of drilling.
Although polarimetric radar sounding data analysis has been implemented to detect horizontally-asymmetric COF for almost 35 half a century (e.g., Hargreaves, 1977), it has not yet been widely implemented, with the majority of radar studies that measure COF variations in ice being conducted within the last fifteen years at coincident ice-coring sites for comparative analysis (Drews et al., 2012;Eisen et al., 2007;Fujita et al., 2006;Jordan et al., 2019Jordan et al., , 2020aLi et al., 2018;Matsuoka et al., 2003Matsuoka et al., , 2009Matsuoka et al., , 2012. Moreover, the majority of previous polarimetric radar studies infer COF at a coarse azimuthal resolution that is limited by the number of observations made along an acquisition plane that rotates around an azimuth centre (Brisbourne et al.,40 2019; Doake et al., 2002Doake et al., , 2003Jordan et al., 2020b;Matsuoka et al., 2003Matsuoka et al., , 2012. Because the maximum azimuthal resolution that can be achieved is subject to human error in measuring the angles between each acquisition plane, there is a coarse limit to the precision of the orientation of fabric asymmetry that can be achieved through this acquisition method. Traditionally, radar studies estimated COF through power-based analyses by investigating the periodicity of birefringent patterns in power anomaly (Fujita et al., 2006;Matsuoka et al., 2012) and phase difference measurements (Brisbourne et al., et al., 2006;Lindbäck et al., 2019;Marsh et al., 2016;Nicholls et al., 2015;Stewart et al., 2019;Sun et al., 2019;Vaňková et al., 2020;Washam et al., 2019), and subsurface imaging .
Instead of using an azimuthal rotational setup, acquisitions can instead be obtained through a combination of four orthogonal antenna orientations, from which the received signal can then be reconstructed at any azimuthal orientation (Fig. 1a) (Fujita et al., 2006;Jordan et al., 2019). This quadrature-(quad-) polarised setup therefore resolves this limitation on azimuthal resolution whilst simultaneously also significantly reducing the field time required to obtain each set of acquisitions. 60 In our study, we use an ApRES and two antennas to acquire a set of quad-polarised measurements at the Western Antarctic Ice Sheet (WAIS) Divide. We apply the polarimetric coherence method to these measurements to present estimates of fabric values that closely align with previous ice core COF measurements at WAIS Divide to a depth of 1500 m at a resolution of 15 m. We show that, using this setup and method, we can estimate the fabric asymmetry at resolutions at or exceeding that from ice cores, and from our results unambiguously identify the principal axis of present and past flow to high angular resolution.

2 Study area
The West Antarctic Ice Sheet (WAIS) divide delineates the surface topographic boundary separating ice flow towards the Ross and Amundsen Sea Embayments (Fig. 1b). A subglacial topographic saddle runs approximately orthogonal to the ice divide in the study region. Ice flow is oriented approximately SWW (250 • ) and NE (45 • ) in the Ross and Amundsen Sea catchments respectively (Conway and Rasmussen, 2009), which is offset from predicted strain configurations, especially in the Amundsen 70 Sea catchment (Matsuoka et al., 2012). The WAIS divide flow boundary is observed to be migrating in the direction towards the Ross Sea at 10 m a −1 , faster than the surface velocity of 3 m a −1 , the migration being attributed to differential flow dynamics in the Ross Sea catchment (Conway and Rasmussen, 2009). Despite this fact, the ice-divide position has likely remained on average within 5 km of its present position throughout the Holocene epoch (Koutnik et al., 2016).

Theory and methods
We primarily follow a combination of three matrix-based methods- Fujita et al. (2006), Dall (2010), andJordan et al. (2019), in 80 which each study builds on the previous-to process the polarimetric measurements, thereby obtaining estimates of anisotropy and fabric strength. Our approach is physically justified through an effective medium model that expresses the bulk dielectric properties of anisotropic polar ice in terms of the birefringence of individual ice crystals (Fujita et al., 2006). First, we follow the framework of Fujita et al. (2006) and Brisbourne et al. (2019) in modelling the expected power anomaly and phase difference (Section 3.2) to explicate the ApRES' response to the underlying COF parameters from data collected in proximity at WAIS  Morlighem et al., 2020) in the WAIS Divide area, as well as GPS-measured surface velocities (black lines, Conway and Rasmussen, 2009) and strain configurations (blue and red arrows, Matsuoka et al., 2012). The locations of the ApRES polarimetry experiment (red dot) and the WAIS Divide Deep Ice Core (WDC, cyan star) are separated by~5 km. The WAIS Divide is delineated as a thick dotted white line, with ice flowing northwards towards the Amundsen Sea and southwards towards the Ross Sea. Location of (b) is shown as a red box in the map inset.
Divide. We then calculate the fabric asymmetry using the polarimetric coherence methods outlined in Jordan et al. (2019Jordan et al. ( , 2020b.

Electromagnetic propagation and COF representation in anisotropic ice
In an anisotropic medium such as polar ice, birefringence and anisotropic scattering are two related, but separate mechanisms that affect the polarisation and azimuthal variation in power of radar returns (Brisbourne et al., 2019). For downward-looking 90 ice-penetrating radar, birefringence occurs as a result of a phase shift between two orthogonally-oriented waves travelling between the surface and the interior of an ice mass, the phase shift manifested in the radar return as characteristic variations with azimuth and depth in power and phase. As a result, birefringence reflects the bulk COF as and is azimuthally asymmetric in the direction of radio wave propagation. On the other hand, anisotropic scattering arises as a consequence of rapid but microscopic continuous depth variations in the orientation of the bulk COF. Therefore, the polarimetric response of radio 95 waves is determined by the bulk (macroscopic) birefringence of the COF (Hargreaves, 1978), of which the area illuminated by 4 https://doi.org/10.5194/tc-2020-264 Preprint. Discussion started: 28 September 2020 c Author(s) 2020. CC BY 4.0 License. the waves is a function of the radar antenna footprint through depth. The birefringence of an individual ice crystal and its COF are related to the bulk dielectric properties of anisotropic polar ice (Appendix of Fujita et al., 2006): where ε (z) is the bulk birefringence tensor, ∆ε = ε − ε ⊥ is the crystal (microscopic) birefringence with ε and ε ⊥ the 100 dielectric permittivities for polarisation planes parallel and perpendicular to the crystallographic (c)-axis. Across the spectrum of ice-penetrating radar frequencies and ice temperatures, ε and ε ⊥ vary within a narrow band of 3.16-3.18 and 3.12-3.14 respectively (Fujita et al., 2000). In this study, following Jordan et al. (2020b), we assign ε = 3.169 and ε ⊥ = 3.134, with ∆ε = 0.035.
The tensor eigenvalue E describes the relative concentration of c-axes aligned with each principal coordinate eigenvector,

105
with E 1 + E 2 + E 3 = 1 and E 3 > E 2 > E 1 following conventional radar notation, which is opposite to conventions normally in ice core studies (E 1 > E 2 > E 3 ). The relative proportions of E can be used to describe different fabric patterns, including . When ice deforms solely by vertical uniaxial compression, such as at the centre of an ice dome, the c-axes rotates towards the vertical and forms a cluster fabric; where lateral tension exists from flow extension, such as at an ice 110 divide, the c-axes orient in a vertical girdle distribution orthogonal to the direction of strain extension (Alley, 1988). Following previous studies (Fujita et al., 2006;Drews et al., 2012;Brisbourne et al., 2019;Jordan et al., 2019Jordan et al., , 2020b, we assume that the E 3 eigenvector is aligned in the vertical direction, and the E 1 and E 2 eigenvectors are parallel to the horizontal plane. The direction of the greatest horizontal c-axis concentration reflects the orientation of horizontal strain extension (Brisbourne et al., 2019;Matsuoka et al., 2012) and corresponds with E 2 in our notation (Jordan et al., 2020b). The E 1 eigenvector is 115 orthogonally oriented to both E 2 and E 3 , and represents the symmetry axis, which is normal to the girdle plane (Brisbourne et al., 2019).
In the horizontal plane, Equation 1 simplifies to ε (z) = ∆ε (E 2 − E 1 ), where the horizontal eigenvalue difference E 2 − E 1 quantifies the horizontal asymmetry of the ice fabric (i.e. strength of the vertical girdle). Equation 1 directly relates both the macro-and microscopic ice fabric anisotropy to dielectric anisotropy, which serves as the basis for the radar processing 120 methods that follow.

Modelling radio-wave signal propagation
The matrix-based formulation calculates the backscatter that is transmitted, reflected, and received at the antennas for each depth step and azimuthal orientation: 125 5 https://doi.org/10.5194/tc-2020-264 Preprint. Discussion started: 28 September 2020 c Author(s) 2020. CC BY 4.0 License.
Equation 2 represents the polarimetric backscatter model described in Eqs. 9-12 of Fujita et al. (2006), which calculates the polarimetric backscatter for each antenna orientation combination as a function of angle (θ), anisotropy (β), and birefringence (ε), through all depth layers z = i to the N th layer. The first term on the right hand side represents the (i) free space propagation (squared to reflect two-way wave travel). Here, j = √ −1 is the imaginary number, and k 0 = 2π/λ 0 (rad m −1 ) the wavenumber in a vacuum with λ 0 the wavelength in a vacuum. Besides the first expression, the second to fourth terms respectively represent 130 three physical processes: (ii) received (upward) propagation; (iii) boundary scattering; and (iv) transmitted (downward) propagation to each boundary depth i. T and Γ respectively represent the transmission and scattering matrices, and R represents the rotation matrix with R its transposition, each detailed respectively in Eqs. 5 and 6, 8, and 10 of Fujita et al. (2006), with the propagation constants required to calculate T defined specifically for the ApRES unit in Brennan et al. (2014). MHz. An ensemble (burst) of 100 chirps were recorded for each polarimetric measurement. Two open-structure antennas (1 140 transmitting and 1 receiving) identical to those described in Nicholls et al. (2015) were used to transmit and receive each burst.

Radar data acquisition
Although not implemented in this study, we note that the four quad-polarimetric measurements can potentially also be obtained simultaneously in one single burst using a multiple-input multiple-output (MIMO) configuration (e.g., Young et al., 2018) with two transmitting and two receiving antennas.

Radar data processing
Data were pre-processed and range-processed following procedures detailed in Stewart et al. (2019). Specifically, for each of the four bursts, the 20 noisiest chirps were culled and the remaining chirps averaged. Each resulting burst-mean was then weighted with a Blackman window, zero-padded, time-shifted to align the phase centre with the start of the signal, and Fourier-155 transformed. The resulting complex-valued spectra (referred to as complex amplitude) store the amplitude and the phase of the signal as the magnitude and the angle of the spectra respectively.
Using the four processed quad-polarised complex amplitudes, the 2 × 2 integrated scattering matrix S (Eq. 2) in the frame of reference can be constructed as (Doake et al., 2003): Equation 3 is often referred to as the Sinclair matrix. From here, we can reconstruct the ApRES received signal S from any transmission angle through the application of an azimuthal (rotational) shift of principal axes at the transmitting and receiving antennas (e.g., Mott, 2006): Note that in Eq. 4 the cross-polarised measurements obtained from the ApRES are geometrically congruent by the Lorentz Reciprocity Theorem (i.e. s vh = s hv ) and therefore the two measurements should be identical, with any deviations being the result of performance variations between the transmitting and receiving antennas and/or error in orienting the antennas along the acquisition geometry axes (Doake et al., 2003).
Because the complex amplitudes retrieved from the ApRES (Eq. 4) are phasors representing the radar return signal, the 170 phase of the signal at a given depth is simply the argument of the complex number. To avoid the typical problems of working with phase-that is, employing phase unwrapping methods for sampled data within [0, 2π]-we calculate the phase difference with respect to azimuth: where the asterisk represents the complex conjugate of its respective phasor.

175
The shift in phase also results in the modulation of received power as a function of azimuth. This can be visualised by calculating the power anomaly from the resulting multipolarisation data (e.g. Eq. 7 in Matsuoka et al., 2003).
The polarimetric coherence and its corresponding phase is computed over a local window via the discrete approximation (Eq. 1 in Dall, 2010): where the superscript stars in Eq. 6 account for the use of the deramped phase stored by the ApRES rather than the original received signal phase (hereafter we do not notate this explicitly) (Jordan et al., 2020b). From Eq. 6b, we can then estimate the with the associated phase error was estimated through the Cramér-Rao Bound, following the methods of Jordan et al. (2019). In Eq. 7b, R and I respectively represent the real and imaginary components of c hhvv . f (ν) represents a reduction parameter for the birefringence of firn with respect to solid ice, with ν the firn density, as detailed in the Appendix of Jordan et al. (2020b).  Figure 2 shows the processed results from the ApRES polarimetry experiment. A pad factor of 2 (equating to a depth resolution of 0.27 m) and an azimuthal resolution of 1 • was used in the phase processing steps to produce co-and cross-polarised profiles of power anomaly and phase difference (Fig. 2b,c,d). A depth window of 15 m was used in the hhvv coherence and phase 195 estimates (Fig. 2e). We evaluated dφ hhvv /dz and estimated its respective error following the suggested procedures in Jordan et al. (2019), with the exception that, in substitution of the FIR filter, we used a 2-D median filter consisting of a 1 • ×5 m matrix moved over the profile, and then a 2-D peak convolution using a Gaussian low-pass filter with the same moving matrix dimensions . We compare the measured results in Fig. 2 with modelled results in Fig. 3 to parse the relative influence of birefringent propagation and anisotropic scattering on the ice column. From here, we estimate fabric asymmetry 200 via Eq. 7a to solve for E 2 −E 1 (Fig. 4). We restrict our observations only to measurements with sufficiently high signal-to-noise ratios (SNR). We also do not make any inferences in the uppermost 20 m due to potential radiation pattern effects.

Results
The depth-azimuth variation of the radar return power anomaly in the co-polarised and cross-polarised measurements are respectively shown in Fig. 2b and Fig. 2c, variations in the co-polarised return signal phase difference in Fig. 2d, and variations in the hhvv phase angle in Fig. 2e. Here, we observe azimuthal variations larger than 10 dB in the observed backscatter 205 power both in the co-polarised and cross-polarised measurements, and variations larger than 6 dB in the hhvv phase measurements. The backscatter variations in the co-polarised power anomaly (Fig. 2b) and hhvv phase measurements (Fig. 2e) show reflectional symmetry at 90 • at all measured depths in our experiment reference frame. The co-polarised phase difference measurements (Fig. 2d) show characteristic 'four-quadrant patterns' formed by azimuthal rotation of the phase-difference sign reversals (Brisbourne et al., 2019) along the same reflection axis. The cross-polarised power anomaly measurements, in 210 contrast, show a 90 • periodicity in the return power difference with a depth-constant azimuthal minima at 0 • and 90 • (Fig. 2c).
By tracing the azimuthal minima in the cross-polarised power anomaly profiles through depth (   Fig. 2b, and the latter using eigenvalues from the WDC that specify the bulk COF (Fitzpatrick et al., 2014). Value ranges for β (white to dark green) and ε (z) (white to dark purple) are [1 5] and 0 1.5 × 10 −2 respectively.
The identified symmetry axis orientation does not azimuthally migrate through the observed depth range to 1500 m, equivalent to a depth age of 7400 years (Sigl et al., 2016). Though we similarly observe no azimuthal migration beyond 1500 m, we do not extend our findings further due to the limited depth samples with sufficient SNR available within this range. Our estimate of the principal axis aligns exactly with the nearest strain configuration (~5 km southwest) as estimated by Matsuoka et al. (2012), and is +14 • from the direction of flow as estimated by Conway and Rasmussen (2009) (Fig 2f). We note that there 220 may be additional unquantified azimuthal error from human inaccuracy in establishing the orientation of the measurement axis during data collection.
The backscatter and phase patterns in the co-polarised measurements, as well as in the hhvv coherence phase (φ hhvv ), all show variation with depth, which indicate changes in either or a combination of birefringence and anisotropy. To better understand what drives these changes, we modelled the azimuth and phase dependence of these three measurements (Fig. 3)  eigenvalues from the ice core analysis suggest a gradual linear transition from an isotropic fabric (E 2 − E 1 ≈ 0.04) near the 230 ice surface to a moderately-strong girdle fabric (E 2 − E 1 ≈ 0.3) at 1500 m depth (Fig. 4). As the principal axis observed in the polarised experiment is depth-invariant, we fixed the principal axes at 90 • following our observations in Fig. 2c. The remaining parameter, anisotropy, was simply estimated as an optimisation problem through identifying the azimuthal distance of nodes (minima) in power anomaly from the principal axis orientation at 100 m depth intervals. At each depth interval, the anisotropy parameter β was estimated to the nearest integer and linearly interpolated to match the model depth step of 1 m.

235
The model outputs for the co-polarised power anomaly and phase difference, as well as the hhvv phase angle cases are shown in Fig. 3 measured phenomenon that arises as a consequence of bulk COF (Fig. 3d).
Several minor differences observed between the measured and modelled results include: (i) a scalar reduction in the overall ranges of power anomaly and phase difference values; (ii) the absence of the upper half of the shallowest phase reversal in the measured co-polarised phase difference and phase angle profiles (Fig. 2d,e); (iii) an offset in pattern repetition (e.g., nodes, quadrants, phase reversals) that increases in depth; and (iv) an absence of noticeable patterns in the shallowest 200 m of the 250 vertical ice column. Notwithstanding these differences, model results in Fig. 3 overall match their counterparts in Fig. 2b, d, and e to a high degree. With the exception of (ii), the corresponding locations and sizes of nodes in the power anomaly, quadrants in the phase difference, and asymptotic zones predicted by the model can all be observed in the measured results at 200-1500 m. Observations below 1500 m depth remain inconclusive, due to the low SNR present at this depth range (Fig. 2a).
Estimates of E 2 − E 1 , a measure of fabric asymmetry, were made by calculating the depth-dependent gradient of the hhvv 255 phase difference along the symmetry axis (Fig. 2e) at each 15-m depth window to obtain dφ hhvv /dz (Eq. 7b), and are shown alongside equivalent E 2 − E 1 measurements from the WDC (Fig. 4). We do not calculate E 2 − E 1 in the uppermost 20 m due to the antennas being subjected to near-surface non-axisymmetric antenna radiation pattern effects. Overall, fabric asymmetry estimates from ApRES match well to WDC ice core estimates, especially between 600-1200 m, where the mismatch between the two depth series averaged less than 0.02. When comparing the two independently-calculated fabric asymmetry datasets 260 at the discrete depths at which the ice core thin sections were extracted from, the resulting correlation was moderately high between 600-1200 m (r 2 = 0.55), and slightly lower over the entirety of the depth series overlap to 1500 m (r 2 = 0.45). Both datasets show a general increase in E 2 − E 1 with increasing depth, with ApRES fabric measurements overall showing more variability than estimates from the WDC along their respective trends. The one exception to this positive correlation occurs towards a vertical plane transverse to flow,but with stronger concentration near the vertical than near the horizontal in that plane,and with increasing vertical concentration with increasing depth. Grain subdivision by polygonization causes only small changes in c-axis orientation, and, as noted below,nucleation and growth of new grains with caxes at large angles to their neighbors does not become dominant until much deeper (see also Budd and Jacka, 1989; Alley, 1992). Lacking independent evidence of the orientation of core sections,this physical understanding has been applied;after rotating all sections to place the core axis in the center of the Schmidt plot,those sections with a clear symmetry plane were rotated so that these planes are

Eigenvalues
Eigenvalues provide a means by which th the fabric can be described using three sca ized eigenvalues have the property that S that S 1 > S 2 > S 3 . Ideally,a totally uniform fa S 1 ⇡ S 2 ⇡ S 3 ⇡ 1/3,a fabric with all c-axes a in a plane would have S 1 = S 2 = 1/2,S 3 = 0, all c-axes pointing in the same direction w S 2 = S 3 = 0. In reality, these values are Previous workers have used eigenvalue me the fabric in a quantitative manner (e.g. DiP 2005;Kennedy and others,2013). Eigenva  towards a vertical plane transverse to flow,but with stronger concentration near the vertical than near the horizontal in that plane,and with increasing vertical concentration with

Eigenvalues
Eigenvalues provide a m Fig. 19. Representative selection of horizontal Schmidt plots showing the evolution of fabric wit rotated so that all plots have the same orientation,which is assumed to be perpendicular to the dire depths (m) are indicated above each plot. Number of points (n) successfully measured by the c-axis below each plot. Fitzpatrick and others: Physi 1194 towards a vertical plane transverse to flow,but with stronger concentration near the vertical than near the horizontal in that plane,and with increasing vertical concentration with increasing depth. Grain subdivision by polygonization causes only small changes in c-axis orientation, and, as noted below,nucleation and growth of new grains with caxes at large angles to their neighbors does not become dominant until much deeper (see also Budd and Jacka, 1989;Alley, 1992). Lacking independent evidence of the orientation of core sections,this physical understanding has been applied;after rotating all sections to place the core axis in the center of the Schmidt plot,those sections with a clear symmetry plane were rotated so that these planes are parallel,making comparisons easier (Fig. 19).

Eigenvalue
Eigenvalues prov the fabric can be ized eigenvalues that S 1 > S 2 > S 3 . I S 1 ⇡ S 2 ⇡ S 3 ⇡ 1/3 in a plane would all c-axes pointin S 2 = S 3 = 0. In re Previous workers the fabric in a qua 2005;Kennedy a that they do not n such as between towards a vertical plane transverse to flow,but with stronger concentration near the vertical than near the horizontal in that plane,and with increasing vertical concentration with increasing depth. Grain subdivision by polygonization causes only small changes in c-axis orientation, and, as noted below,nucleation and growth of new grains with caxes at large angles to their neighbors does not become dominant until much deeper (see also Budd and Jacka, 1989;Alley, 1992). Lacking independent evidence of the orientation of core sections,this physical understanding has been applied;after rotating all sections to place the core axis in the center of the Schmidt plot,those sections with a clear symmetry plane were rotated so that these planes are parallel,making comparisons easier (Fig. 19).

Eigenvalues
Eigenvalues provide a m the fabric can be descri ized eigenvalues have t that S 1 > S 2 > S 3 . Ideally, S 1 ⇡ S 2 ⇡ S 3 ⇡ 1/3,a fab in a plane would have S all c-axes pointing in th S 2 = S 3 = 0. In reality, Previous workers have u the fabric in a quantitativ 2005;Kennedy and othe that they do not necessa such as between multip Fig. 19. Representative selection of horizontal Schmidt plots showing the evolution of fabric wi rotated so that all plots have the same orientation,which is assumed to be perpendicular to the dire depths (m) are indicated above each plot. Number of points (n) successfully measured by the c-axis below each plot.

Fitzpatrick and others: Physi 1194
towards a vertical plane transverse to flow,but with stronger concentration near the vertical than near the horizontal in that plane,and with increasing vertical concentration with increasing depth. Grain subdivision by polygonization causes only small changes in c-axis orientation, and, as noted below,nucleation and growth of new grains with caxes at large angles to their neighbors does not become dominant until much deeper (see also Budd and Jacka, 1989;Alley, 1992). Lacking independent evidence of the orientation of core sections,this physical understanding has been applied;after rotating all sections to place the core axis in the center of the Schmidt plot,those sections with a clear symmetry plane were rotated so that these planes are parallel,making comparisons easier (Fig. 19).

Eigenvalue
Eigenvalues prov the fabric can be ized eigenvalues that S 1 > S 2 > S 3 . I S 1 ⇡ S 2 ⇡ S 3 ⇡ 1/3 in a plane would all c-axes pointin S 2 = S 3 = 0. In re Previous workers the fabric in a qua 2005;Kennedy a that they do not n such as between Fig. 19. Representative selection of horizontal Schmidt plots showing the evolution of fa rotated so that all plots have the same orientation,which is assumed to be perpendicular to depths (m) are indicated above each plot. Number of points (n) successfully measured by th below each plot.

Fitzpatrick and othe 1194
towards a vertical plane transverse to flow,but with stronger concentration near the vertical than near the horizontal in that plane,and with increasing vertical concentration with increasing depth. Grain subdivision by polygonization causes only small changes in c-axis orientation, and, as noted below,nucleation and growth of new grains with caxes at large angles to their neighbors does not become dominant until much deeper (see also Budd and Jacka, 1989;Alley, 1992). Lacking independent evidence of the orientation of core sections,this physical understanding has been applied;after rotating all sections to place the core axis in the center of the Schmidt plot,those sections with a clear symmetry plane were rotated so that these planes are parallel,making comparisons easier (Fig. 19).  to 1500 m, ApRES measurements show a marked increase in variability that, although centred around corresponding depth values in the WDC, varied between 0.04 to 0.42. We similarly observe a sevenfold jump increase in the associated standard error, ranging from values averaging 0.006 at depths of 200-1200 m to 0.04 within the depth range of 1200-1500 m. We do not calculate E 2 − E 1 from the ApRES record beyond 1500 m due to low SNR in this range (Fig. 2a).

Competing influences between anisotropy and birefringence
The azimuthal dependency of backscattered power in ground-penetrating radar is a function of both anisotropic scattering and birefringence (Hargreaves, 1977). Although these two terms are related, they manifest from different electromagnetic phenomena. The birefringent propagation of radio waves arises from differences in dielectric permittivity along two axes perpendicular 275 to the propagation direction, the two axes often referred to as the fast and slow axes or the ordinary and extraordinary axes. On the other hand, anisotropic scattering is a consequence of changes with depth in the anisotropy of permittivity in ice crystals that may not necessarily be related to changes in crystal orientation fabric (COF) (Drews et al., 2012). Both phenomena often occur simultaneously, but our technique focuses on the analysis of the birefringent signals, which provides information on the bulk COF (Brisbourne et al., 2019).

280
Past analyses (Brisbourne et al., 2019;Fujita et al., 2006) have qualitatively inferred COF through identifying the depthperiodicity of both the node-pair minima in the co-polarised power anomaly measurements and the quadrant centres in the copolarised phase difference measurements-in other words, the more closely spaced the nodes and quadrants are in depth, the stronger the fabric asymmetry. The periodicity of these phenomena is a manifestation of the crystal birefringence as identified by the received radio-wave signal (Eq. 1). In this study, we extended the qualitative analyses suggested by these previous 285 studies to obtain quantitative measurements of the fabric asymmetry at WAIS Divide via the polarimetric coherence method (Jordan et al., 2019), and show that the results gleaned from qualitative analyses are compatible with quantitative E 2 − E 1 fabric asymmetry values. The power anomaly model that best matches observed results (Fig. 3a) incorporated a variable anisotropic parameter, which predicted an isotropic medium at depths above 200 m and below 1500 m, and an anisotropic medium between these two 290 bounds with a maximum scattering anomaly of 5 dB observed at 850 m (Fig. 3d). Drews et al. (2012) also observed azimuthal variations in backscatter power anomaly that varied similarly through depth, the variations which they attributed to microscopic (sub-metre) depth transitions in COF. The study also observed the superimposition of elongated bubbles that varied similarly with depth. While the induced polarimetric dependence through these bubbles was calculated to be minimal, where microscopic vertically-varying COF dominates the observed anisotropy, this effect was observed to be amplified at shallower 295 depths. Although we calculate COF in our analysis, we cannot confirm the mechanisms for the observed anisotropy in our results, as the COF values represent not only a bulk-depth average within the calculated depth bin but also a manifestation of the horizontal rather than vertical asymmetry. Interestingly, however, we note that the anisotropic zone parsed from the ApRES measurements (~400-1200 m) coincides almost exactly with the brittle ice zone (650-1300 m, Fitzpatrick et al., 2014), the bubble-clathrate zone (700-1250(700- m, Fitzpatrick et al., 2014, and the zone where bubble shapes were least-equant (400-1500 300 m, Fegyveresi et al., 2016). Although these correlations warrant further investigation, they are all outside the scope of our study. Rather, we focus our analysis on the observed bulk birefringence that arises as a result of the observed COF asymmetry.
Our method applies density-dependent firn correction as suggested by Jordan et al. (2020b), which effectively reduces the value of ∆ε by taking into account the birefringence of firn with respect to solid ice, which in turn increase the fabric asymmetry 13 https://doi.org/10.5194/tc-2020-264 Preprint. Discussion started: 28 September 2020 c Author(s) 2020. CC BY 4.0 License.
proportional to depth. As the quad-polarisation measurements were conducted using a ground-based monostatic antenna setup, we are able to discount the possibility of a tilt angle between the E 3 eigenvector and the direction of radio wave propagation being the source of this observed firn asymmetry (Jordan et al., 2019;Matsuoka et al., 2009). On the other hand, crystal anisotropy in snow (Calonne et al., 2017) and firn (DiPrinzio et al., 2005;Fujita et al., 2009) have previously been observed on multiple occasions, and is thought to be induced by perturbations in climate such as temperature, solar radiation, winds, and 310 deposition, where the combined effects of these variables influence the initial orientation and size of ice crystals (Kennedy et al., 2013), potentially amplifying the resulting ice flow dynamics (Wang et al., 2018). While we cannot at this point conclusively link the observed fabric asymmetry in the uppermost 200 m of the ice column to climate perturbations, it is certainly a plausible explanation especially given increasingly volatile climatic conditions over the Western Antarctic Ice Sheet over the past century that will likely intensify in the near future (Nicolas and Bromwich, 2011;Scott et al., 2019).

Flow history at WAIS Divide
With the addition of our ApRES-derived dataset, there are now three calculations of ice fabric at WAIS Divide. Each dataset was obtained using independent methods, with the other two arising from ice core observations (Fitzpatrick et al., 2014) and sonic logging (Kluskiewicz et al., 2017). A fourth method, observed from shear-wave splitting in seismic surveys, was conducted in the 2018 summer at WAIS Divide and results gleaned from this seismic experiment (Nakata et al., in prep) would likely 320 complement observations made from the corresponding datasets as an additional independent experiment.
We also note that an areal radar polarimetric study was conducted at WAIS Divide by Matsuoka et al. (2012), which quantified the relative orientation of the fabric asymmetry across a 60 km × 150 km study area, but stopped short of calculating the COF structure within their polarimetric measurements. Across the study area, they observed depth-variable azimuthal shifts that varied according to the strain regime at the corresponding age-depth period. At the WDC (their S-W24), however, their 325 results were inconclusive due to the multiple unevenly-spaced azimuthal power maxima at depth in both radar frequency returns. Our results contrast with those of Matsuoka et al. (2012) in that we observe no azimuthal ambiguity in our determination of the symmetry axis down to a depth of 1500 m, which translates to a depth-age of~7400 years BP (before 1950), which encompasses the majority of the Holocene epoch. The depth-age of the record is short because the rates of accumulation over this time period are high in comparison to both historical rates over the same area (Fudge et al., 2013) and present-day rates 330 over other areas across the Antarctic Ice Sheet (Koutnik et al., 2016). Therefore, given that our ApRES measurement location is situated within 20 km from the present-day location of WAIS Divide, our observations of a depth-invariant symmetry axis is consistent with the proposition that the ice divide has likely remained on average within 5 km of its present position throughout the Holocene (Koutnik et al., 2016). 335 Jordan et al. (2019) puts forth the advantage of using polarimetric radar methods to determine the orientation of fabric, especially with regards to the fact that ice core studies can only be conducted in a relative azimuthal reference frame due to the 14 https://doi.org/10.5194/tc-2020-264 Preprint. Discussion started: 28 September 2020 c Author(s) 2020. CC BY 4.0 License. rotational spin of the core during the drilling process. Our study confirms this proposition by establishing the orientation of the symmetry axis through the identification of depth-local minima in the cross-polarised power anomaly measurements (Fig. 2c).

Method comparisons and limitations
Although the cross-polarised power anomaly has a 90 • ambiguity (Li et al., 2018), this ambiguity can be resolved if the cor-340 responding co-polarised measurements detect anisotropic scattering (Brisbourne et al., 2019). We find that the symmetry axis orientation observed in our measurements aligns exactly with the observed present-day surface strain regime rather than the orientation of flow, at 250 • ±7 • (Fig. 2f), consistent with theory relating ice flow and crystal anisotropy (Azuma, 1994).
Of the methods available to quantify depth changes in ice COF, only ice core analyses are currently able to produce a fully three-dimensional set of fabric estimates. Thin (~10 cm) section analyses from ice cores, while providing direct orientation 345 estimates for the majority of grains within each section, are often conducted at depth intervals of tens of metres, with each analysis capturing the local decimetre-wide fabric regime, as is the case for the WAIS Divide deep ice core (WDC) at an average depth interval of 40 metres (Fitzpatrick et al., 2014). In contrast, waveform-based methods average out fabric properties in bulk where, for radar systems, the planar footprint of which the COF is averaged from is dependent on the radiation patterns of the radar antennas and increases linearly through depth, and can be approximated as the radius of the first Fresnel zone 350 (Haynes et al., 2018). As the antennas used in this study have a relatively wide half-power beamwidth (±30 • ), the footprint of the ApRES system would have a radius of approximately 1/4 of the centre wavelength (~1 m) near the surface, but reaches to approximately 1000 m at a depth of 1500 m. Therefore, the bulk COF estimates obtained from ApRES is averaged from a much larger area at depth than near the surface. Along with the SNR, this observed scale-dependence would therefore heavily influence the accuracy and error of results, especially in areas of complex flow and deformation.

355
The COF resolution vary significantly between waveform-based methods. For example, seismic surveys generate high azimuthal resolution at the cost of depth resolution (e.g., Horgan et al., 2011;Brisbourne et al., 2019), while equivalent results from sonic logging techniques reveal the converse (e.g., Gusmeroli et al., 2012;Kluskiewicz et al., 2017). While bulk-averaging over a large amount of crystals induces higher amounts of statistical noise compared to thin-section analysis, they are better able to resolve smaller scale features and discontinuities in the observed fabric (Wilen et al., 2003). Our results using a phase-360 sensitive radio echo sounder, while reconstructed from orthogonal measurements, offer a compromise between azimuthal and depth precision, providing COF estimates comparable to results from the WDC at resolutions comparable to that of seismics in azimuth and sonic logging in depth.
The use of quad-polarised measurements is depth-limited by the signal-to-noise ratio of the cross-polarised terms (s hv and s vh ). Our COF measurements obtained from polarimetric radar, where chirps were coherently summed during pre-processing 365 for each measurement, correlate well with equivalent ice-core measurements made at WDC to a depth of 1500 m, after which the phase is dominated by noise. Given that the relative power of the cross-polarised terms are almost a magnitude lower than that of the co-polarised terms (Fig. 2a), accurate COF calculations using ApRES most likely can not exceed this depth limit with the current specifications.
As this study and those by Jordan et al. (2019Jordan et al. ( , 2020b show, the polarimetric coherence method measures the horizontal 370 asymmetry of the vertical ice column through quantifying the birefringence effects from orthogonally-oriented measurements and relating this to the difference in magnitude between the E 2 and E 1 eigenvectors. As such, an obvious limitation of this 15 https://doi.org/10.5194/tc-2020-264 Preprint. Discussion started: 28 September 2020 c Author(s) 2020. CC BY 4.0 License. method is its inability to discern azimuthally-invariant fabric such as single-cluster crystal distributions or more complex fabrics such as horizontal girdles from each other. However, given the high resolution of the results, ApRES-derived COF are directly complementary to results from sonic logging in that the former calculates the horizontal asymmetry of the fabric through 375 quantification of E 2 −E 1 , whilst the latter quantifies the vertical fabric asymmetry (the strength of the E 3 eigenvector) through P-wave interpretation (Kluskiewicz et al., 2017).
An obvious advantage of seismic and radar surveys, in comparison to ice coring and sonic logging, is that they can be implemented as a much smaller operation in terms of team size, cost, and field time, with ground-based radar surveys an order of magnitude even lower than that of seismics and airborne radar in all three aspects. As a reference, the quad-polarised 380 measurements collected in this study took approximately thirty minutes inclusive of radar and antenna setup, and involved only two persons (co-authors Young and Dawson). While radar surveys can be conducted efficiently in comparison to that of seismics surveys, the latter is able to estimate COF in three dimensions (albeit with a vertically-integrated value) and can distinguish azimuthally-symmetric forms of anisotropy that the former at present cannot. A combination of these two methods can therefore reduce ambiguity in fabric estimation (Brisbourne et al., 2019), especially in areas of complex flow which have 385 the potential to produce complex fabric. Additionally, these methods have the ability to measure COF in a wide variety of flow regimes that may not be suitable as an ice core drilling site. In this respect, the use of one or a combination of geophysical methods has the potential to investigate more dynamic areas that reveal the complexities of ice flow, such as shear margins or grounding zones. In conclusion, as shown in our study, the use of an ApRES in conjunction with the polarimetric coherence method can produce estimates of fabric asymmetry with accuracies comparable to and resolutions exceeding that of thin ice 390 core section analyses.

Conclusions
Using a phase-sensitive radar and two open-structure antennas, we conducted a suite of quadrature-polarimetric measurements within the proximity of the WAIS Divide Ice Core (WDC). Using a combination of the matrix-based backscatter model (Brisbourne et al., 2019;Fujita et al., 2006) and the polarimetric coherence method (Jordan et al., 2019(Jordan et al., , 2020b, we were able to 395 (i) quantify the horizontal asymmetry (E 2 − E 1 ) of the crystal orientation fabric (COF) to a depth of 1500 m; and (ii) unambiguously identify the fabric asymmetry to be depth invariant with the symmetry axis oriented at 250 • ±7 • to the same depth of 1500 m. Our findings in (i) were conducted at an angular and vertical resolution of 1 • and 15 m respectively, exceeding that of ice core-derived measurements made at the WDC (Fitzpatrick et al., 2014). The correlation between these two independent measurements of fabric asymmetry is moderately high over the depth range of 600-1200 m (r 2 = 0.55), and is slightly lower 400 over the entire depth range to 1500 m (r 2 = 0.45). Our findings in (ii) align exactly with the direction of principal strain independently measured by Matsuoka et al. (2012). Our determination of depth-invariant fabric orientation to at least 1500 m, equivalent to 7400 years BP (years before 1950), covers~65% of the Holocene epoch, which suggests that the deformation regime at WAIS Divide has not changed substantially during this period. While ice core-based measurements still represent the only method that is able to estimate COF in three dimensions, the logistics required to conduct polarimetric radar mea-