Anisotropic attenuation in rocks : theory , modelling and lab measurements

S U M M A R Y Anisotropic attenuation affects seismic observations and complicates their interpretations. Its accurate determination is, however, difficult and needs extensive measurements of wavefields in many directions. So far, the traveltime and amplitude decay of waves are usually measured along a sparse grid of propagation directions, and methods for inverting for anisotropic attenuation are not fully developed. In this paper, we present theory allowing a description and parametrization of general triclinic anisotropic attenuation. We focus on a correct recalculation of ray quantities usually measured in lab to phase quantities needed in the inversion. We develop and numerically test an iterative inversion scheme for determining the parameters of anisotropic attenuation. We present a lab facility that allows for measuring anisotropic attenuation using the P-wave ultrasonic sounding of spherical samples in 132 directions distributed regularly over the sphere. The applicability of the proposed inversion method and the performance of the experimental setup are exemplified by determining triclinic anisotropic attenuation of the serpentinite rock from Val Malenco, Northern Italy. The ray velocity and ray attenuation were measured on a spherical sample of the rock with diameter of 45.5 mm at the room temperature and under two pressure levels: 0.1 and 20 MPa. The measurements confirmed that anisotropic attenuation is remarkably sensitive to confining pressure. Since cracks are closing with increasing pressure, attenuation decreases. However, changes in pressure can also induce changes in the directional variation of attenuation and rotation of anisotropy axes. The obtained results for the serpentinite rock sample are unique because they represent the first accurately determined triclinic anisotropic attenuation from lab measurements.


