Deriving micro to macroscale seismic velocities from ice core c-axis orientations

One of the great challenges in glaciology is the ability to estimate the bulk ice anisotropy in ice sheets and glaciers, which is needed to improve our understanding of ice-sheet dynamics. We investigate the effect of crystal anisotropy on seismic velocities in glacier ice and revisit the framework which is based on fabric eigenvalues to derive approximate seismic velocities by exploiting the assumed symmetry. In contrast to previous studies, we calculate the seismic velocities using the exact c-axis angles describing the orientations of the crystal ensemble in an ice-core sample. We apply this approach to fabric data sets 5 from an alpine (KCC) and a polar (EDML) ice core. Our results provide a quantitative evaluation of the earlier approximative eigenvalue framework. For near-vertical incidence our results differ by up to 135 m s−1 for P-wave and 200 m s−1 for S-wave velocity compared to the earlier framework (estimated 1 % difference in average P-wave velocity at bedrock for the short alpine ice core). We quantify the influence of shear-wave splitting at bedrock as 45 m s−1 for the alpine ice core and 59 m s−1 for the polar ice core. At non-vertical incidence we obtain differences of up to 185 m s−1 for P-wave and 280 m s−1 for S-wave 10 velocities. Additionally, our findings highlight the variation in seismic velocity at non-vertical incidence as a function of the horizontal azimuth of the seismic plane, which can be significant in case of non-symmetric orientation distributions and results in a strong azimuth-dependent shear-wave splitting of max. 281 m s−1 at some depths. For a given incidence angle and depth we estimated changes in phase velocity of almost 200 m s−1 for P-wave and more than 200 m s−1 for S-wave and shear-wave splitting under a rotating seismic plane. We assess for the first time the change in seismic anisotropy that can be expected 15 on a short spatial (vertical) scale in a glacier due to a strong variability in crystal-orientation fabric (±50 m s−1 per 10 cm). Our investigation of seismic anisotropy based on ice-core data contributes to advancing the interpretation of seismic data, with respect to extracting bulk information about crystal anisotropy without having to drill an ice core and with special regard to future applications employing ultrasonic sounding.

