Ultrasonic and seismic constraints on crystallographic preferred orientations of the Priestley Glacier shear margin, Antarctica

Abstract

Crystallographic preferred orientations (CPOs) are particularly important in controlling the mechanical properties of glacial shear margins. Logistical and safety considerations often make direct sampling of shear margins difficult, and geophysical measurements are commonly used to constrain the CPOs. We present here the first direct comparison of seismic and ultrasonic data with measured CPOs in a polar shear margin. The measured CPO from ice samples from a $\mathrm{58}$ m deep borehole in the left lateral shear margin of the Priestley Glacier, Antarctica, is dominated by horizontal $c$ axes aligned sub-perpendicularly to flow. A vertical-seismic-profile experiment with hammer shots up to $\mathrm{50}$ m away from the borehole, in four different azimuthal directions, shows velocity anisotropy of both P waves and S waves. Matching P-wave data to the anisotropy corresponding to CPO models defined by horizontally aligned $c$ axes gives two possible solutions for the $c$-axis azimuth, one of which matches the $c$-axis measurements. If both P-wave and S-wave data are used, there is one best fit for the azimuth and intensity of $c$-axis alignment that matches the measurements well. Azimuthal P-wave and S-wave ultrasonic data recorded in the laboratory on the ice core show clear anisotropy of P-wave and S-wave velocities in the horizontal plane that match that predicted from the CPO of the samples. With quality data, azimuthal increments of $\mathrm{30}{}^{\circ }$ or less will constrain well the orientation and intensity of $c$-axis alignment. Our experiments provide a good framework for planning seismic surveys aimed at constraining the anisotropy of shear margins.

1 Introduction

Ice streams and glaciers are localised regions of high ice flow velocity inside otherwise mostly stationary ice masses of Antarctica and Greenland and play a key role in the mass balance of polar ice masses. As a result of high-velocity flow, the margins of these streaming ice bodies undergo strain as they are in contact with stationary ice or rock. Crystallographic preferred orientation (CPO) patterns observed inside glacier margins indicate a very high degree of crystal alignment in the horizontal direction . These results are consistent with observations from shear deformation experiments of ice where $c$-axis maxima are always aligned perpendicularly to the shear plane .

The presence of a CPO results in anisotropic mechanical properties and so influences the viscous behaviour of ice significantly . Shear margins of glaciers therefore can affect the character of ice flow in ice streams due to their distinct mechanical properties . The advection of the shear margins during the flow of ice from land to sea can result in bands of strongly deformed ice that transect ice shelves and can potentially affect the stability of ice shelves. Modelling suggests that remnant CPO resultant in shear margins can be present 10 000 years after advection downstream .

A better understanding of CPO patterns in glacier shear margins is therefore highly desirable to accurately determine their mechanical properties. Ice core drilling, the primary direct information source for CPO, is however rarely performed on fast-flowing ice because of difficulties in access and on-site safety. Geophysical studies, e.g. seismic or radar surveys, provide an alternative way of constraining bulk CPO. Ideally, geophysical work should be combined with drilling to recover ice samples for microstructure analysis and so enable a cross-calibration of CPO constraints.

A continuous ice core of $\mathrm{58}$ m length was drilled and recovered in December 2019 and January 2020 in a lateral shear margin of the Priestley Glacier, located in Victoria Land, Antarctica; see Fig. a. There is no firn layer on the glacier; however a snow cover of approximately $\mathrm{0.4}$ m was present at the time of drilling. At the site the glacier is $\sim \mathrm{1000}$ m thick and $\mathrm{8}$ km wide . Ice flow velocities are measured to increase from $\sim \mathrm{45}$ m yr${}^{-\mathrm{1}}$ close to the glacier margin to $\sim \mathrm{130}$ m yr${}^{-\mathrm{1}}$ towards the centre , resulting in strong shear, with shear strain rates at the site calculated to be $\mathrm{6}\times {\mathrm{10}}^{-\mathrm{10}}$ s${}^{-\mathrm{1}}$ . Core samples were analysed for CPO using electron backscatter diffraction (EBSD) measurements , and a strong horizontal clustering of $c$ axes was observed throughout the entire length of the core. The core orientation was carefully preserved during drilling, which enabled azimuthal orientation of the CPO. The horizontal $c$-axis cluster orientation varies between approximately perpendicular to ice flow to $\mathrm{130}{}^{\circ }$ clockwise (looking down) of the ice flow direction.

Figure 1

(a) Regional map of the Priestley Glacier site. Coordinates are given by the polar stereographic projection with latitude of true scale at $-\mathrm{71}{}^{\circ }$ S. Landsat image courtesy of USGS. (b) Map view illustration of VSP geometry. Black markers show shotpoint locations; the borehole location is shown in blue. (c) Cross-section illustration of VSP geometry showing the surface shotpoints in black and three borehole seismometer positions in blue. Lines connecting sources and receivers indicate seismic ray paths. (d) Vertical component traces of borehole seismometer for shots at $\mathrm{50}$ m offset along the $\mathrm{45}{}^{\circ }$ profile. Example P- and S-wave arrival time picks are shown at $\mathrm{25}$ m depth as vertical blue and red lines.

[Figure omitted. See PDF]

After completion of drilling and core retrieval, the open borehole was used to conduct a vertical-seismic-profile (VSP) experiment to constrain seismic properties of the near-surface glacier ice, with a particular focus on seismic anisotropy. To complete the link between seismic anisotropy of the ice volume around the borehole and CPO measurements from the core, multi-azimuthal ultrasonic velocity measurements were made on core samples in the laboratory. The anisotropy of single ice crystals displays hexagonal symmetry; i.e. wave velocity is only dependent on the angle between the propagation direction and $c$ axis, with the fast P-wave propagation direction parallel to the $c$ axis . Seismic P- and S-wave velocities at different propagation angles relative to the $c$ axis of a single crystal are shown in Fig. .

Figure 2

Seismic phase velocities for a single crystal dependent on angle $\mathit{\theta }$ to the $c$ axis, calculated with values for the elastic constants from adiabatic single-crystal artificial water ice at $-\mathrm{16}$ ${}^{\circ }$C from using the MTEX toolbox . (a) P-wave velocity $v_{\mathrm{P}}$. (b) S-wave velocities $v_{\mathrm{S}}$.

[Figure omitted. See PDF]

This is the first study of the horizontal-cluster CPO type observed in shear margins with seismic methods. Analyses of seismic, ultrasonic and measured CPO datasets are combined to assess the potential of active-source seismic surveys for the constraint of shear margin anisotropy.

2 Analysis of a vertical seismic profile

2.1 Data acquisition

A VSP dataset was recorded at the Priestley drill site using a three-component borehole seismometer (built by ESS Earth Sciences, Victoria, Australia) with a pneumatic clamping system which was installed at depths $z$ between $\mathrm{5}$ and $\mathrm{57}$ m below the glacier surface. The suspension cable, to the centre of the seismometer, remained close to the centre of the hole at the surface for all seismometer depths, indicating true verticality of the borehole. Signals from a hammer-and-plate source were recorded using a Geometrics Geode for a walkaway-VSP geometry along four profiles with different azimuths, where shots were placed at offsets $x$ of $\mathrm{0}$ m (“zero-offset” geometry), $\mathrm{10}$, $\mathrm{30}$ and $\mathrm{50}$ m from the top of the borehole. Profile names (“Flow”, “Perp”, “$\mathrm{45}{}^{\circ }$”, “$-\mathrm{45}{}^{\circ }$”) indicate their orientation relative to the glacier flow direction (see Table ). The survey geometry is illustrated in Fig. b. The depth increments of the borehole seismometer along the four profiles are shown in Table .

