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

Agnès Maurel, Jean-François Mercier, Maurine Montagnat

► To cite this version:

Agnès Maurel, Jean-François Mercier, Maurine Montagnat. Critical investigation of calculation methods for the elastic velocities in anisotropic ice polycrystals. The Cryosphere, 2016, 10, pp.3063 - 3070. 10.5194/tc-10-3063-2016. hal-01418303

HAL Id: hal-01418303

https://inria.hal.science/hal-01418303

Submitted on 24 Mar 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

Agnès Maurel<fnref n="1" />, Jean-François Mercier<fnref n="2" />, and Maurine Montagnat<fnref n="3" />

*Correspondence to:* Agnès Maurel (agnes.maurel@espci.fr)

Received: 25 January 2016 – Published in The Cryosphere Discuss.: 9 February 2016

Revised: 20 May 2016 – Accepted: 9 June 2016 – Published: 16 December 2016

**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 inversion

method based on the effective medium theory provides physically based results and should therefore be favored.

## 1 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-

<footnote n="1">¹Institut Langevin, CNRS, ESPCI ParisTech, 1 rue Jussieu, 75005 Paris, France</footnote>

<footnote n="2">²Poems, CNRS, ENSTA ParisTech, INRIA, 828 boulevard des Maréchaux, 91762 Palaiseau, France</footnote>

<footnote n="3">³LGGE, CNRS, Université Grenoble Alpes, 38041 Grenoble, France.</footnote>
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.

## 2 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.

### 2.1 Wave propagation in a uniform anisotropic medium – the Christoffel equation

In the following, we denote $\mathbf{x} = (x_1, x_2, x_3)$ the spatial coordinates, $\mathbf{u}(\mathbf{x}) = (u_1(\mathbf{x}), u_2(\mathbf{x}), u_3(\mathbf{x}))$ the elastic displacement vector, $\rho$ 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 $\omega$ is described by the wave equation

$$ \rho\omega^2 u_a(\mathbf{x}) + c_{abcd} \frac{\partial^2}{\partial x_b \partial x_c} u_d(\mathbf{x}) = 0, \quad (1) $$

with $a = 1, 2, 3$ and where repeated indices mean summation (Einstein convention). Denoting $\mathbf{k} = \mathbf{k}\mathbf{n}$ the wave vector ($\mathbf{k}$ is the wave number) with $\mathbf{n} = (n_1, n_2, n_3)$ the unitary vector along $\mathbf{k}$, the elastic displacement reads $u_a(\mathbf{x}) = U_a e^{i\mathbf{k}\cdot\mathbf{x}}$ leading to $\rho\omega^2 U_a - k^2 c_{abcd} n_b n_c U_d = 0$ for $a = 1, 2, 3$. This system of equations admits nonzero solutions for $(U_1, U_2, U_3)$ if the determinant of the matrix $\rho\omega^2\delta_{ad} - k^2 c_{abcd} n_b n_c$ vanishes ($\delta_{ab}$ is the Kronecker symbol). One gets the dispersion relation that links $\mathbf{k}$ and $\omega$; equivalently, introducing the phase velocity $V = \omega/k$, we get the classical form of the Christoffel equation

