Critical investigation of calculation methods for the elastic velocities in anisotropic ice polycrystals

Sonic logging measurements in ice core boreholes allow the determination of velocities of the elastic waves used as proxy for the variation of ice polycrystal anisotropy due to the texture (or fabric) evolution with depth. The idea beyond these measurements is to invert the elastic velocities to deduce the ice polycrystal anisotropy; this requires a model which links the elastic velocities to the ice texture. A classical model is based on the effective medium theory; the velocities are associated 5 to elastic waves propagating in a homogeneous medium characterized by an average elasticity tensor. Alternatively, the velocity averaging model uses the average of the velocities in single crystals with different orientations. In this paper, we show that the velocity average method is erroneous in the present context. This is demonstrated for waves propagating in a clustered polycrystal along the c-axis; two different 10 shear wave velocities are obtained while a unique velocity is expected. This unphysical result is not recovered by Bennett (1968) although the velocity averaging method is used; we show that this is due to erroneous expressions of the shear wave velocities in single crystals, used as the starting point in the average process. Because of the weak elastic anisotropy of ice, the inversion from the measured velocities requires the use of an accurate model, models based on effective medium theory 15 should be preferred.

Abstract.Crystallographic texture (or fabric) evolution with depth along ice cores can be evaluated using borehole sonic logging measurements.These measurements provide the velocities of elastic waves that depend on the ice polycrystal anisotropy, and they can further be related to the ice texture.To do so, elastic velocities need to be inverted from a modeling approach that relate elastic velocities to ice texture.So far, two different approaches can be found.A classical model is based on the effective medium theory; the velocities are derived from elastic wave propagation in a homogeneous medium characterized by an average elasticity tensor.Alternatively, a velocity averaging approach was used in the glaciology community that averages the velocities from a given population of single crystals with different orientations.
In this paper, we show that the velocity averaging method is erroneous in the present context.This is demonstrated for the case of waves propagating along the clustering direction of a highly textured polycrystal, characterized by crystallographic c axes oriented along a single maximum (cluster).In this case, two different shear wave velocities are obtained while a unique velocity is theoretically expected.While making use of this velocity averaging method, reference work by Bennett (1968) does not end with such an unphysical result.We show that this is due to the use of erroneous expressions for the shear wave velocities in a single crystal, as the starting point of the averaging process.
Because of the weak elastic anisotropy of ice single crystal, the inversion of the measured velocities requires accurate modeling approaches.We demonstrate here that the inver-sion method based on the effective medium theory provides physically based results and should therefore be favored.

Introduction
Wave propagation in glaciology is mostly regarded in the context of seismic waves; see, e.g., Kohnen (1974) and Blankenship et al. (1987).More recently, a new interest has emerged with the idea of using in situ wave velocity measurements along ice boreholes in order to evaluate the texture (fabric) anisotropy and its evolution with depth along deep ice cores.Such measurements are based on sonic logging that evaluates the travel time of elastic waves propagating in the ice over a short distance (typically few meters).A classical sonic logger was recently used during a measurement campaign at EPICA Dome C (East Antarctica) (Gusmeroli et al., 2012) and revealed the sensitivity of the elastic velocities on the degree of anisotropy of ice polycrystals characterized by a cluster-type texture (with c axis orientations clustered around the vertical direction) that varies with depth.The measured velocity changes are small, and the inversion procedure to access the texture estimation is dependent on the level of uncertainties in the measurements.Strong motivation therefore arises for the development of accurate modeling for the inversion of the elastic velocities into the local ice anisotropy.
In a previous paper (Maurel et al., 2015) we applied a classical model for wave propagation in polycrystals to some of the most specific textures that are found in ice recovered from deep ice cores, namely clusters with a vertical trans-Published by Copernicus Publications on behalf of the European Geosciences Union.
verse isotropy (VTI), as described before, and girdles with a horizontal transverse isotropy (HTI).The latter corresponds to crystallographic c axes mostly aligned in a vertical plan, which form in situations where a horizontal tension dominates the stress field.The model used relies on the definition of an "effective medium" characterized by an averaged elasticity tensor defined at the scale of many grains (the polycrystal).Maurel et al. (2015) presented a comparison with Bennett's predictions (Bennett, 1968), which are based on the velocity averaging method, and illustrated the weak differences in the velocities predicted by the two models for VTI clustered textures.Contrary to the effective medium theory, the velocity averaging method does not rely on a rigorous mathematical formalism.Thus, beyond the relative agreement observed for clustered textures, the correctness of the velocity averaging method can be questioned.This is the goal of the present paper.We will show that the velocity averaging method is based on an erroneous fundamental assumption.To do so, we perform a demonstration for the cases of cluster and girdle ice textures with VTI.It is found that a wave propagating along the vertical axis is associated with two different shear velocities, which is unphysical since the polarizations of the shear waves are in the plane of isotropy.In Bennett (1968) a unique expression of the shear velocity was obtained, starting with modified expressions of the two shear wave velocities in ice single crystals.The modification consisted in attributing symmetrical weights to both velocities in order to get the same value after the velocity averaging.It resulted in a velocity value close to the harmonic mean of the two unphysical shear velocities derived from the direct use of the velocity averaging method.
Considering the low elastic anisotropy of the ice single crystal, the accuracy of the inversion procedure from the measured velocities to ice texture is essential.In particular, a model relying on a rigorous mathematical formalism is required, in order to describe the case of complex textures, beyond the simple case of VTI clusters.With regard to this requirement, we will show that the model based on the effective medium theory appears more reliable.