Table 1

Multi-azimuth VSP survey parameters.

Shot profile orientation	Seismometer depth
(profile name)	increment (m)
Zero offset	$\mathrm{1}$
Flow-parallel (Flow)	$\mathrm{1}$
Flow-perpendicular (Perp)	$\mathrm{2.5}$
$\mathrm{45}{}^{\circ }$ to flow	$\mathrm{5}$
$-\mathrm{45}{}^{\circ }$ to flow	$\mathrm{5}$

Geode data recording for each shot in the field was initiated by a Geometrics switch trigger taped to the sledgehammer handle. Quality control of traces in the field found that this trigger type produces inconsistent zero times; i.e. repeat shots for a given source–receiver pair exhibit different arrival times. To enable the determination of absolute velocities, the recorded signal from surface geophones, collocated at the shotpoints, was used to define shot times. The Geode was set to record $\mathrm{100}$ ms before the hammer switch trigger signal to enable recording of the full first-arrival signal on the shotpoint surface geophone. Manual picking of the first-arrival time on the shotpoint surface geophone trace gives the time of hammer impact and so the true source time. Total recording length was $\mathrm{2}$ s with a sampling rate of $\mathrm{8000}$ Hz. Shot and geophone locations were cleared of all snow cover to ensure direct contact of source and receivers with the hard ice surface. Several repeat shots were recorded at all shotpoints along a given profile for one seismometer depth before the instrument was lowered to the next depth. After all shots for all seismometer depths were completed for a profile, surface geophones were moved to the next profile position.

Polarisation patterns indicate a ringing effect of the pneumatic borehole seismometer, where phase arrivals are followed by a tail of mono-frequent oscillations (see traces in Fig. d). The oscillations are distributed along all three seismometer components, even after separation of P- and S-wave signals through component rotation into ray coordinates . This noise in the phase arrival signals impedes investigating polarisation patterns, such as S-wave splitting, as a constraint on seismic anisotropy .

2.2 Multi-azimuth VSP travel times

P- and S-wave first-arrival signals were recorded with a high signal-to-noise ratio, and phase arrival times can be clearly identified. Picking of seismic phase arrivals is performed manually for each shot to determine travel times: one first-arrival pick is made on the surface geophone trace at the shot location in addition to picks of P- and S-wave arrivals on the borehole seismometer traces. The P- and S-wave travel time $t$ is then calculated as the difference between the picked arrival time on the borehole seismometer and the picked source time on the surface geophone.

Observed P-wave travel times along the seismic profiles are presented in Fig. a–d. The zero-offset P-wave travel times in Fig. a can be approximated by a constant velocity model ($v_{\mathrm{P}}=\mathrm{3850}$ m s${}^{-\mathrm{1}}$, solid line). This highlights that there is little velocity variation due to heterogeneity in the shallow glacier ice.

Differences in travel times between the profiles become apparent for shots at offsets $\ge \mathrm{30}$ m and below seismometer depths of $\sim \mathrm{25}$ m. The Flow and $-\mathrm{45}{}^{\circ }$ profiles show consistently earlier P-wave arrivals than the Perp and $\mathrm{45}{}^{\circ }$ profiles.

Observed S-wave travel times are presented in Fig. f–h for shots from the four walkaway VSPs. Signal-to-noise ratios from zero-offset shots were found to be insufficient to allow picking of S-wave arrivals, likely a consequence of the radiation patterns of a plate seismic source which produces no S-wave energy in the vertical direction. The first indication of S-wave arrivals was picked consistently; travel times are therefore interpreted to be representative of the fast S-wave velocity $v_{\mathrm{S}\mathrm{1}}$. Systematic travel time differences between profiles are observed particularly for shots from a $\mathrm{50}$ m offset, where the $-\mathrm{45}{}^{\circ }$ profile shows the earliest arrival times over almost the entire depth range. Only for seismometer depths below $\mathrm{50}$ m are travel times from the Perp profile measured at the same time as or earlier than travel times from along the $-\mathrm{45}{}^{\circ }$ profile. The Flow and $\mathrm{45}{}^{\circ }$ profiles show similar travel times that are slower than the travel times along the other profiles.

Figure 3

VSP travel times observed along the four shot profiles. (a) P-wave travel times for shots from $\mathrm{0}$ m offset. (b) P-wave travel times for shots from $\mathrm{10}$ m offset. (c) P-wave travel times for shots from $\mathrm{30}$ m offset. (d) P-wave travel times for shots from $\mathrm{50}$ m offset. (e) No S-wave travel times were picked from $\mathrm{0}$ m offset shots due to a poor signal-to-noise ratio. (f) S-wave travel times for shots from $\mathrm{10}$ m offset. (g) S-wave travel times for shots from $\mathrm{30}$ m offset. (h) S-wave travel times for shots from $\mathrm{50}$ m offset.

[Figure omitted. See PDF]

2.3 Velocity calculation

Seismic velocities are calculated from the presented travel times for different depths $z$, offsets $x$ and azimuths in Fig. :

1$${\displaystyle }v={\displaystyle \frac{\sqrt{x^{\mathrm{2}}+z^{\mathrm{2}}}}{t}}.$$ The calculation is based on the assumption of straight ray paths between the shot and receiver, which is regarded as a valid approximation for a site location on hard ice without snow cover. The seismic waves travel entirely in ice, and no large velocity gradients that could result in significant ray path bending are therefore expected . The zero-offset travel times in Fig. a exhibit no apparent indications of velocity gradients since they can be approximated by a constant vertical P-wave velocity. We therefore assume the ice volume sampled by the VSP experiment to be characterised by a homogeneous, constant-CPO layer. The incidence angle $\mathit{\theta }$ from the vertical direction is defined as 2$${\displaystyle }\mathit{\theta }=\arctan \left({\displaystyle \frac{x}{z}}\right).$$

Velocity uncertainties $\mathrm{\Delta }v$ are calculated from uncertainty estimates for P-wave travel time $\mathrm{\Delta }{t}_{\mathrm{P}}=\mathrm{0.125}$ ms, S-wave travel time $\mathrm{\Delta }{t}_{\mathrm{S}}=\mathrm{1}$ ms, offset $\mathrm{\Delta }x=\mathrm{0.2}$ m and depth $\mathrm{\Delta }z=\mathrm{0.1}$ m. 3$$\mathrm{\Delta }v={\displaystyle \frac{x\mathrm{\Delta }x+z\mathrm{\Delta }z}{t\cdot \sqrt{x^{\mathrm{2}}+z^{\mathrm{2}}}}}+{\displaystyle \frac{\mathrm{\Delta }t\cdot \sqrt{x^{\mathrm{2}}+z^{\mathrm{2}}}}{t^{\mathrm{2}}}}$$ Observations with relative uncertainty $\mathrm{\Delta }v/v>\mathrm{0.05}$ are discarded from further analysis.

3 Analysis of ultrasonic experiments

Ultrasonic experiments were performed inside a freezer at the temperature $-\mathrm{23}\pm \mathrm{2}$ ${}^{\circ }$C on a subset of samples from the continuous ice core. Samples were lathed to the shape of a highly regular cylinder. Travel time measurements were made perpendicularly to the core axis in multiple azimuths to measure ultrasonic velocities in the horizontal direction of the glacier ice. Olympus ultrasonic P- or S-wave transducers were spring-loaded against the cylinder surface, on opposing sides of the cylinder (Fig. a and b). S-wave transducers were used for excitation and recording of S waves with polarisations in the vertical (parallel to the long core axis) and horizontal (perpendicular to the long core axis) direction. Coupling of transducers to the ice core samples was ensured through the use of synthetic high-performance low-temperature grease.