velocities. Additionally, our findings highlight the variation in seismic velocity at non-vertical incidence as a function of the horizontal azimuth of the seismic plane, which can be significant in case of non-symmetric orientation distributions and results in a strong azimuth-dependent shear-wave splitting of max. 281 m s −1 at some depths. For a given incidence angle and depth we estimated changes in phase velocity of almost 200 m s −1 for P-wave and more than 200 m s −1 for S-wave and shear-wave splitting under a rotating seismic plane. We assess for the first time the change in seismic anisotropy that can be expected 15 on a short spatial (vertical) scale in a glacier due to a strong variability in crystal-orientation fabric (± 50 m s −1 per 10 cm).
Our investigation of seismic anisotropy based on ice-core data contributes to advancing the interpretation of seismic data, with respect to extracting bulk information about crystal anisotropy without having to drill an ice core and with special regard to future applications employing ultrasonic sounding.
One of the most important goals for glaciological research is the establishment of a thorough understanding of ice dynamics, in which internal deformation plays a crucial role. This deformation is predominantly evident and described on a macro-scale (O(km)) as most observations rely on remote sensing or ice sheet modeling. However, it is necessary to connect the bulk behaviour with the governing processes on the micro-scale (O(µm)) to be able to develop a comprehensive understanding of the 5 deformation mechanisms (Weikusat et al., 2017). The fundamental deformation mechanisms on the atomic scale are driven by the external stress field and lead to the alignment of single ice crystals in preferential directions (previous works are reviewed in Faria et al., 2014). Due to the intrinsic anisotropy of each ice crystal an anisotropic bulk medium is formed as a result of the crystal-preferred orientation (CPO, also known as lattice-preferred orientation, LPO, and crystal-orientation fabric, COF). The anisotropy is evident in elastic, plastic and electromagnetic properties of the ice and the respective parameters can be connected 10 to each other. The plastic anisotropy can have a considerable influence on the bulk deformation rate and vice versa. Hence, it is desirable to incorporate the development of anisotropy in flow models (Pettit et al., 2007;Seddik et al., 2008;Martin et al., 2009;Graham et al., 2018).
Currently, the development and extent of fabric anisotropy is mainly investigated by laboratory measurements on ice core samples which provide one-dimensional data (along the core axis, z-axis in Fig. 1) that only partially cover the length of the core. 15 However, geophysical evidence of crystal-orientation fabric can also be obtained by exploiting the elastic anisotropy which influences the propagation of seismic waves in ice (Blankenship and Bentley, 1987;Smith et al., 2017). Seismic waves propagate from a seismic source to the seismic receivers on the glacier surface within the seismic plane (Fig. 1). This is the vertical plane underneath the horizontal seismic profile, which runs along the surface of the glacier, and may also contain the vertical ice core axis, along which fabric information is collected. Seismic reflections occur due to sudden changes of fabric (Horgan 20 et al., 2008(Horgan 20 et al., , 2011Hofstede et al., 2013) and offer the chance to obtain spatially distributed information on the COF structure in various depths of the ice column, the acquisition of which would be unfeasible using direct sampling via ice-coring.
The motivation of this study is to improve the interpretation of seismic data by connecting the micro-and the macro-scale using the elastic properties of ice. Early work to this end were accomplished by Bennett (1968); Bentley (1972) ;Blankenship and Bentley (1987) and more recent approaches include Gusmeroli et al. (2012), Maurel et al. (2015), and Diez and Eisen 25 (2015). Specifically, the starting point of this paper is the study by , who establish a connection between the commonly reported fabric parameter of second-order orientation tensor eigenvalues and the elasticity tensor describing the polycrystalline medium to calculate seismic velocities from ice-core fabric data. (We refer to the method of  as the ev framework.) They illustrate the proposed calculation framework on fabric data from the polar ice core EDML (European Project for Ice Coring in Antarctica in Dronning Maud Land). 30 Our main objective is to present an improved method for the estimation of the bulk elasticity tensor, and to use this to (i) evaluate the use of the ev framework, and (ii) demonstrate the effect of a real fabric on seismic velocities. Our study concentrates on the forward calculation of seismic velocities from the c-axis orientation distribution. The application of an inverse method will likely always require some simplification to symmetries, for which we now can quantify involved uncertainties -a required component of the covariance matrices for inverse methods.
We first present experimental measurements, theoretical basis and mathematical algorithm of our new framework (cx). We apply this framework to two ice cores and further investigate how fabric variability at the submetre scale is reflected in theoretical seismic interval velocities. Finally, we assess the effect of asymmetrical fabric distributions on seismic velocity and explore the potential of azimuth-dependent seismic surveys. For our analysis of seismic velocities we use fabric data from the polar ice core EDML and the alpine ice core KCC. The EDML ice core was drilled as part of EPICA (European Project for Ice Coring in Antarctica) between 2001 and 2006 at Kohnen Station, Antarctica, and reaches to a depth of 2774 m (Oerter et al., 2009;Weikusat et al., 2017). The KCC ice core  (Bohleber et al., 2018). KCC is 72 m long with the firn-ice-transition at a depth of about 36 m and a borehole temperature between −13.6°C in 13 m depth and −12.4°C at the bedrock, measured in 2014 (pers. comm. M. Hoelzle, University of Fribourg, 2014). Both KCC and EDML were stored at a minimum of −18°C during transport and at −30°C during processing. 15 Vertical and horizontal thin sections of the ice cores were prepared and measured using polarised light microscopy (e.g. Wilson et al., 2003;Peternell et al., 2009) with an automatic fabric analyser from Russell-Head Instruments (models G20 and G50 in case of EDML and G50 for KCC). For each identified ice crystal in the thin section the measurement provides the orientation of the crystallographic c-axis c by two spherical coordinates, azimuth ϑ in the interval (0, 2π) and colatitude ϕ in the interval (0, π/2), with respect to the vertical ice-core axis that we define to coincide with the z-axis of the global coordinate system 20 ( Fig. 1). The c-axis is expressed in the global cartesian coordinate system {x, y, z} as a vector with unit length: c(ϑ, ϕ) = (sin(ϕ) cos(ϑ), sin(ϕ) sin(ϑ), cos(ϕ)) (1) Figure 1 illustrates the geometric relation between the angles ϑ, ϕ for describing the c-axis orientations from an ice core sample and the setup of a seismic survey profile across an ice core. The incidence of a seismic wave on an ice sample is determined by the angle of incidence ψ and the azimuth angle ϑ s of the seismic plane. 25 The EDML fabric data (data sets: Weikusat et al., 2013a, b, c, d) are presented in detail in Weikusat et al. (2017). The total data set used in this study comprises 154 samples between 104 and 2563 m depth with a coarse sampling interval and 40 additional vertical section samples that were measured continuously in several intervals between 2359 and 2380 m. These high resolution measurements were performed with the G50 instrument. The KCC fabric data (data set: Kerch et al., 2016a, b) consist of 85 vertical thin sections covering 11 % of the entire ice core. 30 Eigenvalues λ i (i = 1, 2, 3) of the second-order orientation tensor a ij are usually calculated from the c-axis distribution within a thin section sample and can be grain-, area-or volume-weighted to describe the fabric (Woodcock, 1977;Durand et al., 2006;Figure 1. The global coordinate system {x, y, z}, typically corresponding to East, North and up, used for the description of a c-axis with unit length: spherical coordinates ϑ and ϕ specify the orientation of the c-axis. For each grain the c-axis can be expressed in its local cartesian coordinate system {p, q, r} by (0, 0, 1); p,q are not shown here. The equatorial plane corresponds to the horizontal thin section plane. The z-axis is assumed to be parallel to the ice core axis. The orientation of a hypothetical seismic plane (grey) is defined by the seismic azimuth angle ϑs, with the angle ψ of an incident seismic wave (dashed line) with unit vector normal to the plane wavefront n. The star symbolises the seismic source and the triangles represent a line of seismic receivers. Mainprice et al., 2011). They describe the type and strength of anisotropy in the crystal ensemble visible in the thin section (e.g. Cuffey and Paterson, 2010). Typically, different types of fabric are identified  by the relation of the three eigenvalues with λ 1 ≤ λ 2 ≤ λ 3 and λ i = 1. By this classification the crystal anisotropy of the bulk can be described in a convenient way, if a unimodal distribution can be assumed, and can be associated with different deformation regimes (e.g. Weikusat et al., 2017).

Seismic wave propagation in anisotropic ice
In a glacier, the fabric anisotropy also introduces an anisotropy of the elastic properties of the material. This elastic anisotropy results in a seismic anisotropy, which means the propagation of seismic waves is influenced by the fabric anisotropy. To study this connection, theoretical velocities can be calculated if the fabric anisotropy is known.
The mathematical background for the calculation of seismic phase velocities from the elastic properties in anisotropic ice can 10 be found in many publications (e.g. Tsvankin, 2001;. For convenience the essential concepts are shortly repeated in the following. Here we focus on phase velocities and group velocities are not considered.
For an anisotropic elastic medium -ice behaves elastically for the propagation of seismic waves -stress and strain are linearly connected following the generalised Hooke's law: where c mnop is the elasticity tensor, a fourth-order tensor which describes the medium's elastic properties. The inverse relation uses the compliance tensor s mnop . Due to the symmetry of strain and stress tensor and thermodynamic considerations (Aki and Richards, 2002), the 81 unknowns of the elasticity tensor reduce to 21 independent components for general anisotropy.
The elasticity tensor can then be expressed in a simplified manner, known as Voigt notation (Voigt, 1910), where pairs of indices from the fourth-order tensor are replaced by single indices. The resulting elasticity tensor in Voigt notation C ij (i, j = 5 1, 2, 3, 4, 5, 6) is a symmetric second-order tensor. In the case of monocrystalline ice I h , the components of the elasticity tensor were measured in the laboratory by means of Brillouin spectroscopy (Gammon et al., 1983). There are five independent components due to the hexagonal crystal symmetry. Several sets of values for the elastic moduli have been found by different authors, as summarised in . Here, the monocrystal elasticity tensor C m by Gammon et al. (1983), as measured on samples of artificial ice at −16°C, is used for all calculations: To apply this description to the study of large ice sheets and glaciers, we have to consider the elastic properties of the polycrystal. The understanding of the elastic behaviour of a monocrystal can be used together with the fabric description to estimate the elastic properties of the polycrystal. Different theoretical models have been developed for the estimation of the elasticity tensor of an anisotropic polycrystal, usually making use of fabric symmetries (e.g. Nanthikesan and Sunder, 1994;Maurel 15 et al., 2015) or by calculating orientation density functions (ODF, Mainprice et al., 2011). In this context some authors refer to the polycrystal as the "effective medium" (Maurel et al., 2015).
For the calculation of the polycrystal elastic properties from anisotropic monocrystal properties the concept of Voigt-Reussbounds is often used (Nanthikesan and Sunder, 1994;Bons and Cox, 1994;Helgerud et al., 2009;Vaughan et al., 2016), and we also use this approach here. Voigt-Reuss bounds provide estimates of the upper and lower limits for the elastic moduli 20 of the polycrystal, as was shown by Hill (1952), with the Reuss bound exceeding the Voigt bound. Nanthikesan and Sunder (1994) find that the difference of the Voigt-Reuss-bounds for the elastic moduli of polycrystalline ice does not exceed 4.2 % and conclude that either of the bounds or an average is a good approximation.
Once the elastic properties for the polycrystal are known, the Christoffel equation provides the relationship to calculate seismic velocities. For a linearly elastic, arbitrarily anisotropic homogeneous medium the wave equation is solved by a harmonic 25 steady-state plane wave and we obtain the Christoffel equation: with the density ρ, the polarisation vector U , the phase velocity v ph , the unit vector normal to the plane wavefront n, and the Kronecker delta δ mo . G mo = c mnop n n n p is the positive definite, thus symmetric Christoffel matrix. The vector n indicates the 5 direction of wave propagation and depends on the angle of incidence ψ, which is measured from the vertical, and, if applicable, the azimuth angle ϑ s between the vertical plane of incidence and the azimuthal orientation of the ice core with respect to the geocoordinates ( Fig. 1): Eq.
(3) constitutes an eigenvalue problem for G mo . The real and positive eigenvalues are identified with the phase velocities v p , v sh , v sv for P-wave, SH-wave and SV-wave respectively. Different solutions are proposed, depending on the form of the elasticity tensor. The solution used in this study for an arbitrarily anisotropic medium is outlined in section 2.4.
Instead of interval velocities often the root mean square (RMS) velocity v rms is considered, which gives the velocity of the homogeneous half-space equivalent to the stack of N horizontal layers (i) to this depth: 15 with the two-way travel time (TWT) t 0 of a seismic wave that travels vertically through a single layer with the corresponding interval velocity v. For a layered anisotropic medium a reliable depth-conversion from travel times is only feasible if the RMS velocity for zero-offset can be deduced (Diez et al., 2014).
In-situ temperature and density are essential when comparing seismic velocities. However, as this study is focused on the comparison of calculation frameworks that use the same elastic moduli, a temperature correction is not applied.

Recap: Eigenvalue framework
Diez and Eisen (2015) presented a framework for calculating seismic velocities from COF data in form of eigenvalues, which we briefly recapture here. In the following, this framework is referred to as ev framework and associated variables are indicated with ev . 25 The ev framework can be summarised in three steps:

From eigenvalues to seismic interval velocities
1. The fabric data in the standard parameterisation of second-order orientation tensor eigenvalues are sorted into three fabric classes (cone, thick girdle, partial girdle), where each is defined by one or two opening angles χ, φ, symmetrical with respect to the vertical, and enveloping the c-axis distribution of a sample.
2. The opening angles characterising the fabric of each sample are used to integrate the elasticity tensor of the monocrystals, Eq. (2), to obtain the elasticity tensor of the polycrystal, which exhibits an orthorhombic symmetry with respect to the vertical.
3. From the polycrystal elasticity tensor the approximative solutions to the Christoffel equation (3) for the orthorhombic case provided by Daley and Krebes (2004) are applied to obtain seismic interval phase velocities v ev p , v ev sh , v ev sv , which 5 can be used for comparison with measured seismic data. Voigt calculation is used following the argument that Reuss and Voigt bounds are close enough.

Uncertainty of ev framework
The advantages of this approach are : • Eigenvalues are a standard parameter for expressing the strength of fabric and can be directly used for the ev framework 10 without additional information about the particular measurement of thin sections from an ice core.
• By assuming an orthorhombic symmetry the solution to the Christoffel equation can be readily found. No information on the azimuthal orientation of the ice core (relative to any seismic measurements on a glacier) is needed, although this could be considered to improve the results in case of girdle fabric.
However, some uncertainty is inherent in the framework: 15 • The eigenvalues of the second-order orientation tensor do not constitute a complete and unambiguous description of the fabric. Specifically, they do not provide information on preferential orientations with regard to the coordinate system.
To get a rough idea about the orientation of the fabric the eigenvectors would have to be used in addition, an approach seldomly followed. Instead, the orientation of the eigenvector to the largest eigenvalue is typically assumed to correspond to the vertical, which may in fact not be the case and could introduce an unknown uncertainty.

20
• By assuming an orthorhombic symmetry while using opening angles to describe the c-axis distribution any information on asymmetric fabric (with respect to the vertical) is dismissed and approximation errors are introduced for more asymmetric c-axes distributions.
• Fabric data from ice cores indicate that transitions between fabric classes usually occur gradually, and sudden changes are only expected to occur due to changes in impurity content or deformation regime (Montagnat et al., 2014;Weikusat 25 et al., 2017). However, the classification into fabric groups based on threshold values for the eigenvalues can introduce artificial discontinuities in the calculated velocity profile.
We will provide a quantitative estimate of the uncertainty of the ev framework in the following sections.

C-axes framework
In this study we aim to provide a quantitative estimate of the error introduced by the approximation of the ev framework and to assess the potential of the hitherto neglected information for future analyses. For that purpose the exact angle information of each individual c-axis is used in the following to derive the elasticity tensor C p of the polycrystal. Then, the phase velocities in an arbitrarily anisotropic medium are calculated. We refer to the new approach as cx framework. Both frameworks are 5 summarised in Figure 2.

Calculating the elasticity tensor for discrete crystal ensemble
If not indicated otherwise, elasticity/compliance tensors and velocities are calculated for the effective medium, which, in this study, is typically represented by a thin section comprising a number of grains (N g ) of the order of a hundred to a thousand.
Specifically, for EDML the number of grains is mostly between 100 and 650 grains per sample with the exception of some large-grained samples from the depth interval 2300 − 2380 m with less than 100 grains. For KCC the number of grains is between 155 and 1707 grains per sample and there are more than 250 grains per sample in the lowest 5 m of the ice core.
A data set of COF measurements from an ice core is considered that gives pairs of angles determining the c-axis of each grain c(ϑ, ϕ) in a grain ensemble per thin section. We apply the following steps to obtain the effective elasticity tensor for a thin section sample: 5 1. Transformation of the monocrystal elasticity tensor: Considering the monocrystal elasticity tensor C m,k , given by Eq.
(2), in the k-th grain's local coordinate frame {p, q, r} with c = (0, 0, 1), a transformation (indicated by t ) to the global coordinate frame {x, y, z} by using the angles ϑ, ϕ is necessary: with rotation matrix R C as given by Eq. (A3) and R C its transpose matrix. C t m,k is unlikely to have vertical transversely 10 isotropic (VTI) symmetry, as most c-axes in a real fabric do not coincide with the z-axis, but will lie obliquely in the {x, y, z} coordinate frame.
2. Grain area weighting: If grain size information is available, each transformed monocrystal elasticity tensor C t m,k is multiplied by the grain cross-section area (A k ) fraction f k = A k / k A k . Otherwise, it is multiplied by 1/N g for an equal contribution of all grains to the effective medium elasticity tensor. 15 3. Discrete summation over the transformed monocrystal elasticity tensor for all grains to obtain the polycrystal elasticity tensor C p in the global coordinate frame: The obtained elasticity tensor C p is very likely to have only non-zero components and describes an arbitrarily anisotropic medium. 20 Derivation via the compliance tensor: To determine the Voigt-Reuss bounds as introduced above, the polycrystal elasticity tensor is also calculated via the compliance tensor S m . To accomplish this, the monocrystal elasticity tensor is inverted: Steps 1 to 3 are then applied accordingly using Eq. (A4) to derive the compliance tensor of the polycrystal S p , which is then again inverted to provide C R p , where R denotes Reuss.
2.4.2 Deriving seismic interval phase velocities for an arbitrarily anisotropic medium 25 The phase velocities v ph (ψ, ϑ s ) are obtained from the polycrystal elasticity tensor C p by applying the analytical solution to find the eigenvalues v ph of the Christoffel matrix for an arbitrarily anisotropic medium, following Tsvankin (2001). The algorithm is presented in Appendix B and variables associated with the cx framework are annotated by superscript cx . Thus, the velocities are calculated for any fabric, incidence angle ψ, and azimuthal orientation ϑ s of the seismic plane. . a) P-wave velocity v cx p for cone fabric from randomly generated c-axes. b) Difference in P-wave velocity between the two frameworks for cone fabric. Blue color shows where the ev framework obtains higher velocities than the cx framework. Red shades indicate the opposite. They differ by −50 to 20 m s −1 (corresponding to ≤ 1.5 %).