Classical results on polycrystal effective anisotropy and wave propagation in anisotropic media
In this section, we recall classical results (i) on the propagation of elastic waves in a crystal or in a polycrystal characterized by a uniform elasticity tensor and (ii) on the elasticity tensor of polycrystals, considered at the scale of many grains as an equivalent "single" crystal.This allows us to introduce the notations that will be used further in the text and to clarify in a self-consistent way some properties that will be needed.

Wave propagation in a uniform anisotropic medium -the Christoffel equation
In the following, we denote x = (x 1 , x 2 , x 3 ) the spatial coordinates, u(x) = (u 1 (x), u 2 (x), u 3 (x)) the elastic displacement vector, ρ the constant mass density and c abcd the elasticity tensor being uniform in space.The case of a uniform c abcd corresponds either to a single crystal or to a polycrystal thought as an effective medium; in this latter case, c abcd characterizes the texture of the polycrystal.The propagation of monochromatic waves of frequency ω is described by the wave equation with a = 1, 2, 3 and where repeated indices mean summation (Einstein convention).Denoting k = kn the wave vector (k is the wave number) with n = (n 1 , n 2 , n 3 ) the unitary vector along k, the elastic displacement reads (2) The Christoffel equation admits in general three solutions, V = V P , V SH , V SV , that are the elastic velocities of the longitudinal wave and the two transverse waves (P, SH and SV stand for pressure wave, shear horizontal and shear vertical waves).It is important to stress at this point that because the three values of V 2 are the eigenvalues of the matrix c abcd n b n c /ρ, they do not depend on the particular frame (e 1 , e 2 , e 3 ) used to express c abcd .On the contrary, the eigenvectors (U 1 , U 2 , U 3 ) associated with the eigenvalues obviously depend on the frame where they are expressed.