Travel time measurements across the ice core were made in azimuthal increments of $\mathrm{10}{}^{\circ }$ covering the full core diameter, resulting in $N=\mathrm{36}$ measurements per transducer arrangement and core. A fiducial line that was made in the field to provide geographic reference of the core orientation served as the $\mathrm{0}{}^{\circ }$ reference on individual samples. The fiducial line was made perpendicularly to the glacier flow direction on the surface . New fiducial lines were started wherever a core break could not be fitted together, and the relative orientations of the lines were reconstructed using the core CPO . Measurements are on core sections $\mathrm{003}$ from a depth of $\sim \mathrm{2.5}$ m (diameter $d=\mathrm{100.6}\pm \mathrm{0.3}$ mm), $\mathrm{007}$ from a depth of $\sim \mathrm{6.0}$ m (diameter $d=\mathrm{99.9}\pm \mathrm{0.3}$ mm) and $\mathrm{010}$ from a depth of $\sim \mathrm{8.5}$ m (diameter $d=\mathrm{100.3}\pm \mathrm{0.3}$ mm). The in situ temperature at the sample depths was observed to be between $\approx -\mathrm{7}$ and $\approx -\mathrm{15}$ ${}^{\circ }$C . Ultrasonic measurements made at these warmer temperatures did not result in any measurable effect on seismic anisotropy. Since warmer temperatures caused problems with transducer coupling to the ice cores, the initially described freezer temperature setting at $-\mathrm{23}\pm \mathrm{2}$ ${}^{\circ }$C was used for the ultrasonic experiments.

The ultrasonic source signal pulse was created by a JSR Ultrasonics DPR300 pulser unit and shows a dominant frequency $f\approx \mathrm{1}$ MHz, resulting in a dominant wavelength of $\mathit{\lambda }\approx \mathrm{3.8}$ mm in the ice. Recording of the source and receiver signal was performed on separate oscilloscope channels by a PicoScope digital oscilloscope with a sampling interval $t_{\mathrm{S}}=\mathrm{0.8}$ ns. Signals were recorded directly by the oscilloscope without the use of amplifiers. The source signal was used to trigger signal recording by exceedance of an amplitude threshold. A signal length of $\mathrm{100}$ $\mathrm{\mu }\mathrm{s}$ after triggering is recorded.

3.1 Processing

At each azimuth $\mathrm{64}$ individual waveforms were recorded. The mean amplitude and linear trend are removed from these individual traces before they are stacked to increase the signal-to-noise ratio. The stacked traces are tapered and filtered with a bandpass filter with corner frequencies $f_{\mathrm{1}}$ $=$ 50 KHz and $f_{\mathrm{2}}=\mathrm{5}$ MHz and normalised for plotting.

Waveforms recorded on sample $\mathrm{007}$ are shown in Fig. c–e. The core orientation angle in these plots is relative to the fiducial marker on the core sample. Recorded waveforms using the P-wave transducers are shown in Fig. c. A rotational symmetry with periodicity of $\mathrm{180}{}^{\circ }$ is present, which confirms a consistent and reliable excitation and recording of P-wave signals for all azimuths around the core.

Waveforms recorded using S-wave transducers set for vertical vibration are shown in Fig. d, indicating a clear periodicity of $\mathrm{180}{}^{\circ }$ and are predominantly indicative of the slow S-wave velocity $v_{\mathrm{S}\mathrm{2}}$.

Waveforms from S-wave transducers set for horizontal vibration are shown in Fig. e. Here, four local maxima in velocity are observed for the full range of azimuths. This arrangement of transducers enables the measurement of fast S-wave velocity $v_{\mathrm{S}\mathrm{1}}$ at most azimuths.

Figure 4

(a) Side view of ultrasonic rig with slots for different transducer orientations highlighted. (b) Map view of ultrasonic rig. Transducers are set for vertical vibration. (c) Multi-azimuth signal traces using P-wave transducers on sample $\mathrm{007}$. (d) Multi-azimuth signal traces using S-wave transducers set for vertical vibration on sample $\mathrm{007}$. (e) Multi-azimuth signal traces using S-wave transducers set for horizontal vibration on sample $\mathrm{007}$.

[Figure omitted. See PDF]

3.2 P-wave and S-wave ultrasonic velocity results

Arrival times $t$ of the direct P- and S-wave phases are picked by hand on the stacked traces and used to calculate ultrasonic velocities $v=d/t$. Successful measurements of $v_{\mathrm{P}}$, $v_{\mathrm{S}\mathrm{1}}$ and $v_{\mathrm{S}\mathrm{2}}$ are made at all azimuths and samples, which constitute an unprecedented characterisation of ultrasonic velocity anisotropy in high detail. Velocity uncertainties are calculated from diameter uncertainty $\mathrm{\Delta }d=\mathrm{0.3}$ mm and picking uncertainty $\mathrm{\Delta }t=\mathrm{0.1}$ $\mathrm{\mu }\mathrm{s}$.

Ultrasonic P- and S-wave velocities of the three core sections studied are shown in Fig. a–c. The velocities are presented in the kinematic reference frame, where the angles relative to the fiducial lines on individual samples are transformed to give the orientation relative to glacier flow at the site . The azimuth given in the plots presents the angle from the flow-parallel VSP in Fig. b in a clockwise direction.

The largest P- and S-wave velocities are observed in sample $\mathrm{007}$; sample $\mathrm{003}$ shows the lowest velocities, and intermediate velocities are measured in sample $\mathrm{010}$. A variation in seismic velocities with the azimuth is present in all samples. Shapes of the velocity curves are very similar for the three individual samples; however the observed velocities of $\mathrm{003}$ are shifted by an angle of $\approx \mathrm{20}{}^{\circ }$ relative to the curves of 007 and 010.

Ultrasonic $v_{\mathrm{P}}$ measurements were also made along the core axis, sampling the vertical direction at the temperature $-\mathrm{23}\pm \mathrm{2}$ ${}^{\circ }$C. Sample velocities of $v_{\mathrm{P},\mathrm{vert},\mathrm{003}}$ $=$ 3838 $\pm$ 20 m s${}^{-\mathrm{1}}$, $v_{\mathrm{P},\mathrm{vert},\mathrm{007}}$ $=$ 3822 $\pm$ 15 m s${}^{-\mathrm{1}}$ and $v_{\mathrm{P},\mathrm{vert},\mathrm{010}}=\mathrm{3842}\pm \mathrm{18}$ m s${}^{-\mathrm{1}}$ again highlight near-constant velocity with depth implicit in the zero-offset VSP P-wave travel times (Fig. a). The range of vertical $v_{\mathrm{P}}$ is shown in Fig. a. The observation that vertical $v_{\mathrm{P}}$ is faster than horizontal $v_{\mathrm{P}}$ for most azimuths is inconsistent with the expected velocities for a horizontal-cluster CPO, and no definite explanation for this observation is available at this point. Potential additional influences on velocities could be given by populations of “oddly” oriented grains that show orientations in the horizontal plane but outside of the $c$-axis cluster. They could therefore reduce the maximum $v_{\mathrm{P}}$ along the cluster axis but not affect $v_{\mathrm{P}}$ in the vertical direction. Anisotropy related to preferred bubble shape or aligned cracks, both observed in the uppermost $\approx \mathrm{10}$ m of the ice core , could also play a role.

An attempted measurement of $v_{\mathrm{S}}$ in the vertical direction failed because no clear ultrasonic S-wave signals that allowed travel time picking could be recorded in this geometry.

