Propagation of elastic waves through textured polycrystals: application to ice

The propagation of elastic waves in polycrystals is revisited, with an emphasis on configurations relevant to the study of ice. Randomly oriented hexagonal single crystals are considered with specific, non-uniform, probability distributions for their major axis. Three typical textures or fabrics (i.e. preferred grain orientations) are studied in detail: one cluster fabric and two girdle fabrics, as found in ice recovered from deep ice cores. After computing the averaged elasticity tensor for the considered textures, wave propagation is studied using a wave equation with elastic constants c=〈c〉+δc that are equal to an average plus deviations, presumed small, from that average. This allows for the use of the Voigt average in the wave equation, and velocities are obtained solving the appropriate Christoffel equation. The velocity for vertical propagation, as appropriate to interpret sonic logging measurements, is analysed in more details. Our formulae are shown to be accurate at the 0.5% level and they provide a rationale for previous empirical fits to wave propagation velocities with a quantitative agreement at the 0.07–0.7% level. We conclude that, within the formalism presented here, it is appropriate to use, with confidence, velocity measurements to characterize ice fabrics.

the axes of tensional stress (or equivalently the a-axes in the basal plane tend to be aligned with the tensional axis). Among the many possible fabrics, girdle fabrics (defined below) are characteristic of ice encountered in convergent flow regions, where the vertical axis corresponds to the compression due to gravity and with a tensional axis due to a given ice flow direction [19]. This is the case for ice along a ridge, a few hundred metres from a geographical dome; here, the dissymmetry in the surface slopes along and perpendicular to the ridge can induce a horizontal tension component in the strain field [20].
Until recently, this fabric was mostly determined by recovering the individual grain orientations on an ice thin section (dimensions of about 10 × 10 cm 2 or 100-500 grains) and this method is certainly the most precise when a statistically sufficient number of grains can be found in the thin section. An alternative approach was initiated in the 1970s by Bennett [21] (see also [22, ch. 6]) and Bentley [23]: the use of borehole sonic measurements to determine in situ depthcontinuous fabric. This method has been tested at Dome C, East Antarctica [24], and it reveals the sensitivity of both the longitudinal and transverse waves with depth, and thus, possibly, the variation of ice fabrics with depth. However, uncertainties appear in the inversion process as used in [24], which are partially due to the measurements themselves but which are also due to the model of the sonic velocities that is employed (this latter point will be discussed in this study). Thus, there is a need to provide accurate models for these media to answer the question of whether or not sonic logging is able to discriminate the different fabrics and, within a given fabric, to determine the degree of anisotropy. This is the goal of this paper, where the propagation of elastic waves in polycrystals with cluster and girdle fabrics is studied. Note that the presence of stress may affect the wave speeds (acoustoelastic effect) but it is excluded from consideration in our study.
Section 2 presents the definition and derivation of the second-order orientation tensor, whose maximum eigenvalue is a common characteristic of the degree of anisotropy of ice polycrystals [17]. This is done in order to get closed forms of the orientation tensors by defining simple orientation distribution functions (ODFs) of the c-axes for cluster and girdle fabrics. The realism of these ODFs with respect to existing physical modelling of the fabrics (see [25,26]) will be discussed elsewhere. In §3, we derive the Voigt matrices for the considered fabrics, starting by averaging the tensors of a single crystal expressed in an arbitrary coordinate frame. The expected effective anisotropies are obtained: hexagonal with vertical transverse isotropy (VTI) for clusters, orthorhombic for partial girdles and hexagonal with horizontal transverse isotropy (HTI) for thick girdles. Direct application of these calculations consists in setting the determinant of the Christoffel matrix equal to zero and this is presented in §4 for an arbitrary direction of the wave propagation. Except in the case of the partial girdle, closed forms for the velocities for the P-wave and for the S-waves are obtained. Because our calculations are primary motivated by their application in glaciology, we collect in §5 elements of discussion in this particular context. This concerns notably the inspection of the accuracy of our approximation, which assumes that the local deviation δc ijkl of the elastic stiffness is small with respect to its average value. In the Introduction and throughout the paper, we use interchangeably the term anisotropy for the anisotropy (in its purest form) of the single crystals and for the effective anisotropy characteristic of the fabric, both resulting from a homogenization process (at the scale of many atoms and at the scale of many grains). Finally, when used, the values of the elastic constants for ice single crystals are taken from [21] and reported below ice single crystal