Elasticity tensor of the polycrystal resulting from anisotropy from grain to grain
In the previous section, we have considered the uniform elasticity tensor of a polycrystal.This elasticity tensor is derived from the characteristics of the grains which compose the polycrystal.For ice, each grain is composed of the same single crystal, with hexagonal symmetry being characterized by its elasticity tensor c 0 ij kl , written in terms of C 0 I J in the Voigt's notation (Voigt, 1928).The superscript 0 refers to an elasticity tensor expressed in the frame of its principal axes, with the c axis being oriented along e 3 , hereafter referred as the vertical axis.The Voigt's matrix C 0 (with elements C 0 I J , I = 1, . . .6 and J = 1, . . .6) reads with the standard notations For an arbitrary direction of the c axis (Fig. 1) the elasticity tensor c abcd is deduced from c 0 ij kl following with R the rotation matrix (with elements R ij ): c abcd depends on (θ , ϕ), which are the usual angles in spherical coordinates.
The anisotropy at the macroscopic scale (at the scale of many grains) results from the many (or few) possible orientations of the c axis in each grain.This distribution of c axes is defined by a probability distribution function p(θ , ϕ) satisfying d p(θ, ϕ) = 1, with d = sin θ dθ dϕ. (8) The elasticity tensor c eff ij kl and the associated Voigt's matrix C eff I J of the polycrystal can be calculated by means of the average of the elasticity tensors of the grains following with the same index convention, Eq. ( 4), between the elasticity tensor c eff ij kl and the Voigt matrix C eff I J .For VTI textures, the direction of the axis of anisotropy is given by the effective c axis, denoted ĉeff .
For the numerical application, we will use the elastic constants as used in Bennett (1968): 3 Elastic wave velocities in the velocity averaging method and velocities of the effective medium for VTI textures In this section, we compare the elastic velocities obtained from the velocity averaging method, as used in Bennett (1968), Gusmeroli et al. (2012) and Vélez et al. (2016), and the elastic velocities obtained from the effective medium theory.The two approaches are as follows.
1. Velocities from the velocity averaging method: first, solve the Christoffel equation for given (θ , ϕ) using Then compute the average of the slowness, from which 2. Velocities of the effective medium: first, compute the effective elasticity tensor using Then, solve the Christoffel equation: Below, we report examples of VTI textures for which the velocity averaging method leads to inconsistent results.In such media and for a wave propagating along the vertical direction e 3 , the shear displacements are in the plane (e 1 , e 2 ), that is, in the plane where the response of the polycrystal is isotropic.Thus, there is a unique shear velocity in this particular case (V SV = V SH ).

Velocities from the velocity averaging method
In this method, we first derive the velocities in a grain, and this is done without lack of generality for a wave propagating along e 3 .The Christoffel equation, Eq. ( 2) with n b = δ b3 (same for n c ), simplifies to where the C I J is derived from c ij kl in Eq. ( 6), for a c axis given by Eq. ( 5).
We have used the notations cϕ ≡ cos ϕ, sϕ ≡ sin ϕ, sθ ≡ sin θ and cθ ≡ cos θ.The determinant, Eq. ( 14), can be calculated and we get the roots ρV 2 : The first root corresponds to a pure shear wave polarized in a direction perpendicular to both k and ĉ, referred as the SH wave.The two roots of the second equation are associated with the so-called quasi-shear and quasi-longitudinal waves, being coupled.The directions of polarization of the three waves are orthogonal (because the determinant is associated with a symmetric matrix) but the quasi-longitudinal wave is in general not along e 3 and the quasi-shear wave is in general not in the (e 1 , e 2 ) plane.More explicitly, the three velocities read with These are the expressions of the velocities in a single grain with a c axis forming an angle θ with k.
The second step in the velocity averaging method can be applied: for V = V SH , V SV , V P taken from Eq. ( 17).These last averages depend further on P θ 0 (θ ).

Velocities of the effective medium
With the probability function given by Eq. ( 13), the effective medium is characterized by a Voigt matrix: The Voigt matrix C I J (θ , ϕ) has 21 nonzero coefficients.Averaging C I J (θ , ϕ) over ϕ ∈ [0, 2π ] causes 15 of them to vanish, and the resulting Voigt matrix has a VTI symmetry, as expected; see Eqs. ( A2)-(A3) in Maurel et al. (2015).The remaining integrations over θ depend on P θ 0 (θ ).
Next, the velocities of the elastic waves propagating along the e 3 axis can be derived by solving the Christoffel equation, Eq. ( 2) with C I J → C eff I J .We get and we report below the intermediate results on C 33 and C 44 after ϕ averaging, from Eq. ( 15): Afterwards, The resulting effective velocities read The longitudinal wave has a velocity associated with the polarization along e 3 ; more importantly for the present demonstration, the two transverse waves have polarizations in the (e 1 , e 2 ) plane and they are associated with the same velocity.