Figure 5

Ultrasonic multi-azimuth velocities for different core sections relative to the macroscopic glacier flow direction ($\mathrm{0}{}^{\circ }$). (a) P-wave velocity $v_{\mathrm{P}}$. The grey-shaded box shows the range of measured vertical ultrasonic $v_{\mathrm{P}}$. (b) Fast S-wave velocity $v_{\mathrm{S}\mathrm{1}}$. (c) Slow S-wave velocity $v_{\mathrm{S}\mathrm{2}}$.

[Figure omitted. See PDF]

4 CPO modelling

The observed high degree of seismic anisotropy in VSP seismic and multi-azimuth ultrasonic data is consistent with the observation of strong CPO in the retrieved core samples from the site. EBSD measurements on core samples constrain the CPO to be characterised by a strong clustering of $c$ axes in the horizontal plane . The availability of constraints from sample microscopic analysis and geophysical data at the Priestley Glacier site is therefore suitable for a calibration of seismic properties to the known CPO. The observed P- and S-wave velocity anisotropies from VSP and multi-azimuth ultrasonic observations are compared to models of polycrystal elasticity connected to CPOs comprising horizontally clustered $c$ axes. In general, when interpreting geophysical data, no CPO constraints are available and a range of different CPO geometries have to be tested to match observations . We skip this step and focus entirely on the influence of survey geometry and suitability of P- and S-wave velocity information for constraining ice CPO.

4.1 Horizontal-cluster CPO

CPO models are created using the MTEX toolbox for MATLAB , consisting of $\mathrm{1000}$ clustered $c$ axes with the cluster symmetry axis in the horizontal direction. A geometrical description of this horizontal-cluster CPO forward model can be given by a diffused cluster with the symmetry axis in the horizontal plane. The model parameters of this cluster geometry are defined by the azimuth $\mathit{\phi }$ of the cluster symmetry axis and the cluster width angle $\mathit{\alpha }$. An illustration of this CPO type is presented in Fig.  in an upper-hemisphere stereographic projection.

Figure 6

Upper-hemisphere stereographic projection $c$-axis orientations of an example horizontal-cluster CPO model with illustration of model parameters: cluster orientation $\mathit{\phi }$ and opening angle $\mathit{\alpha }$.

[Figure omitted. See PDF]

Seismic properties of individual crystals are characterised by the elasticity tensor $C$ in Eq. () of synthetic ice at $-\mathrm{16}$ ${}^{\circ }$C by and density $\mathit{\rho }=\mathrm{919.1}$ kg m${}^{-\mathrm{3}}$. For any given set of CPO model parameters $\mathit{\phi }$ and $\mathit{\alpha }$, synthetic seismic velocities associated with the CPO are calculated after forming the Voigt–Reuss–Hill average elastic properties . 4$$\begin{array}{rl}& C=\\ & \left[\begin{array}{cccccc}\mathrm{13.913}& \mathrm{7.026}& \mathrm{5.801}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{7.026}& \mathrm{13.913}& \mathrm{5.801}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{5.801}& \mathrm{5.801}& \mathrm{15.059}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{3.011}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{3.011}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{3.4435}\end{array}\right]\text{ GPa}\end{array}$$

Model parameters $\mathit{\phi }$ and $\mathit{\alpha }$ are systematically varied in regular step sizes, and all combinations within the individual parameter ranges presented in Table  are explored. Synthetic seismic properties are compared to observations for each model realisation.

Table 2

Horizontal-cluster CPO model parameters.

Parameter	Start value	End value	Step size
Cluster orientation $\mathit{\phi }$	$\mathrm{0}{}^{\circ }$	$\mathrm{180}{}^{\circ }$	$\mathrm{5}{}^{\circ }$
Opening angle $\mathit{\alpha }$	$\mathrm{0.1}{}^{\circ }$	$\mathrm{40.1}{}^{\circ }$	$\mathrm{1}{}^{\circ }$

4.2 Model misfit calculation

The ability of a CPO model to explain measured seismic anisotropy is assessed by introduction of a misfit between synthetic forward-modelled seismic properties and observations. We found that CPO-predicted seismic velocities are consistently faster than observed velocities, an effect which can be attributed to the absence of grain boundary effects or air bubbles in synthetic CPOs. The elasticity tensor by was derived using oscillations in gigahertz (GHz) frequencies; therefore dispersion is another potential factor to introduce differences between modelled and observed seismic or ultrasonic velocities with hertz (Hz) to kilohertz (kHz) frequencies. A normalisation is introduced to address these described effects on absolute velocities, and we focus entirely on modelled and observed anisotropic velocity variations rather than absolute velocities, studying the variation in velocities relative to the mean velocity. This definition of seismic anisotropy $\mathit{\delta }v$ is given in Eq. (), where $\overline{v}$ is the mean wave velocity of all observations.

5$$\mathit{\delta }v=\left({\displaystyle \frac{v}{\overline{v}}}-\mathrm{1}\right)\cdot \mathrm{100}$$

We consider the elasticity data from to provide high accuracy of relative velocity variation. Therefore we have chosen to rely upon these data, and no other published elasticities of ice are investigated to match the observed anisotropy as has been done in previous studies which have compared absolute measured and synthetic velocities . A misfit ${\mathit{\chi }}_{\mathrm{raw}}$ is calculated by Eq. () between measured $\left(\mathit{\delta }v\right)$ and synthetic $\left(\mathit{\delta }{v}_M\right)$ velocity variations. $N$ is the number of velocity observations. 6$${\mathit{\chi }}_{\mathrm{raw}}\left(v\right)=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N{\left(\mathit{\delta }{v}_i-\mathit{\delta }{v}_{M,i}\right)}^{\mathrm{2}}}$$

Misfits ${\mathit{\chi }}_{\mathrm{raw}}\left(v\right)$ are calculated individually for the $v_{\mathrm{P}}$, $v_{\mathrm{S}\mathrm{1}}$ and $v_{\mathrm{S}\mathrm{2}}$ observations and all synthetic CPO models. The $\mathit{\chi }\left(v_{\mathrm{P}}\right)$, $\mathit{\chi }\left(v_{\mathrm{S}\mathrm{1}}\right)$ and $\mathit{\chi }\left(v_{\mathrm{S}\mathrm{2}}\right)$ misfits are then calculated through division by the largest misfit of all models. 7$$\mathit{\chi }\left(v\right)={\mathit{\chi }}_{\mathrm{raw}}\left(v\right)/\text{max}\left({\mathit{\chi }}_{\mathrm{raw}}\right(v\left)\right)$$ In addition, misfit averages are defined as $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\mathit{\chi }\left(v_{\mathrm{P}}\right)+\mathit{\chi }\left(v_{\mathrm{S}\mathrm{1}}\right)=\\ \text{8}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{{\mathit{\chi }}_{\mathrm{raw}}\left(v_{\mathrm{P}}\right)}{\text{max}\left({\mathit{\chi }}_{\mathrm{raw}}\right(v_{\mathrm{P}}\left)\right)}}+{\displaystyle \frac{{\mathit{\chi }}_{\mathrm{raw}}\left(v_{\mathrm{S}\mathrm{1}}\right)}{\text{max}\left({\mathit{\chi }}_{\mathrm{raw}}\right(v_{\mathrm{S}\mathrm{1}}\left)\right)}}\right),{\displaystyle }& {\displaystyle }\mathit{\chi }\left(v_{\mathrm{P}}\right)+\mathit{\chi }\left(v_{\mathrm{S}\mathrm{1}}\right)+\mathit{\chi }\left(v_{\mathrm{S}\mathrm{2}}\right)=\\ \text{9}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\left({\displaystyle \frac{{\mathit{\chi }}_{\mathrm{raw}}\left(v_{\mathrm{P}}\right)}{\text{max}\left({\mathit{\chi }}_{\mathrm{raw}}\right(v_{\mathrm{P}}\left)\right)}}+{\displaystyle \frac{{\mathit{\chi }}_{\mathrm{raw}}\left(v_{\mathrm{S}\mathrm{1}}\right)}{\text{max}\left({\mathit{\chi }}_{\mathrm{raw}}\right(v_{\mathrm{S}\mathrm{1}}\left)\right)}}+{\displaystyle \frac{{\mathit{\chi }}_{\mathrm{raw}}\left(v_{\mathrm{S}\mathrm{2}}\right)}{\text{max}\left({\mathit{\chi }}_{\mathrm{raw}}\right(v_{\mathrm{S}\mathrm{2}}\left)\right)}}\right).\end{array}$$