I N T RO D U C T I O N
Many rocks of the lithosphere and uppermost mantle are seismically anisotropic and attenuating (Babuška & Cara 1991;Savage 1999;Romanowicz 2003;Fouch & Rondenay 2006).These two phenomena are closely related and jointly affect the seismic waves (Burton 2007;Carcione 2014).They cause that the propagation velocity of seismic waves and dissipation of energy are directionally dependent and the signal is dispersive.Anisotropy and dispersion complicates modelling of wave propagation as well as interpretation of seismic observations.
Anisotropic velocity and attenuation are efficiently described using a model of viscoelastic anisotropy.The viscoelastic parameters are complex-valued and frequency dependent (Auld 1973;Carcione 2014).Their real and imaginary parts describe elastic and attenuation anisotropy, respectively.The use of complex algebra allows for generalizing the theory developed for elastic anisotropy to viscoelastic anisotropy.The equations for waves in viscoelastic media are formally the same as in elastic media except for being complex.Implementing complex algebra into equations is mathematically straightforward, but still some care is needed for understanding properly the physical meaning of all complex-valued quantities in the equations.
Anisotropic attenuation of rocks and other materials has been modelled theoretically (Mukerji & Mavko 1994;Carcione 2000;Chapman 2003;Wenzlau et al. 2010) and also measured in laboratory (Johnston & Toksöz 1980;Kim et al. 1983;Kern et al. 1997;Stanchits et al. 2003;Zhu et al. 2007).So far, most experiments were based on measuring elastic properties and attenuation of rocks by ultrasonic waves propagating in a limited range of directions of cubic, block or cylindrical rock samples (Hosten et al. 1987;Kern et al. 1997;Stanchits et al. 2003).Therefore, they allowed for rough estimating anisotropic attenuation of high symmetry only.The elastic anisotropy has also been measured on spherical rock samples (Pros et al. 1998(Pros et al. , 2003;;Lokajíček et al. 2014;Svitek et al. 2014).Such arrangement is more convenient and robust because it allows for measuring rock properties in a rather high number of uniformly distributed directions, which is necessary for determining parameters of general anisotropy.On the other hand, the inversion of measurements on spherical samples is more involved.The standard schemes are based on the inversion of phase quantities, which are related to propagation of plane waves.However, in lab experiments on spherical samples, waves with curved wave fronts are generated and the ray quantities are measured and interpreted.Similarly, field experiments provide usually measurements of ray rather than phase quantities.
In this paper, we extend the approach originally developed for studying general elastic anisotropy on spherical samples by incorporating measurements of anisotropic attenuation.We extend theory of inversion for viscoelastic anisotropy presented by Vavryčuk (2015) to be applicable to measurements in lab or field experiments.Specifically, we focus on the following unsolved issues: (1) how to calculate ray attenuation from amplitudes of signals propagated in anisotropic rocks, (2) how to eliminate effects of elastic anisotropy in amplitudes of waves and (3) how to construct the complex phase velocity surface from the complex energy velocity surface, which is a necessary step before applying an inversion scheme for parameters of viscoelastic triclinic anisotropy.The accuracy and robustness of the developed approach is numerically tested.It is shown that the inversion is applicable to weak as well as strong anisotropy.Finally, we demonstrate the determination of parameters of anisotropic attenuation on measurements of the serpentinite rock from Val Malenco, Northern Italy studied by Kern et al. (1997).

A N I S O T RO P I C V I S C O E L A S T I C M E D I A
In this section, basic formulae for waves propagating in anisotropic viscoelastic media are shortly reviewed.In the formulae, real and imaginary parts of complex-valued quantities are denoted by superscripts R and I.The complex-conjugate quantity is denoted by an asterisk.The magnitude of complex-valued vector v is √ v • v.If any complex-valued vector has a real-valued direction, it is called homogeneous.If its direction is complex, it is called inhomogeneous.In formulae, the Einstein summation convention is used for repeated subscripts.

Viscloelastic parameters
A viscoelastic anisotropic medium is defined by density-normalized stiffness parameters a i jkl = c i jkl /ρ which are, in general, frequency-dependent and complex-valued.The real and imaginary parts of a ijkl , define elastic and viscous properties of the medium.The Christoffel tensor jk jk (p) = a ijkl p i p l , (2) is frequency-dependent and complex.Similarly as in elastic media, the eigenvalue of the Christoffel tensor and its derivative are expressed ( Červený 2001) where p is the slowness vector p = n/c, c is the complex phase velocity, n is the complex slowness direction and v is the complex energy velocity vector.The polarization vector g is calculated as the eigenvector of the Christoffel tensor and normalized as g • g = 1.A more common normalization in complex algebra g • g * = 1 is not used because it leads to inconsistencies between equations in elastic and viscoelastic media.Vectors v and p are related by the equation, expressing their polar reciprocity similarly as in elastic media (Helbig 1994).
The above equations are formally the same for elastic and viscoelastic media.They just differ in whether the wave quantities are real or complex.The wave quantities are real in elastic media but generally complex in viscoelastic media.

Phase and ray quantities
The phase quantities describe propagation of plane waves and are obtained by decomposing the complex slowness vector p: where V phase , A phase and D phase are the real phase velocity, phase attenuation and phase inhomogeneity.Vectors n || and n ⊥ are real, mutually perpendicular unit vectors, n || is normal to the wave front and n ⊥ lies in the wave front.The phase inhomogeneity D phase depends on the boundary conditions.If D phase is known, the phase velocity V phase and attenuation A phase are uniquely determined from the eikonal equation (3).If D phase is zero, slowness vector p is homogeneous, its direction n is real and coincides with n || .Velocity V phase and attenuation A phase are then called intrinsic.The ray quantities describe propagation of energy along a ray.They are obtained by decomposing the complex energy velocity vector where V ray , A ray and D ray are the real ray velocity, ray attenuation and ray inhomogeneity.Vectors N || and N ⊥ are real, mutually perpendicular unit vectors.The ray inhomogeneity D ray depends on the boundary conditions.For point sources situated in homogeneous media, D ray is zero (see Vavryčuk 2007a).In this case, energy velocity vector v is homogeneous, its direction N is real and coincides with N || (see appendix F of Carcione & Ursin 2016).Velocity V ray and attenuation A ray are then called intrinsic, similarly as for the phase quantities.
The intrinsic phase and ray quantities are related to the complex phase velocity c and complex energy velocity v as follows: hence where The intrinsic complex phase velocity c is calculated straightforwardly using the eikonal equation ( 3), but calculating the intrinsic complex energy velocity v using eq.( 4) is more involved.It requires computing the so-called 'stationary' slowness vector p 0 introduced in ray theory (Vavryčuk 2007a).The stationary slowness vector is generally inhomogeneous being uniquely constrained by the boundary conditions.It can be calculated by complex ray tracing (Kravtsov et al. 1999;Hanyga & Seredynska 2000;Vavryčuk 2008aVavryčuk , 2010Vavryčuk , 2012) ) and corresponds to a ray connecting the source of waves with a receiver.While the source and receiver are points in real space, the ray is a curve in complex space.The problem is simplified in homogeneous media, where the ray becomes a straight line in real space.To determine the stationary slowness vector in homogeneous media, either a system of polynomial equations for unknown components of p 0 or an inverse problem for p 0 must be solved for each ray coming out of the source (for details, see Vavryčuk 2007a,b).