Orientation distribution function and second-order orientation tensor
It is usual to characterize the anisotropy of ice polycrystals using the second-order orientation tensor a defined by a ij = c i c j , (2.1) where . denotes ensemble averages. That is, an average over all possible realizations of the directions of the c-axis [27]. In this section, we define the possible orientations of the c-axis using the ODF associated with the angles θ and ϕ, the co-latitude and longitudinal angles for cluster and girdle fabrics (figure 1); within girdle fabrics, we have restricted our study to the cases of partial and girdle fabrics (e.g. [28]). These three fabrics are relatively simple as they are defined by a single parameter in the ODF (the parameter is θ 0 for the cluster fabric and for the partial girdle and it is ξ 0 for the thick girdle; figure 2). Obviously, more complex fabrics may involve two or more parameters [26,29].
The largest eigenvalue of a, denoted a 1 in the following, will be used further to characterize in a unified way the strength of the anisotropy of the resulting structure (limiting cases being a 1 = 1/3 for a fully isotropic structure and a 1 = 1 for a fabric with a single maximum; that is, when all the c-axes have the same direction) [27]. From figure 1, one can anticipate that cluster and thick girdle fabrics correspond to a resulting structure with hexagonal symmetry, whereas the partial girdle corresponds to a resulting structure with orthorhombic symmetry. In the case of the clusters, the structure is invariant by rotation along the vertical axis e 3 , thus it is isotropic in the transverse plane (often referred to as VTI for vertical transverse isotropy). In the case of a thick girdle, the structure is isotropic in a plane perpendicular to a horizontal axis (e 1 in figure 2c, referred to as HTI for horizontal transverse isotropy).

(a) Cluster fabric
Clustered fabrics, figure 1a, correspond to c-axes being oriented within a cone about the vertical axis e 3 and we denote θ 0 the opening angle this cone makes with the vertical. This configuration is naturally described by the spherical angles (ϕ, θ ) (figure 1), with the c-axis defined by with ϕ ∈ [0, 2π ] and θ ∈ [0, θ 0 ]. The probability distribution function p(ϕ, θ) is, as in [22,27], with H 1 (θ ) = 1, for 0 ≤ θ ≤ θ 0 and zero otherwise. For the spherical angles, the average of any quantity A is defined by It follows, using equations (2.1)-(2.4), that the second orientation tensor a is given as (denoting (sθ, cθ ) for (sin θ , cos θ ), respectively, (cϕ, sϕ) for (cos ϕ, sin ϕ))    in agreement with [22,27]. We have ordered the eigenvalues a 1 > a 2 = a 3 , and a 2 = a 3 is obtained because the directions along e 1 and e 2 are equivalent for the VTI symmetry. The cluster fabrics cover all the degrees of anisotropy, from full isotropy with a 1 = 1/3 (θ 0 = π/2, corresponding to all possible orientations of the c-axes) to the maximum anisotropy, with a 1 = 1 (θ 0 = 0, corresponding to all the c-axes being aligned with e 3 , often referred to as a single-maximum fabric).

(b) Partial girdle
This configuration is two dimensional (figure 2b), and the polar angle θ is sufficient to describe the orientation of the c-axis for θ ∈ [0, θ 0 ] (the ± signs correspond to ϕ = π/2(+) and 3π/2(−)). The problem being two dimensional in (e 2 , e 3 ), the projection onto (e 1 , e 2 ) is not needed and the average is defined with θ being considered as a polar angle. For any quantity A, we now have and we have kept ϕ in the average as the plane of rotation of the c-axes has to be defined by a delta-function to fix the ϕ-values, with  Figure 3. Parametrization angles (φ,θ ) used for the thick girdle and the correspondences between the coordinate frames (ê 1 ,ê 2 ,ê 3 ) and (e 1 , e 2 , e 3 ).
(c) Evolution of the anisotropy with the a i Figure 4 shows the variation of the a i , i = 1, 2, 3, as a function of θ 0 for the cluster fabric and as a function of θ 0 and ξ 0 for the partial and thick girdles, respectively. Inspecting the value of a 1 allows us to characterize the degree of anisotropy of a polycrystalline structure in a unified way, that is, independently of the considered fabric.