Framework comparison for cone fabrics
The frameworks (ev and cx) were compared by applying them to cone fabric for all cone angles 0 • ≤ φ ≤ 90 • , thus excluding any effects from asymmetric fabric. We generated artificial fabric with 1000 c-axes, randomly distributed in a solid (cone) angle, in 1°steps, and calculated the respective eigenvalues. Figure 3a shows the theoretical P-wave velocity distribution v cx p (ψ, φ) for all cone and incidence angles as calculated with the cx framework and Figure 3b gives the difference between the 5 calculation results from both frameworks. As is to be expected the maximum velocity is found for a seismic wave at vertical incidence on a narrow single maximum fabric. The strongest velocity deviation between the frameworks is found for cone angles of approximately 50 − 60°at vertical incidence, where the ev velocity exceeds the cx velocity by approx. 50 m s −1 (deviation of 1.5 %).
3 From ice core fabric to seismic velocities -case studies 10 We apply the cx framework, outlined in section 2.4, to two fabric data sets from ice cores EDML and KCC, respectively. Thus, we investigate the potential of the new framework with respect to the earlier established ev framework, which was already where high resolution measurements where taken (Fig. 5). c) relates the difference ∆vp0 = v cx p0 − v ev p0 to the fabric classes that are indicated by shading (dark gray: cone, light gray: thick girdle, white: partial girdle). The shading extends for each data point to half the distance to the neighbouring data points. d) shows the seismic RMS velocities resulting from the interval velocities integrated from the surface (without taking density, temperature and firn into account); S-wave velocities refer to the upper x-axis and P-wave velocities to the lower x-axis. applied to fabric data from EDML . We use the same EDML data set (c-axis angles and grain-weighted eigenvalues), complemented by 40 additional thin sections measured more recently. The threshold eigenvalues for classifying the EDML fabric within the ev framework are as follows: girdle fabric is given if λ 2 ≥ 0.2 and λ 1 ≤ 0.1, with thick girdle fabric for 0.05 < λ 1 ≤ 0.1 and partial girdle for λ 1 ≤ 0.05; cone fabric is identified otherwise. The threshold values for classifying the KCC fabric are chosen such that only cone fabric is recognised by the algorithm, i.e. the threshold for girdle fabric is set 5 to λ 2 ≥ 0.4 and λ 1 ≤ 0.1; cone fabric is identified otherwise. This is justified by the evaluation of stereographic projections of the c-axis distributions which shows that cone fabric is dominant in all KCC samples, although some tendencies towards girdle is observed in deeper samples, and artificial discontinuities are prevented. KCC eigenvalues are area-weighted as grain size information is available from automatic image processing. The results obtained from the cx framework are considered to be more accurate for the purpose of comparing the frameworks in the following case studies. If not stated otherwise, all velocities are interval velocites, i.e. the seismic wave velocity within a layer, for which an elasticity tensor is calculated based on the fabric in this layer. The velocity differences between the frameworks for the two case studies are summarised in Table 2.