Asymptotic Green's function
Radiation of waves generated by a point source is usually calculated using the Green's function.The asymptotic Green's function in homogeneous anisotropic viscoelastic media reads (Vavryčuk 2007a, his eq. 18): or where and p 0 is the stationary slowness vector (see Vavryčuk 2007a,b).Quantity K = K 1 K 2 is the Gaussian curvature of the slowness surface, K 1 and K 2 are the principal curvatures and ϕ 1 and ϕ 2 are their phase angles.All quantities dependent on p in eqs ( 13) and ( 14) are taken at stationary point p 0 .Position vector x = r N, distance r, ray vector N, frequency ω, phase angles ϕ 1 and ϕ 2 and density ρ are real-valued; but polarization vector g, Gaussian curvature K, principal curvatures K 1 and K 2 , energy velocity v and slowness vector p 0 are complex-valued.The most efficient way how to calculate the Gaussian curvature K in eqs ( 13) and ( 14) is by computing the determinant of the wave metric tensor (see Vavryčuk 2003, his eqs 10 and 15).

Phase and ray measurements
The parameters of anisotropic attenuating media can be determined from a directionally dependent phase or ray velocity and attenuation.
The phase quantities are measured if the source of waves has a finite dimension and can generate a planar wave front.If the source is much smaller than the distance between the source and the receiver (the so-called point-like source), the wave front is curved and the ray quantities are measured.Since inverting phase quantities is computationally much simpler than inverting ray quantities, measurements of the phase velocity and attenuations are preferable.From the phase quantities, we can directly construct the complex phase velocity and then invert it for parameters of viscoelastic anisotropy.However, devising an experiment with planar wave fronts excited and propagating in a medium in many directions is complicated.Instead, we often measure just ray quantities on wave fronts generated by a point-like source.In this case, the problem is more involved and we have to first recalculate the ray quantities to the phase quantities before inverting them for parameters of viscoelastic anisotropy.

Ray velocity and attenuation
In lab or field experiments, propagation time of a signal along a ray and its amplitude are usually measured.Determining ray velocity from propagation time is straightforward but determining ray attenuation from amplitudes of the signal is more involved.The attenuation must be calculated from amplitude decay of a signal obtained as the difference between amplitudes measured at two points of a ray and corrected for geometrical spreading.
If the amplitude decay is measured along the source-receiver distance, the amplitudes must be corrected for the radiation pattern of a source (see eq. 13) Assuming the emitter and receiver of waves oriented along ray direction N, the radiation function takes the following form: The radiation function R in eq. ( 16) is complex and directionally dependent in anisotropic viscoelastic media.The imaginary part of R introduces phase shifts of the signal but the amplitude of the signal is influenced negligibly.If we focus on processing of amplitudes only, the radiation function R for attenuation anisotropy can be well approximated by that in elastic anisotropy.In this way, we can separate the effects of elastic anisotropy and attenuation and determine ray attenuation in the following way: (1) We determine elastic anisotropy of a rock from measurements of the ray velocity (Svitek et al. 2014).
(2) We calculate the radiation pattern and geometrical spreading of elastic waves propagating in the anisotropic rock.
(3) We correct the measured amplitude decay of the signal for these two factors and normalize the corrected value to a unit ray length.