Elasticity tensors of single crystals and textured polycrystals
In this section, we derive the elasticity tensor c ijkl of a single crystal when expressed in an arbitrary coordinate frame, with the c-axis being parametrized by the two spherical angles (θ, ϕ). The result is further used to derive the elasticity tensors for the cluster and girdle fabrics, describing average macroscopic properties of the polycrystalline materials. We use the elasticity tensor in Voigt's notation [30], C IJ with the standard correspondences   Note a non-conventional form of the rotation matrix when compared with the rotation matrix that takes Cartesian to spherical coordinates; indeed, in this usual notation, the vector e 1 is transformed into the radial vector, which is the first vector in spherical coordinates, but for convenience it is the third vector e 0 3 in the present case. This produces a c-axis oriented along the vertical direction e 3 for ϕ = θ = 0.
When expressed in the local frame (e 0 1 , e 0 2 , e 0 3 ), the single crystal with hexagonal symmetry is described by the elasticity tensor c 0 ijkl , written as C 0 in Voigt's notation, In a frame (e 1 , e 2 , e 3 ), the elasticity tensor c abcd is deduced from c 0 ijkl using according to the transformation laws for tensors (see, for example, [31] for the calculation of anisotropic properties from texture data using an open source package). The corresponding relation between the Voigt matrices C 0 and C is in general involved, and the use of the Bond matrix is often preferred [32]. We do not use the Bond matrix formalism but, owing to the hexagonal symmetry in C 0 , we use a reasonably tractable expression, similar to the one used in [33], These terms are calculated and we obtain, for the diagonal terms, C 44 = (A + C − 2F)sθ 2 cθ 2 sϕ 2 + L[sϕ 2 (cθ 2 − sθ 2 ) 2 + cθ 2 cϕ 2 ] + Nsθ 2 cϕ 2 , (b) Elasticity tensors of textured polycrystals: cluster and girdle fabrics Next, we derive the averaged elasticity tensors that are needed for the description of the effective macroscopic properties of polycrystals with cluster or girdle fabrics. This is done by averaging the Voigt matrices of single crystals (equations (3.6) and (3.7)) using the ODF in equations (2.3), (2.9) and (2.15).

(i) The clustered fabric
The averages of the Voigt matrix with elements C ij in equations (3.6) and (3.7) are averaged according to with the ODF for the clustered fabric in equation (2.3). It is shown in appendix A that the resulting Voigt matrix is associated with a structure with VTI, as expected (1) The single maximum fabrics for θ 0 = 0 (X = 3, Y = 2) lead to C 11 = A, C 33 = C, C 44 = L, C 66 = N and C 13 = F; the effective medium is, as expected, the same as a homogeneous medium with a unique orientation of the c-axis along e 3 , in agreement with the expression of the Voigt matrix C 0 in equation (3.3). (2) For θ 0 = π/2 (X = 1, Y = 0), the average is done for c-axes varying randomly over all directions, resulting in an effective isotropy; we obtain effective isotropy and C 13 iso = C 11 iso − 2 C 44 iso , resulting in effective Lamé coefficients These elastic constants in the isotropic case agree with previous derivations, referred to as the Markham derivation [2] (see also [33, eqn (4.b)]).

(ii) Partial girdle fabrics
The average of the Voigt matrix is carried out according to the two-dimensional configuration with p being defined in equation (2.9). Calculations are straightforward (see details in appendix B) and the resulting Voigt matrix is associated with a structure with orthorhombic symmetry ⎛ where we have defined S 2 ≡ sinc2θ 0 and S 4 ≡ sinc4θ 0 , (3.17) and sinc x ≡ sin x/x (and sinc 0 = 1). As previously noted, partial girdles always have isotropic structure. For θ 0 = 0 (S 2 = S 4 = 1), we recover the hexagonal symmetry of a single crystal with the c-axis along e 3 (equation (3.3)). Also, for θ 0 = π/2 (S 2 = S 4 = 0), the effective medium has hexagonal symmetry with HTI corresponding to a perfect girdle, with C 23 (3.18)