Example 1: the VTI girdle with a single zenith angle θ 0
The first configuration is shown in Fig. 2. All the grains within the polycrystal have the same angle θ = θ 0 but different ϕ randomly distributed in [0, 2π ], where The velocities obtained from the velocity averaging method, Sect.3.1.1,are obtained with V in Eq. ( 17) and Eqs. ( 18) and ( 24), leading to  2a and b.The velocities V P and V S are reported as functions of θ 0 , from the effective medium theory (black curves) and from the velocity averaging method (plain red lines).The dotted red lines show the shear wave velocity from Bennett's expression, close to the harmonic mean of the two unphysical shear wave velocities in plain red line. with These velocities are reported in Fig. 4a (red curves) and they will be discussed later.We now derive the two shear and transverse velocities of the effective medium, Sect.3.1.2,which are given by Eq. ( 23) with C eff I J (θ 0 ) = dθ C I J ϕ (θ ) sin θ P θ 0 (θ ), and Eqs. ( 21) and ( 24).We get and these velocities are reported in black lines in Fig. 4a.

Example 2: the VTI clustered texture with opening angle θ 0
This texture is shown in Fig. 3.It corresponds to where H θ 0 is the rectangular function, equaling unity for 0 ≤ θ ≤ θ 0 and zero otherwise.The velocities (V av SH , V av SV , V av P ) in the velocity averaging method, Sect.3.1.1,cannot be calculated analytically in this case.The averages on θ in Eq. ( 18) are performed numerically (with Eqs. 17 and 27), and they are shown in red lines in Fig. 4b.The velocities of the effective medium (V eff S and V eff P ; Sect.3.1.2)are calculated as previously, using Eqs.( 21) to ( 23), with Eq. ( 27), and we get with X ≡ 1 + cos θ 0 + cos 2 θ 0 and Y ≡ cos 3 θ 0 + cos 4 θ 0 .These results are shown in black lines in Fig. 4b.Let us now discuss Fig. 4. Black lines correspond to the velocities V eff P and V eff S and as previously said, a unique velocity V eff S is found from the effective medium theory, by construction.Red lines correspond to the velocities V av P and (V av SH , V av SV ) from the velocity averaging method.The discrepancy is incidental for the longitudinal wave velocity and does not allow to discriminate between the two methods.The important result is that we find two distinct shear wave velocities using the velocity averaging method, which is unphysical in the present case.
4 Comment on Bennett's derivation of the velocities in the usual clustered texture Bennett (1968) presented a calculation of the elastic velocities in the VTI cluster (Fig. 3), and he obtained a unique expression of the shear velocity.This result is in contradiction with the calculations performed in the previous section and reported in plain red lines in Fig. 4. To understand this discrepancy, we analyze the derivation proposed in Bennett (1968).
An additional angle φ, not represented in Fig. 5a, between the planes (k, e 3 ) and (k, ĉ) is used.Next, an ensemble of relations between the different angles is derived, among which cos θ = sin σ cos ϕ sin θ + cos σ cos θ and These relations are correct; for instance if σ = 0, then θ = θ and φ = ϕ, as expected.Next, it is said that the slowness of the wave depends on φ following This is postulated by Bennett as being an intuitive approximation.Since φ can vary while θ remains constant (see Eq. 30), the above equations pretend that the velocity of the shear waves depends on more than just θ, which Eq. ( 29) shows as incorrect.More explicitly, using σ = 0 (thus, θ = θ, φ = ϕ) in Eq. ( 31) and using Eqs.( 29)-( 30), we get similarly to Eqs. ( 5)-(15) in Bennett (1968).These expressions for the velocities in a single crystal contain an unphysical dependence vs. ϕ and they are not in agreement with the original expressions, Eq. ( 29), they derive from.The modification performed in Eq. ( 31), when compared with the original expressions, consists of attributing symmetrical weights to both S 1 and S 2 in order to get an identical value for the two shear wave velocities, after velocity averaging.Bennett's results are reported in dotted red lines in Fig. 4. As expected, the shear wave velocity appears to be close to the harmonic mean of the two unphysical velocities (resulting from the velocity averaging calculation).The velocity of the P wave, being not modified by artificial weights, remains the same.
It is difficult to anticipate the consequences of such an approach in the case of other textures, since the weights used to deduce the shear wave velocity in Eq. ( 31) were chosen while anticipating a result expected after velocity averaging.The intuition that one can have for simple texture configuration may become hazardous when complex textures are dealt with.