4.3 CPO model fit using VSP velocities

The model misfit of the VSP data is calculated using $N_{\mathrm{P}}=\mathrm{1354}$ and $N_{\mathrm{S}\mathrm{1}}=\mathrm{1036}$ velocity observations. Figure a shows the model misfits based on P-wave velocities $\mathit{\chi }\left(v_{\mathrm{P}}\right)$. The best-fitting model is indicated by the red dot. This model CPO is shown in Fig. d. Figure b and e show S-wave misfits $\mathit{\chi }\left(v_{\mathrm{S}\mathrm{1}}\right)$ and the corresponding best-fitting CPO model. The misfit average $\mathit{\chi }\left(v_{\mathrm{P}}\right)+\mathit{\chi }(v_{\mathrm{S}\mathrm{1}}$) is presented in Fig. c, with its best-fitting CPO model presented in Fig. f. Parameters of the best-fitting CPO models based on VSP observations are provided in Table . Parameter uncertainties are taken from the range of models that produce the minimum $\mathrm{1}$ % of misfits.

Figure 7

(a) Misfit surface considering only P-wave velocities. (b) Misfit surface considering only S-wave velocities. (c) Misfit surface considering P- and S-wave velocities. (d) Best-fitting CPO model result considering only P-wave velocities. Contour lines indicate $c$-axis density in an upper-hemisphere stereographic projection. Labels indicate the orientation of the Flow and Perp profile shots in the VSP survey (compare Fig. b). (e) Best-fitting CPO model result considering only S-wave velocities. (f) Best-fitting CPO model result considering P- and S-wave velocities.

[Figure omitted. See PDF]

Table 3

Best-fitting CPO model parameters and uncertainties informed by VSP data.

Misfit	Cluster orientation	Cluster width	Misfit
type	$\mathit{\phi }$ (${}^{\circ }$)	$\mathit{\alpha }$ (${}^{\circ }$)	[1 % misfit]
$\mathit{\chi }\left(v_{\mathrm{P}}\right)$	$\mathrm{65}\begin{array}{c}+\mathrm{5}\\ -\mathrm{0}\end{array}$	$\mathrm{2.1}\begin{array}{c}+\mathrm{7}\\ -\mathrm{2}\end{array}$	0.69 [0.70]
$\mathit{\chi }\left(v_{\mathrm{S}\mathrm{1}}\right)$	$\mathrm{110}\begin{array}{c}+\mathrm{5}\\ -\mathrm{5}\end{array}$	$\mathrm{22.1}\begin{array}{c}+\mathrm{4}\\ -\mathrm{3}\end{array}$	0.58 [0.59]
$\sum \mathit{\chi }$	$\mathrm{115}\begin{array}{c}+\mathrm{0}\\ -\mathrm{5}\end{array}$	$\mathrm{18.1}\begin{array}{c}+\mathrm{6}\\ -\mathrm{2}\end{array}$	0.66 [0.66]

Figure shows the observed P-wave anisotropy $\mathit{\delta }{v}_{\mathrm{P}}$ (symbols) along the four survey profiles together with predicted anisotropy from the CPO models presented in Table  (dashed lines). The observed variation in velocities with incidence angle is consistent across all seismometer depths, strengthening the assumption that no change in CPO with depth is present within the investigated ice volume. Modelled anisotropies generally match observations within uncertainties along the Flow, Perp and $-\mathrm{45}{}^{\circ }$ profiles. It is notable that the three different models exhibit very similar seismic properties along these profiles, although the CPO model informed by $\mathit{\chi }\left(v_{\mathrm{P}}\right)$ in Fig. d shows a fundamentally different geometry than the CPO models shown in Fig. e and f. The $-\mathrm{45}{}^{\circ }$ profile shows a larger difference between the models with a steady decrease in $\mathit{\delta }{v}_{\mathrm{P}}$ predicted by the models informed by $\mathit{\chi }\left(v_{\mathrm{S}\mathrm{1}}\right)$ and $\mathit{\chi }\left(v_{\mathrm{P}}\right)+\mathit{\chi }\left(v_{\mathrm{S}\mathrm{1}}\right)$, while the $\mathit{\chi }\left(v_{\mathrm{P}}\right)$ CPO model suggests a decrease in $\mathit{\delta }{v}_{\mathrm{P}}$ to incidence angles of $\sim \mathrm{40}{}^{\circ }$, followed by increasing $\mathit{\delta }{v}_{\mathrm{P}}$ towards larger incidence angles. The difference between data and models is largest along this profile; however the observed $\mathit{\delta }{v}_{\mathrm{P}}$ values support the steady decrease in P-wave velocity between the vertical ($\mathrm{0}{}^{\circ }$) and horizontal incidence ($\mathrm{90}{}^{\circ }$) that is suggested by the $\mathit{\chi }\left(v_{\mathrm{S}\mathrm{1}}\right)$ and $\mathit{\chi }\left(v_{\mathrm{P}}\right)+\mathit{\chi }\left(v_{\mathrm{S}\mathrm{1}}\right)$ CPO models.

Figure 8

Multi-azimuth P-wave observed velocity anisotropy (symbols) and model results (dashed lines) along the different profiles. Symbol colour shows sensor depths. (a) Flow profile; (b) Perp profile; (c) $\mathrm{45}{}^{\circ }$ profile; (d) $-\mathrm{45}{}^{\circ }$ profile.

[Figure omitted. See PDF]

Fast S-wave anisotropy $\mathit{\delta }{v}_{\mathrm{S}\mathrm{1}}$ observations are shown in Fig.  alongside the modelled anisotropy. Measured velocities exhibit a variation with incidence angle which is consistent across all seismometer depths. Again, observations and models are all in general agreement on the Flow and Perp profile, albeit with an observed decrease in $\mathit{\delta }{v}_{\mathrm{S}\mathrm{1}}$ on the Perp profile for incidence angles $\ge \mathrm{60}{}^{\circ }$ which is not satisfyingly matched by any model. The most significant differences are seen along the $\mathrm{45}{}^{\circ }$ and $-\mathrm{45}{}^{\circ }$ profile: here only the models informed by $\mathit{\chi }\left(v_{\mathrm{S}\mathrm{1}}\right)$ and $\mathit{\chi }\left(v_{\mathrm{P}}\right)+\mathit{\chi }\left(v_{\mathrm{S}\mathrm{1}}\right)$ match observations, and the model informed by $\mathit{\chi }\left(v_{\mathrm{P}}\right)$ is strongly mismatched.

Figure 9