(iii) The thick girdle
The elasticity tensor for the thick girdle is derived using the same trick that was previously used for the derivation of the second-order orientation tensor in §2b(i). It allows for a straightforward use of the expression of the elasticity tensor derived in §3a for single crystals. Namely, It is shown in appendix C that the resulting Voigt matrix is associated with a structure with HTI symmetry (with respect to e 1 ), ⎛ (4) Thick girdles (with a 1 between 1/3 and 1/2) have in general a lower degree of anisotropy than partial girdles (a 1 between 1/2 and 1). Nevertheless, the two structures coincide when realizing a perfect girdle (partial with θ 0 = π/2 and thick with ξ 0 = 0); in that case, the effective parameters in equation (3.22) (with sξ 0 = 0) coincide with those in equation (3.18). (5) The thick girdle becomes isotropic ξ 0 = π/2 (sξ 0 = 1), and we indeed recover the effective isotropic parameters (equation (3.12)).

Wave propagation in a textured polycrystal
In this section, we derive the velocities of elastic waves propagating in polycrystals (with cluster or girdle fabrics). Using the Voigt average, the calculation is straightforward using the expressions of the averaged elasticity tensors established in the previous section. Velocities are given for an arbitrary direction of the wave propagation, afterwards we will focus on the case of a propagation along the vertical direction e 3 ( §5). This is because our motivation comes from application to sonic logging measurements in ice cores, as previously commented. As a warm up, the case of propagation in single crystals is first considered; this also allows for a validation of our expressions in equations (3.6) and (3.7). Indeed, these are given as a function of the longitudinal angle ϕ of the c-axis; with a wave propagation along the vertical direction, the problem becomes invariant by rotation around e 3 so that the wave velocities have to be found independent of ϕ. Then, the propagation in polycrystals using ensemble averages is briefly recalled and the wave velocities are derived.
(a) The case of propagation in a single crystal (i) Exact expressions of the velocities for single crystals The propagation of monochromatic waves of frequency ω in single crystals is described by the wave equation In our case, the problem will be solved considering the wave propagating along the vertical e 3 direction for a c-axis having arbitrary direction (thus, without loss of generality at this stage). This approach slightly differs from previous ones where the elasticity tensor is written in the crystal frame (e 0 1 , e 0 2 , e 0 3 associated with the principal directions of anisotropy of the crystal) but where an arbitrary direction of the wave propagation is considered (see [34]; see also [35], which corrects a misprint in the former reference). Obviously, the two approaches are equivalent, and, when propagation in a single crystal only is needed, this latter approach is simpler.
Waves propagating along e 3 correspond to u a = U a e ikx 3 , with k the wavenumber, and equation which admits a non-zero solution for (U 1 , U 2 , U 3 ) if the discriminant of the matrix vanishes, leading to a dispersion relation D(ω, k) = 0. The dispersion relation admits in general three solutions for the frequencey ω as a function of wavenumber k. They correspond to three eigenvectors: one longitudinal wave and two transverse waves. One obtains the dispersion relation in the form of a Christoffel equation with C, the Voigt matrix, associated with the elasticity tensor c abcd , in equations (3.6) and (3.7). Obviously, for a single crystal, the use of ϕ is useless as the problem is invariant by rotation around e 3 and this has been done in [34]. We have checked that the solutions of the dispersion relation do not depend on ϕ (this has been done numerically by solving the eigenvalue problem associated with equation   Inboth(a,b), red curves refer to the exact values (equations (4.6)), and green curves refer to Thomsen's approximation (equations (4.7)).
(ii) Approximate expressions: the case of an ice single crystal The variations of the wave velocities as a function of θ (equations (4.6)) are shown in figure 5 (with values from (1.1)).
As ice single crystals are only weakly anisotropic, approximate expressions for the velocities are often used. In 1986, Thomsen [35] remarked that a common approximation is to neglect the anisotropy because 'the mathematical equations for anisotropic wave propagation are algebraically daunting' and he proposed a simplified version, which is commonly used in the context of geophysics. They are

(b) Effective elastic velocities in textured polycrystals
When an ensemble of grains with different anisotropy directions is considered, it is difficult and not useful to describe wave propagation through a single, particular, distribution. Rather, as discussed in the Introduction, ensemble averages can be considered, leading to a description of the properties averaged over all possible realizations of the disorder. In this paper, we shall write the elastic constants as a sum of an average plus deviations from that average: c abcd = c abcd + δc abcd , so the wave propagation in a given realization reads where δc abcd is the local deviation of the elasticity tensor from its average value, and it is assumed to be small; that is, if c and δc are the typical magnitudes of c abcd and δc abcd , we define the small parameter = δc/c. Equation (4.9) is the wave equation through a medium with average properties given by the effective stiffness tensor c abcd and submitted to a 'potential' V ad = (∂/∂x b )δc abcd (∂/∂x c ) which reflects the fluctuations of the real medium with respect to the averaged one. If 1, the wave propagation is mainly determined by the simplified equation and any particular realization of the propagation will be close to this average, up to 2 . Indeed, because δc abcd = 0 by construction, the potential V ad has no contribution at first order in . At second and higher orders, it produces a hierarchy of corrections, including wave attenuation and backscattering (e.g. [1,15]). The successive iterations correspond to the successive approximations of the Dyson equation. Below, the zeroth-order approximation is considered and the accuracy of this zeroth-order will be discussed in more detail in the case of ice polycrystals in §5c. The elasticity tensor c ijkl has been expressed in the previous section using a frame where the direction e 1 , e 2 or e 3 is an axis of symmetry. To get the dispersion relation in a general case, one has to consider now an arbitrary direction of the wave propagation given by the wave vector k = (k 1 , k 2 , k 3 ) (and we defined k using the co-latitude and longitudinal angles Θ and Φ; figure 6). Our Voigt matrices for the considered cluster and girdle fabrics are at least orthorhombic (equation (3.15)), and we restrict the dispersion relations to this symmetry (which includes the hexagonal symmetry). Within this symmetry and looking for a wave u a = U a e i(k 1 x 1 +k 2 x 2 +k 3 x 3 ) , the dispersion relation takes the form, Depending on the fabric, we will derive in the following sections the corresponding velocities owing to the Voigt matrices previously calculated (equations (3.10), (3.16) and (3.22)).

(i) Cluster fabric
The cluster fabric has hexagonal symmetry, being isotropic in the (e 1 , e 2 ) plane. It is thus sufficient to define the angle Θ that k forms with e 3 . With k 1 = 0, k 2 = k sin Θ and k 3 = k cos Θ, equation (with sΘ, cΘ denoting sin Θ, cos Θ), and with the Voigt matrix defined in equation (3.10), with C 55 = C 44 . We deduce the velocities (the plane of incidence being (e 2 , e 3 )) (4.13) Obviously, these expressions are analogous to the expressions found for single crystals in equations (4.6) as the hexagonal symmetry is the same. Also, the limiting cases are the direct consequences of the limiting cases for the Voigt matrix: for θ 0 = π/2 in equations (3.10) (leading to equations (3.12)), the polycrystal is isotropic with shear and compressional velocities simplifying to isotropic case (4.14) These expressions of the velocities in the isotropic case do not coincide with the leading order in Thomsen's approximation [35] (equations (4.7)), with β = δ = γ 0 (leading to V Th P √ C/ρ and V Th S √ L/ρ). This is because Thomsen defines the anisotropy of polycrystals with respect to the isotropy of a single crystal (a 1 = 1), while the isotropy leading to the velocities above refers to an isotropic polycrystal with a 1 = 1/3; this will be discussed further in §5a. The dependence of the largest S-wave velocity and of the P-wave velocities on the k-direction of propagation is illustrated in figure 7 for different a 1 values. For a value of a 1 between 1 and 1/3, the range of the Θ-dependent velocities decreases to a single value in the isotropic case. For a 1 = 1, we recover the angular dependence of single crystals (

(ii) Girdle fabrics
The case of partial girdles does not lead to substantial simplifications of the dispersion relation, as one has to account in general for the three components of k. The dispersion relation has to be solved for a given Θ and Φ, the co-latitude and longitudinal angles of k in (e 1 , e 2 , e 3 ), and no tractable expression can be given. Closed forms of the velocities will be considered in §5 for waves propagating along the vertical direction. Nevertheless, it is straightforward to solve the dispersion relation, for instance by determining the eigenvalues of the Voigt matrix divided by ρ, with V 2 = ω 2 /k 2 . This has been done numerically and the results for a 1 = 0.8 and 0.5 (perfect girdle with symmetry by rotation around e 1 ) are reported in figure 8.
For a thick girdle, the symmetry is hexagonal with isotropic behaviour in the (e 2 , e 3 ) plane. Using the notations in figure 6b, we can use directly the result of equations (4.16), with the Voigt matrixĈ. Next, using the correspondences in equations (3.20), we get, for k 2 = 0, k 1 = k cosΘ and k 3 = k sinΘ (that is, forΦ = π/2), the dispersion relation where we have used that C 33 = C 22 , C 66 = C 55 and C 13 = C 12 (equation (3.21)). We deduce the velocities ρV 2 SH = C 55 cos 2Θ + C 44 sin 2Θ , and

Comments and comparison with previous work on ice polycrystals
In this section, we focus more specifically on the application to ice fabric characterization using a sonic log of the deep boreholes, as tested during January 2011 at Dome C, East Antarctica, by Gusmeroli et al. [24]. To the best of our knowledge, this is only the second in situ measurement campaign, the first one being the deep drill hole at Byrd Station, Antarctica, during the 1969-1970 field seasons by Bentley [23]. In between, laboratory measurements of the velocity have been performed using core samples: in the Greenland Ice Sheet Project II (GISP2) deep core [36,37] and in the Dye 3, Greenland, deep ice core [38,39]. Except in [24], these studies have confirmed qualitatively, but not quantitatively, the sensitivity of the elastic waves to the degree of anisotropy of cluster fabrics (which are the dominant fabrics at Dye 3, GISP2 and Dome C). The exception in [24] is the quantitative inspection that is proposed, by means of an inversion procedure and Then, an inverse procedure has to be used to obtain information on the fabric in the interrogated zone. This procedure raises several questions that we will address in this section. They are (i) a comparison of our expressions of the velocities with the semi-empirical determination of the elastic velocities, as proposed by Bennett [21] for clustered fabrics because of their averred successful comparisons with experiments. We did not find an expression of the velocities for girdle fabrics in the previous literature that could be used for comparison with our present result. (ii) A discussion on the variation of the velocities from one fabric to another; and within a given fabric, how the velocities change with the degree of anisotropy. It is, of course, necessary that the velocity change that a purported fabric or anisotropy induces is quantified and that this change is larger than the accuracy of the experimental data, if it is to be used as initial data in an inverse problem. (iii) In connection with the previous question, the modelling presented in this paper has to be sufficiently accurate if it is to be used in the inverse process. Thus, the accuracy of the zeroth-order Dyson equation that we employ is addressed. As a prerequisite, note that the sonic measurements are performed at about 20 KHz, so the wavelength in ice is typically 10-20 cm, much larger than the grain size (1 mm-1 cm).

(a) Velocities for vertical wave propagation in an ice polycrystal
As discussed in §1, we have chosen the e 3 axis to coincide with the vertical direction, a convenient choice to analyse sonic logging measurements of antarctic ice. As a first step, we report below the expression of the velocities for propagation of the wave along the vertical e 3 axis. In this case, the expressions from equations (4.11) simplify to V P = √ C 33 /ρ, V S1 = √ C 44 /ρ and V S2 = √ C 55 /ρ (to see this, take k 1 = k 2 = 0 in equations (4.11); it leads to the determinant of a diagonal matrix). Here, we use V S1 and V S2 with V S1 the largest velocity (when different, the two S-waves cannot be identified a priori and the term of 'quasi-' S-wave is preferred, e.g. [40]). For cluster fabrics, from equations (3.10) (or using Θ = 0 in equations (4.16) and For partial girdle fabrics (P.G.), with C 33 , C 44 and C 55 given by equations (3.16), we get P.G.
For the thick girdle (T.G.), with C 33 , C 44 and C 55 being given in equations (3.22), we have T.G.
(5.3) Figure 9 shows the variation of the velocities as a function of the degree of anisotropy of the polycrystal, measured with a 1 (defined in equations (2.6), (2.11) and (2.17)). For comparison, we also present the velocities for the isotropic single crystal, which is the leading order in the Thomsen's approximation (equations (4.7) with β = δ = γ 0), and the velocities that correspond to an ice polycrystal with randomly oriented c-axes (equations (4.14)). As expected, these two limits correspond to velocities obtained at a 1 = 1 (single maximum) and a 1 = 1/3 (isotropic polycrystal), respectively. It is often tempting to neglect the anisotropy of an ice polycrystal and we report here the relative errors in the velocities when using this approximation-when using single-crystal isotropy (Thomsen's leading order): for the cluster 3.5% and 5.7% (P-and S-waves) and for the girdles 3.8%, 7.7% and 3.7% (for P-, SV-and SH-waves, respectively); when using polycrystal isotropy: for the cluster 1.9% and 2.5% (P-and S-waves) and for the girdles 1.4%, 2.3% and 3.7% (for P-, SV-and SH-waves, respectively); thus, polycrystal isotropy is a better approximation.
Next, we come back to the range of variation of the velocities with respect to the degree of anisotropy of the fabrics. We report the relative variation 2(V max − V min )/(V max + V min ) for each fabric for clusters: S-9.0%, P-5.5%. for girdle: S 1 -11.7%, S 2 -7.7%, P-5.2%.
These ranges of variation are small, because an ice single crystal is weakly anisotropic; with the same definitions of the velocity relative variation (for ice single crystal: SV-17.8%, SH-6.1%, P-7.0%). Thus, in all the steps of the measurements and of the inversion process, the errors have to be smaller than, say, 1%.

(b) Comparison with previous results (i) Remark on the importance of the ϕ-average
As discussed in the Introduction, a common approach to derive the effective velocities in a polycrystal consists in averaging the 'time of flight', that is, and, in most of the studies on clusters (as used in [24]), the average is done on the colatitude θ only, that is, the angle between the c-axis and the direction e 3 of the wave propagation, and not on the longitude ϕ. At first sight, this approach appears reasonable. However, there are two x 0 (q 0 = p/2) pitfalls. (i) The times of flight are calculated for each S-or P-wave separately, assuming that the wave energies originally distributed into S-and P-waves remain unchanged; in other words, the re-repartition of the energy between SH-, SV-and P-waves at each grain boundary by mode conversion is disregarded. (ii) More seriously in our opinion, it is assumed that an ensemble of grains with the same θ but different ϕ-value behaves as a homogeneous medium, which is clearly not the case. Although this crystallographic arrangement is artificial, it allows us to exemplify the error made when ignoring the ϕ-dependence. Let us first remark that the structure resulting from such an arrangement has an anisotropy with hexagonal symmetry (effective c-axis along e 3 ), thus the two shear velocities are the same for waves propagating along the e 3 -axis. From appendix A, the velocities for grains with the same θ and different ϕ-values are equivalent to an HTI structure (with symmetry around e 3 ), and, from equations (3.7) and (A 3), we get all c with same θ and the S-velocities are correctly found to be the same. On the contrary, if the average is performed on θ only (and the average in this case consists in selecting a unique value of θ by a delta function), equations (5.5) predict that the ensemble of grains behaves as a single crystal with SH-and SVwave velocities being given by equations (4.6), and this cannot be the case. Incidentally, note that not only the V S (θ ) but also V P (θ ) is false when using equations (5.5) (see V P in equations (4.6) and (5.6)).

(ii) Bennett's prediction
Alternatively to the slowness average, Bennett [21] provides a prediction for V P and V S based on acoustic measurements. The details of the method are given in [21, ch. 5] and in [22, ch. 6] and they are not discussed here. These semi-empirical derivations of the velocities have been shown to be in agreement with experiments so it makes sense to use them as a validation. For a propagation along the e 3 -axis, Bennett's velocities read where (X, Y) are the same as used in equations (3.11), Inspecting equations (5.1) and (5.7), we get the correspondences between our velocities and Bennett's ones. The results of (5.1) are obtained from (5.8) by 9) and the results of (5.7) are obtained from (5.8) by from (5.7) (5.10) Table 1 reports the values of the coefficients defined in equations (5.8), from equations (5.9) using Bennett's values of the elastic constants (A, L, N, F) and density ρ given by (1.1), and from equations (5.10) using the values of (a i , b i ) given above. The agreement is excellent, although less good for the S-wave than for the P-wave, and this will be analysed in more detail in forthcoming work. The resulting agreement on the velocities is 0.07% for the P-wave and 0.7% for the S-wave, without any adjustment (figure 10) (red and black curves). For completeness, we also report in figure 10 the results obtained from slowness averaging as used in [24], omitting the ϕ-average (equations (5.5)) (green curves). This latter case leads to two different S-wave velocities, which is unphysical as the wave propagates along the symmetry axis. Incidentally, there is a slightly more notable disagreement with Bennett [21], with 0.9% for both the P-and the S-wave velocities (when compared with the highest S-velocity).

