the Creative Commons Attribution 4.0 License.

Special issue: Advances in seismic imaging across the scales

**Research article**
04 Nov 2019

**Research article** | 04 Nov 2019

# Anisotropic P-wave travel-time tomography implementing Thomsen's weak approximation in TOMO3D

Adrià Meléndez Clara Estela Jiménez Valentí Sallarès and César R. Ranero

^{1},

^{1},

^{1},

^{2}

**Adrià Meléndez et al.**Adrià Meléndez Clara Estela Jiménez Valentí Sallarès and César R. Ranero

^{1},

^{1},

^{1},

^{2}

^{1}Barcelona Center for Subsurface Imaging, Institut de Ciències del Mar (CSIC), Barcelona, 08003, Spain^{2}Barcelona Center for Subsurface Imaging, ICREA at Institut de Ciències del Mar (CSIC), Barcelona, 08003, Spain

^{1}Barcelona Center for Subsurface Imaging, Institut de Ciències del Mar (CSIC), Barcelona, 08003, Spain^{2}Barcelona Center for Subsurface Imaging, ICREA at Institut de Ciències del Mar (CSIC), Barcelona, 08003, Spain

**Correspondence**: Adrià Meléndez (melendez@icm.csic.es)

**Correspondence**: Adrià Meléndez (melendez@icm.csic.es)

Received: 28 Feb 2019 – Discussion started: 13 Mar 2019 – Revised: 06 Sep 2019 – Accepted: 18 Sep 2019 – Published: 04 Nov 2019

We present the implementation of Thomsen's weak
anisotropy approximation for vertical transverse isotropy (VTI) media within TOMO3D, our code for 2-D and
3-D joint refraction and reflection travel-time tomographic inversion. In
addition to the inversion of seismic P-wave velocity and reflector depth,
the code can now retrieve models of Thomsen's parameters (*δ* and
*ε*). Here, we test this new implementation following four
different strategies on a canonical synthetic experiment in ideal conditions
with the purpose of estimating the maximum capabilities and potential weak
points of our modeling tool and strategies. First, we study the sensitivity
of travel times to the presence of a 25 % anomaly in each of the
parameters. Next, we invert for two combinations of parameters (*v*, *δ*, *ε* and *v*, *δ*, *v*^{⊥}), following two inversion strategies,
simultaneous and sequential, and compare the results to study their
performance and discuss their advantages and disadvantages. Simultaneous
inversion is the preferred strategy and the parameter combination (*v*, *δ*, *ε*)
produces the best overall results. The only advantage of
the parameter combination (*v*, *δ*, *v*^{⊥}) is a better recovery of the
magnitude of *v*. In each case, we derive the fourth parameter from the equation
relating *ε*, *v*^{⊥} and *v*. Recovery of *v*, *ε* and
*v*^{⊥} is satisfactory, whereas *δ* proves to be impossible to
recover even in the most favorable scenario. However, this does not hinder
the recovery of the other parameters, and we show that it is still possible
to obtain a rough approximation of the *δ* distribution in the medium by
sampling a reasonable range of homogeneous initial *δ* models and
averaging the final *δ* models that are satisfactory in terms of data
fit.

- Article
(9285 KB) -
Supplement
(911 KB) - BibTeX
- EndNote

An isotropic velocity field is the rare exception in the Earth subsurface. Anisotropy is a multiscale phenomenon, and its causes are diverse. In the crust, it can be produced by the preferred orientation of mineral grains or their crystal axes (Schulte-Pelkum and Mahan, 2014; Almqvist and Mainprice, 2017), the alignment of cracks and fracture networks and the presence of fluids (Crampin, 1981; Maultzsch et al., 2003; Yousef and Angus, 2016), or the bedding of layers much thinner than the wavelength used to explore them (Backus, 1962; Johnston and Christensen, 1995; Sayers, 2005). In the mantle, anisotropy is related to the alignment of olivine crystals due to mantle flow (Nicolas and Christensen, 1987; Montagner et al., 2007), aligned melt inclusions (Holtzman et al., 2003; Kendall et al., 2005), large-scale deformation (Vinnik et al., 1992; Vauchez et al., 2000) and pre-existing lithospheric fabric (Kendall et al., 2006), among others. Anisotropy has proven an informative physical property in the understanding of the Earth's interior (Ismaïl and Mainprice, 1998; Long and Becker, 2010), most particularly in continental rifts (Eilon et al., 2016), mid-ocean ridges (Dunn et al., 2001) and subduction zones (Long and Silver, 2008).

The theory of anisotropic wave propagation has been described in several publications (Kraut, 1963; Babuska and Cara, 1991). Numerous formulations of varying complexity have been proposed to approximate anisotropy depending on its magnitude and the symmetry conditions of the medium (Nye, 1957). Overall, 21 elastic stiffness parameters define the most general anisotropic medium with the lowest symmetry conditions, whereas for the highest symmetry not equivalent to isotropy, only five parameters are needed (Almqvist and Mainprice, 2017). Regarding the strength of anisotropy, in view of the overall success of isotropic methods in studying the Earth's subsurface, and of the experimental evidence and sample measurements available, it is admitted that the anisotropy is generally weak (Thomsen, 1986). Specifically, anisotropy is considered weak when Thomsen's parameters are much smaller than 1, i.e., for a ∼20 % or smaller velocity variation with angle. Precisely Thomsen (1986) presented the formulation for the transverse isotropy symmetry on weakly anisotropic media, which is the reference that we follow in this work. Thomsen's parameters are by far the most common and convenient combinations of stiffness tensor elements used in seismic anisotropy modeling (Tsvankin, 1996; Thomsen and Anderson, 2015). Applications of simpler approximations exist, as is the case of elliptical symmetry (Song et al., 1998; Giroux and Gloaguen, 2012), as well as others that assume the most general anisotropic model (Zhou and Greenhalgh, 2008).

The objectives of this work are (1) presenting the anisotropic version of
TOMO3D (Meléndez et al., 2015b) for the study of vertical transverse isotropy (VTI) weakly anisotropic
media in terms of Thomsen's parameters (*δ* and *ε*)
using P-wave arrival times and (2) comparing several parameterizations and
inversion strategies under optimal and equal conditions for all parameters
with the purpose of defining an upper limit of the code's capabilities, an
ideal but generalizable estimation of the code's performance, as well as
highlighting its potential weaknesses. Moreover, the development of this
code is motivated by the need to combine wide-angle and near-vertical
travel-time picks in field data applications, as we plan to do with the
trench-parallel 2-D profile in Sallarès et al. (2013), which is affected
by a ∼15 % anisotropy judging from the mismatch in the
interplate boundary locations obtained separately from near-vertical and
wide-angle data. The P-wave velocity, *δ* and *ε* models
obtained from the modeling of field data would be useful and geologically
informative by themselves, but they would also serve as initial models in
anisotropic full waveform inversion (FWI).