Complex phase velocity
Having calculated ray velocity V ray and ray attenuation A ray , we determine the complex energy velocity v for a set of ray directions N using eq.( 8) and construct the complex energy velocity surface v = v(N).Decomposing the energy velocity vector v and the slowness vector p into their real and imaginary parts, the polar reciprocity relation ( 5 Taking into account that vector v is homogeneous and vector p is inhomogeneous where N is the real ray direction, and s R and s I are real directions of vectors p R and p I , respectively, eqs ( 17) and ( 18) imply and consequently Finally, complex slowness vector p is and its direction The above equations can be used for calculating the complex slowness surface p = p(n) from the complex energy velocity surface v = v(N).First, we calculate real unit vectors s R and s I as the normals to the real part v R = v R (N) and the imaginary part are real, we can use standard formulae of differential geometry (Lipschutz 1969).Then, we calculate p R and p I using eqs ( 22) and (23), and subsequently the complex slowness vectors p and its direction n using eqs ( 24) and ( 25).
If the procedure is applied to a sufficiently dense set of ray directions N, we can construct the whole slowness surface p = p(n) and subsequently the complex phase velocity surface c = c(n): and invert for parameters of viscoelastic anisotropy.

Inversion scheme
Determination of complex viscoelastic parameters a i jkl from the complex phase velocity surface c = c(n) is a non-linear inverse problem which can be solved using perturbation theory and iterations (Mensch & Rasolofosaon 1997;Vavryčuk 2008bVavryčuk , 2015;;Svitek et al. 2014).In perturbation theory, we assume that the anisotropic medium defined by unknown parameters a i jkl can be obtained by a small perturbation of a known reference medium where a 0 i jkl defines the viscoelastic reference medium and a i jkl its viscoelastic perturbation.Under this assumption, the cubic equation for the square of the phase velocity c 2 can be linearized as follows (Pšenčík & Vavryčuk 2002;Svitek et al. 2014): where c 0 and g 0 define the complex phase velocity and complex polarization vector in the reference medium, c 2 is the misfit between the squared phase velocity calculated from measurements and the velocity in the reference medium.The reference medium can be isotropic or anisotropic, elastic or viscoelastic.From perturbations a i jkl , we obtain parameters a i jkl using eq.( 28).This new medium serves as the reference medium in the next iteration.The iterations are repeated until perturbations a i jkl are negligibly small.The number of iterations depends on anisotropy strength: the stronger the viscoelastic anisotropy is, the more iterations are needed.For weak anisotropy (WA), the first iteration is usually sufficient for getting accurate results.

Weak-anisotropy-attenuation parameters
Determination of the complete set of 21 viscoelastic parameters requires measurements of the velocity and attenuation for the P and S waves.In the case of P-wave mesurements only, six parameters related to the S-wave propagation (called the S-wave related parameters) a 44 , a 55 , a 66 , a 45 , a 46 and a 56 cannot be well resolved (see Svitek et al. 2014).Under weak velocity and attenuation anisotropy, the S-wave related parameters cannot be determined at all.For this reason, in analogy to WA parameters (see e.g.Pšenčík & Gajewski 1998;Farra & Pšenčík 2003, 2008) we define the so-called weak anisotropy-attenuation (WAA) parameters and modify the inversion scheme to invert for 15 WAA parameters only (see Vavryčuk 2009).Since the reference medium is isotropic, eq. ( 28) is modified for the P wave as follows: and further rewritten in the form where the complex WAA parameters are defined as and Quantities α 0 and Q 0 are the mean real P-wave velocity and P-wave Q-factor in the reference medium.
Knowing the P-wave complex phase velocity c for a set of complex slowness directions n, eq. ( 30) can be solved for unknown complex WAA parameters using the generalized inversion (Menke 1989).The robustness and accuracy of this inversion will be demonstrated on synthetic data as well as lab measurements in the next sections.

N U M E R I C A L M O D E L L I N G
In this section, the inversion for viscoelastic anisotropy is tested on an example of the P wave propagating in a model of an orthorhombic viscoelastic medium (see Tables 1 and 2).The velocity anisotropy corresponds to the Torre Alfina xenolith with parameters taken from Pera et al. (2003, their table 3).The attenuation anisotropy is not based on measurements but it is synthetic.The mean P-wave phase velocity of the model is 7.93 km s −1 and the P-wave velocity anisotropy is 12.4 per cent.The mean P-wave phase Q-factor is 87.7, and the P-wave Q-factor anisotropy is 84.0 per cent.The density of the Torre Alfina xenolith is 3.31 g cm −3 .
The intrinsic phase and ray velocities, attenuations and Q-factors (see Fig. 1) for the synthetic model are calculated using equations from Sections 2.1 and 2.2.The quantities are evaluated in a regular grid of directions with the incidence ranging from −90 • to 90 • and with azimuth from 0 • to 180 • .The grid step is alternatively 5 • , 10 • , 15 • , 20 • and 30 • .Hence, in total of 1260, 306, 132, 72 and 30 independent measurements describe the directional variation of the velocity and attenuation.The ray velocity and attenuation values are used for calculating the complex phase velocity (see Section 3.3) which is inverted for the WAA parameters of the medium.The inversion is performed with noise-free and noisy data.Noise is random with a uniform distribution.The noise level is different for the real and imaginary parts of the complex phase velocity c reflecting that the propagation velocity is usually measured with much higher accuracy than attenuation.The noise level for the real part of c varies from 0 to ±3 per cent with step of 0.2 per cent, and the level for the imaginary part of c varies from 0 to ±15 per cent with step of 1 per cent.To get statistically robust results, the random noise is generated 100 times and the results of the inversion are averaged over all noise realizations.Since the convergence of the iterative process is quite fast (see Fig. 2), we performed the inversion with six iterations only.
The inversion is run for the velocity and attenuation measured along phase normals (Fig. 3) or along rays (Fig. 4).Fig. 3 shows the maximum errors of the predicted phase velocity V phase and phase attenuation A phase for various grid steps.The maximum value is calculated over all directions and it is further averaged over noise realizations.As expected, the higher the number of input measurements, Table 1.Elastic and attenuation parameters of the numerical model.Notes: The two-index Voigt notation A R K L and Q K L , K, L = 1, . . ., 6 is used for density-normalized elastic parameters a R i jkl and quality-factor parameters q i jkl = −a I i jkl /a R i jkl (no summation over repeated indexes).The elastic and attenuation parameters not shown in the table are identically zero.    1.

Elastic parameters
the lower the sensitivity of the results to noise.The attenuation errors are significantly higher than the velocity errors because they comprise errors produced by eliminating effects of the inaccurately determined elastic Green's function.Fig. 4 shows similar results as Fig. 3 but for the inversion of ray quantities.The attenuation errors in the inversion of noise-free ray data are due to the recalculation of ray to phase quantities before the inversion of the complex phase velocity c.In the case of a denser grid of ray measurements, this error decreases.For an increasing noise level, the role of the errors due to the recalculation of ray to phase quantities is suppressed.The upper and lower panels in Fig. 5 display the errors in the WAA parameters inverted using the phase and ray quantities for step of 15 • in the grid of directions.The numbers identifying the WAA parameters are defined in Table 4 (first two columns).Since the velocity anisotropy is much weaker than the attenuation anisotropy, the real parts of the WAA parameters are more sensitive to noise and thus less accurate than their imaginary parts.Nevertheless, the accuracy of predicted quantities is satisfactory even when inverting data with the highest noise level (Fig. 6).Again, the accuracy of the results is worse for the inversion of ray quantities (Figs 1 and  6) because of the necessity to recalculate the ray to phase quantities before the inversion.

Setup of the experiment
The determination of viscoelastic anisotropy is exemplified using an experimental facility that allows a measurement of the P-wave ray velocities and P-wave amplitudes on spherical rock samples.A diameter of the analysed rock sample is 45.5 mm but the equipment is capable of measuring also spherical samples of other diameters using specially designed arms.These arms enable insertion of a setting ring to adapt a different diameter of individual spherical samples (see Fig. 7).The sample was measured at room temperature being exposed to confining pressures of 0.1 and 20 MPa.Two piezoceramic sensors had a plain contact surface that resulted in a point contact between the sensor and the sphere.The ultrasonic signals were excited and recorded by a piezoceramics sensor pair with a resonant frequency of 2 MHz.The equipment allows for an ultrasonic sounding of rock samples in 132 independent directions with a pair of the P-wave sensors where both transmitter and receiver are polarized along the radial direction.The waveforms of ultrasonic signals were recorded using an A/D convertor with the dynamic range of 8 bits and with the sampling frequency of 100 MHz.Each waveform was recorded 10 times and then averaged in order to reduce noise.The sensitivity of the convertor was set up to record individual signals without any distortion and with the maximum possible dynamics.During the measurement, the point contact between the spherical sample and the transducers was established by a metal spring that acted by constant force at every direction.Transform oil was used as pressure medium.

Rock sample
We measured a sample of serpentinite from Val Malenco, Northern Italy, which comes from the block of serpentinite with the Kiel sample number 987 (Kern et al. 2015).At ambient conditions, the bulk density of the serpentinite derived from the volume and the mass of the sample sphere is 2.71 g cm −3 .The constituent minerals are 83 vol.% antigorite, 13 vol.%olivine, 2 vol.% magnetite and 2 vol.% chromite.The average grain size is about 200 μm.The porosity is rather low; the volume compaction of 1.63 per cent represents mainly the crack porosity.The sample displays a pronounced lattice preferred orientation of antigorite, that has been documented by 13.2 and 13.6 m.r.d.maxima for (001) pole figures in Kern et al. (2015), values that are consistent with other foliated serpentinites.To protect the pore space against the oil, the spherical surface of the sample was covered by a thin layer of an epoxy resin.
The P-and S-wave velocity and attenuation anisotropy of this rock has also been studied by Kern et al. (1997).In that study, the authors used a cubic sample with axes of the cube referring to macroscopic fabric coordinates: x parallel to the lineation, y normal to the lineation within the foliation plane and z normal to the foliation plane.The velocity and attenuation has been measured as a function of pressure up to 600 MPa and temperature up to 600 • C. The most prominent changes have been detected for pressure ranging from 0 to 100 MPa related to the closure of grain boundary cracks.While the mean velocity slightly increased, the mean attenuation remarkably decreased with increasing pressure (see Kern et al. 1997, their fig. 9a).Obviously, this dependence must be visible also in the measurements performed on the spherical sample in this paper.

Data processing
The ultrasonic measurements provide waveforms in 132 independent directions (see Fig. 8).The waveforms consist of a direct P wave followed by scattered coda waves produced by multiple reflections and conversions at the surface of the sample.As a consequence, the complete waveforms are complex and difficult to interpret (Fig. 8a).Nevertheless, the analysis of the direct P wave is rather simple because the secondary reflections arrive at times sufficient for separation of the direct P wave.The P-wave onset is usually noise-free and well distinguished and its determination can be performed using an automatic picking algorithm (Allen 1982;Sedlák et al. 2009;Svitek et al. 2010).The amplitude of the P wave is defined as the amplitude of the first maximum (or minimum) after the P-wave arrival and its automatic measurement is reliable with an error less than few samples (see Fig. 9).
Since the piezoceramics pair of sensors used for emitting and receiving the ultrasonic waves are planar, they touch the spherical sample at an area of about 2 × 2 mm or less.Since the size of the contact is much smaller than the radius of the sample, the contact between the sensors and the rock sample is considered to be pointlike.Consequently, the wave front of excited waves is not planar but curved and the velocity and amplitude of the signals are not measured along a phase normal but along a ray.Completing measurements of the traveltime and amplitude of the ultrasonic signal propagating through the sample in all directions, we construct a directional dependence of the ray velocity and ray amplitude for both pressure levels (see Fig. 10).As expected, the average velocity is slightly higher for confining pressure of 20 MPa because of a closure of cracks (see Table 3).This effect is even more pronounced in amplitudes indicating distinctly lower attenuation at this pressure level.The directional patterns of velocities and amplitudes are, however, similar for both pressure levels.Only, the amplitude variations seem to be slightly rotated relative to each other.

Results
The determination of parameters of viscoelastic anisotropy from measured velocities and amplitudes is performed in the following steps.First, the elastic anisotropy must be calculated from the directional variation of ray velocities.The procedure is quite analogous to that described in Section 3.4, but all quantities are real in the inversion.Second, the radiation pattern of a point source in elastic anisotropy is calculated (see Section 3.2).Third, the directionally dependent amplitudes of the signal emitted by the source and the amplitudes recorded at the receiver (Figs 8 and 9) are used for calculating the corrected amplitude decay along a ray and the ray attenuation (see Section 3.2).Fourth, an additional effect of a thin layer of an epoxy resin on the surface of the sample must be taken into account.The layer of the epoxy resin causes a uniform diminution of amplitudes in all directions and thus affects the scale of the determined ray attenuation and Q-factor.Therefore, the scale was adjusted by measurements of Kern et al. (1997) carried out for the same rock in the directions of the symmetry axes.The resultant ray attenuation and Q-factor patterns (see Fig. 11 and Table 3)  confirm a remarkable decrease of attenuation for the pressure of 20 MPa compared to the atmospheric pressure.However, the attenuation anisotropy and the Q-factor anisotropy increase.The increase of anisotropy strength indicates that the effects of aligned structure elements and inhomogeneities are more pronounced under confining pressure.Finally, the phase velocity and attenuation were calculated from the ray quantities according to Section 3.3 (see Fig. 12) and inverted for the parameters of viscoelastic anisotropy (see Table 4).
The predicted phase velocity, attenuation and Q-factor are displayed in Fig. 13.The predicted quantities are very similar to the input phase data indicating that the inversion was successful.The basic pattern of the predicted quantities at two pressure levels is also similar with some minor differences in attenuation and Q-factor.In particular, we observe a clockwise rotation of attenuation and Q-factors patterns by about 30 • .Since the velocity pattern is unchanged for both pressure levels, the difference in attenuation  Notes: For the definition and meaning of the quantities, see Table 2.
and Q-factors could not be produced by a misorientation of the sample in the two measurements, but it should be a real stress-induced effect.This indicates that anelastic properties of rocks are more sensitive to confining pressure than their elastic properties.

D I S C U S S I O N A N D C O N C L U S I O N S
An accurate and correct determination of anisotropic attenuation of rocks in laboratory is not an easy task.It needs an experimental facility that allows for measurement of the traveltime and amplitude decay of waves propagating along a dense grid of directions.The determination of viscoelastic parameters is simpler if we use measurements of plane waves than of waves generated by a point source.For plane waves propagating through the rock sample, we measure the traveltime and amplitude decay along the phase normal.Subsequently, we calculate the real phase velocity and phase attenuation, construct the complex phase velocity and invert for viscoelastic anisotropy.
For waves generated by point sources, we measure the traveltime and amplitudes along a ray.In order to calculate the amplitude decay and subsequently the ray attenuation, we have to measure wave amplitudes at two points along a ray.This can be realized by measuring amplitudes for two samples of a different size or by measuring amplitudes just for one sample provided we know the amplitude of the wave emitted at the source.In this case, however, the amplitude at the source must be corrected for geometrical spreading and for the radiation pattern of a point source in anisotropic media.Knowing the ray velocity and ray attenuation, we can construct the complex energy velocity.This velocity is further recalculated to the complex phase velocity using the polar reciprocity relation.After that we can invert for the viscoelastic anisotropy.
Using measurements of P waves we cannot determine the complete tensor of viscoelastic parameters.The six of the 21 viscoelastic parameters related to propagation of the S waves cannot be determined.The inversion for viscoelastic anisotropy is non-linear, but it can be linearized.The linearization and elimination of the S-wave related parameters can conveniently be performed using the so-called WAA parametrization.Since the inversion is run in iterations, it is applicable even to strong attenuation anisotropy.
Numerical tests on synthetic data proved that the inversion for the WAA parameters worked well.The convergence of the iterative   process was quite fast.The accuracy of the retrieved parameters was higher if the phase velocity and attenuation were inverted.Utilizing ray velocity and attenuation worked slightly worse because of the necessity of recalculating the ray to phase quantities before the inversion.In both inversions, the accuracy of the individual WAA parameters significantly varied.The lowest accuracy was achieved for parameters η x , η y and η z reflecting a rather low sensitivity of the P-wave velocity and attenuation to them.The applicability of the proposed methods was exemplified by determining anisotropic attenuation of a serpentinite rock from Val Malenco, Northern Italy.The ray velocity and ray attenuation were measured on a spherical sample of the rock with the diameter of 45.5 mm in 132 independent directions at the room temperature and under two pressure levels: 0.1 and 20 MPa.First, elastic anisotropy of the rock sample was determined and then attenuation anisotropy was inverted for.To our knowledge, the results obtained for the serpentinite rock sample represent the first accurately determined triclinic anisotropic attenuation from lab measurements.The measurements confirmed that anisotropic attenuation is remarkably sensitive to confining pressure.Since cracks are closing with increasing pressure, attenuation decreases.If cracks are preferentially oriented, the crack closure can diminish also the velocity anisotropy.However, changes in pressure acting on a rock with preferentially oriented crack systems can also induce changes in the directional variation of attenuation and a rotation of anisotropy axes.
We applied the proposed method to determination of anisotropic attenuation of a dry rock sample from lab measurements, but the applicability of the approach is broader.The approach can be applied to measurements of saturated rock samples to study fluid-driven attenuation and fluid flow in rocks.Attenuation effects similar to those measured in the lab will probably be observed also in seismic field experiments.Since the Earth's crust and the mantle are anisotropic and attenuating, future detailed seismological tomographic studies will include effects of anisotropic attenuation.Such studies can advantageously apply the model of viscoelastic anisotropy and follow the inversion schemes presented in this paper after adapting them to a different configuration and scale of the experiment and to a different frequency range of studied waves.

A C K N O W L E D G E M E N T S
We thank Jörg Renner, Georg Dresen, José Carcione and two anonymous reviewers for their helpful comments, and Prof H. Kern for providing us with the serpentinite sample.The technical help of the P-wave phase velocity anisotropy, attenuation anisotropy and Q-factor anisotropy.The anisotropy is calculated as a = 200 (U MAX − U MIN )/ (U MAX + U MIN ), where U MAX and U MIN are the maximum and minimum values of the respective quantity.

Figure 1 .
Figure 1.Stereographic plots of the P-wave velocities, attenuations and Q-factors in the model of the Torre Alfina xenolith.Left: phase intrinsic quantities and right: ray intrinsic quantities.For parameters of the model, see Table 1.

Figure 2 .
Figure2.The relative error between the true and predicted complex phase velocity as a function of the number of iterations of the inversion.The error is calculated as the maximum over all directions.The inversion is run for (a) the phase quantities and (b) the ray quantities.Blue/red: the real/imaginary part of the complex phase velocity.The higher errors in (b) than in (a) are due to the recalculation of ray to phase quantities in a sparse grid of directions (step of 15 • ) before the inversion.

Figure 3 .
Figure 3. Errors of the predicted (a) phase velocity and (b) attenuation using the inversion of phase quantities.The errors for several alternative grid steps (5 • , 10 • , 15 • , 20 • and 30 • ) are displayed as a function of noise in input data.The real part of noise (0-3 per cent) is five times lower than its imaginary part (0-15 per cent) expressing higher uncertainties in measurements of signal amplitudes than in traveltimes.

Figure 4 .
Figure 4. Errors of the predicted (a) phase velocity and (b) attenuation using the inversion of ray quantities.The errors for several alternative grid steps (5 • , 10 • , 15 • , 20 • and 30 • ) are displayed as a function of noise in input data.The real part of noise (0-3 per cent) is five times lower than its imaginary part (0-15 per cent) expressing higher uncertainties in measurements of signal amplitudes than in traveltimes.

Figure 5 .
Figure 5. Errors of the real (left) and imaginary (right) parts of the retrieved weak-anisotropy-attenuation (WAA) parameters.The errors are displayed as a function of noise in input data.The inversion is run for (a) the phase quantities and (b) the ray quantities in a grid of directions with step of 15 • .

Figure 6 .
Figure 6.Stereographic plots of the predicted P-wave phase velocities, attenuations and Q-factors in the model of the Torre Alfina xenolith retrieved by the inversion in a grid of directions with step of 15 • .Left: inversion of the phase quantities and right: inversion of the ray quantities.The noise level was up to ±3 and ±15 per cent in the real and imaginary parts of the complex phase velocity, respectively.

Figure 7 .
Figure 7.The experimental setup and measurement geometry shown in the horizontal plane; TP and RP are the radial transmitter and receiver, respectively (mostly generating and sensitive to the P waves).The axis of rotation is vertical.The red colour marks the setting ring for adapting a different diameter of individual spherical samples.

Figure 8 .
Figure 8. Ultrasonic signals observed for the serpentinite sample at confining pressure of 0.1 MPa.The waveforms are recorded in all directions at two different timescales.The origin of the timescale corresponds to the excitation time of the signal.

Figure 9 .
Figure 9. Example of one trace taken from the set in Fig. 8 showing how the amplitude of the direct P wave was measured.

Figure 10 .
Figure 10.Directional variations of ray velocity and ray amplitude measured for the serpentinite sample at confining pressure of (a) 0.1 MPa and (b) 20 MPa.The amplitudes in (a) and (b) are normalized to the maximum amplitude measured at 0.1 MPa.

Figure 11 .
Figure 11.The ray attenuation and ray quality factor for the serpentinite sample at confining pressure of (a) 0.1 MPa and (b) 20 MPa.The attenuation and quality factor are determined from the ray amplitudes shown in Fig. 10 and from the radiation pattern of the point source in the anisotropic serpentinite rock (see Section 3.2).

Figure 12 .
Figure 12.The phase velocity, attenuation and quality factor for the serpentinite sample at confining pressure of (a) 0.1 MPa and (b) 20 MPa.The phase quantities were calculated from the measured ray quantities and served as the input data for the inversion for parameters of viscoelastic anisotropy.

Figure 13 .
Figure 13.The predicted phase velocity, attenuation and quality factor for the serpentinite sample at confining pressure of (a) 0.1 MPa and (b) 20 MPa.The phase quantities were calculated using the weak-anisotropy-attenuation (WAA) parameters retrieved by the inversion.

Table 2 .
P-wave velocity and attenuation anisotropy of the synthetic model.
Notes: V phase , Āphase and Qphase are the mean P-wave phase velocity, attenuation and Q-factor; a phase V

Table 3 .
P-wave velocity and attenuation anisotropy of the serpentinite rock sample.

Table 4 .
WAA parameters of the serpentinite rock sample.