(c) Inspecting the accuracy of our prediction
In our study, c ijkl corresponds to anisotropic, or textured, effective media, with orthorhombic or hexagonal symmetry. For ice, δc is a small correction as ice single-crystal anisotropy is weak,   .7)) and as used in [24], omitting the ϕ-average. For V S , this latter calculation leads to two (unphysical) velocities (green solid and dotted lines).  (4.9), is justified but the results are accurate up to 2 = (δc/c) 2 only. It is of importance to get an estimate of 2 as the accuracy of the prediction on the velocities has to be sufficient to resolve the typical variations of the velocities. This is reported in figures 11 and 12. First, we report typical fluctuations of the C ij parameters (four terms are shown in figure 11), for a cluster fabric as a function of a 1 , that is, for decreasing cone angle θ 0 ; each point in these plots corresponds to a value of C ij of one grain with an orientation randomly chosen in [0, θ 0 ]. The dispersion of the C ij indicates how far the actual grains are from their effective medium. For all a 1 values, the effective medium has hexagonal symmetry, thus vanishing C 14 , C 15 and they indeed fluctuate with zero mean. For the non-zero terms of the Voigt matrix (C 12 and C 13 are considered), fluctuations occur, with mean values being our C 12 and C 13 ; the dotted line shows the value C iso (equations (3.12) with C iso 12 = C iso 13 ), which appears to be a good estimate of the mean for low a 1 values. Next, we can get an estimate of 2 , the accuracy of our zeroth-order approximation. This is done by calculating numerically the with the averages defined in equations (2.4), (2.8) and (2.14). Results are reported in figure 12a; the error appears to be less than 0.5% for any degree of anisotropy a 1 and for both fabrics. We also report (figure 12b) the error iso , which defines the accuracy of the results if effective isotropy is assumed, (5.12) and the error is here of first order, as C ij − C iso ij does not vanish, except for a 1 = 1/3. As expected, the error is significantly higher and it increases with the degree of anisotropy a 1 up to about 10%.   .10) and (b) at first-order iso when considering affective isotropy. In red, for clusters and in green for girdles.

Concluding remarks
We have derived closed forms of the elasticity tensors and of the elastic wave velocities for ice textured polycrystals with cluster and girdle fabrics. The prediction is accurate up to the second order in the local deviation of the elasticity parameters from their average values. This is justified for weak anisotropic single crystals or for fabrics with concentrated c-axes, that is, in general a system with a low variability in the elasticity tensor of the grains. Motivated by measurements in deep ice cores, we have inspected this small variability case in more detail, revealing an accuracy better than 0.5%. This suggests that velocity measurements can confidently be used to characterize ice fabrics. In a forthcoming publication, we will address the specifics related to a more complete inspection of this problem. More generally, extensions of this study are at least twofold. (i) As few theoretical studies of the effects of texture on backscattering and attenuation have been conducted, it should be interesting to iterate the Dyson equation using the potential defined in equation (4.9) to quantify such effects. (ii) We have used simple ODFs which are based on the observation of the textures; other ODFs, as proposed in [27,41], use two or three tuneable parameters, being related to the fabric evolution under external stresses (due to the gravity and to possible ice flows). Being more involved, it does not allow for analytical calculations, but it should be used to quantify the sensitivity of the velocities on the form of the ODF.
Ethics statement.