In the following section, we describe the anisotropy formulation and the modifications implemented on our 3-D joint refraction and reflection travel-time tomography code TOMO3D to incorporate the inversion of Thomsen's parameters. Next, in Sect. 3, we present the synthetic tests performed and their results, including accuracy and sensitivity analyses and synthetic inversions. In Sect. 4, these results are discussed and interpreted in terms of the ability of the code to retrieve both the velocity field and Thomsen's anisotropy parameters. Finally, in the last section, we summarize the main conclusions of this work.

The first part of this section is a general overview of the treatment of anisotropy within the field of seismic inversion, while the second one describes the implementation of Thomsen's weak VTI anisotropy formulation in the 3-D joint refraction and reflection travel-time tomography code TOMO3D.

## 2.1 Anisotropy in seismic inversion methods

Anisotropy was first incorporated to seismic inversion methods in travel-time tomography with the development of the linearized perturbation theory (Cervený, 1982; Cervený and Jech, 1982; Jech and Psencik, 1989). Previously, the approach to deal with anisotropy was to approximately remove its estimated effect to then apply an isotropic method (e.g., McCann et al., 1989). Linearized perturbation theory was first implemented in anisotropic travel-time tomography by Chapman and Pratt (1992) and Pratt and Chapman (1992), assuming the weak anisotropy approximation, which allowed them to use isotropic ray tracing and approximate anisotropy effects as being caused by small perturbations of the isotropic system. The initial development of anisotropic ray tracing is attributed to Cervený (1972). Methods for anisotropic ray tracing and travel-time computation depend on the symmetry assumptions made regarding the medium. The most common of those is rotational symmetry around a vertical pole. This formulation is known as VTI and also polar anisotropy (e.g., Rüger and Alkhalifah, 1996; Alkhalifah, 2002), and it is the simplest geologically applicable case: it reproduces the symmetry exhibited by minerals in sedimentary rocks and that produced by parallel cracks or fine layering. Furthermore, it significantly simplifies the mathematical formulae since anisotropy is defined by only five parameters, which contributes to a greater computational efficiency. The generalization of VTI to a tilted symmetry axis is the so-called tilted transverse isotropy (TTI). Some authors argue that it is not possible to distinguish TTI from VTI in real experimental cases without a priori information (Bakulin et al., 2009). Assuming the most general anisotropic media has also become rather usual, in particular with the improvement of computational resources, allowing for a more detailed and complex reconstruction of the subsurface physical properties (e.g., Zhou and Greenhalgh, 2005), although successful field data applications are yet to be achieved to the best of our knowledge. Regarding the inversion process, the main difficulty arises from the trade-off between velocity heterogeneity and anisotropy (e.g., Bezada et al., 2014). To the best of our knowledge, Stewart (1988) was the first to propose an inversion algorithm, specifically for the recovery of Thomsen's parameters in a weakly anisotropic VTI medium. Other authors have produced inversion algorithms for different formulations such as azimuthal anisotropy (e.g., Eberhart-Phillips and Henderson, 2004; Dunn et al., 2005) or a 3-D TTI medium (e.g., Zhou and Greenhalgh, 2008). Concerning FWI, anisotropy in active data is typically modeled following Thomsen's parameters and the VTI and/or TTI approximation for the medium. The first anisotropic wave propagators appeared during the 1980s and 1990s (e.g., Helbig, 1983; Alkhalifah, 1998) and new improvements on this matter continue today (e.g., Fowler et al., 2010; Duveneck and Bakker, 2011). When performing anisotropic FWI, both 2-D and 3-D, some authors choose to invert only for the velocity field, fixing the initial anisotropy models throughout the inversion because it simplifies the process (e.g., Prieux et al., 2011; Warner et al., 2013). However, other works have explored the feasibility of multiparameter inversions, that is, using different combinations of velocity and anisotropy parameters and of inversion strategies (e.g., Gholami et al., 2013a, b; Alkhalifah and Plessix, 2014).

## 2.2 Anisotropy in TOMO3D: Thomsen's weak VTI anisotropy formulation

We adapted TOMO3D (Meléndez et al., 2015b) to perform anisotropic ray tracing and travel-time calculations, as well as inversion of Thomsen's parameters for P-wave data (Meléndez et al., 2015a). In TOMO3D, the forward problem solver is parallelized to simultaneously trace rays for multiple sources and receivers, and it uses an hybrid ray-tracing algorithm that combines the graph or shortest path method (Moser, 1991) and the bending refinement method (Moser et al., 1992). The inverse problem is solved sequentially using the LSQR algorithm (Paige and Saunders, 1982). Velocity models are discretized as 3-D orthogonal and vertically sheared grids that can account for topography and/or bathymetry. Velocity values are assigned to the grid nodes, and the velocity field is built by trilinear interpolation within each cell. Apart from first-arrival travel times, the code allows for the inversion of reflection travel times to obtain the geometry of major geological boundaries associated with impedance contrasts that produce strong seismic energy reflections in the data recordings. Such reflecting interfaces are modeled as 2-D grids independent of the velocity grid. The code is also prepared to extract information from the water-layer multiple of refracted and reflected seismic phases (Meléndez et al., 2014). A detailed description of the code can be found in Meléndez (2014).

Our anisotropy formulation is based on Thomsen (1986) and specifically in the following weakly anisotropic velocity equation for the P-wave velocity:

where *v*^{a} is the anisotropic velocity, *v* is the velocity along the
symmetry axis (*α*_{0} in Thomsen, 1986), *θ* is the angle with
respect to the symmetry axis, and *δ* and *ε* are Thomsen's
parameters controlling the anisotropic P-wave propagation. Studying the
cases of *θ*=0 and $\mathit{\theta}=\mathit{\pi}/\mathrm{2}$, the meaning of
*ε* becomes clear: it is the relative difference between the
velocities along and across the symmetry axis that we refer to as parallel
and perpendicular velocities, respectively.

According to Thomsen (1986), the meaning of *δ* is far from intuitive,
but the author states that it is associated with the near-vertical anisotropic
response and shows that it relates *v* and the normal move-out velocity
(*V*_{NMO}). *V*_{NMO} models are a byproduct of multichannel seismic
reflection data processing, a mathematical construct that involves the
assumptions of a stratified media with constant velocity layers and of small
spread, i.e., near-vertical propagation. It does not seem wise to try
estimating it by other means, less so if the data and modeling used do not
necessarily fulfill the assumptions for the normal move-out correction that
define *V*_{NMO}. Moreover, we want to combine travel times from as many
types of seismic data sets as possible, notably from multichannel reflection
(near-vertical propagation) and wide-angle (subhorizontal propagation)
experiments, in order to have the best polar coverage, and with that, the
best recovery of *v* and *ε* (or *v*^{⊥}). Thus, we do not
consider *V*_{NMO} to be a useful parameter in describing the general
anisotropic VTI media for our modeling method, and we did not implement
parameterizations (*v*,*V*_{NMO},*ε* and *v*,${V}_{\mathrm{NMO}},{v}^{\perp}$). We do
think, however, that *V*_{NMO} models can help in the building of initial
*δ* models, just as the comparison between near-vertical propagation
and subhorizontal data can provide an initial estimation of *ε*.
In conclusion, we only considered Eq. (2), and we implemented two
parameterizations of the medium: (*v*, *δ*, *ε*) and (*v*, *δ*, *v*^{⊥}). From here on, for simplicity, we will refer to them as P[*ε*] and P[*v*^{⊥}], respectively.