Multi-azimuth S-wave observed velocity anisotropy (symbols) and model results (dashed lines) along the different profiles. Symbol colour shows sensor depths. (a) Flow profile; (b) Perp profile; (c) $\mathrm{45}{}^{\circ }$ profile; (d) $-\mathrm{45}{}^{\circ }$ profile.

[Figure omitted. See PDF]

4.4 CPO model fit using ultrasonic velocities

Misfits $\mathit{\chi }$ of the horizontal-cluster CPO are calculated from the multi-azimuth observations of $v_{\mathrm{P}}$, $v_{\mathrm{S}\mathrm{1}}$ and $v_{\mathrm{S}\mathrm{2}}$ on the core samples $\mathrm{003}$, $\mathrm{007}$ and $\mathrm{010}$.

Figure shows misfits $\mathit{\chi }$ of the different velocities for sample $\mathrm{003}$, where the red mark indicates the best-fitting model. The individual misfits in Fig. a–c and the sum of all misfits in Fig. d show best-fitting models with a consistent cluster orientation and a small scatter in reconstructed cluster opening angles. The shape of the areas of lowest misfit shows that cluster orientation can be constrained with higher confidence than cluster width. An overview of CPO model parameters informed by the ultrasonic measurements is provided in Table . The misfit informed by $v_{\mathrm{S}\mathrm{1}}$ in Fig. b shows a local minimum at an orientation that is $\sim \mathrm{90}{}^{\circ }$ offset from the best-fitting model. This reflects the shape of the measured $v_{\mathrm{S}\mathrm{1}}$ variation with the azimuth, which shows approximately $\mathrm{90}{}^{\circ }$ periodicity, whereas $v_{\mathrm{P}}$ and $v_{\mathrm{S}\mathrm{2}}$ show $\mathrm{180}{}^{\circ }$ periodicity (see Fig. ). The sum of all misfits (Fig. d) is regarded as providing the most reliable constraint on model fit. Here, the local minimum present in $\mathit{\chi }\left(v_{\mathrm{S}\mathrm{1}}\right)$ is attenuated, which resolves the ambiguity in the CPO model constraint.

Figure 10

Misfit surfaces of seismic velocities showing fit to model parameters of the horizontal-cluster CPO for sample $\mathrm{003}$: (a) $\mathit{\chi }\left(v_{\mathrm{P}}\right)$, (b) $\mathit{\chi }\left(v_{\mathrm{S}\mathrm{1}}\right)$, (c) $\mathit{\chi }\left(v_{\mathrm{S}\mathrm{2}}\right)$ and (d) $\mathit{\chi }\left(v_{\mathrm{P}}\right)+\mathit{\chi }\left(v_{\mathrm{S}\mathrm{1}}\right)+\mathit{\chi }\left(v_{\mathrm{S}\mathrm{2}}\right)$.

[Figure omitted. See PDF]

Table 4

Best-fitting CPO model parameters and uncertainties informed by ultrasonic data. The notation $\sum \mathit{\chi }=\mathit{\chi }\left(v_{\mathrm{P}}\right)+\mathit{\chi }\left(v_{\mathrm{S}\mathrm{1}}\right)+\mathit{\chi }\left(v_{\mathrm{S}\mathrm{2}}\right)$ is introduced.

Misfit type (sample)	Cluster orientation $\mathit{\phi }$ (${}^{\circ }$)	Cluster width $\mathit{\alpha }$ (${}^{\circ }$)	Misfit
$\mathit{\chi }\left(v_{\mathrm{P}}\right)$ ($\mathrm{003}$)	100 [100, 105]	13.1 [9.1, 18.1]	0.13 [0.17]
$\mathit{\chi }\left(v_{\mathrm{S}\mathrm{1}}\right)$ ($\mathrm{003}$)	105 [100, 105]	19.1 [16.1, 23.1]	0.18 [0.20]
$\mathit{\chi }\left(v_{\mathrm{S}\mathrm{2}}\right)$ ($\mathrm{003}$)	105 [95, 105]	16.1 [12.1, 18.1]	0.31 [0.33]
$\sum \mathit{\chi }$ ($\mathrm{003}$)	105 [100, 105]	17.1 [13.1, 20.1]	0.22 [0.24]
$\sum \mathit{\chi }$ ($\mathrm{007}$)	125 [120, 130]	14.1 [10.1, 19.1]	0.16 [0.22]
$\sum \mathit{\chi }$ ($\mathrm{010}$)	130 [125, 135]	14.1 [8.1, 18.1]	0.16 [0.22]

The best-fitting models of the three core samples are identified by the minimum in the misfit sum $\mathit{\chi }\left(v_{\mathrm{P}}\right)+\mathit{\chi }\left(v_{\mathrm{S}\mathrm{1}}\right)+\mathit{\chi }\left(v_{\mathrm{S}\mathrm{2}}\right)$. The $c$-axis distribution of these CPO models is shown in upper-hemisphere plots in Fig.  alongside the measured CPO of the samples.

Figure 11

Measured (black) and modelled (red) CPO geometries. $c$-axis orientations in an upper-hemisphere stereographic projection are shown. The VSP orientations of the Flow and Perp profiles are labelled. (a) Sample $\mathrm{003}$; (b) sample $\mathrm{007}$; (c) sample $\mathrm{010}$.

[Figure omitted. See PDF]

5 Discussion

5.1 Horizontal-cluster reconstruction

Both the VSP data and the ultrasonic data are best matched by CPOs with cluster azimuths between $\mathrm{105}$ and $\mathrm{130}{}^{\circ }$ from the macroscopic ice flow direction at the site. This is consistent with CPO measurements made on core samples using EBSD analysis . The agreement in model results between ultrasonic and seismic data bridges the two different scales at which these measurements are made. Ultrasonic and EBSD measurements that are made in ice core samples are shown to be representative of anisotropy within the macroscopic volume that is sampled by the VSP survey.

The inter-sample variation in CPO cluster orientations inside the shallow ice at the site is also clearly visible in the ultrasonic velocity variation: Fig.  shows that the ultrasonic velocities in the three samples follow similar patterns, which are slightly offset between the individual samples. Different fast seismic directions in the studied samples are apparent in Fig. , especially between sample $\mathrm{003}$ and samples 007 and 010. Consequently the CPO model of $\mathrm{003}$ shows a different cluster azimuth than the models of samples $\mathrm{007}$ and $\mathrm{010}$ in Fig. .

Figure  also shows that the $c$-axis maxima modelled from seismic anisotropies for samples $\mathrm{007}$ and $\mathrm{010}$ are rotated around the vertical axis relative to CPO measurements. The measured and modelled CPOs of sample $\mathrm{003}$ show excellent agreement; however for samples $\mathrm{007}$ and $\mathrm{010}$ the modelled CPO clusters are slightly rotated in a clockwise direction relative to the measured CPO geometry. Using the EBSD-measured $c$-axis orientations, a forward model of velocity anisotropy is calculated with MTEX. A comparison of observed and forward-modelled velocity variations is given for sample $\mathrm{003}$ in Fig. , which highlights that the constrained CPO model explains measured and CPO-predicted seismic anisotropy to a high degree.

Figure 12

Observed, modelled and CPO-predicted anisotropy for sample $\mathrm{003}$. (a) $\mathit{\delta }{v}_{\mathrm{P}}$; (b) $\mathit{\delta }{v}_{\mathrm{S}\mathrm{1}}$; (c) $\mathit{\delta }{v}_{\mathrm{S}\mathrm{2}}$.

[Figure omitted. See PDF]