Seismic interval velocity for vertical incidence
We now assess the velocity difference between the eigenvalue and the cx framework v p0 at vertical incidence of a seismic wave (Figure 1, ψ = 0°) as indicated by subscript 0 , with a focus on the effect of the ev framework fabric classification.

5
Vertical incidence refers to the direction parallel to the ice core axis, which we assume to be normal to the glacier surface. As the seismic P-wave velocity for vertical incidence is invariant under azimuthal rotation of the seismic plane of the core, it is possible to assess the uncertainty introduced by using the eigenvalues. We mainly show results for the P-wave velocity, but included the S-wave velocity in our assessment of RMS velocities.
Vertical incidence at EDML 10 The evolution of the fabric of the EDML ice core becomes apparent when assessing the eigenvalues (Fig. 4a) and is discussed in detail in Weikusat et al. (2017). In the following, observations are made for the comparison of the velocitites from the EDML ice core.
The general trend of the velocities of the two frameworks is in good agreement (Fig. 4b). However, a systematic difference can be observed (Fig. 4c): for cone fabric the P-wave velocity is overestimated by the ev framework, for girdle fabric the P-wave velocity is underestimated.
In the upper 1785 m the velocity from the cx framework clearly exhibits a higher variability, as quantified by the standard 5 deviation s(v p0 ) ( Table 1). Below that depth, there is less variation in the velocity of the cx framework. The higher variability in the eigenvalue velocity below 1785 m is due to several transitions between fabric classes in the depth interval of 1800 to 2035 m; notably the prominent peak at 1802 m is an example of this. In the lower part of the core at 2306 m a sudden weakening of the fabric anisotropy is evident in the results of both frameworks. The velocity is, however, underestimated by the ev framework by 48 m s −1 due to a switch from cone to girdle fabric classification. RMS velocities were calculated from the 10 interval velocities for P-and S-wave (Fig. 4d) using Eq. (5) in order to assess the cumulated effect of the velocity deviation. In the anisotropic case the zero-offset RMS velocities are needed for the depth conversion in classical reflection seismic profiles (Diez et al., 2014). For EDML, the P-wave RMS velocities v p0,rms for the two methods converge towards the bedrock as a serendipitous result of the systematic under-and overestimation described before. The S-wave RMS velocities v s0,rms show a similar shear-wave splitting of 67 m s −1 and 59 m s −1 at the bedrock but the cx velocities also show a small split in the upper 15 750 m of the ice core, where the ev framework assumes a VTI fabric with no resulting shear-wave splitting. Figure 5 is a close-up of the shaded depth in Fig. 4b and shows high resolution measurements completed since  providing ten data points per metre (filled diamonds). The new data exhibit a strong submetre-scale variability in fabric strength (Weikusat et al., 2017) which has not been regularly observed in ice-core fabric data so far (Fig. 5a). Only in recent ice core projects have fabric measurements covered continuous intervals (ongoing measurements at the site of the East Greenland Vertical incidence at KCC 25 We show the results of the velocity calculations for vertical incidence from the KCC fabric data in Fig. 6. The cone angle calculated from the eigenvalues varies between 10-30°for each depth interval (Fig. 6a). A detailed discussion of the fabric data is beyond the scope of this paper. The P-wave interval velocities calculated with both frameworks (Fig. 6b) increase with depth as a stronger anisotropic single-maximum fabric evolves and show high variability between adjacent 10 cm long samples. The difference in P-wave velocity between the two frameworks is shown in Fig. 6c. The ev framework overestimates 30 the P-wave velocity on average by 46 m s −1 . Hence, differences between the frameworks are similar for the KCC ice core as for the cone fabric regions of the EDML ice core. In the bottom layer the largest difference in P-wave velocity is −135 m s −1 , due to the strong single maximum inclined to the vertical. The change in c-axes velocity δv cx p0 of each 10 cm sample to the previous within a continuous measurement interval can exceed 40 m s −1 (Fig. 6d). For the estimate of P-wave and S-wave

Seismic interval velocities for non-vertical incidence
During typical seismic profile surveys the seismic wave will have an inclined angle of incidence with respect to the vertical, normal to the glacier surface. The velocities will vary depending on the incidence angle if the medium is anisotropic and this will affect the recorded travel times . For a single maximum (or cone) fabric that is symmetric around 10 the vertical this angle dependency is invariant under rotation of the seismic plane. However, the symmetry axis of a fabric described by a set of eigenvalues could also be inclined with respect to the vertical depending on the deformation regime on site. The recorded travel times will then depend on the direction of the seismic profile on the glacier surface. This, in turn, can be exploited to acquire additional information on the anisotropy. As the cx framework does not restrict the description of the crystal anisotropy of the effective medium to a symmetry with respect to the vertical, the variation of seismic velocities under 15 a rotating seismic plane can be studied.
In the following we assess how the seismic velocities will change when the ice core fabric data and the seismic plane of incidence are rotated with respect to each other. The zero orientation (ϑ s = 0) is not associated with any specific geographical direction. The largest uncertainty for this assessment originates in the difficulty to identify the ice core's orientation during drilling. Although it is usually tried to fit the consecutive ice core pieces during logging to maintain the correct relative orien-20 tation it is not guaranteed that there are no discontinuities. We use the term difference to refer to the comparison of differently Figure 6. Comparison of zero-offset velocities calculated from KCC fabric data via ev and cx framework. a) shows cone angles derived from eigenvalues following  and Schmidt diagrams illustrating the distribution of c-axes in the upper and the lower part of the core (projection of c-axes onto the horizontal ice core plane). b) presents the calculated P-wave velocities for the two frameworks for all thin sections (symbols) and the average velocities for each continuously sampled depth interval (lines). c) shows the difference ∆ vp0 = v cx p0 −v ev p0 per sample and per interval average. d) illustrates the velocity change δv cx p0 = v cx p0 (di+1) − v cx p0 (di) between subsequent 10 cm sections at depths di. e) shows the RMS velocities which were calculated from the averaged velocities for layers centered around the measurement intervals (alternating shading); S-wave velocities refer to the upper x-axis and P-wave velocities to the lower x-axis. calculated velocities, while change is used here for the azimuth-dependent observations. In this section we focus on, and begin with, the results of the alpine ice core KCC to demonstrate the relevance of the cx framework for asymmetric fabric.

Non-vertical incidence at KCC
The change of the P-wave velocity with increasing angle of incidence and rotated seismic plane as calculated with the azimuth sensitive cx framework is displayed in Fig. 7. The seismic plane is rotated around the ice core axis in steps of ∆ϑ s = 45 • .

5
Several core pieces were presumably rotated relative to the majority of all ice core pieces during processing to optimise the aliquot cutting. The rotation was estimated and the data is corrected accordingly before applying the cx framework algorithm.
The influence of the asymmetry of the anisotropic fabric in the deeper part of KCC appears very clear. For some layers a spread of velocities of up to 120 m s −1 is observed for a given angle of incidence when considering different seismic plane azimuth angles.   The difference between the framework velocities v cx p (ψ) − v ev p (ψ) is shown in Fig. 8, for ϑ s = 0. v ev p (ψ) is invariant under the rotation of the seismic plane in case of cone fabric. Thus, only the cx velocity is changing with rotation. The change from v cx p (ψ, ϑ s = 0 • ) to v cx p (ψ, ϑ s ) is shown for ϑ s > 0 • . The difference in P-wave velocity when comparing the calculation frameworks reaches up to ± 190 m s −1 for the bottom layer and ± 50 − 100 m s −1 for most depths below 48 m ice depth for various incidence angles and seismic plane azimuth angles.

5
Although the slower S-waves are not routinely acquired during seismic imaging in polar environments, they provide a better resolution and are of special interest for the study of the elastic properties of ice from traditional seismic reflection profiles (Picotti et al., 2015). In particular, the splitting of the shear waves can provide unique information about the anisotropy of the medium (Anandakrishnan et al., 1994;Smith et al., 2017). In the case of the evidently asymmetric fabric of the KCC ice core we observe S-wave splitting of well above 200 m s −1 in the lower half of the ice core with a maximum value of 281 m s −1 . The 10 strength of the shear-wave splitting for a particular seismic incidence angle changes when rotating the seismic plane. Figure 9 shows the difference between SV-and SH-wave velocities for non-vertical incidence (for ϑ s = 0 • ) and investigates how the difference between the S-wave modes changes when rotating the seismic plane. The initial difference v cx sv − v cx sh at ϑ s = 0 • is low for small angles except for the bottom samples. It reaches more than 200 m s −1 for angles > 40°. For specific azimuth angles the change in shear-wave splitting reaches about 200 m s −1 for many depths below 48 m ice depth for incidence angles   Table 2. Where the core orientation is insufficiently known and corrected with respect to neighbouring core pieces, vertical variation in velocity in dependence of the incidence angle may not be a true variation.

Non-vertical incidence at EDML
For the EDML core, no information on the core pieces' azimuth angle relative to the ice sheet or to each other is provided.
However, it is assumed that no sudden short-scale change in the flow regime can occur. Thus, abrupt offsets in girdle orientation must be caused by the unnoticed rotation of core pieces. Prior to the application of the cx framework any misorientation needs to be corrected, or at least highlighted, to avoid misinterpretation of the results from applying the cx framework for seismic 5 velocity calculation considering non-vertical incidence angles. For the EDML data set the orientation of several single thin sections was corrected according to the girdle orientation of the neighbouring thin sections. For instance, a sharp change of girdle direction of about 45°in 1705 m (Weikusat et al., 2017) could not be corrected and has to be kept in mind when looking at the velocity calculation results for non-vertical incidences.
As the ev framework does not aim to include the orientation of the girdle, the velocity in girdle fabric is considered as invariant 10 under the rotation as well. We therefore only assess the change in P-wave velocity v cx p as calculated with the cx framework (Fig. 10). The respective figures for S-wave velocities can be found in Kerch (2016). The highest seismic P-wave velocities (∼ 4028 m s −1 ) calculated with the cx framework for non-vertical incidence are found deeper than 2000 m, where the fabric anisotropy is strong, for incidence angles below 20°. Seismic P-wave velocities are underestimated by the ev framework by max. 131 m s −1 and overestimated by max. 84 m s −1 . The difference is only small (± 20 m s −1 ) for cone fabric in the upper 15 part (0 − 800 m). The highest change is apparent for the lower part of the girdle fabric, below the earlier mentioned sudden rotation of the dominant azimuth direction, and for cone fabric in the deep part of the core. There, the change in interval velocity can exceed 100 m s −1 for some seismic azimuth planes as compared to the defined 0°-plane.

Evaluation of the cx framework
The cx framework provides a refined approach for the use of fabric information to obtain seismic velocities in ice. By including 5 all the c-axis observations, instead of using eigenvalue representation, we keep information that is lost with the ev framework and we avoid the approximation to the true c-axis distribution by deriving opening angles. We average on the crystal scale to obtain the full elasticity tensor for the polycrystalline ice. This is, to our knowledge, the first time this approach has been applied to actual ice-core fabric data. Recent work from Vaughan et al. (2017) presents P-wave velocities from cryo-EBSD data on artificial ice using the MTEX toolbox (Mainprice et al., 2011). 10 By using the fabric data from thin sections we acknowledge the uncertainty which arises from sampling with a relatively small sample size. We use less than 1 % of the ice core EDML and 11 % of the ice KCC to infer the fabric development in the ice cores.
There is currently no comprehensive data available to investigate the sampling effect on real ice. As we are concerned with the comparison of theoretical seismic velocities calculated from the same fabric data, we assume that the sampling uncertainty can be neglected. For the comparison with measured seismic data the uncertainty needs to be considered, and appropriate density 15 and temperature corrections are required.
The observed variation in eigenvalues in the EDML ice core (Fig. 4a) can partly be attributed to a systematic deviation between horizontal and vertical thin sections which is a bias produced by the older fabric analyser model G20 (Weikusat et al., 2017).
Both the short-scale variation in the high resolution intervals in the EDML ice core and in the KCC ice core are not affected by the instrument bias. The cx framework propagates this systematic variability more that the ev framework, with a higher standard 20 deviation for the EDML depth interval 0 − 1785 m (Table 1), illustrating the higher sensitivity to small fabric differences. The measurement of c-axes from thin sections with the instrument and the subsequent automatic image processing, which provides the c-axis angles as an average per grain, contribute to a smaller extent to the overall uncertainty at a level which is difficult to quantify. However, the processing routine (Eichler, 2013) has proven to provide robust results with respect to minor changes in the procedure and algorithm. 25 The currently employed algorithms for the calculation of seismic velocities in ice polycrystals on the crystal scale (including this study) do not consider any possible effects on the grain boundaries. For laboratory measurements the difference in stress on a polycrystalline ice sample as compared with in-situ conditions can affect the degree to which grains are bonding and, thus, the elasticity (Helgerud et al., 2009). Processes like grain boundary sliding are currently explored in the context of deformation mechanism on the micro-scale (pers. comm. E.-J. Kuiper, Univ. of Utrecht, 2017) but can also influence the elastic behaviour 30 of ice (Elvin, 1996). These issues should be addressed for future applications employing ultrasonic methods for the estimation of elastic properties of ice.
The lack of knowledge about the dispersion of seismic waves in ice introduces an unknown uncertainty to the calculation based on a monocrystal elasticity tensor that was measured in the laboratory by means of ultrasonic sounding. Again, for the application of ultrasonic methods, which operate in the same frequency range, this uncertainty can be neglected. The connection of fabric and seismic velocities on the crystal scale we present here complements this advancing field of study.
We have shown in section 2.4.3 that the ev and cx frameworks differ slightly in the case of vertically symmetric cone fabric for vertical incidence and large cone angles. This type of fabric can commonly be expected in the shallower depth of any 5 glacier where vertical compression is dominant. We conclude that the observed deviation in the vertical P-wave velocity profile (EDML) between the ev and the cx velocity for cone fabric could partly be attributed to this inherent difference between the frameworks.
In the case of asymmetric c-axis distributions, as observed in the KCC ice core, we obtain large differences between the interval velocities of the two frameworks, resulting in a detectable difference between the RMS velocities at the bedrock which 10 is relevant for the depth conversion. We can confirm the assessment of Voigt-Reuss bounds to lie below 1 % (for P-wave) in our study.
An advantage of the cx framework is the lack of a need for the fabric classification, thus eliminating artificial discontinuities.
In synthetic seismograms derived from the modelled velocities, such artefacts could result in artificial reflectors and, thus, lead to false interpretations. The example of high resolution sampling in the EDML ice core demonstrates the importance of this 15 advance, allowing us to separate the true high variability in seismic velocities from the artificially enhanced variability. This finding could, however, be used to tune the threshold values for the fabric classification in the ev framework.
Potentially, our framework can be used in principle for the development of inverse methods to derive the fabric distribution from seismic velocities. Following experience from other fields of active seismology, this would, first, most likely require comprehensive data sets suitable for full-waveform inversion not yet available for glaciological applications, and, second, 20 some simplifying assumptions on the distribution of crystal fabric, e.g. in terms of considered symmetries. The framework we presented allows to quantify the potential effect of simplifying assumptions and could help to more accurately specify covariance matrices, thus enabling the quantification of uncertainties coming along with the results produced by application of an inverse method.
Azimuth-sensitive seismic velocities 25 The cx framework we developed and employed in this study takes into account the asymmetry of anisotropic fabric, with respect to the vertical. This is especially relevant for glacial environments with a complex flow pattern, for example in sloping mountain glaciers, fast-flowing polar outlet glaciers (Hofstede et al., 2018) and ice streams (Smith et al., 2017). For such sites the approximation of the fabric by opening angles centered around the vertical can deviate much more from the reality than for sites that are located in the vicinity of an ice divide. It becomes evident from the presented KCC case study that the 30 azimuthal change of the fabric and the resulting velocities are not negligible. On the contrary, the velocities calculated with the cx framework for non-vertical incidence angles from an arbitrary seismic azimuth can change strongly for both P-and S-waves and the associated shear-wave splitting. If the velocity depth profile changes continuously, as is illustrated in Fig. 7, 8 and 9, this should, in principle, be resolved in seismic surface profile data from different seismic azimuth directions, providing information about the (asymmetric) crystal anisotropy evolution with depth.
A requirement of the cx framework is the dependency on accurate core orientation information, i.e. the orientation of the fabric distribution in the equatorial plane has to be known for the consecutive fabric samples. To date, the retrieval of ice cores with known azimuth remains a challenge. Hence, the uncertainty in the calculation of seismic velocities is much larger in the vertical direction than under azimuthal rotation. On the other hand, analysing seismic data with azimuthal resolution around 5 an ice core drilling site could provide the information to improve the reconstruction of the core orientation. The appearance of a non-symmetric fabric might also be induced by inclined drilling. Ideally, to be able to link calculated and measured seismic velocities a possible inclination of the ice core with respect to the vertical and to the horizontal seismic profile should be considered.
Rapid velocity changes over short vertical distances 10 We use COF measurements on a submetre scale for our analysis of seismic velocities. The results suggest the existence of closely spaced reflective surfaces for elastic seismic waves (and also radar waves). The relevance of the presented analysis for real seismic data is based on the major assumption of a laterally extended and coherent fabric layering on the scale of the first Fresnel zone (Drews et al., 2012). Although fabric layering is regularly observed in the KCC ice core and also in the continuously sampled depth intervals in EDML, it is still unclear how representative these short-scale variations are for 15 both the close vicinity and a larger region in a glacier. However, evidence has been presented for abrupt COF changes as a frequent cause of seismic reflectivity (Horgan et al., 2011). Other studies do not observe such a high reflectivity due to COF but identify a high degree of gradually evolving fabric anisotropy (Picotti et al., 2015) or single strong reflections associated with transitions in fabric classes, e.g. from cone to girdle . The coherence of thin layers with distinct fabric will largely depend on the unresolved question of how they evolve exactly. If the short-scale fabric stratigraphy is largely governed 20 by local conditions and heterogeneous small-scale deformation, possibly resulting in "layer roughness" (Drews et al., 2009), no coherent structure is to be expected . In this case, it should be challenged, how representative the elastic properties derived from thin sections are, and the question arises, how incoherent short-scale fabric changes alter the rheological properties of the bulk. It can be hypothesised that under the increasing influence of large-scale shear deformation in the deeper regions of the glacier coherent fabric layers might develop. Accordingly, more seismic reflectivity should be expected in depth 25 and from more dynamic settings, as proposed by Horgan et al. (2011). Eisen et al. (2007) show that transitions in COF in the deep ice can be followed with radio-echo sounding over longer horizontal distances of several kilometres. However, variations in seismic velocity on short vertical scales cannot be resolved with conventional surface-based seismic techniques with large wavelengths of the order of 10 m, depending on the source of the seismic waves and the sounding depth. Still, Hofstede et al. (2013) obtain numerous laterally continous reflections at Halvfarryggen, Antarctica. They suggest that closely spaced layers 30 ("stacks") of varying fabric, possibly as have been observed in this study, are the major cause for the reflections.
Far more fabric data than is currently sampled in ice core studies, is required to pursue this hypothesis in the future. To this end, ultrasonic methods can be applied in ice core boreholes (Bentley, 1972;Gusmeroli et al., 2012) to infer crystalorientation fabric in situ. Although the interpretation of these data is not straightforward (Maurel et al., 2015), it is currently the only technique that is capable of a continuous fabric measurement. However, a sonic pulse samples the volume around the borehole (∼ 2 m 3 , Gusmeroli et al., 2012), which means the method is not azimuth-sensitive. While it cannot provide the two-dimensional microstructure nor exact and highly resolved fabric information, it can help to bridge the gap between laboratory-based interval fabric measurements and large-scale seismic data.
Following the findings of our study we recommend for seismic data acquisition in the field to (1) consider polarimetric survey 5 setups (with two or even more crosslines) with both reflection and wideangle measurements, and to (2) focus on accurate travel time recordings at high source frequencies. This should be supported by 3-component vertical seismic profiling where boreholes are available. Also, S-waves should be acquired as they provide useful information on crystal anisotropy due to shear-wave splitting. On the crystal scale, we suggest to include the investigation of the possible influence of variations in grain size for the elastic wave propagation in polycrystalline ice, which is currently not considered for theoretical calculations, to 10 complement recent work on the temperature dependency of elastic properties (Vaughan et al., 2016). Ongoing microstructure studies on both alpine and polar ice provide indications of considerable vertical short-scale variability in grain topology. Recent laser ultrasound measurements on ice have provided first high-resolution data (Mikesell et al., 2017) and promise further advances towards understanding and efficiently measuring the elastic properties of polycrystalline ice on the crystal scale.

15
The presented cx framework contributes to the understanding of the propagation of seismic velocities in glacial ice by deriving bulk elastic properties on the crystal scale. Based on anisotropic fabric from two ice cores, we showed that the fabric classification scheme in the ev framework can mask the true velocity variability by producing artificially enhanced peaks in the velocity profile. By applying the cx framework we extract the velocity variability that is caused by the actual fabric variability.
The velocity difference between the cx and ev frameworks is larger for the alpine than for the polar core. This suggests that 20 the ev framework provides a good enough approximation for the polar site, situated on an ice divide, for the current degree of seismic resolution and interpretation of physical properties, not considering the artificial discontinuities, but is not adequate for the alpine site.
We found that the azimuthal change in P-wave velocity and shear-wave splitting can be as large as ∼ 200 m s −1 . We conclude that the possibility of an azimuthal asymmetry of the fabric distribution should be considered when planning seismic surveys 25 (e.g. polarimetric profiles around a drilling site) as well as for the calculation of seismic velocities from fabric data. This also offers an opportunity to further constrain azimuthal ice-core orientation.
The results of our study demonstrate for the first time that a short-scale variability in anisotropic fabric as observed in these polar and alpine ice cores causes a corresponding high short-scale variability in seismic interval velocities. Current laboratory fabric measurements from an ice core drilled on an ice stream also show early indications of a high fabric variability and 30 unexpected fabric types (pers. comm. J. Eichler, 2017), offering an ideal target for extending this study to an environment with another deformation regime. Based on the presented evidence in this study the next steps should include the investigation of how a succession of short-scale fabric layers could induce englacial reflections as has been reported and hypothesised in earlier studies (Horgan et al., 2011;Hofstede et al., 2013).
As conventional surface-based seismic surveys are not likely to resolve this short-scale variability, ultrasonic techniques for borehole and laboratory studies could be the solution to both issues of lost core orientation and low resolution. For this emerging field of applications, we offer further insight into what to expect from crystal-orientation fabric anisotropy in ice. Equally, our results can provide context for data collected with frozen-in seismometers in boreholes, where evidence for shear-wave 5 splitting on non-vertical ray paths was found (pers. comm. David Prior, 2018). Lastly, we want to highlight that while the depth scale of the KCC ice core differs from that of the EDML ice core by a factor of 1/35, the presented case study is another example (Eisen et al., 2003;Diez et al., 2014) of the importance of mid-latitude high-altitude glaciers as in-situ laboratories to study fundamental processes in glaciers.
Data availability. The fabric and eigenvalue data sets for the ice cores KCC (Kerch et al., 2016a, b) and EDML (Weikusat et al., 2013a, b, 10 c, d;Weikusat et al., 2018) are published in the open-access database PANGAEA ® and available upon request.

Appendix A: Tensor transformation
A fourth-order tensor rotation is expressed as: The general rotation matrix in three dimensions is given by the cosines between the axes of local {p, q, r} and global {x, y, z} coordinate frame: cos(x, p) cos(x, q) cos(x, r) cos(y, p) cos(y, q) cos(y, r) cos(z, p) cos(z, q) cos(z, r) It is possible to express both rotations in a single rotation matrix (as done by Maurel et al., 2015, section 3).
20 Table 2. Summary of the results from the seismic velocity comparison between frameworks. Values are calculated depth-profile average (with standard deviation) and/or extreme (±) interval velocity differences (other than RMS) for incidence angles of 0 − 70°. Negative values indicate smaller velocities from the cx framework relative to the ev framework. The reported extreme values can be influenced by outliers from the general trend. (*) Example of how to read the table: "For a specific seismic plane azimuth ϑs, an incidence angle ψ and a specific interval at the KCC site the SV-wave velocity as calculated with the cx framework is found to be 279 m s −1 higher than is calculated with the ev framework which is the maximum difference for any combination of ϑs, ψ and depth."