The linearized inverse problem matrix equation including anisotropy parameters for a refraction-only case is as follows:

Smoothing (**L**) and damping (**D**)
constraints for *δ* and *ε* parameters follow the same
formulation described in Meléndez (2014) for velocity parameters. The
kernels (**G**) have been modified to account for anisotropy.
The linearized and discretized equation that relates the travel-time residual
of the *n*th refracted pick to changes in the model parameters is written as

In each of the three terms, the first summation corresponds to the model
cells that are illuminated by the *n*th ray path. A cell is considered
illuminated if it contains a ray path segment. The second summation is over
the eight nodes of each of those cells. In the third term, *ε* is
replaced by *v*^{⊥} when using this alternative parameterization.
$\mathrm{\Delta}{t}_{n}^{\mathrm{0}}$ is the travel-time residual for the *n*th refracted pick,
*u*^{a}=1/*v*^{a} is the anisotropic slowness, $u=\mathrm{1}/v$ is the along-axis
slowness, Δ*u*_{i}, Δ*δ*_{j} and Δ*ε*_{k} (or $\mathrm{\Delta}{v}_{k}^{\perp}$) are the
parameter perturbations for each illuminated cell in their respective grids,
and the *r*_{m} factors are the weights that distribute these perturbations
among the eight nodes of each illuminated cell according to the trilinear
interpolation used to define the four fields (*u*, *δ*, *ε*
and *v*^{⊥}).

In order to build the kernel matrices, we need to compute two partial
derivatives for each model parameter. The first-order partial derivative of
travel time with respect to the anisotropic slowness is the ray path segment
*s*_{i} within each cell that is covered in a given travel time at a given
slowness:

From Eq. (1), the first-order partial derivatives of the anisotropic slowness
with respect to the model parameters (*u*, *δ*, *ε* and
*v*^{⊥}) are as follows:

We have performed a number of tests using canonical synthetic models made of an anomaly centered in a uniform background with two main objectives: (1) checking that the newly implemented anisotropic travel-time tomography method works properly and (2) providing a quantitative measure of the potential recovery of anisotropy based on P-wave travel times alone. All data and files used for these synthetic tests are available at the digital CSIC repository (Meléndez et al., 2019). First, we run a sensitivity test to assess the effect that a variation in each model parameter has in the synthetic travel times, and we calibrated the code by comparing the synthetic data that it generates to analytically calculated data. Next, we performed a number of synthetic inversion tests considering both possible parameterizations and inversion strategies. These tests are conducted under ideal and equal conditions for all parameters with the purpose of obtaining an upper limit but widely applicable estimation of the code's performance and detecting any potential weak points.

The models in all these tests are cubes with 5 km long edges. The background
model of all four parameters is set to a constant value; i.e., *v*, *δ*,
*ε* and *v*^{⊥} background models are homogeneous. Note that
the *z*-axis positive direction points downwards. Grid spacing is 0.125 km for
the four parameters in all three dimensions, so that differences in model
discretization do not influence the test results. The volume of the anomaly
is determined by the 3*σ* region of a 3-D Gaussian function centered
in the cube setting 3*σ*=0.5 km. The values of *v*, *δ*,
*ε* and *v*^{⊥} within this volume are homogeneously increased,
resulting in a discretized representation of a spherical anomalous body.

## 3.1 Sensitivity

We define the sensitivity of a parameter as the difference between the first-arrival travel times with and without the anomaly in that parameter. Figure 2 shows synthetic and analytic sensitivities for two selected meridians. We express sensitivity both as normalized and as relative travel-time difference. In its normalized form, travel-time differences are divided by the greatest of these differences among all parameters, i.e., normalized to 1, whereas in its relative form, these differences are given with respect to the travel time without the anomaly and multiplied by 100. Note that the analytic response does not change between meridians, i.e., given the symmetry of the models and of the anisotropic formulation, sensitivity is independent of the azimuth angle.

The acquisition configuration consists of 482 diametrically opposed source–receiver pairs. Each receiver records exclusively the first-arrival travel time from its corresponding source for a total of 482 travel times (Fig. 1). Sources and receivers are located at 2.5 km of the center of the model, at the locus defined by the surface of the sphere inscribed in the cube and placed at the crossing points of 32 meridians with 15 parallels and at each pole.

According to the definition of sensitivity, alternately for *v*, *δ* and
*ε*, the said anomaly was added at the center of the cube
representing a 25 % increase on the background value, while the models
for the rest of parameters in the parameterization remained homogeneous.
Table 1 summarizes the background and anomaly values for all parameters in
each sensitivity test, along with their equivalence in the alternative
parameterization. Since *v* is related to *ε* and *v*^{⊥} through
Eq. (2), its sensitivity pattern changes depending on the parameterization
used (Fig. 2). Indeed, a 25 % increase in *v* with respect to equivalent
background models in P[*ε*], and P[*v*^{⊥}] (Table 1) yields
different sensitivities because of the different parameters involved
(*ε* or *v*^{⊥}) in the representation of the medium.
Contrarily, the *δ* sensitivity pattern is independent of the
parameterization used. *ε* and *v*^{⊥} sensitivities are only
defined in their respective parameterizations. However, a 25 % increase in
*v*^{⊥} yields an equivalent *ε* of 0.45, which is greater than
the ∼0.2 limit for weak anisotropy approximation. Thus,
instead of establishing the comparison with *v*^{⊥} sensitivity based on
a proportional anomaly increment and measuring its effect on travel times, we
based it on an equal travel-time change; i.e., the same change in travel time
requires a change of 25 % in *ε* but only a ∼3.4 % change in *v*^{⊥} (Table 1), indicating that data are
∼7 times more sensitive to *v*^{⊥} than to *ε*
changes.

*v* sensitivity in P[*ε*] is the highest for all angles, 4.5 %
to 5 % (Fig. 2), and it can be shown that expressed in its relative form,
the analytic solution is a constant 5 % (see mathematical proof in the
Supplement). In its normalized form, it follows the same
sinusoidal pattern as in P[*v*^{⊥}] (Fig. 2e), and both have equal maxima
in the directions parallel to the symmetry axis. In both cases, minima are
found in the directions perpendicular to the symmetry axis, but in
P[*v*^{⊥}] they reach down to 0, whereas in P[*ε*] the value is
∼0.85. As expected from Eq. (1), *ε* sensitivity
goes to 0 % in the directions parallel to the symmetry axis and has its
maxima (∼0.8 % or ∼0.15) for the polar
angles perpendicular to it. *v*^{⊥} sensitivity would follow the same
angular dependence but with maxima of the order of magnitude of *v*
sensitivities. Finally, *δ* sensitivity is, at its maxima, more than
1 order of magnitude smaller than *v* sensitivity, around 0.25 % or 0.05.
These sensitivity results indicate that we can generally expect similar
recoveries for *v* and *v*^{⊥}, better than for *ε*, and that
retrieving *δ* might prove complicated. Keep in mind that these
sensitivities for *v*, *δ* and *ε* are produced by an anomaly
that represents a 25 % increase with respect to the background value, and
that an anomaly in *v*^{⊥} would produce the same sensitivity pattern as
*ε* with only a ∼3.4 % increment.

The differences in synthetic sensitivities between meridians arise from the
discretization of the model space in a Cartesian system of coordinates. Such
approximation inevitably defines privileged directions for ray tracing and
consequently produces differences in synthetic travel times. The mismatch
between synthetic and analytic sensitivities (Fig. 2) occurs because the
discretization used cannot represent the surface of a perfect sphere. These
effects are most notable in the *v* sensitivity, in P[*ε*] and to a
lesser extent in P[*v*^{⊥}], precisely because it is the most sensitive
parameter, and thus the errors in the representation of a sphere and the
existence of privileged directions have a much larger influence on the
calculated travel times. Figure S1 shows how refining the *v* model reduces the
relative travel-time error (Fig. S1a in the Supplement) and generates a more accurate
sensitivity pattern (Fig. S1b). In a real case study, one can always refine
the grid spacing of a particular parameter to achieve better accuracy, but
here we wish to test the performance of the code in the modeling and
recovery of each parameter under the same conditions, i.e., equivalent
anomalies and identical model discretization.

## 3.2 Accuracy

For the four simulations in the sensitivity analysis, we compared the
synthetic travel times obtained with our code to the analytic solution to
quantify the accuracy of the code's performance. The comparison along the
two selected meridians at 0 rad and *π*∕4 rad azimuths is displayed in
Fig. 3 expressed as relative travel-time error, i.e., the difference between
synthetic and analytic travel times relative to the latter in percentage.
Table 2 contains the mean of the relative travel-time errors and their
respective mean deviations for these four tests and along the two selected
meridians, as well as the overall values for each of them.

Comparison of Fig. 2a–b with the values in Table 2 and Fig. 3 indicates that the forward calculation of travel times is accurate enough with respect to the travel-time residuals expected for the selected anomalies. Sensitivity is at least 5 times and up to 2 orders of magnitude greater than travel-time accuracy errors depending on the parameters, with the exception of angles for which sensitivity tends to zero. Figure S2 illustrates that the code is able to reproduce nearly identical accuracies using both parameterizations. Furthermore, given that we are using the same synthetic travel times that we used for the sensitivity analysis, this also implies that the sensitivity patterns obtained with the alternative parameterization would be virtually equal to those in Fig. 2.

## 3.3 Inversion results

For the inversion tests, we considered a synthetic medium defined by the
anomaly models of all four parameters. Here, we refer to these models as
target models, and the goal of the inversion is to retrieve the
heterogeneity in each of them. These tests are conducted for the two
parameterizations of the anisotropic medium described in Sect. 2:
P[*ε*] and P[*v*^{⊥}]. Note that, in order to perform the
inversion tests on equivalent cases for both parameterizations, the
heterogeneity in *v*^{⊥} is calculated with Eq. (2) considering the 25 %
anomalies in *v* and *ε*, which yield a ∼29.3 % anomaly in *v*^{⊥} (Table 3). If not indicated otherwise, we use
background models as initial models. Finally, we study the potential
recovery of *δ* because inverting this particular parameter proves
notoriously difficult due to its low sensitivity (Fig. 2).

The synthetic data set is made of 114 sources, each recorded at 113 receivers for a total 12 882 first-arrival travel times. For the acquisition geometry, again sources and receivers are located at the surface defined by the sphere inscribed in the cube (Fig. 1). The 114 positions at the surface of this sphere are shared by sources and receivers, and each receiver records all sources, except for the one source located at its same position.

For both parameterizations, we compared two inversion strategies: simultaneously inverting for all parameters and a two-step sequential inversion. First, in Figs. 5 and 6, we show the best results for the simultaneous inversion strategy. For each parameterization, we derived the fourth parameter by applying Eq. (2).

Table 4 shows several statistical measures to quantify the quality of these
inversion results. As a measure of data fit improvement, we provide the root mean square (rms)
of travel-time residuals for the first and last iterations. As a measure of
model recovery or fit, for each parameter, we calculate the mean relative
difference for the background area between the inverted model and either the
target or the initial one, as they are identical in this area, as well as
for the anomaly area comparing the inverted model to both the target and the
initial models. In the case of a perfect recovery, mean relative difference
for the background area would be 0 %, whereas for the anomaly area, it
would be 0 % when using the target model as a reference and 25 %
(∼29.3 % for *v*^{⊥}) when comparing to the initial
model.

In an attempt to improve the recovery of *δ*, we repeated these two
tests for different values of the smoothing constraints, but it proved
impossible. Correlation lengths tested for all four parameters include 0.25 km
(twice the grid spacing), 0.5 and 1 km. The weights of the smoothing
submatrices for each parameter, *λ* in Eq. (3), were varied between 1
and 100, with intermediate values of 2, 5, 10, 20, 30 and 60. For
successful inversions, very similar results for the other parameters were
obtained regardless of the final *δ* model. In other words, the low
sensitivity of *δ* makes it extremely hard, if not impossible, to
recover this parameter from travel-time data, but for this same reason it has
little or no influence on the recovery of *v* and *ε* or *v*^{⊥}.

The two-step sequential inversion strategy was also tested for both
parameterizations, P[*ε*] and P[*v*^{⊥}]. For the first step, we
tested two options: (a) inverting for *v* while fixing *δ* and
*ε* or *v*^{⊥} and (b) fixing only *δ*. In the second
step, we used the inverted models from step 1 as initial models and tested
three options: (c) inverting for all three parameters, (d) fixing only
*δ*, when following option (a) in step 1, and (e) fixing *v* and/or
*ε*or *v*^{⊥}, when following option (b) in step 1. Figure 4
shows a flowchart describing all the sequential inversion combinations
tested. Again, smoothing constraints were varied for similar correlation
lengths and submatrix weights as detailed for the simultaneous inversion
case.

As indicated in Fig. 4, in the case of P[*ε*], the best result
(Fig. 7) was obtained inverting for *v* and *ε* while fixing
*δ* in the first step and fixing only *ε* in the second
step, whereas for P[*v*^{⊥}], the best combination for the two-step
inversion (Fig. 8) was fixing only *δ* in step 1 and inverting for
the three parameters in step 2. Tables 5 and 6 summarize the statistical
quantification of data and model fit for each parameterization.

## 3.4 Modeling *δ*

Observing that good results for *v* and *ε* or *v*^{⊥} are
achieved regardless of the result in *δ* and knowing that the
sensitivity of *δ* is notably smaller than that of the other
parameters, we explored a strategy to have an estimate of this parameter.
First, as a reference, we considered an unrealistically optimal scenario in
which the real *v* and *ε* or *v*^{⊥} models are known to us.
Figures 9 and 10 show the resulting *δ* achieved by repeating inversions
in Figs. 5 and 6 with *v* and *ε* or *v*^{⊥} target models as
initial models. Table 7 summarizes the travel-time residuals' rms and the mean
relative difference for each parameter in these inversions. Again, these two
tests were repeated for ranges of smoothing constraints in all four
parameters, as described for the cases in Figs. 5 and 6. Table 7 and Figs. 9
and 10 correspond to the best results obtained, which indicate that the
recovery of *δ* is, at best, extremely complicated due to the limited
sensitivity of travel-time data to changes in this parameter.

Next, we decided to try neglecting *δ* in Eq. (1), and we repeated a
number of inversions, such as the ones displayed in Figs. 5 and 6, following

The purpose of these tests was checking whether it was possible to invert
*v* and *ε* or *v*^{⊥} with data generated following Eq. (1)
using the approximation in Eq. (11), given that the influence of *δ* on
the results for other parameters is rather small, that *δ* cannot be
accurately retrieved from travel time alone and that it has the smallest
sensitivity. To do so, a homogeneous model of *δ*=0 was fixed
throughout the inversions. These tests were unsuccessful, with noticeably
poorer results than when considering a dependence on *δ* (Table 8).
However, they were useful in proving that even if a detailed *δ* model
is not necessary to successfully retrieve the other parameters, at least a
rough approximation of the *δ* field is needed to recover the other
parameters, e.g., the background *δ* model that we used as the initial
model in inversions displayed in Figs. 5 and 6.

Finally, we tested whether it would be possible to obtain at least this
rough approximation of *δ* in the medium, valid in the sense that it
allows for the successful recovery of the rest of parameters, using any a
priori information available such as compilations of anisotropy measurements
(e.g., Thomsen, 1986; Almqvist and Mainprice, 2017). Once again, we repeated
inversions from Figs. 5 and 6 (initial *δ*=0.16), now using
different homogeneous initial models for *δ* within a range of
possible values from 0.1 to 0.24. Table 9 contains the initial and final rms
of travel-time residuals, as well as *δ* mean values for the inverted
model along with the corresponding mean deviations. It is straightforward to
note that, for a central subrange of the tested initial *δ* values,
final rms values are an order of magnitude smaller, a few tenths of
a millisecond compared to the few milliseconds outside this subrange.
Specifically, very similar results to those in Figs. 5 and 6 in terms of
travel-time, residuals' rms are produced by initial *δ* values between
0.13 and 0.22 for P[*ε*], and between 0.12 and 0.22 for
P[*v*^{⊥}]. The narrowing of the initial *δ* distribution to a
smaller subrange of mean *δ* values for the inverted models is
indicative of a good general convergence trend.

The rough estimate of the *δ* field could be built, for instance, as
the average of the mean *δ* values for the inverted models in the
central subranges defined by the change in magnitude of the final rms of
travel-time residuals. One final inversion could be run using a homogeneous
initial *δ* model with this average value, with the additional option
of fixing it and inverting only for the other two parameters. As mentioned
in Sect. 2.2, potentially more detailed initial *δ* models could
be obtained from the normal move-out correction of near-vertical reflection
seismic data.

## 3.5 Discussion

We have tested two parameterizations of the VTI anisotropic media,
P[*ε*] (*v*, *δ*, *ε*) and P[*v*^{⊥}] (*v*, *δ*, *v*^{⊥}), and two inversion strategies, simultaneously inverting for all three
parameters and a two-step sequential process fixing some of the parameters
in each step. We consider three criteria for evaluating and comparing the
quality of the inversion results obtained following the four possible
combinations of strategies and parameterizations: visual inspection of the
results, as well as travel-time data and model fits. For both inversion
strategies and both parameterizations, the recovery of the parameters is
positively correlated with their respective sensitivities; the more
sensitive parameters are systematically better recovered.

For the simultaneous inversion, both parameterizations were able to produce
acceptable final results (Figs. 5 and 6). According to our tests,
P[*ε*] provides the best outcome, specifically because data and
model fits (Table 4) are better for this option but, more importantly,
because of the difference in the quality of the recovery of *ε*.
However, visual comparison of the recovery of *v* as well as the *v* model fit in
Table 4 indicates that P[*v*^{⊥}] yields a slightly better result
regarding the magnitude of the anomaly for this parameter, as is
particularly evident at its center. In P[*ε*], given the
disparity in sensitivities between *v* and *ε*, the former might be
prone to accounting for data misfit that actually corresponds to the latter,
whilst in P[*v*^{⊥}] this disparity is less acute. The boundaries of the
anomalies for both velocities are still better retrieved with
P[*ε*]. Also, the recovery of *δ*, even though it is far
from acceptable, is significantly better in the case of P[*ε*].
This might be explained because the disparity of *δ* sensitivity is
greater with respect to the other two parameters in P[*v*^{⊥}] than in
P[*ε*].

In general, sequential inversion is a more complex process that requires
more human intervention and fine tuning in each step. In addition, fixing
some of the parameters in the first step may result in the inverted
parameters artificially accounting for part of the data misfit that is
actually related to the fixed ones. This can more easily lead convergence
into a local minimum, and it might be impossible to correct this tendency in
the second step. For this inversion strategy, it is also P[*ε*]
that produces the best results (Figs. 7 and 8). Data fit and the model fit of
*ε* are slightly better than for the simultaneous inversion of
this parameterization (Tables 4 and 5), whereas model fits and the visual
aspect of both velocities are almost identical to those obtained by
simultaneously inverting all three parameters (Figs. 5 and 7). Visually, it
is difficult to decide whether the recovery of *ε* is better or
not than for the simultaneous inversion (Figs. 5 and 7). As for *δ*,
recovery is unsuccessful and artifacts appear in the background area of the
model but, according to both its model fit and its visual aspect, it is
notably better than for the simultaneous inversion. As in the case of
simultaneous inversion, the only advantage of using P[*v*^{⊥}] instead of
P[*ε*] is that it yields a better recovery of the anomaly
magnitude of *v* (Tables 4 and 6). The results for both velocities are
virtually identical to those obtained by simultaneous inversion of this
parameterization (Figs. 6 and 8). *δ* and *ε* are not
properly retrieved, but the results are significantly better than for the
simultaneous inversion of this same parameterization.

*δ* has been shown to be by far the most complicated parameter to
retrieve because of the low sensitivity of travel time to its variation (Fig. 2).
Even when excellent *v* and *ε* or *v*^{⊥} models are
available, i.e., the target models for these parameters in our synthetic
tests, the recovery of *δ* is limited at best (Figs. 9 and 10 and Table 7).
However, and for the same reason, poor recoveries of *δ* do not
affect the recovery of the other two parameters, meaning that a detailed
*δ* model is not necessary to satisfactorily retrieve *v* and
*ε* or *v*^{⊥} (Figs. 5–10). Still, our inversion tests also
proved that neglecting *δ* in Eqs. (1) and (3) is not an option, the
accuracy in the recovery of the other parameters resulting severely affected
(Table 8). Thus, even if a detailed inversion of *δ* is, at the very
least, hard to achieve, and it is not needed for a successful result in the
other parameters, some sort of simple, even homogeneous, initial *δ*
model with a value or values about the average *δ* in the medium is
necessary for a good recovery of the other parameters.

We showed that given some a priori information on the range of possible
*δ* values in the medium, it should be possible to create the
necessary initial *δ* model. In order to illustrate this, we chose a
range of *δ* values for the initial model and reran the inversions in
Figs. 5 and 6. The results indicate that for any homogeneous initial
*δ* model in a certain subrange close to the actual average *δ*
value of the medium, the results for *v* and *ε* or *v*^{⊥} are
satisfactory and virtually identical to those of the original inversions
(Table 9). This subrange is easily defined by looking at the final rms of
travel-time residuals, which experiences a notorious change of 1 order of
magnitude. Any model within this subrange works similarly well as an initial
*δ* model. Alternatively, a possibly more robust selection of the
constant value for a homogeneous initial *δ* model might be the mean
(or also the median or the mode) of the mean *δ* values for the
inverted models in this subrange. It is worth noting that, whereas for the
purpose of this work we used the same discretization for all parameters, in
a real case study it would probably be recommendable to use a coarser
discretization for *δ* than for the other parameters and in general a
finer discretization for the more sensitive parameters (Fig. S1). Indeed, a
heterogeneity of a given spatial scale and relative variation will produce a
greater effect on data for a parameter of greater sensitivity. Thus, for a
parameter of greater sensitivity, it will be easier for the code to identify
smaller heterogeneities both in scale and variation, which will require a
finer grid.

We have successfully implemented and tested a new anisotropic travel-time
tomography code. For this implementation, we had to modify both the forward
problem and the inversion algorithms of the TOMO3D code (Meléndez et
al., 2015b). The forward problem was adapted to compute the velocities
observed by rays considering Eq. (1) for the weak VTI anisotropy formulation
in Thomsen (1986). The inversion solver was extended to include the
*δ*, *ε* and *v*^{⊥} kernels in the linearized forward
problem matrix equation, as well as smoothing and damping matrices for these
parameters defined following the same scheme as for velocity in the
isotropic code (Eqs. 3–10).

Regarding the synthetic tests, we first determined the sensitivity of
travel-time data to changes in each of the parameters defining anisotropy in
the medium (*v*, *δ*, *ε*, *v*^{⊥}) (Fig. 2), and we
checked the proper performance of the code by comparing the synthetic
travel times produced with their respective analytic solutions (Fig. 3).
Next, we performed canonical inversion tests to compare two possible media
parameterization, P[*ε*] and P[*v*^{⊥}], and two possible
inversion strategies, simultaneous and sequential. According to our tests,
both parameterizations have their strengths: P[*ε*] produces the
best overall result in the sense that all parameters are acceptably
recovered, with the exception of *δ*, and trade-off between parameters
is lower, but P[*v*^{⊥}] yields the best result for the magnitude of the
anomaly in *v*. Regarding the inversion strategy, simultaneous inversion is
more straightforward and involves less human intervention, and given that
both strategies yield similar results, it would be our first choice.
Sequential inversion is always a more complex process that can be shown to
work in a synthetic case because the target models are available, but in
field data applications the complexity would most likely be unmanageable.
These tests were conducted under ideal conditions, and thus the conclusions
provided by their results are an upper limit but generalizable estimation of
the strengths and weaknesses in the code's performance.

An acceptable recovery of *δ* turned out to be impossible due to the
small sensitivity of travel times to this parameter, but we verified that it
cannot simply be neglected in the equations. Whereas the recovery of the
other parameters is not significantly affected by that of *δ*, a rough
estimate of the average *δ* value in the medium is necessary and
sufficient to generate a homogeneous initial model that allows for
satisfactory inversion results in these other parameters. We also proved
that it is possible to obtain it, provided that some a priori knowledge on
*δ* values in the medium is available to define a range of plausible
values, such as field or laboratory measurements.

The anisotropic version of TOMO3D is available only for academic purposes on our group website at http://www.barcelona-csi.cmima.csic.es/software-development (Meléndes et al., 2015a).

All data and files used in the synthetic tests presented in this article are available at the digital CSIC repository. The corresponding DOI is https://doi.org/10.20350/digitalCSIC/8957 (Meléndes et al., 2019).

The supplement related to this article is available online at: https://doi.org/10.5194/se-10-1857-2019-supplement.

The formulation of the overarching research goals of this work is a product of discussion among the four co-authors. AM was in charge of software development and data curation, analysis, visualization and validation. AM also prepared the manuscript. The methodology for the synthetic tests was designed by AM, CEJ and VS. CRR was responsible for the acquisition of financial support.

The authors declare that they have no conflict of interest.

This article is part of the special issue “Advances in seismic imaging across the scales”. It is a result of the 18th International Symposium on Deep Seismic Profiling of the Continents and their Margins, Cracow, Poland, 17–22 June 2018.

We thank all our fellows at the B-CSI for their contribution to this work. We also wish to thank Marko Riedel for his constructive comments at the 18th edition of the biennial International Symposium on Deep Seismic Profiling of the Continents and their Margins (SEISMIX 2018). The comments of the anonymous referees have also contributed significantly to the improvement of our work.

Adrià Meléndez and Clara Estela Jiménez are funded by Respol through the SOUND collaboration project with CSIC, and the work in this paper was conducted at the Grup de Recerca de la Generalitat de Catalunya 2009SGR146: Barcelona Center for Subsurface Imaging (B-CSI).

This paper was edited by Michal Malinowski and reviewed by two anonymous referees.

Alkhalifah, T.: Acoustic approximations for processing in transversely isotropic media, Geophysics, 63, 623–631, https://doi.org/10.1190/1.1444361, 1998.

Alkhalifah, T.: Travel Time computation with the linearized eikonal equation for anisotropic media, Geophys. Prospect., 50, 373–382, https://doi.org/10.1046/j.1365-2478.2002.00322.x, 2002.

Alkhalifah, T. and Plessix, R. É.: A recipe for practical full-waveform inversion in anisotropic media: An analytical parameter resolution study, Geophysics, 79, R91–R101, https://doi.org/10.1190/geo2013-0366.1, 2014.

Almqvist B. S. G. and Mainprice, D.: Seismic properties and anisotropy of the continental crust: Predictions based on mineral texture and rock microstructure, Rev. Geophys., 55, 367–433, https://doi.org/10.1002/2016RG000552, 2017.

Babuska, V. and Cara, M.: Seismic Anisotropy in the Earth, 1st edition, Kluwer Academic Publishers, Dordrecht, the Netherlands, 1991.

Backus G. E.: Long-wave elastic anisotropy produced by horizontal layering, J. Geophys. Res., 67, 4427–4440, https://doi.org/10.1029/JZ067i011p04427, 1962.

Bakulin, A., Woodward, M., Osypov, K., Nichols, D., and Zdraveva, O.: Can we distinguish TTI and VTI media?, in: SEG Technical Program Expanded Abstracts 2009 of the 79th SEG Annual Meeting, Houston, USA, 25–30 October 2009, 226–230, https://doi.org/10.1190/1.3255313, 2009.

Bezada, M. J., Faccenda, M., Toomey, D. R., and Humphreys, E.: Why ignoring anisotropy when imaging subduction zones could be a bad idea, AGU Fall Meeting, San Francisco, USA, 15–19 December 2014, S41D-04, 2014.

Cervený, V.: Seismic rays and ray intensities in inhomogeneous anisotropic media, Geophys. J. R. Astr. Soc., 29, 1–13, https://doi.org/10.1111/j.1365-246X.1972.tb06147.x, 1972.

Cervený, V.: Direct and inverse kinematic problems for inhomogeneous anisotropic media – linearized approach, Contr. Geophys. Inst. Slov. Aca. Sci., 13, 127–133, 1982.

Cervený, V. and Jech, J.: Linearized solutions of kinematic problems of seismic body waves in inhomogeneous slightly anisotropic media, J. Geophys., 51, 96–104, 1982.

Chapman C. H. and Pratt, R. G.: Traveltime tomography in anisotropic media – I. Theory, Geophys. J. Int., 109, 1–19, https://doi.org/10.1111/j.1365-246X.1992.tb00075.x, 1992.

Crampin, S.: A review of wave motion in anisotropic and cracked elastic-media, Wave Motion, 3, 343–391, https://doi.org/10.1016/0165-2125(81)90026-3, 1981.

Dunn, R. A., Toomey, D. R., Detrick, R. S., and Wilcock, W. S. D.: Continuous Mantle Melt Supply Beneath an Overlapping Spreading Center on the East Pacific Rise, Science, 291, 1955–1958, https://doi.org/10.1126/science.1057683, 2001.

Dunn, R. A., Lekic, V., Detrick, R. S., and Toomey, D. R.: Three-dimensional seismic structure of the Mid-Atlantic Ridge (35^{∘} N): Evidence for focused melt supply and lower crustal dike injection, J. Geophys. Res., 110, B09101, 1–17, https://doi.org/10.1029/2004JB003473, 2005.

Duveneck, E. and Bakker, P. M.: Stable P-wave modeling for reverse-time migration in tilted TI media, Geophysics, 76, S65–S75, https://doi.org/10.1190/1.3533964, 2011.

Eberhart-Phillips, D. and Henderson, M.: Including anisotropy in 3-D velocity inversion and application to Marlborough, New Zealand, Geophys. J. Int., 156, 237–254, https://doi.org/10.1111/j.1365-246X.2003.02044.x, 2004.

Eilon, Z., Abers, G. A., and Gaherty, J. B.: A joint inversion for shear velocity and anisotropy: the Woodlark Rift, Papua New Guinea, Geophys. J. Int., 206, 807–824, https://doi.org/10.1093/gji/ggw177, 2016.

Fowler, P. J., Du, X., and Fletcher, R. P.: Coupled equations for reverse-time migration in transversely isotropic media, Geophysics, 75, S11–S22, https://doi.org/10.1190/1.3294572, 2010.

Gholami, Y., Brossier, R., Operto, S., Ribodetti, A., and Virieux, J.: Which parameterization is suitable for acoustic vertical transverse isotropic full waveform inversion? Part 1: Sensitivity and trade-off analysis, Geophysics, 78, R81–R105, https://doi.org/10.1190/geo2012-0204.1, 2013a.

Gholami, Y., Brossier, R., Operto, S., Prieux, V., Ribodetti, A., and Virieux, J.: Which parameterization is suitable for acoustic vertical transverse isotropic full waveform inversion? Part 2: Synthetic and real data case studies from Valhall, Geophysics, 78, R107–R124, https://doi.org/10.1190/geo2012-0203.1, 2013b.

Giroux, B. and Gloaguen, E.: Geostatistical traveltime tomography in elliptically anisotropic media, Geophys. Prospect., 60, 1133–1149, https://doi.org/10.1111/j.1365-2478.2011.01047.x, 2012.

Helbig, K.: Elliptical anisotropy – Its significance and meaning, Geophysics, 48, 825–832, https://doi.org/10.1190/1.1441514, 1983.

Holtzman, B. K., Kohlstedt, D. L., Zimmerman, M. E., Heidelbach, F., Hiraga, T., and Hustoft, J.: Melt segregation and strain partitioning: implications for seismic anisotropy and mantle flow, Science, 301, 1227–1230, https://doi.org/10.1126/science.1087132, 2003.

Ismaïl, W. B. and Mainprice, D.: An olivine fabric database: an overview of upper mantle fabrics and seismic anisotropy, Tectonophysics, 296, 145–157, https://doi.org/10.1016/S0040-1951(98)00141-3, 1998.

Jech, J. and Psencik, I.: First-order perturbation method for anisotropic media, Geophys. J. Int., 99, 369–376, https://doi.org/10.1111/j.1365-246X.1989.tb01694.x, 1989.

Johnston, J. E. and Christensen, N. I.: Seismic anisotropy of shales, J. Geophys. Res., 100, 5991–6003, https://doi.org/10.1029/95JB00031, 1995.

Kendall, J. M., Stuart, G. W., Ebinger, C. J., Bastow, I. D., and Keir, D.:. Magma-assisted rifting in Ethiopia, Nature, 433, 146–148, https://doi.org/10.1038/nature03161, 2005.

Kendall, J. M., Pilidou, S., Keir, D., Bastow, I. D., Stuart, G. W., and Ayele, A.: Mantle upwellings, melt migration and the rifting of Africa: insights from seismic anisotropy, Geol. Soc. Lond. Spec. Publ., 259, 55–72, https://doi.org/10.1144/GSL.SP.2006.259.01.06, 2006.

Kraut, E. A.: Advances in the Theory of Anisotropic Elastic Wave Propagation, Rev. Geophys., 1, 401–448, https://doi.org/10.1029/RG001i003p00401, 1963.

Long, M. D. and Silver, P. G.: The Subduction Zone Flow Field from Seismic Anisotropy: A Global View, Science, 319, 315–318, https://doi.org/10.1126/science.1150809, 2008.

Long, M. D. and Becker, T. W.: Mantle dynamics and seismic anisotropy, Earth Planet. Sc. Lett., 297, 341–354, https://doi.org/10.1016/j.epsl.2010.06.036, 2010.

Maultzsch, S., Chapman, M., Liu, E., and Li, X. Y.: Modelling frequency-dependent seismic anisotropy in fluid-saturated rock with aligned fractures: implication of fracture size estimation from anisotropic measurements, Geophys. Prospect., 51, 381–392, https://doi.org/10.1046/j.1365-2478.2003.00386.x, 2003.

McCann, C., Assefa, S., Sothcott, J., McCann, D. M., and Jackson, P. D.: In-situ borehole measurements of compressional and shear wave attenuation in Oxford clay, Sci. Drill., 1, 11–20, 1989.

Meléndez, A.: Development of a New Parallel Code for 3-D Joint Refraction and Reflection Travel-Time Tomography of Wide-Angle Seismic Data – Synthetic and Real Data Applications to the Study of Subduction Zones, PhD thesis, Universitat de Barcelona, Barcelona, Spain, 171 pp., http://hdl.handle.net/2445/65200 (last access: 23 October 2019), 2014.

Meléndez, A., Sallarès, V., Ranero, C. R., and Kormann, J.: Origin of water layer multiple phases with anomalously high amplitude in near-seafloor wide-angle seismic recordings, Geophys. J. Int., 196, 243–252, https://doi.org/10.1093/gji/ggt391, 2014.

Meléndez, A., Korenaga, J., and Miniussi, A.: TOMO3D, Barcelona Center for Subsurface Imaging, available at: http://www.barcelona-csi.cmima.csic.es/software-development, last access: on 23 October 2019, 2015a.

Meléndez, A., Korenaga, J., Sallarès, V., Miniussi, A., and Ranero, C. R.: TOMO3D: 3-D joint refraction and reflection traveltime tomography parallel code for active-source seismic data – synthetic test, Geophys. J. Int., 203, 158–174, https://doi.org/10.1093/gji/ggv292, 2015b.

Meléndez, A., Jiménez-Tejero, C. E., Sallarès, V., and Ranero, C. R.: Synthetic data sets for the testing of the anisotropic version of TOMO3D, Digital.CSIC, https://doi.org/10.20350/digitalCSIC/8957, last access: on 23 October 2019.

Montagner, J. P., Marty, B., and Stutzmann, E.: Mantle upwellings and convective instabilities revealed by seismic tomography and helium isotope geochemistry beneath eastern Africa, Geophys. Res. Lett., 34, L21303, https://doi.org/10.1029/2007GL031098, 2007.

Moser, T. J.: Shortest path calculation of seismic rays, Geophysics, 56, 59–67, https://doi.org/10.1190/1.1442958, 1991.

Moser, T. J., Nolet, G., and Snieder, R.: Ray bending revisited, B. Seismol. Soc. Am., 82, 259–288, 1992.

Nicolas, A. and Christensen, N. I.: Formation of Anisotropy in Upper Mantle Peridotites – A Review, in: Composition, Structure and Dynamics of the Lithosphere-Asthenosphere System, edited by: Fuchs, K. and Froidevaux, C., American Geophysical Union, USA, 111–123, 1987.

Nye, J. F.: Physical Properties of Crystals: their representation by tensors and matrices, 1st Edn., Oxford University Press, Oxford, England, 1957.

Paige, C. C. and Saunders, M. A.: LSQR: An algorithm for sparse linear equations and sparse least squares, ACM Transactions on Mathematical Software (TOMS), 8, 43–71, 1982.

Pratt, R. G. and Chapman, C. H.: Traveltime tomography in anisotropic media – II. Application, Geophys. J. Int., 109, 20–37, https://doi.org/10.1111/j.1365-246X.1992.tb00076.x, 1992.

Prieux, V., Brossier, R., Gholami, Y., Operto, S., Virieux, J., Barkved, O. I., and Kommedal, J. H.: On the footprint of anisotropy on isotropic full waveform inversion: the Valhall case study, Geophys. J. Int., 187, 1495–1515, https://doi.org/10.1111/j.1365-246X.2011.05209.x, 2011.

Rüger, A. and Alkhalifah, T.: Efficient two-dimensional anisotropic ray tracing, in: Seismic Anisotropy, edited by: Erling, F., Rune, M. H., and Jaswant, S. R., Soc. Expl. Geophys., 556–600, https://doi.org/10.1190/1.9781560802693.ch18, 1996.

Sallarès, V., Meléndez, A., Prada, M., Ranero, C. R., McIntosh, K., and Grevemeyer, I.: Overriding plate structure of the Nicaragua convergent margin: Relationship to the seismogenic zone of the 1992 tsunami earthquake, Geochem. Geophy. Geosy., 14, 3436–3461, https://doi.org/10.1002/ggge.20214, 2013.

Sayers, C. M.: Seismic anisotropy of shales, Geophys. Prospect., 53, 667–676, https://doi.org/10.1111/j.1365-2478.2005.00495.x, 2005.

Schulte-Pelkum, V. and Mahan, K. H.: A method for mapping crustal deformation and anisotropy with receiver functions and first results from USArray. Earth Planet. Sc. Lett., 402, 221–233, https://doi.org/10.1016/j.epsl.2014.01.050, 2014.

Song, L., Zhang, S., Liu, H., Chun, S., and Song, Z.: A formalism for acoustical traveltime tomography in heterogenous anisotropic media, in: Review of Progress in Quantitative Nondestructive Evaluation, edited by: Thompson, D. O. and Chimenti, D. E., Springer, Boston, USA, 1529–1536, https://doi.org/10.1007/978-1-4615-5339-7_198, 1998.

Stewart, R. R.: An algebraic reconstruction technique for weakly anisotropic velocity, Geophysics, 53, 1613–1615, https://doi.org/10.1190/1.1442445, 1988.

Thomsen, L.: Weak elastic anisotropy, Geophysics, 51, 1954–1966, https://doi.org/10.1190/1.1442051, 1986.

Thomsen, L. and Anderson, D. L.: Weak elastic anisotropy in global seismology, Geol. Soc. Spec. Pap., 514, 39–50, https://doi.org/10.1130/2015.2514(04), 2015.

Tsvankin, I.: P-wave signatures and notation for transversely isotropic media: An overview, Geophysics, 61, 467–483, https://doi.org/10.1190/1.1443974, 1996.

Vauchez, A., Tommasi, A., Barruol, G., and Maumus, J.: Upper mantle deformation and seismic anisotropy in continental rifts, Phys. Chem. Earth A, 25, 111–117, https://doi.org/10.1016/S1464-1895(00)00019-3, 2000.

Vinnik, L. P., Makeyeva, L. I., Milev, A., and Usenko, A. Y.: Global patterns of azimuthal anisotropy and deformations in the continental mantle, Geophys. J. Int., 111, 433–447, https://doi.org/10.1111/j.1365-246X.1992.tb02102.x, 1992.

Warner, M., Ratcliffe, A., Nangoo, T., Morgan, J., Umpleby, A., Shah, N., Vinje, V., Štekl, I., Guasch, L., Win, C., Conroy, G., and Bertrand, A.: Anisotropic 3D full-waveform inversion, Geophysics, 78, R59–R80, https://doi.org/10.1190/geo2012-0338.1, 2013.

Yousef B. M. and Angus, D. A.: When do fractured media become seismically anisotropic? Some implications on quantifying fracture properties, Earth Planet. Sc. Lett., 444, 150–159, https://doi.org/10.1016/j.epsl.2016.03.040, 2016.

Zhou, B. and Greenhalgh, S.: “Shortest path” ray tracing for most general 2D/3D anisotropic media, J. Geophys. Eng., 2, 54–63, https://doi.org/10.1088/1742-2132/2/1/008, 2005.

Zhou, B. and Greenhalgh, S.: Non-linear traveltime inversion for 3-D seismic tomography in strongly anisotropic media, Geophys. J. Int., 172, 383–394, https://doi.org/10.1111/j.1365-246X.2007.03649.x, 2008.