The observed, modelled and forward-modelled velocity anisotropy for sample $\mathrm{007}$ is shown in Fig. . The observed anisotropy and model anisotropy are in excellent agreement, but both curves are clearly offset by $\sim \mathrm{10}{}^{\circ }$ from the predicted velocity variations from measured CPO in this sample.

Figure 13

Observed, modelled and CPO-predicted anisotropy for sample $\mathrm{007}$. (a) $\mathit{\delta }{v}_{\mathrm{P}}$; (b) $\mathit{\delta }{v}_{\mathrm{S}\mathrm{1}}$; (c) $\mathit{\delta }{v}_{\mathrm{S}\mathrm{2}}$.

[Figure omitted. See PDF]

This azimuthal difference between EBSD-measured CPO and CPO models based on ultrasonic data is best explained by small-scale variation within the samples. There is a potential azimuthal error of up to about $\mathrm{3}{}^{\circ }$ in cutting and mounting samples for CPO measurement, and this error will primarily be a rotation around the core axis. Probably more important is the fact that the locations of CPO and ultrasonic measurements for a sample do not coincide. As small-scale rotations of the $c$-axis maximum, around a vertical axis, are observed in the core , a difference in the sample position of a few centimetres along the core axis could give a few degrees rotation of the $c$-axis maxima. Some sample locations, through the depth of the core, have a population of oddly oriented grains, in addition to a main $c$-axis maximum , that would give rise to a rotation around a vertical axis of sample-average acoustic anisotropy of $\sim \mathrm{10}{}^{\circ }$: areas with and without this orientation population can occur within the same core section (Fig. 9 in ).

The seismic VSP data record wavelengths that are larger than the scale of individual samples and therefore do not resolve a variation in cluster orientation at this scale. The CPO model informed by VSP velocities provides an averaging of any given variation in CPO properties over the entire sampled depth range. The CPO model found to best explain the VSP data is characterised by a slightly larger cone opening angle compared to results of the individual ultrasonic measurements and therefore exhibits slightly weaker anisotropy.

5.2 The role of the studied seismic phase (P or S wave) and survey geometry on CPO model ambiguity

The VSP seismic and ultrasonic datasets presented in this study have fundamentally different acquisition geometries which ultimately determine the observed velocity variation due to CPO. The constraint of CPO models from seismic anisotropy is consequently highly sensitive to the sampling geometry.

The acquisition geometry of the VSP survey will have a critical control on the ability to distinguish different CPO patterns. For example the Flow and Perp lines show equally good P-wave (Fig. a, b) and S-wave (Fig. a, b) fits to the $\mathit{\chi }\left(v_{\mathrm{P}}\right)$ model, which has the $c$-axis cluster azimuth at $\mathrm{65}{}^{\circ }$ (Fig. a, Table ), and to the $\mathit{\chi }\left(v_{\mathrm{S}\mathrm{1}}\right)$ and $\sum \mathit{\chi }$ models, which have the $c$-axis cluster azimuths of $\mathrm{110}{}^{\circ }$ to $\mathrm{115}{}^{\circ }$ (Fig. b, c; Table ). The diagonal lines show clear differences in the predictions of the $\mathit{\chi }\left(v_{\mathrm{P}}\right)$ model and the other models in both P-wave (Fig. c, d) and S-wave (Fig. c, d) velocity variations with the incidence angle. Modelling CPOs based on P-wave travel times alone does not give a good fit to the EBSD data on diagonal lines: using S-wave travel times or S-wave travel times combined with P-wave travel times gives a much better fit to the EBSD measurements.

The ultrasonic data offer a dense sampling of velocities along azimuths in the horizontal direction. The given CPO type with horizontal $c$-axis clusters exhibits the largest magnitude of velocity variation in the horizontal plane; therefore this sampling geometry results in high sensitivity to CPO parameters. The individual reconstructed CPO models using $v_{\mathrm{P}}$, $v_{\mathrm{S}\mathrm{1}}$ or $v_{\mathrm{S}\mathrm{2}}$ anisotropy are found to be consistent and also generally agree with measured CPO in core samples (see Figs.  and ). Periodicity of seismic anisotropy can result in ambiguities: $\mathit{\delta }{v}_{\mathrm{P}}$ and $\mathit{\delta }{v}_{\mathrm{S}\mathrm{2}}$ exhibit $\mathrm{180}{}^{\circ }$ periodicity and reconstruct the CPO cluster's measured orientation without ambiguity. The cluster orientation reconstructed by $\mathit{\chi }\left(v_{\mathrm{S}\mathrm{1}}\right)$ is also found to agree with measured orientation; however a misfit local minimum, offset by a $\mathrm{90}{}^{\circ }$ azimuth to the global misfit minimum, reduces confidence in the results. We suggest that it is important to consider the information from velocity measurements of all phases that are available to avoid a local minimum. In this study this is achieved by forming the sum of individual misfits $\mathit{\chi }$ and selecting the model associated with the misfit minimum in this case.

The influence of azimuthal sampling on model results is investigated by sub-sampling of the multi-azimuth ultrasonic dataset. Figure  presents model parameter results for sample $\mathrm{003}$ where, instead of $\mathrm{10}{}^{\circ }$ azimuth spacing between measurements, increments of $\mathrm{20}$, $\mathrm{30}$, $\mathrm{60}$ and $\mathrm{90}{}^{\circ }$ are used to inform CPO models. For each chosen new sampling interval, all possible downsampled datasets are considered. For example, an azimuth spacing of $\mathrm{20}{}^{\circ }$ allows two downsampled datasets to be created: one, where the first measurement is taken at $\mathrm{0}{}^{\circ }$ azimuth and another one, where the first measurement is at $\mathrm{10}{}^{\circ }$.

Figure 14

Parameter results using downsampled velocity measurements with azimuthal spacing of $\mathrm{20}$, $\mathrm{30}$, $\mathrm{60}$ and $\mathrm{90}{}^{\circ }$ (indicated by symbol size) for sample $\mathrm{003}$: (a) $\mathit{\chi }\left(v_{\mathrm{P}}\right)$, (b) $\mathit{\chi }\left(v_{\mathrm{S}\mathrm{1}}\right)$, (c) $\mathit{\chi }\left(v_{\mathrm{S}\mathrm{2}}\right)$ and (d) $\mathit{\chi }\left(v_{\mathrm{P}}\right)+\mathit{\chi }\left(v_{\mathrm{S}\mathrm{1}}\right)+\mathit{\chi }\left(v_{\mathrm{S}\mathrm{2}}\right)$.

[Figure omitted. See PDF]

Models informed by $\mathit{\chi }\left(v_{\mathrm{P}}\right)$ (Fig. a) confirm the result of the full dataset if a coarser sampling of velocities with azimuth increments of $\mathrm{20}$, $\mathrm{30}$ or $\mathrm{60}{}^{\circ }$ is chosen. Models informed by multi-azimuth measurements with $\mathrm{90}{}^{\circ }$ spacing show a large spread of CPO parameter results. Models informed by $\mathit{\chi }\left(v_{\mathrm{S}\mathrm{1}}\right)$ and $\mathit{\chi }\left(v_{\mathrm{S}\mathrm{2}}\right)$ (Fig. b and c) agree for spacing of $\mathrm{20}{}^{\circ }$ and $\mathrm{30}{}^{\circ }$, with larger spacing resulting in a wide scatter of model parameters. At $\mathrm{60}{}^{\circ }$ sampling, $\mathit{\chi }\left(v_{\mathrm{S}\mathrm{1}}\right)$ reconstructs three models at the EBSD-observed cluster orientation and three models at the $\mathrm{90}{}^{\circ }$ offset orientation associated with the misfit local minimum in Fig. b. The approach of using the sum of all misfits (Fig. d) reduces the spread of the model parameters found and finds realistic models for $\mathrm{20}$ and $\mathrm{30}{}^{\circ }$ spacing. For $\mathrm{60}{}^{\circ }$ spacing, the cluster orientation angle is correctly found but the cluster width is poorly constrained. This is a clear improvement relative to the reconstructed model parameters in Fig. b and c for $\mathrm{60}{}^{\circ }$ spacing. The results for $\mathrm{90}{}^{\circ }$ spacing are widely scattered and not in agreement with the model parameters that were found using the full dataset. The discrepancy in model parameter results found by the different $\mathrm{90}{}^{\circ }$ sampling geometry relates to the fact that there is more than one cluster orientation that will give the same ratios of orthogonal velocities. Thus the solution is very dependent on the measurement errors in the velocity dataset. At $\mathrm{90}{}^{\circ }$ sampling the errors in $v_{\mathrm{S}\mathrm{1}}$ will be particularly influential as $v_{\mathrm{S}\mathrm{1}}$ velocities should be very similar for orthogonal azimuths (see predictions in Figs.  and ).

