Matlab Seismic Anisotropy Toolkit - using the functions

The first group of functions generate elasticity matrices from input arguments (see the table below). MS_load and MS_load_list read data from input files while MS_elasticDB extracts elasticity information from a built-in database. Other functions create elasticity matrices from summary information and a knowledge of the symmetry.

In the most general case 21 elastic constants are required to represent minerals and materials with a triclinic symmetry, but systems with a higher symmetry (e.g. orthorhombic, hexagonal or cubic) require fewer elastic constants. In the simplest case, that of isotropic materials, this reduces to two parameters. The limitations of what can usually be resolved using seismology mean that higher symmetry classes are often assumed. A class which is commonly invoked is transverse isotropy (TI, also known as radial or polar anisotropy). In a TI medium properties only vary with angle from the symmetry axis. This can arise in a number of ways and is exhibited by: crystalline materials with a hexagonal symmetry, isotropic thin layers of different stiffness, aligned inclusions or cracks which are uniformly oriented. TI requires five independent constants to parameterise. There are a number of commonly used choices of parameter sets (depending on the area of literature); MSAT has functions to work with several of these.

It is also possible to model the effective anisotropic elasticity caused by particular rock fabrics; an effect that can be important even if all the components of the rock are isotropic on a short length scale compared to the seismic wavelength. Several effective medium models are supported by the MS_effective_medum function. These include the case of an elastically isotropic medium contains aligned ellipsoidal inclusions filled with an elastically isotropic material of a different stiffness, the case of layering of isotropic materials of different stiffness, and the effect of aligned cracks.

A second group of functions include elasticity matrices as input and output. These functions, listed below, are used to manipulate or transform an elasticity matrix.

One important manipulation is tensor rotation. MS_rot3, MS_rotEuler and MS_rotR all rotate an elasticity matrix (the functions differ in the way the rotation is specified: in all cases a rotation matrix is constructed and MS_rotR is used to perform the actual manipulation). To understand the approach recall that a vector, $\mathbf{v}$ is rotated from to a new orientation, $\mathbf{v}^{(R)}$, by multiplication by a rotation matrix: $\mathbf{v}^{(R)}=\mathbf{g}\mathbf{v}$, where the rotation matrix, $\mathbf{g}$, is a $3 \times 3$ orthogonal matrix with determinant 1 and elements that represent the cosines of angles of the rotation. This transformation be also be written element-wise as: $v_i^{(R)}$ = $g_{ij}v_j$ (with an implicit summation over $j=1,3$). This formula can be extended for higher order Cartesian tensors. For example, the second order stress tensor is rotated: $\sigma_{ij}^{(R)}$ = $g_{ik}g_{jl}\sigma_{kl}$ (with two implicit summations) and the elastic stiffness tensor rotated: $C_{ijkl}^{(R)}$ = $g_{im}g_{jn}g_{kp}g_{lq}C_{mnpq}$ (with four summations and avoiding Voigt notation). For second order tensors the transformation can also be written in the form of matrix multiplication by the rotation matrix and its transpose: $\mathbf{\sigma}^{(R)} = \mathbf{g} \mathbf{\sigma} \mathbf{g}^{T}$, but this shortcut is not available for higher order tensors. However, it is possible (in fact, highly desirable for performance reasons) to perform rotations on elastic constants in Voigt form without requiring to transform to the $3\times3\times3\times3$ tensor form. In order to do this we use the basis change formula from Bower (2009; Section 3.2.11):

$ \mathrm{C}^{(R)}=\mathrm{K}\,\mathrm{C}\,\mathrm{K}^T $

where $\mathrm{C}^{(R)}$ is the rotated version of $\mathrm{C}$. In this equation $\mathrm{K}$ is a $6\times6$ matrix derived from combinations of the normal is its matrix transpose.

MS_axes also rotates an input matrix but by an amount designed to allow the anisotropic nature of the matrix to be analysed. Such analysis makes use of two other functions, MS_decomp and MS_norms. Taken together, these three functions allow the symmetry of elasticity matrices to be expressed following the approach of Browaeys and Chevrot (2004). Example of this approach are provided elsewhere in the documentation.

Another type of manipulation concerns the averaging of elastic constants. In many fields we do not deal with single crystals. We need, therefore, a method of estimating the aggregate elastic properties of a polycrystalline material containing a range of crystal types and orientations. While much more rigorous treatments exist for many purposes simple averaging schemes suffice. These either assume constant strain (Voigt averaging) or stress (Reuss averaging) throughout the sample, and essentially amount to taking the arithmetic average of each element of the stiffness or compliance tensor, respectively. These two values place an upper (Voigt) and lower (Reuss) bound on the true value. It was observed by Hill (1952) that experimental values were often close to the average of these two bounds, a value which has become known as the Voigt-Reuss-Hill (VRH) average:

$ C^{\mathrm{VRH}}=\frac{1}{2}\left(C^{\mathrm{Voigt}}+C^{\mathrm{Reuss}}\right) $

where

$ C^{\mathrm{Voigt}}=\sum_{i=1}^{N}{v_iC_i}\;\;\; \mathrm{and}\;\;\;C^{\mathrm{Reuss}}=\left(\sum_{i=1}^{N}{v_iS_i}\right)^{-1}. $

Here, $v_i$ are the volume fractions of the $N$ individual components, with associated elasticities ($C_i$) and compliances ($S_i$). This is very commonly employed. It should be noted that it is an empirical relation with no strict theoretical foundation, albeit a useful one. MS_VRH calculates the Voigt, Reuss, and Voigt-Reuss-Hill average elasticity of an aggregate of materials, each described by their own elasticity matrix and volume fraction, while MS_polyaverage takes this averaging one step further and estimates the isotopic elasticity matrix if the components are randomly oriented.

The third group of functions allow elasticity to be analysed or displayed. They take elasticity matrices as input arguments and either return derived quantities or cause graphics to be produced. Examples include MS_anisotropy, which implements some of the attempts to capture the degree of elastic anisotropy of a material as a single value and MS_poisson which calculates Poisson’s ratio for a general anisotropic material as a function of direction.

For seismology, one of the most important ways to analyse the anisotropic elasticity of a medium is in terms of seismic phase velocities. The propagation of motion through an elastic medium is governed by the elastodynamics equation:

$ \rho\left(\frac{\partial^2u_i}{\partial t^2}\right) = C_{ijkl}\left(\frac{\partial^2u_l}{\partial x_jx_k}\right) $

where $\rho$ is the density, $u_i$ is the displacement, $x_j$ is position and $t$ is time. By substituting into this the displacement associated with a monochromatic plane wave as a function of time, we obtain the well-known Christoffel equation:

$ \left(C_{ijkl}n_jn_l-\rho V^2\delta_{ik}\right)p_{k}=0 $

where $n_j$ is the propagation unit vector, $V$ are the phase velocities, $p_k$ are the polarisation unit vectors and $\delta_{ik}$ is the Kronecker delta function. Solutions to this equation for a given direction ($n_j$) yield the three orthogonal waves with different phase velocities. One of these will be close to the wavefront normal (designated quasi-P or qP) and two will be near perpendicular (quasi-S1, quasi-S2; qS1 and qS2). In practice the quasi- prefix is often omitted. The reader is directed to Mainprice (2007) for a more complete treatment of the subject. MSAT follows the methodology of Mainprice (1990) to calculate phase velocities and this is implemented in MS_phasevels. The function accepts a (set of) propagation direction(s) and returns the three seismic phase velocities along with their polarisations. This function is used by MS_plot and MS_sphere, which plot the resulting seismic anisotropy on a pole figures and a three-dimensional spherical representation, respectively.

Alongside the functions designed to analyse aspects of elasticity, MSAT also contains functions to deal with simple aspects of modelling one of the most common seismological observations related to elastic anisotropy, shear wave splitting. As discussed above, seismic waves can propagate at three velocities through an anisotropy body and two of these velocities correspond to shear waves. If a shear wave passes through an anisotropic region it is split into two shear waves with particle motion polarised at 90 degrees to each other (and approximately perpendicular to the direction of wave propagation. As the two shear waves have different velocities they will emerge from the anisotropic region at different times and with particle motion direction imposed by the elastic anisotropy of the medium. This process is known as shear wave splitting (SWS) and its observation in seismic data is considered strong evidence for anisotropy somewhere along the ray path.

In MSAT we only implement functions to handle the simplest aspects of SWS. Time series of particle displacements are described by three arrays we call time, T00 and T90. T00(i) and T90(i) relate to the particle position of two perpendicular traces at time time(i). The function MS_make_trace generates these assuming a first derivative Gaussian shape for a given dominant frequency making space in the arrays for a user selected maximum amount of splitting. MS_phasevels can calculate splitting parameters from an elasticity matrix and this can be used with MS_split_trace to apply splitting parameters. The effective splitting generated by a wave passing through multiple layers can be found by combing splitting operators using MS_effective_spltting_N. This either implements an analytic approach to calculating the resulting splitting parameters or uses MS_split_trace multiple times and then makes a measurement of the resulting effective splitting using MS_measure_trace_splitting. MS_splitting_misfit allows two splitting operators to be compared (allowing the best fitting elasticity to be found, see the Splitting misfit example) and MS_plot_trace can be used to visualise time, T00 and T90. For complex applications these three arrays can be converted to and from SAC records using the functions available from MSAC.

Finally, range of miscellaneous and utility functions are provided to, for example, convert between different ways of representing elasticity or check that a particular Matlab array is a valid elasticity matrix.