$$ \mathrm{Det} \left[ \rho V^2 \delta_{ad} - c_{abcd} n_b n_c \right] = 0. \quad (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 / \rho$, 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.

### 2.2 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_{ijkl}$, written in terms of $C^0_{IJ}$ 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_{IJ}$, $I = 1, \dots, 6$ and $J = 1, \dots, 6$) reads
$$ \mathbf{C}^0 = \begin{pmatrix} A & A-2N & F & 0 & 0 & 0 \\ A-2N & A & F & 0 & 0 & 0 \\ F & F & C & 0 & 0 & 0 \\ 0 & 0 & 0 & L & 0 & 0 \\ 0 & 0 & 0 & 0 & L & 0 \\ 0 & 0 & 0 & 0 & 0 & N \end{pmatrix}, \quad (3) $$

with the standard notations

$$ \left\{ \begin{array}{l} c_{ijk\ell} \rightarrow \mathbf{C}_{IJ}, \\ \text{for } (i, j) \rightarrow I, (k, l) \rightarrow J, \\ \text{and } (1, 1) \rightarrow 1, (2, 2) \rightarrow 2, (3, 3) \rightarrow 3, \\ (3, 2), (2, 3) \rightarrow 4, (3, 1), (1, 3) \rightarrow 5, (1, 2), (2, 1) \rightarrow 6 \end{array} \right. . \quad (4) $$

For an arbitrary direction of the c axis (Fig. 1)

$$ \hat{\mathbf{c}} = (\sin\theta \cos\varphi, \sin\theta \sin\varphi, \cos\theta), \quad (5) $$

the elasticity tensor $c_{abcd}$ is deduced from $c_{ijkl}^0$ following

$$ c_{abcd} = R_{ia}R_{jb}R_{kc}R_{ld}c_{ijk\ell}^0, \quad (6) $$

with **R** the rotation matrix (with elements $R_{ij}$):

$$ \mathbf{R} \equiv \begin{pmatrix} \cos\theta \cos\varphi & \cos\theta \sin\varphi & -\sin\theta \\ -\sin\varphi & \cos\varphi & 0 \\ \sin\theta \cos\varphi & \sin\theta \sin\varphi & \cos\theta \end{pmatrix}. \quad (7) $$

$c_{abcd}$ depends on $(\theta, \varphi)$, 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(\theta, \varphi)$ satisfying

$$ \int d\Omega p(\theta, \varphi) = 1, \quad \text{with } d\Omega = \sin\theta d\theta d\varphi. \quad (8) $$

The elasticity tensor $c_{ijkl}^{\text{eff}}$ and the associated Voigt's matrix $\mathbf{C}_{IJ}^{\text{eff}}$ of the polycrystal can be calculated by means of the average of the elasticity tensors of the grains following

$$ c_{ijkl}^{\text{eff}} = \int d\Omega p(\varphi, \theta) c_{ijkl}, \quad \mathbf{C}_{IJ}^{\text{eff}} = \int d\Omega p(\varphi, \theta) \mathbf{C}_{IJ}, \quad (9) $$

with the same index convention, Eq. (4), between the elasticity tensor $c_{ijkl}^{\text{eff}}$ and the Voigt matrix $\mathbf{C}_{IJ}^{\text{eff}}$. For VTI textures, the direction of the axis of anisotropy is given by the effective c axis, denoted $\hat{\mathbf{c}}^{\text{eff}}$.

For the numerical application, we will use the elastic constants as used in Bennett (1968):

$$ \text{Ice single crystal} \quad \begin{cases} A = 14.06 \times 10^9 \, \text{Nm}^{-2}, & C = 15.24 \times 10^9 \, \text{Nm}^{-2}, \\ L = 3.06 \times 10^9 \, \text{Nm}^{-2}, & \\ N = 3.455 \times 10^9 \, \text{Nm}^{-2}, & F = 5.88 \times 10^9 \, \text{Nm}^{-2}, \\ \rho = 917 \, \text{kg m}^{-3} \end{cases}. \quad (10) $$

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Spherical angles $(\varphi, \theta)$ used for the orientation of the c axis, ($\hat{\mathbf{c}}$ is a unitary vector).</figcaption></figure>

### 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 $(\theta, \varphi)$ using

$$ V(\theta, \varphi) = V_P, V_{SV}, V_{SH} \text{ are the roots of} $$

$$ \text{Det} \left[ \rho V^2 \delta_{il} - c_{ijkl} n_j n_k \right]. $$

Then compute the average of the slowness, from which

$$ V^{\text{av}} = \left[ \int d\Omega p(\theta, \varphi) V^{-1}(\theta, \varphi) \right]^{-1}. \quad (11) $$

2. Velocities of the effective medium: first, compute the effective elasticity tensor using

$$ c_{ijkl}^{\text{eff}} = \int d\Omega p(\theta, \varphi) c_{ijkl}. $$

Then, solve the Christoffel equation:

$$ V = V_P^{\text{eff}}, V_{SV}^{\text{eff}}, V_{SH}^{\text{eff}} \text{ are the roots of} $$

$$ \text{Det} \left[ \rho V^2 \delta_{il} - c_{ijkl}^{\text{eff}} n_j n_k \right]. \quad (12) $$

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 $\mathbf{e}_3$, the shear displacements are in the plane $(\mathbf{e}_1, \mathbf{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}$).
<figure><img src="image_4.png" /><figcaption>**Figure 2.** First configuration of a VTI girdle with a single zenith angle $\theta_0$. (a) shows a typical c axis $\hat{c} = (\sin\theta_0 \cos\varphi, \sin\theta_0 \sin\varphi, \cos\theta_0)$ within one grain, with a constant $\theta_0$, and $\varphi$ varies randomly from grain to grain. (b) shows the resulting VTI texture of the polycrystal at the macroscopic scale. $\hat{c}^{\text{eff}}$ is the effective c axis.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Second configuration of the usual VTI clustered texture. (a) shows a typical c axis $\hat{c} = (\sin\theta \cos\varphi, \sin\theta \sin\varphi, \cos\theta)$ within one grain, where $\theta$ varies in $[0, \theta_0]$ and $\varphi$ varies randomly from grain to grain. (b) shows the resulting clustered, VTI, texture of the polycrystal at the macroscopic scale. $\hat{c}^{\text{eff}}$ is the effective c axis.</figcaption></figure>

Two examples of VTI structures will be presented and the propagation along the vertical direction is considered. The first structure is artificial, with a c axis associated with a unique $\theta = \theta_0$ value in $[0, \pi/2]$ (and $\varphi \in [0, 2\pi]$), but it provides explicit expressions of the velocities in both approaches. It could be seen as a girdle with VTI (Azuma and Goto-Azuma, 1996) with a single zenith angle $\theta_0$ (Fig. 2). The second example corresponds to a cone representative of clustered textures measured along ice cores (Gusmeroli et al., 2012; Diez and Eisen, 2015), with $\theta \in [0, \theta_0]$ (and $\varphi \in [0, 2\pi]$), and is studied numerically (Fig. 3).

## 3.1 Wave propagation along the vertical axis in a polycrystal with VTI

We consider a wave propagating along the vertical axis, thus $k = (0, 0, k)$, in a polycrystal with VTI (Figs. 2 and 3). The VTI structure provides that $\varphi \in [0, 2\pi]$ and is associated with a probability distribution independent of $\varphi$. Thus, the distribution of c axes is given by probability distribution functions of the form

$$p(\theta, \varphi) = \frac{P_{\theta_0}(\theta)}{2\pi}, \quad \text{with} \quad \int_{0}^{\pi/2} d\theta \sin\theta P_{\theta_0}(\theta) = 1. \qquad (13)$$

### 3.1.1 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 = \delta_{b3}$ (same for $n_c$), simplifies to

$$\left| \begin{array}{ccc} \rho V^2 - C_{55} & -C_{45} & -C_{35} \\ -C_{45} & \rho V^2 - C_{44} & -C_{34} \\ -C_{35} & -C_{34} & \rho V^2 - C_{33} \end{array} \right| = 0, \qquad (14)$$

where the $C_{ijl}$ is derived from $c_{ijkl}$ in Eq. (6), for a c axis given by Eq. (5).
$$
\left\{
\begin{aligned}
C_{33} &= As\theta^4 + 2(2L + F)s\theta^2 c\theta^2 + Cc\theta^4 \\
C_{44} &= (A + C - 2F)s\theta^2 c\theta^2 s\varphi^2 + L \left[s\varphi^2(c\theta^2 - s\theta^2)\right]^2 \\
&\quad + c\theta^2 c\varphi^2\right] + Ns\theta^2 c\varphi^2 \\
C_{55} &= (A + C - 2F)s\theta^2 c\theta^2 c\varphi^2 + L \left[c\varphi^2(c\theta^2 - s\theta^2)\right]^2 \\
&\quad + c\theta^2 s\varphi^2\right] + Ns\theta^2 s\varphi^2 \\
C_{34} &= -s\varphi s\theta c\theta \left[As\theta^2 - Cc\theta^2 + (2L + F)(1 - 2s\theta^2)\right] \\
C_{35} &= -c\varphi s\theta c\theta \left[As\theta^2 - Cc\theta^2 + (2L + F)(1 - 2s\theta^2)\right] \\
C_{45} &= s\varphi c\varphi s\theta^2 \left[(A + C - 2F)c\theta^2 + L(1 - 4c\theta^2) - N\right]
\end{aligned}
\right.
\quad (15)
$$

We have used the notations $c\varphi \equiv \cos\varphi$, $s\varphi \equiv \sin\varphi$, $s\theta \equiv \sin\theta$ and $c\theta \equiv \cos\theta$. The determinant, Eq. (14), can be calculated and we get the roots $\rho V^2$:

$$
\left\{
\begin{aligned}
\rho V^2 &= (Lc\theta^2 + Ns\theta^2), \\
(\rho V^2)^2 &= (As\theta^2 + Cc\theta^2 + L)\rho V^2 - F^2 s\theta^2 c\theta^2 \\
&\quad + ACs\theta^2 c\theta^2 + ALs\theta^4 + CLc\theta^4 \\
&\quad - 2FLs\theta^2 c\theta^2 = 0
\end{aligned}
\right.
\quad (16)
$$

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**<sub>3</sub> and the quasi-shear wave is in general not in the (e<sub>1</sub>, e<sub>2</sub>) plane. More explicitly, the three velocities read

$$
\left\{
\begin{aligned}
\rho V_{\text{SH}}^2 &= Lc\theta^2 + Ns\theta^2, \\
\rho V_{\text{SV}}^2 &= \frac{1}{2} \left[ C + L + (A-C)s\theta^2 - \sqrt{D} \right], \\
\rho V_{\text{P}}^2 &= \frac{1}{2} \left[ C + L + (A-C)s\theta^2 + \sqrt{D} \right],
\end{aligned}
\right.
\quad (17)
$$

with

$$
D \equiv [As\theta^2 - Cc\theta^2] [As\theta^2 - Cc\theta^2 + 2L(c\theta^2 - s\theta^2)] + 4s\theta^2 c\theta^2 (F^2 + 2FL) + L^2.
$$

These are the expressions of the velocities in a single grain with a c axis forming an angle $\theta$ with **k**.

The second step in the velocity averaging method can be applied:

$$
\begin{aligned}
\frac{1}{V^{\text{av}}(\theta_0)} &= \int_0^{2\pi} d\varphi \int_0^{\pi/2} d\theta \sin\theta p(\theta, \varphi) V^{-1}(\theta, \varphi) \\
&= \int_0^{\pi/2} d\theta \sin\theta P_{\theta_0}(\theta) V^{-1}(\theta),
\end{aligned}
\quad (18)
$$

for $V = V_{SH}$, $V_{SV}$, $V_{P}$ taken from Eq. (17). These last aver-
ages depend further on $P_{\theta_0}(\theta)$.

**3.1.2 Velocities of the effective medium**

With the probability function given by Eq. (13), the effective
medium is characterized by a Voigt matrix:

$$
\mathbf{C}_{IJ}^{\text{eff}}(\theta_0) = \int_{0}^{2\pi} \frac{d\varphi}{2\pi} \int_{0}^{\pi/2} d\theta \mathbf{C}_{IJ}(\theta, \varphi) \sin\theta P_{\theta_0}(\theta). \quad (19)
$$

The Voigt matrix $C_{IJ}(\theta, \varphi)$ has 21 nonzero coefficients. Averaging $C_{IJ}(\theta, \varphi)$ over $\varphi \in [0, 2\pi]$ 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 $\theta$ depend on $P_{\theta_0}(\theta)$.

Next, the velocities of the elastic waves propagating along
the **e**<sub>3</sub> axis can be derived by solving the Christoffel equation,
Eq. (2) with **C**<sub>**IJ**</sub> → **C**<sub>**IJ**</sub><sup>eff</sup>. We get

$$
\left|
\begin{array}{ccc}
\rho V^2 - C_{44}^{\text{eff}}(\theta_0) & 0 & 0 \\
0 & \rho V^2 - C_{44}^{\text{eff}}(\theta_0) & 0 \\
0 & 0 & \rho V^2 - C_{33}^{\text{eff}}(\theta_0)
\end{array}
\right| = 0, \quad (20)
$$

and we report below the intermediate results on $C_{33}$ and $C_{44}$ after $\varphi$ averaging, from Eq. (15):

$$
\begin{cases}
\langle C_{33} \rangle_{\varphi}(\theta) = As\theta^4 + 2(2L + F)s\theta^2 c\theta^2 + Cc\theta^4, \\
\langle C_{44} \rangle_{\varphi}(\theta) = Lc\theta^2 + Ns\theta^2.
\end{cases}
\quad (21)
$$

Afterwards,

$$
C_{IJ}^{\text{eff}}(\theta_0) = \int d\theta \langle C_{IJ} \rangle_{\varphi}(\theta) \sin\theta P_{\theta_0}(\theta). \quad (22)
$$

The resulting effective velocities read

$$
\begin{cases}
V_{\mathrm{P}}^{\mathrm{eff}} = \sqrt{C_{33}^{\mathrm{eff}}(\theta_0) / \rho}, \\
V_{\mathrm{S}}^{\mathrm{eff}} = \sqrt{C_{44}^{\mathrm{eff}}(\theta_0) / \rho}.
\end{cases}
\qquad (23)
$$

The longitudinal wave has a velocity associated with the po-
larization along **e**<sub>3</sub>; more importantly for the present demon-
stration, the two transverse waves have polarizations in the
(**e**<sub>1</sub>, **e**<sub>2</sub>) plane and they are associated with the same velocity.

**3.2 Examples of the velocity averaging method for polycrystals with VTI**

**3.2.1 Example 1: the VTI girdle with a single zenith angle θ₀**

The first configuration is shown in Fig. 2. All the grains
within the polycrystal have the same angle $\theta = \theta_0$ but dif-
ferent $\varphi$ randomly distributed in $[0, 2\pi]$, where

$$
P_{\theta_0}(\theta) = \frac{\delta(\theta - \theta_0)}{\sin\theta_0}. \qquad (24)
$$

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
<figure><img src="image_6.png" /><figcaption>**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 $\theta_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.</figcaption></figure>

$$ 
\left\{
\begin{aligned}
V_{SH}^{av} &= \sqrt{\frac{Lc\theta_0^2 + Ns\theta_0^2}{\rho}}, \\
V_{SV}^{av} &= \sqrt{\frac{C + L + (A - C)s\theta_0^2 - \sqrt{D}}{2\rho}}, \\
V_{P}^{av} &= \sqrt{\frac{C + L + (A - C)s\theta_0^2 + \sqrt{D}}{2\rho}},
\end{aligned}
\right.
\qquad (25) $$

with

$$ D = [As\theta_0^2 - Cc\theta_0^2][As\theta_0^2 - Cc\theta_0^2 + 2L(c\theta_0^2 - s\theta_0^2)] + 4s\theta_0^2c\theta_0^2(F^2 + 2FL) + L^2. $$

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_{IJ}^{eff}(\theta_0) = \int d\theta \langle C_{IJ} \rangle_{\varphi}(\theta) \sin\theta P_{\theta_0}(\theta)$, and Eqs. (21) and (24). We get

$$ 
\left\{
\begin{aligned}
V_{P}^{\text{eff}} &= \sqrt{\frac{1}{\rho} \left[ As\theta_0^4 + 2(2L + F)s\theta_0^2c\theta_0^2 + Cc\theta_0^4 \right]}, \\
V_{S}^{\text{eff}} &= \sqrt{\frac{1}{2\rho} \left[ (A + C - 2F)s\theta_0^2c\theta_0^2 + L(4s\theta_0^4 - 5s\theta_0^2 + 2) + Ns\theta_0^2 \right]}.
\end{aligned}
\right.
\qquad (26) $$

and these velocities are reported in black lines in Fig. 4a.

### 3.2.2 Example 2: the VTI clustered texture with opening angle $\theta_0$

This texture is shown in Fig. 3. It corresponds to

$$ P_{\theta_0}(\theta) = \frac{H_{\theta_0}(\theta)}{1 - \cos\theta_0}, \qquad (27) $$

where $H_{\theta_0}$ is the rectangular function, equating unity for $0 \le \theta \le \theta_0$ and zero otherwise. The velocities ($V_{SH}^{av}$, $V_{SV}^{av}$, $V_{P}^{av}$) in the velocity averaging method, Sect. 3.1.1, cannot

be calculated analytically in this case. The averages on $\theta$ 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_S^{eff}$ and $V_P^{eff}$; Sect. 3.1.2) are calculated as previously, using Eqs. (21) to (23), with Eq. (27), and we get

$$ \text{Cluster} \left\{
\begin{aligned}
V_S^{\text{eff}} &= \sqrt{\frac{(L+N)}{2\rho} + \frac{[2(A+C)-4F-3L-5N]}{30\rho} \times \frac{10\rho}{[(A+C)-2(2L+F)]Y}} \\
V_P^{\text{eff}} &= \sqrt{\frac{A}{\rho} + \frac{[-7A+3C+4(2L+F)]}{15\rho} \times \frac{10\rho}{[A+C-2(2L+F)]Y}}
\end{aligned}
\right. \qquad (28) $$

with $X \equiv 1 + \cos\theta_0 + \cos^2\theta_0$ and $Y \equiv \cos^3\theta_0 + \cos^4\theta_0$. These results are shown in black lines in Fig. 4b.

Let us now discuss Fig. 4. Black lines correspond to the velocities $V_P^{eff}$ and $V_S^{eff}$ and as previously said, a unique velocity $V_S^{eff}$ is found from the effective medium theory, by construction. Red lines correspond to the velocities $V_P^{av}$ and ($V_{SH}^{av}$, $V_{SV}^{av}$) 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).
## 4.1 Bennett's calculations in Bennett (1968)

Bennett starts with the slowness in a single crystal

$$S = 1/V,$$

given by

$$ \begin{cases} S_1 = S_P \approx a_1 - b_1 \cos 4\theta - c_1 \cos 2\theta, \\ S_2 = S_{SV} \approx a_2 + b_2 \cos 4\theta, \\ S_3 = S_{SH} \approx a_3 + b_3 \cos 2\theta, \end{cases} \quad (29) $$

where $(a_1, a_2, a_3, b_2, b_3)$ are constant for ice. The above expressions are close to the Thomsen approximations (Thomsen, 1986). They correspond to approximate expressions of the inverses of the velocities for the single crystals as given in Eq. (17). In the above expressions, $\theta$ is the angle $(k, \hat{c})$ (Fig. 5a). Bennett used a different frame where neither $k$ nor $\hat{c}$ are along the direction $e_3$. In this frame, he defined new multiple angles (Fig. 5b): $\sigma \equiv (k, e_3)$, $\theta = (\hat{c}, e_3)$ and $\hat{\theta} \equiv (k, \hat{c})$.

An additional angle $\hat{\phi}$, not represented in Fig. 5a, between the planes $(k, e_3)$ and $(k, \hat{c})$ is used. Next, an ensemble of relations between the different angles is derived, among which $\cos \hat{\theta} = \sin \sigma \cos \varphi \sin \theta + \cos \sigma \cos \theta$ and

$$ \begin{cases} \sin^2 \hat{\phi} = \frac{\sin^2 \varphi \sin^2 \theta}{\sin^2 \hat{\theta}}, \\ \cos^2 \hat{\phi} = \frac{(\cos \sigma \cos \varphi \sin \theta - \sin \sigma \cos \theta)^2}{\sin^2 \hat{\theta}}. \end{cases} \quad (30) $$

These relations are correct; for instance if $\sigma = 0$, then $\hat{\theta} = \theta$ and $\hat{\phi} = \varphi$, as expected. Next, it is said that the slowness of the wave depends on $\hat{\phi}$ following

$$ \begin{cases} S_P = S_1, \\ S_{SV} = S_2 \cos^2 \hat{\phi} + S_3 \sin^2 \hat{\phi}, \\ S_{SH} = S_2 \sin^2 \hat{\phi} + S_3 \cos^2 \hat{\phi}. \end{cases} \quad (31) $$

This is postulated by Bennett as being an intuitive approximation. Since $\hat{\phi}$ can vary while $\hat{\theta}$ remains constant (see Eq. 30), the above equations pretend that the velocity of the shear waves depends on more than just $\hat{\theta}$, which Eq. (29) shows as incorrect. More explicitly, using $\sigma = 0$ (thus, $\hat{\theta} = \theta$, $\hat{\phi} = \varphi$) in Eq. (31) and using Eqs. (29)–(30), we get

$$ \begin{cases} S_{SV} = (a_2 + b_2) - 8b_2 \cos^2 \varphi \cos^2 \theta \sin^2 \theta \\ -2b_3 \sin^2 \varphi \sin^2 \theta, \\ S_{SH} = (a_2 + b_2) - 8b_2 \sin^2 \varphi \cos^2 \theta \sin^2 \theta \\ -2b_3 \cos^2 \varphi \sin^2 \theta, \end{cases} \quad (32) $$

similarly to Eqs. (5)–(15) in Bennett (1968). These expressions for the velocities in a single crystal contain an unphysical dependence vs. $\varphi$ 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

<figure><img src="image_7.png" /><figcaption>**Figure 5.** (a) System of angles used in Bennett's calculations. $\hat{c}$ is given by $(\varphi, \theta)$ and $k$ is given by $\sigma$, being otherwise in the $(e_1, e_3)$ plane. The extra angle $\hat{\theta}$ denotes the angle between $k$ and $\hat{c}$, thus $\cos \hat{\theta} = \sin \sigma \cos \varphi \sin \theta + \cos \sigma \cos \theta$. (b) Particular case of the Bennett configuration, for $\sigma = 0$; in this case, $\hat{\theta} = \theta$.</figcaption></figure>

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.

## 5 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
single crystal is weak, and thus the anisotropy of ice poly-
crystals is even weaker. On the other hand, ice along bore-
holes 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 (Gus-
meroli 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 de-
pendences of these properties with pressure and temperature
reaches precisions of the order of a few percents. This in-
crease in the precision renders the inverse problem accessi-
ble, 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 can-
didate and should be favored.

*Acknowledgements.* We are thankful to A. Diez and an anonymous referee for helping to improve the paper. Maurine Montagnat is supported by the French CNRS (Centre National de la Recherche Scientifique) and the institutes INSIS (Institut des Sciences de l'Ingénierie et des Systemes) and INSU (Institut National des Sciences de l'Univers).

Edited by: O. Eisen

Reviewed by: A. Diez and one anonymous referee

References

Azuma, N. and Goto-Azuma, K.: An anisotropic flow law for ice-sheet ice and its implications, Ann. Glaciol., 23, 202–208, 1996.

Bennett, H. F.: An investigation into velocity anisotropy through measurements of ultrasonic wave velocities in snow and ice cores from Greenland and Antarctica, PhD thesis, University of Wisconsin, Madison, 1968.

Blankenship, D. D., Bentley, C. R., Rooney, S. T., and Alley, R. B.: Till beneath Ice Stream B: 1. Properties derived from seismic travel times, J. Geophys. Res.-Solid Ea., 92, 8903–8911, 1987.

Diez, A. and Eisen, O.: Seismic wave propagation in anisotropic ice – Part 1: Elasticity tensor and derived quantities from ice-core properties, The Cryosphere, 9, 367–384, doi:10.5194/tc-9-367-2015, 2015.

Gusmeroli, A., Pettit, E. C., Kennedy, J. H., and Ritz, C.: The crystal fabric of ice from full?waveform borehole sonic logging, J. Geophys. Res.-Ea. Surf., 117, F03021, doi:10.1029/2012JF002343, 2012.

Kohnen, H.: The temperature dependence of seismic waves in ice, J. Glaciol., 13, 144–147, 1974.

Maurel, A., Lund, F., and Montagnat, M.: Propagation of elastic waves through textured polycrystals: application to ice, P. Roy. Soc. Lond. A, 471, 2177, doi:10.1098/rspa.2014.0988, 2015.

Thomsen, G. W.: Weak elastic anisotropy, Geophysics, 15, 1954–1966, 1986.

Vélez, J. A., Tsoflias, G. P., Black, R. A., Van der Veen, C. J., and Anandakrishnan, S.: Distribution of preferred ice crystal orientation determined from seismic anisotropy: Evidence from Jakobshavn Isbræ and the North Greenland Eemian Ice Drilling facility, Greenland, Geophysics, 81, WA111–WA118, 2016.

Voigt, W.: Lehrbuch der Kristallphysik, reprint of the 1st Edn., Teubner, Leipzig, 1928.