Our downsampling analysis shows that in an ideal survey geometry, which is in this case given by the sampling of a horizontal-cluster CPO by horizontal velocity measurements, a realistic CPO model can be created from sampling in up to $\mathrm{30}{}^{\circ }$ steps regardless of the seismic phase studied. Observed scatter in model parameter results shows that the set of actual sampled azimuths becomes critical for larger sampling intervals: the datasets with $\mathrm{20}$ and $\mathrm{30}{}^{\circ }$ sampling are much less sensitive to the choice of the first sampled azimuth compared to the datasets with $\mathrm{60}$ and $\mathrm{90}{}^{\circ }$ spacing. The consideration of the full range of seismic velocity information ($\mathit{\chi }\left(v_{\mathrm{P}}\right)+\mathit{\chi }\left(v_{\mathrm{S}\mathrm{1}}\right)+\mathit{\chi }\left(v_{\mathrm{S}\mathrm{2}}\right)$) can reduce scatter and aid in reconstructing a realistic CPO model from data with wider azimuthal spacing.

The VSP survey CPO modelling presents ambiguity in cluster orientation if only P-wave velocities are considered, as shown in Fig. a. The best-fitting model in this case does not match the CPO observations from the site, orienting the cluster at an azimuth that is offset by $\sim \mathrm{45}{}^{\circ }$ in the anticlockwise direction from observed cluster orientations. In Fig.  the uncertainties in $\mathit{\delta }{v}_{\mathrm{P}}$ observations are large compared to the P-wave anisotropy of the best-fitting model.

The difficulty in reconstructing a realistic CPO model from P-wave velocities in the VSP dataset is a consequence of poor azimuthal sampling. The VSP survey is characterised by coverage of ray paths from a range of incidence angles from multiple azimuths. For the given horizontal-cluster CPO this geometry is clearly not as sensitive to anisotropy as the dense azimuthal sampling of horizontal velocities, which characterises the ultrasonic measurements. The inclusion of P- and S-wave phase information mitigates this shortcoming of the VSP data by resolving the ambiguity in cluster orientation and identifying a CPO model which is in agreement with measured CPO. The study of all available seismic phases should therefore become standard in seismic CPO constraints in ice, rather than the commonly encountered focus on P-wave velocities.

For the horizontal-cluster CPO, the variation in seismic velocities with incidence angle is highly dependent on the azimuth. Therefore, VSP data might be unable to constrain this CPO if azimuths with strong variation are not sampled. The difficulty in finding a correct CPO model from VSP P-wave velocities could be a consequence of this problem, highlighted in Fig.  by large error bars relative to the overall $v_{\mathrm{P}}$ variation with incidence angle. A greater azimuthal sampling of VSP data is required to improve the CPO model constraint. Furthermore, larger shot offsets should be targeted to sample near-horizontal ray paths for the study of horizontal-cluster CPOs, since the largest amplitude of the anisotropy signal is in the horizontal plane.

6 Conclusions

We have conducted a vertical-seismic-profile (VSP) experiment and laboratory ultrasonic experiments aimed at measuring the seismic anisotropy of ice from the lateral shear margin of the Priestley Glacier, Antarctica, and at linking these data to seismic anisotropy model predictions based on measured crystallographic preferred orientations (CPOs) in EBSD data.

Horizontal-cluster CPO models are informed by P-wave and S-wave velocity anisotropy data from a $\sim \mathrm{50}$ m scale, four-azimuth, walkaway-VSP experiment. A best-fitting CPO model informed by VSP P-wave data alone is found to exhibit a $c$-axis cluster azimuth, which is inconsistent with measured $c$-axis orientations and also overestimates the strength of $c$-axis clustering. Best-fitting CPO models based on VSP S-wave data or a combination of P-wave and S-wave data are found to give a good match to the direct CPO measurements.

Azimuthal ultrasonic P-wave and S-wave velocity measurements, made in $\mathrm{10}{}^{\circ }$ increments, on core samples from the borehole show the pattern of anisotropy with considerable detail. The anisotropy pattern matches the pattern predicted from the CPO in the shallowest sample $\mathrm{003}$. Differences in CPO-predicted and measured anisotropy for deeper samples $\mathrm{007}$ and $\mathrm{010}$ are an indication of intra-sample CPO variability. Inter-sample variability is also apparent from ultrasonic velocity measurements: the anisotropy patterns in different samples are rotated through small angles relative to each other around the core axis, a pattern that matches direct CPO measurements.

The ultrasonic data have been downsampled to larger azimuthal increments ($\mathrm{20}$, $\mathrm{30}$, $\mathrm{60}$ and $\mathrm{90}{}^{\circ }$) to explore how well these lower-resolution data constrain the CPO responsible for velocity anisotropy. Increments of up to $\mathrm{30}{}^{\circ }$ constrain both the azimuthal orientation and the intensity of horizontal $c$-axis alignment well. Increments of $\mathrm{60}{}^{\circ }$ constrain orientation. Increments of $\mathrm{90}{}^{\circ }$ do not provide useful constraints. The design of field seismic surveys should be designed to address ambiguity as result of sampling geometry by considering two main points. First, a dense sampling of propagation angles must be realised by the acquisition geometry. Second, sources and receivers must enable the derivation of both P-wave and S-wave anisotropy.

Data availability

VSP seismic data, ultrasonic data and travel time picks are available at 10.17608/k6.auckland.17108639 .

Author contributions

DJP and CLH led the Priestley Glacier project. DJP and FL designed the seismic and ultrasonic experiments. FL collected the multi-azimuth ultrasonic data. HS collected the axial ultrasonic data. SF wrote MTEX code used to simulate CPOs. All authors except CLH and SF were involved in seismic data collection in the field. FL wrote the manuscript in collaboration with DJP. All authors edited the manuscript.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acknowledgements

We would like to acknowledge the logistics support from Antarctica New Zealand and field support from staff at Scott Base, Mario Zucchelli and Jang Bogo research stations. We would also like to thank Brent Pooley for building the ultrasonic rig. We would like to acknowledge Adam Booth for the editorial work and the reviewer Chao Qi and the anonymous reviewer for their insightful comments.

Financial support

This research has been supported by the Marsden Fund (grant no. UOO1716), the Korea Polar Research Institute (grant no. PE22430) and the New Zealand Antarctic Research Institute (Early Career Researcher Seed Grant, grant no. NZARI 2020-1-3).

Review statement

This paper was edited by Adam Booth and reviewed by Chao Qi and one anonymous referee.