Conclusions
In this paper, we have proposed a critical analysis of the velocity averaging method that has been used recently in the post treatment of the wave velocities deduced from borehole sonic logging measurements along ice cores.We have illustrated the error made in the fundamental assumption postulated at the basis of this method.This critical analysis is performed here in the case of simple VTI textures for which the error leads to an unphysical result and thus allows for a clear demonstration.For VTI textures, Bennett circumvented this problem by postulating weighted forms of the shear wave velocities (Bennett, 1968).The method is clever but it requires us to be able to anticipate good guesses for the selected weights, which may become difficult in the case of complex crystallographic textures.
The sonic logging measurements are sought as a proxy of the ice polycrystal anisotropy.Beyond the case of clustered VTI textures, complex textures could be characterized by such methods (or at least transitions between various textures with depth).On the one hand, the anisotropy of the ice www.the-cryosphere.net/10/3063/2016/The Cryosphere, 10, 3063-3070, 2016 single crystal is weak, and thus the anisotropy of ice polycrystals is even weaker.On the other hand, ice along boreholes is reasonably free of impurities when compared to most of the materials studied by similar sonic measurements.The recent sonic logging campaign performed at EPICA Dome C has revealed the feasibility of such measurements (Gusmeroli et al., 2012).Nowadays, commercial loggers measure precisely the velocity.In parallel, the characterization of the physical properties of the ice single crystal and of the dependences of these properties with pressure and temperature reaches precisions of the order of a few percents.This increase in the precision renders the inverse problem accessible, from the measured velocities to the ice texture, but the inversion model still requires a high level of accuracy.With respect to this requirement, a rigorous model as the model based on the effective medium theory appears as a good candidate and should be favored.
2, 3.This system of equations admits nonzero solutions for (U 1 , U 2 , U 3 ) if the determinant of the matrix ρω 2 δ ad − k 2 c abcd n b n c vanishes (δ ab is the Kronecker symbol).One gets the dispersion relation that links k and ω; equivalently, introducing the phase velocity V = ω/k, we get the classical form of the Christoffel equation Det ρV 2 δ ad − c abcd n b n c = 0.

Figure 2 .
Figure 2. First configuration of a VTI girdle with a single zenith angle θ 0 .(a) shows a typical c axis ĉ = (sin θ 0 cos ϕ, sin θ 0 sin ϕ, cos θ 0 ) within one grain, with a constant θ 0 , and ϕ varies randomly from grain to grain.(b) shows the resulting VTI texture of the polycrystal at the macroscopic scale.ĉeff is the effective c axis.

Figure 3 .
Figure 3. Second configuration of the usual VTI clustered texture.(a) shows a typical c axis ĉ = (sin θ cos ϕ, sin θ sin ϕ, cos θ ) within one grain, where θ varies in [0, θ 0 ] and ϕ varies randomly from grain to grain.(b) shows the resulting clustered, VTI, texture of the polycrystal at the macroscopic scale.ĉeff is the effective c axis.

Figure 4 .
Figure 4. Illustration of the inconsistency of the velocity averaging method for the two VTI textures shown in Fig.2a and b.The velocities V P and V S are reported as functions of θ 0 , from the effective medium theory (black curves) and from the velocity averaging method (plain red lines).The dotted red lines show the shear wave velocity from Bennett's expression, close to the harmonic mean of the two unphysical shear wave velocities in plain red line.

Figure 5 .
Figure 5. (a) System of angles used in Bennett's calculations.ĉ is given by (ϕ, θ ) and k is given by σ , being otherwise in the (e 1 , e 3 ) plane.The extra angle θ denotes the angle between k and ĉ, thus cos θ = sin σ cos ϕ sin θ + cos σ cos θ .(b) Particular case of the Bennett configuration, for σ = 0; in this case, θ = θ .