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 propagation of sound in polycrystals has long been the subject of studies (for a review, see ). Polycrystals are inhomogeneous media, the inhomogeneity being in the size and/or in the crystallographic orientations of the grains. To describe the propagation of elastic waves through a large number of grains, averaged wave properties have to be determined and, until the 1980s, these properties were determined by considering a spatial average. This average was obtained from the inverse of the velocity following the intuitive argument that the time of flight is the correct ‘intensive’ variable to consider: the time of flight needed to travel along a given distance is equal to the sum of the times needed to go through subsequent grains on the wave path [2,3]. However, this intuitive approach neglects the mode conversion between longitudinal and transverse waves at the grain boundaries, and it does not guarantee that the obtained velocities respect the symmetries of the polycrystal. Alternatively, velocities have been derived using averaged elasticity parameters in the static limit, that is, interrogating the average of the stress to strain relation, rather than the wave equation (the wave velocity is a sub-result of the primary average). Thus, the square of the velocity (Voigt average) or the square of the inverse of the velocity (Reuss average) have been considered, and it is commonly accepted that these two average values bound the observed or measured values. Most of these studies considered the case of randomly oriented grains, with a uniform distribution of orientations, that is, an isotropic fabric in which preferred orientations were disregarded, and they used an average for the longitudinal wave and for the transverse waves independently, thus again neglecting mode conversion . The influence of microstructure on the macroscopic properties of random heterogeneous materials has been surveyed by Torquato .
Starting with Karal and Keller in the 1960s [6,7], stochastic methods have been developed in which the spatial averages are replaced by ensemble averages . These statistical approaches simplify the problem in a mathematically rigorous formalism and it is widely accepted that the assumption of ergodicity in the stochastic process ensures that the two averages, spatial average and ensemble average, are equivalent . The basic idea behind ergodicity is that, for a sufficiently long wave propagation distance, all the realizations of the averaged quantity (for instance, all possible crystallographic orientations) will be interrogated by the wave. Note that, if this is not the case, a spatial average is questionable as well. Indeed, in that case, the apparent velocity along a short propagation distance will be dependent on the particular realization, that is, dependent on the states (size, crystallographic orientations, etc.) of a small number of grains. Rather, spatial or ensemble averages have to refer to the most probable realization, to which any particular single realization is believed to be close. Nowadays, the statistical approach is the method used in most studies. An important advantage is that a perturbative method can be properly iterated to extract information on the effective medium, such as the attenuation and the backscattering, which appear as a second-order effect, that cannot be obtained by simple-minded averaging. The successive iterations are often referred to as the successive Born approximations [1,10,11] or the successive approximations of the Dyson equation [12–15].
Textured polycrystals are structures which have crystallized with a preferred orientation or which have been loaded by a non-uniform stress after formation. This latter case concerns notably ice polycrystals, as recovered in ice sheets along deep ice cores, which is the motivation of this study. We consider hexagonal single crystals, the predominant form of ice found on Earth. As usual, we denote by c the long axis of the unit cell and by a the short axis on the plane perpendicular to c. In the central part of deep ice sheets, ice is deformed at very low strain rates (10−10 to 10−12 s−1) induced by progressive accumulation of snow layers at the surface. In the simplest case (for instance, under perfect domes), ice is loaded by uniaxial compression due to gravity, and the fabric evolves from nearly isotropic at the surface towards a cluster fabric, where the c-axes are randomly distributed within a cone with vertical axis. The cone angle decreases with depth or ice age, until a strong single maximum close to the bedrock is reached [16,17]. At the scale of thousands of years, ice may also flow along flow lines that depend on a deformation field that is mainly influenced by the surface slope and the bedrock properties, and also because a geographical dome can shift in position . Under these conditions, ice fabrics result from a complex stress field, the c-axes rotating towards the axes of compressive stress and away from 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 . 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 .
Until recently, this fabric was mostly determined by recovering the individual grain orientations on an ice thin section (dimensions of about 10×10 cm2 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  (see also , ch. 6) and Bentley : the use of borehole sonic measurements to determine in situ depth-continuous fabric. This method has been tested at Dome C, East Antarctica , 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 , 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 . 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 δcijkl 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  and reported below 1.1where (A,C,L,N,F) are the elastic parameters of the elasticity tensor (see equation (3.3)).
2. 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 2.1where 〈.〉 denotes ensemble averages. That is, an average over all possible realizations of the directions of the c-axis . 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. ). 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 a1 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 a1=1/3 for a fully isotropic structure and a1=1 for a fabric with a single maximum; that is, when all the c-axes have the same direction) . 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 e3, 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 (e1 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 e3 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 2.2with φ∈[0,2π] and θ∈[0,θ0]. The probability distribution function p(φ,θ) is, as in [22,27], 2.3with H1(θ)=1, for 0≤θ≤θ0 and zero otherwise. For the spherical angles, the average of any quantity A is defined by 2.4
It follows, using equations (2.1)–(2.4), that the second orientation tensor a is given as (denoting (sθ, cθ) for (), respectively, (cφ,sφ) for ) 2.5with 2.6in agreement with [22,27]. We have ordered the eigenvalues a1>a2=a3, and a2=a3 is obtained because the directions along e1 and e2 are equivalent for the VTI symmetry. The cluster fabrics cover all the degrees of anisotropy, from full isotropy with a1=1/3 (θ0=π/2, corresponding to all possible orientations of the c-axes) to the maximum anisotropy, with a1=1 (θ0=0, corresponding to all the c-axes being aligned with e3, 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 2.7for θ∈[0,θ0] (the ± signs correspond to φ=π/2(+) and 3π/2(−)). The problem being two dimensional in (e2,e3), the projection onto (e1,e2) is not needed and the average is defined with θ being considered as a polar angle. For any quantity A, we now have 2.8and 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 2.9and H2(θ)=1 for θ in [0,θ0], zero otherwise. Note that we could use θ∈[−θ0,θ0] and φ=π/2 but the calculations shown below are easier using the form above. From equations (2.7) and (2.8), and summing the two contributions of the c-axes for φ=π/2 and 3π/2 (with projections on e2), a is diagonal, with 2.10and 2.11where we have defined , and sinc0=1. As we could have anticipated, the partial girdle is always anisotropic, with a1 varying between (perfect girdle, for θ0=π/2) and a1=1 (single maximum fabric for θ0=0).
(i) Thick girdle
The thick girdle is parametrized by two angles but, as defined, θ and φ are not adapted to deal with it in a simple way: indeed, as defined, θ∈[π/2−ξ0;π/2+ξ0] and φ∈[0;2π] describe a thick girdle with VTI symmetry (symmetry by rotation with respect to e3); the thick girdle has HTI symmetry, with a symmetry by rotation with respect to e1. To avoid additional calculations in the forthcoming calculations of the Voigt matrix (see §3), we use the following ‘trick’. (i) We denote ( the co-latitude and longitudinal angles in a frame , as previously (figure 3); (ii) the ODF is defined to produce a thick girdle with rotational invariance around , thus and ; and (iii) the characteristics of the thick girdle with rotational invariance around e1 are simply deduced by defining 2.12
Now, the c-axis is parametrized by and in equation (2.2) when expressed in terms of . When expressed in terms of (e1,e2,e3), it is 2.13with and . Averages are defined by 2.14with the simple ODF 2.15and for , and zero otherwise. The c-axis is given by equations (2.13), and, using equations (2.14) and (2.15), the a tensor is again diagonal 2.16with 2.17
As for the clustered fabrics, two of the eigenvalues are equal. This reflects the equivalence between two of the principal directions of the structure (e2 and e3 for the thick girdle). As expected, the thick girdles prolongate the range of anisotropy of the partial girdles, with a1 varying from a1=1/2 (perfect partial girdle for ξ0=0) to a1=1/3 (perfect isotropy for ξ0=π/2).
(c) Evolution of the anisotropy with the ai
Figure 4 shows the variation of the ai, 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 a1 allows us to characterize the degree of anisotropy of a polycrystalline structure in a unified way, that is, independently of the considered fabric.
3. Elasticity tensors of single crystals and textured polycrystals
In this section, we derive the elasticity tensor cijkl 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 , CIJ with the standard correspondences cijkℓ→CIJ, for (i,j)→I, (k,l)→J and (1,1)→1, (2,2)→2, (3,3)→3, (3,2),(2,3)→4, (3,1),(1,3)→5 and (1,2),(2,1)→6.
(a) The elasticity tensor of a single crystal
The elasticity tensor is derived for an arbitrary direction of the c-axis (equation (2.2)) in (e1,e2,e3) (figure 1). The crystal frame is defined by the principal directions of anisotropy, , being in the basal plane and . In the reference frame (e1,e2,e3), we have 3.1and we define the rotation matrix R 3.2Note 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 e1 is transformed into the radial vector, which is the first vector in spherical coordinates, but for convenience it is the third vector in the present case. This produces a c-axis oriented along the vertical direction e3 for φ=θ=0.
When expressed in the local frame , the single crystal with hexagonal symmetry is described by the elasticity tensor , written as C0 in Voigt's notation, 3.3
In a frame (e1,e2,e3), the elasticity tensor cabcd is deduced from using 3.4according to the transformation laws for tensors (see, for example,  for the calculation of anisotropic properties from texture data using an open source package). The corresponding relation between the Voigt matrices C0 and C is in general involved, and the use of the Bond matrix is often preferred . We do not use the Bond matrix formalism but, owing to the hexagonal symmetry in C0, we use a reasonably tractable expression, similar to the one used in , 3.5These terms are calculated and we obtain, for the diagonal terms, 3.6and the off-diagonal terms 3.7
(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 Cij in equations (3.6) and (3.7) are averaged according to 3.8with 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 3.9with the elastic constants 3.10with 3.11
(1) The single maximum fabrics for θ0=0 (X=3, Y =2) lead to 〈C11〉=A, 〈C33〉=C, 〈C44〉=L, 〈C66〉=N and 〈C13〉=F; the effective medium is, as expected, the same as a homogeneous medium with a unique orientation of the c-axis along e3, in agreement with the expression of the Voigt matrix C0 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 3.12and 〈C13〉iso=〈C11〉iso−2〈C44〉iso, resulting in effective Lamé coefficients 3.13These elastic constants in the isotropic case agree with previous derivations, referred to as the Markham derivation  (see also , eqn (4.b)).
(ii) Partial girdle fabrics
The average of the Voigt matrix is carried out according to the two-dimensional configuration 3.14with 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 3.15with 3.16where we have defined 3.17and (and sinc 0=1). As previously noted, partial girdles always have isotropic structure. For θ0=0 (S2=S4=1), we recover the hexagonal symmetry of a single crystal with the c-axis along e3 (equation (3.3)). Also, for θ0=π/2 (S2=S4=0), the effective medium has hexagonal symmetry with HTI corresponding to a perfect girdle, with 〈C23〉PG=〈C22〉PG−2〈C44〉PG and 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,
(3) The Voigt matrix Cij expressed in the material frame (e1,e2,e3) is rearranged according to correspondences , and , leading to 3.20
It is shown in appendix C that the resulting Voigt matrix is associated with a structure with HTI symmetry (with respect to e1), 3.21with 3.22where 3.23
(4) Thick girdles (with a1 between 1/3 and 1/2) have in general a lower degree of anisotropy than partial girdles (a1 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)).
4. 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 e3 (§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 e3 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 4.1In our case, the problem will be solved considering the wave propagating along the vertical e3 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 ( associated with the principal directions of anisotropy of the crystal) but where an arbitrary direction of the wave propagation is considered (see ; see also , 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 e3 correspond to ua=Ua eikx3, with k the wavenumber, and equation (4.1) simplifies to 4.2which admits a non-zero solution for (U1,U2,U3) 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 4.3with C, the Voigt matrix, associated with the elasticity tensor cabcd, 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 e3 and this has been done in . 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 (4.3) for various φ-values and the eigenvalues are found to be independent of φ, as expected). For simplicity, we report below the result for φ=0, where the coefficients Cii, i=3,4,5 and C34, C45, C35 needed in the dispersion relation (equation (4.3)), have the simpler form (from equations (3.6) and (3.7)) 4.4and C34=C45=0. The discriminant simplifies to 4.5The wavenumber k2=ρω2/C44 corresponds to a non-zero eigenmode (0,U2,0) associated with the transverse SH-wave perpendicular to the propagation plane (e3,c). The two other wavenumbers satisfy a second-order equation (vanishing second term in equation (4.5)), and they correspond to the longitudinal wave P and the transverse wave SV that are coupled when propagating. The associated velocities are 4.6in agreement with , eqns A(10) and A(11).
(ii) Approximate expressions: the case of an ice single crystal
As ice single crystals are only weakly anisotropic, approximate expressions for the velocities are often used. In 1986, Thomsen  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 4.7with 4.8being small parameters (β∼4%, γ∼7% and δ∼20%), and the reference velocities being and . These approximate velocities are also reported in figure 5 for comparison with the exact ones. The errors in Thomsen's approximations are indeed small: 0.5% for an SP-wave, 0.3% for SH-waves and 0.4% for P-waves.
(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: cabcd=〈cabcd〉+δcabcd, so the wave propagation in a given realization reads 4.9where δcabcd 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 〈cabcd〉 and δcabcd, 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 〈cabcd〉 and submitted to a ‘potential’ Vad=(∂/∂xb)δcabcd(∂/∂xc) 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 4.10and any particular realization of the propagation will be close to this average, up to ϵ2. Indeed, because 〈δcabcd〉=0 by construction, the potential Vad 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 〈cijkl〉 has been expressed in the previous section using a frame where the direction e1,e2 or e3 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=(k1,k2,k3) (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 ua=Ua ei(k1x1+k2x2+k3x3), the dispersion relation takes the form, 4.11where we have defined and Bijk≡[〈Cij〉+〈Ckk〉]. 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)).
(iii) Cluster fabric
The cluster fabric has hexagonal symmetry, being isotropic in the (e1,e2) plane. It is thus sufficient to define the angle Θ that k forms with e3. With , equation (4.11) simplifies to 4.12(with sΘ, cΘ denoting , ), and with the Voigt matrix defined in equation (3.10), with 〈C55〉=〈C44〉. We deduce the velocities (the plane of incidence being (e2,e3)) 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 4.14These expressions of the velocities in the isotropic case do not coincide with the leading order in Thomsen's approximation  (equations (4.7)), with β=δ=γ≃0 (leading to and ). This is because Thomsen defines the anisotropy of polycrystals with respect to the isotropy of a single crystal (a1=1), while the isotropy leading to the velocities above refers to an isotropic polycrystal with a1=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 a1 values. For a value of a1 between 1 and 1/3, the range of the Θ-dependent velocities decreases to a single value in the isotropic case. For a1=1, we recover the angular dependence of single crystals (figure 5), with VS∈[1831.8,2177.5] m s−1 and VP∈[3775.7,4044.5] m s−1. For the isotropic case, a1=1/3, we get VisoS=1956.1 m s−1 and VisoP=3847.7 m s−1.
(iv) 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 (e1,e2,e3), 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 V2=ω2/k2. This has been done numerically and the results for a1=0.8 and 0.5 (perfect girdle with symmetry by rotation around e1) are reported in figure 8.
For a thick girdle, the symmetry is hexagonal with isotropic behaviour in the (e2,e3) 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 k2=0, and (that is, for ), the dispersion relation 4.15
where we have used that 〈C33〉=〈C22〉, 〈C66〉=〈C55〉 and 〈C13〉=〈C12〉 (equation (3.21)). We deduce the velocities 4.16with D≡[〈C22〉sΘ2−〈C11〉cΘ2][〈C22〉sΘ2−〈C11〉cΘ2+2〈C55〉(cΘ2−sΘ2)]+4sΘ2cΘ2(〈C12〉2+2〈C12〉〈C55〉)+〈C55〉2.
5. 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. . 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 . 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 , 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  is the quantitative inspection that is proposed, by means of an inversion procedure and by means of a comparison between the in situ measurements of the velocities and values of the velocities deduced from the measures of c-axis orientations in thin sections of core samples. The logs as used in in situ measurements provide the time of flight between two receivers located on a vertical logger which can be used to obtain the S- and P-wave velocities. 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  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 e3 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 e3 axis. In this case, the expressions from equations (4.11) simplify to , and (to see this, take k1=k2=0 in equations (4.11); it leads to the determinant of a diagonal matrix). Here, we use VS1 and VS2 with VS1 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. ). For cluster fabrics, from equations (3.10) (or using Θ=0 in equations (4.16) and ), we get 5.1with, as previously, and .
Figure 9 shows the variation of the velocities as a function of the degree of anisotropy of the polycrystal, measured with a1 (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 a1=1 (single maximum) and a1=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(Vmax−Vmin)/(Vmax+Vmin) for each fabric 5.4These 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, 5.5and, in most of the studies on clusters (as used in ), the average is done on the colatitude θ only, that is, the angle between the c-axis and the direction e3 of the wave propagation, and not on the longitude φ. At first sight, this approach appears reasonable. However, there are two 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 e3), thus the two shear velocities are the same for waves propagating along the e3-axis. From appendix A, the velocities for grains with the same θ and different φ-values are equivalent to an HTI structure (with symmetry around e3), and, from equations (3.7) and (A 3), we get 5.6and 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 SV-wave velocities being given by equations (4.6), and this cannot be the case. Incidentally, note that not only the VS(θ) but also VP(θ) is false when using equations (5.5) (see VP in equations (4.6) and (5.6)).
(ii) Bennett's prediction
Alternatively to the slowness average, Bennett  provides a prediction for VP and VS based on acoustic measurements. The details of the method are given in , ch. 5 and in , 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 e3-axis, Bennett's velocities read 5.7where (X,Y) are the same as used in equations (3.11), and the ai and bi are constants given by Bennett's fits: a1=256.28, b1=5.92, c1=5.08, a3=531.40, b2=45.37 and b3=15.94 (in μs m−1). Bennett's fits are usually reported in the literature with a different form. For convenience, we have rewritten them in the form of equations (5.7): 5.8Inspecting 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 5.9and the results of (5.7) are obtained from (5.8) by 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 (ai,bi) 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 , 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 , 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, 〈cijkl〉 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, as previously commented. Thus, the zeroth order of the Dyson equation, which we have used by solving equation (4.10) instead of equation (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 Cij parameters (four terms are shown in figure 11), for a cluster fabric as a function of a1, that is, for decreasing cone angle θ0; each point in these plots corresponds to a value of Cij of one grain with an orientation randomly chosen in [0,θ0]. The dispersion of the Cij indicates how far the actual grains are from their effective medium. For all a1 values, the effective medium has hexagonal symmetry, thus vanishing 〈C14〉, 〈C15〉 and they indeed fluctuate with zero mean. For the non-zero terms of the Voigt matrix (C12 and C13 are considered), fluctuations occur, with mean values being our 〈C12〉 and 〈C13〉; the dotted line shows the value Ciso (equations (3.12) with ), which appears to be a good estimate of the mean for low a1 values. Next, we can get an estimate of ϵ2, the accuracy of our zeroth-order approximation. This is done by calculating numerically the error ϵ2 5.11with 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 a1 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.12and the error is here of first order, as does not vanish, except for a1=1/3. As expected, the error is significantly higher and it increases with the degree of anisotropy a1 up to about 10%.
6. 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 [41,27], 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.
In this study, we have not conducted any experiments on human beings or animals, nor used any collection of human or animal data.
This work does not have any experimental data.
F.L. acknowledges the support of FONDECYT grant no. 1130382. This work benefited from support of the French institutes INSIS and INSU of CNRS. M.M. is part of Labex(ANR10 LABX56). A.M. is grateful for the support of LABEX WIFI (Laboratory of Excellence within the French Programme ‘Investments for the Future’) under references ANR-10-LABX-24 and ANR-10-IDEX-0001-02 PSL*.
M.M. provided the knowledge and bibliography for the ice core context of fabric development. A.M. and F.L. conceived the mathematical models and performed the calculations. A.M., F.L. and M.M. interpreted the results and compared them with the previous literature.
Conflict of interests
We have no competing interests.
A.M. thanks Benoit Froelich for enlightening discussions.
Appendix A. Effective stiffness tensor for a cluster fabric
(a) The average over φ
This first average over φ corresponds to a c-axis rotating about the vertical axis e3 with a constant angle θ. To omit multiple notations, we note Cij(θ) is the resulting average (〈Cij〉 is the final result). With φ∈[0,2π], we thus have A 1Using 〈cφ2〉=〈sφ2〉=1/2, 〈cφ4〉=〈sφ4〉=3/8 and 〈sφ2cφ2〉=1/8 (the odd power vanishes on average), we get, from equations (3.6) and (3.7), the following form of the Voigt matrix: A 2with the non-zero coefficients for this medium with transversely isotropic symmetry A 3
(b) The average over θ
Next, we average the coefficients in equations (A 3) over θ. In agreement with equations (2.3) and (2.4), we use A 4with Cij(θ) defined in equations (A 1) and (A 3). We use 〈cθ2〉=X/3, 〈sθ2cθ2〉=(2X−3Y)/15 and 〈cθ4〉=(X+Y)/5 with and . We get A 5
Appendix B. Effective stiffness tensor for a partial girdle fabric
First, the elastic Voigt matrix is considered for φ=π/2,3π/2 (δ-functions in equation (2.9)) to make the configuration symmetric, so we consider B 1We get a Voigt matrix which has orthorhombic symmetry B 2with non-zero terms B 3We now average over θ∈[0;θ0], according to equation (2.9) (the problem in θ is two dimensional), thus B 4We use 〈cθ2〉=(1−sinc2θ0)/2 and 〈cθ4〉=(3+4sinc2θ0+sinc4θ0)/8, with (and sinc0=1). The resulting Voigt matrix is, with S2≡sinc2θ0, and S4≡sinc4θ0, B 5
(a) Comparison with 
In their paper , Nanthikesan and Shyam Sunder considered a configuration equivalent to the partial girdle. For some reason, they took the c-axes in the horizontal plane with opening angle ψ0, and we take their results with , and (N refers to Nanthikesan and Shyam Sunder's notations). Their results on the mean elastic constants are denoted by (eqn 15(b) in , not reported here for brevity) and are easily compared with our own. Owing to the correspondences between the coefficients (ψ0,α,β,b1,b2,b3) (defined after eqns (15b) in ) and our coefficients S2 and S4 B 6It is now sufficient to rearrange the Voigt matrix according to the changes of axes for a=1,2,3,4,5,6 (resp. b) corresponding to i=2,3,1,5,6,4 (resp. a). We get , , , , , , , and , in agreement with eqns (15b) in  and with our equation (3.16).
Appendix C. Effective stiffness tensor for a thick girdle fabric
We use the notations of figure 3, and in for the ODF in equation (2.15). The derivation of is the same as in appendix A, resulting in equations (A 2) and (A 3) C 1with VTI symmetry in this frame (). The next average is now done on , according to equations (2.14) and (2.15), with C 2with Cij(ξ) given by equations (B 3). We use , (thus, and ) leading to , and C 3and with . It is now sufficient to use the correspondences resulting from equation (2.12) to HTI symmetry with respect to e1 (see equations (3.21) and (3.22)), C 4
- Received December 20, 2014.
- Accepted March 18, 2015.
© 2015 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.