# Seismological Research Letters

- © Seismological Society of America

## ABSTRACT

Shear‐wave splitting analysis is commonly used to infer deep anisotropic structure. For simple cases, obtained delay times and fast‐axis orientations are averaged from reliable results to define anisotropy beneath recording seismic stations. However, splitting parameters show systematic variations with back azimuth in the presence of complex anisotropy and cannot be represented by average time delay and fast‐axis orientation. Previous researchers had identified anisotropic complexities at different tectonic settings and applied various approaches to model them. Most commonly, such complexities are modeled using multiple anisotropic layers with priori constraints from geologic data. In this study, a graphical user interface called M‐Split is developed to easily process and model multilayered anisotropy with capabilities to properly address the inherited nonuniqueness. M‐Split program runs user‐defined grid searches through the model parameter space for two‐layered anisotropy using formulation of Silver and Savage (1994) and creates sensitivity contour plots to locate the local maxima and analyze all possible models with parameter trade‐offs. To minimize model ambiguity and identify the robust model parameters, various misfit calculation procedures are also developed and embedded in M‐Split and can be used depending on the quality of the observations and their back azimuthal coverage. M‐Split is an open‐source program and can be extended by users for additional capabilities or for other applications.

## INTRODUCTION

Shear‐wave splitting analysis of teleseismic data is widely used to image the anisotropy in the Earth’s interior. Splitting observations are used to infer crust and mantle fabric beneath recording stations related to past and current deformations. Normally, shear‐wave splitting analyses are performed under the assumption of a simple model composed of just one horizontal layer with a horizontal axis of symmetry. But in some complex cases, the splitting parameters show strong back‐azimuthal dependency. There are three different interpretations for variation of splitting parameters with respect to back azimuth: (1) two‐layered horizontal anisotropy (Savage and Silver, 1993; Silver and Savage, 1994) or three or more anisotropic layers (e.g., Yang *et al.*, 2014) can generate back‐azimuthal dependencies with 90° of periodicity, (2) laterally varying anisotropy (Alsina and Snieder, 1995), and (3) inclined symmetry axis in the lithosphere (Babuška *et al.*, 1984, 1993; Plomerová *et al.*, 1996; Šílený and Plomerová, 1996).

Through time, numerous codes are developed to model splitting observations using simple one‐layered anisotropy. Three techniques are commonly used to estimate the splitting parameters (i.e., fast polarization direction and delay time). These techniques are minimum energy (SC; Silver and Chan, 1991), rotation correlation (RC; Bowman and Ando, 1987), and eigenvalue (EU; Silver and Chan, 1991). Recently, Wüstefeld *et al.* (2007) developed the SplitLab program to implement all three techniques along with a workflow setup to identify null measurement. An interactive forward‐modeling tool that has preliminary insight about the two‐layered anisotropy is also set in the SplitLab program. Unfortunately, this tool has neither proper grid search nor sensitivity analysis capabilities and thus is not commonly used. Until now, forward codes were developed and used in individual cases and there was no comprehensive tool covering all those methods and also offering tools to analyze multilayered anisotropy. We introduce our user friendly graphical program M‐Split that uses the same forward formulas introduced by Silver and Savage (1994) and enhanced misfit calculation algorithms for different observation conditions and limitations. Analysis of complex anisotropy suffers from the presence of high level of nonuniqueness in the forward problem; hence, in this program a sensitivity analysis tool is also included to find all possible solutions in model space by plotting all local extremums.

## METHODOLOGY

M‐Split performs a grid search over the model parameter space that contains in this case fast polarization directions and delay times of two horizontal layers. It is also worth to mention that in grid‐search procedure; there is a trade‐off between the accuracy of the final output model and processing time. If a very fine grid search is carried out, the number of models that will be tested and therefore the time of the calculations will be highly increased but the final result will be more accurate. After performing the grid search, formulas introduced by Silver and Savage (1994) are used to calculate the theoretical splitting parameters for each set of model parameters. Finally, the theoretical/calculated splitting parameters are compared with the observations and the model with the most resemblance to the observations is extracted. To compare the observations with calculated parameters, a misfit function calculated from residuals similar to Fontaine *et al.* (2007) is utilized, and the model with the smallest residual is extracted as the best‐fitting solution.

To test two horizontal layered anisotropy assumptions properly, the same procedure, including the grid search, is also performed for one‐layered models. This enabled us to assess how much better a two‐layered model fits the observations compared with the one‐layered simple model, and it is quantified by calculating the coefficient of determination *R*^{2} (Walker *et al.*, 2004) and the (Walker *et al.*, 2005). In M‐Split, various misfit calculation algorithms are developed and utilized which, depending on the quality and back‐azimuthal coverage of the observations, can be used for better results. The first algorithm is the unweighted method that gives the same weight to all observations. In other words, all observations will have the same effect in the misfit calculation and therefore have the same role in the final model extraction. If there is no measurement representing the observation quality or if the calculated errors are not reliable, this method will be appropriate. As shown in equation (1), in misfit calculation, both phi and *dt* segments are normalized by the corresponding mean value of all observations: (1)in which , , *δt*_{obs,i} and *δt*_{calc,i} are *i*th observed and calculated fast polarization direction and observed and calculated delay times, respectively, and and mean_{δt} are mathematical mean values of fast polarization directions and delay times for all observations, respectively.

Second algorithm is the weighted method that gives weights according to measurement errors so that observations with small errors will have greater effect in the misfit calculation and therefore in the extracted final model. This method is proper for the observations that have solitary points with small errors that are a critical part of the model diagram and therefore it is important that these points have more effect on the final calculated model (equation 2). According to the equation, each part is normalized by its own error. The errors are half of the difference of minimum and maximum values of that parameter (i.e., −errPhiXX and +errPhiXX and −errdtXX and +errdtXX, in which XX stands for the method used for splitting parameter estimation) for that observation; in other words, they are the average of differences between a given value from its minimum and maximum (see equation 2): (2)in which , , *δt*_{obs,i} and *δt*_{calc,i} are *i*th observed and calculated fast polarization direction and observed and calculated delay times, respectively. and *E*_{δt,i} are errors on fast polarization direction and delay time for *i*th observation.

Finally, the third method is a Band‐Fit method that is developed to reduce the effect of individual observations with small errors that dominate the final model. This method is proper for the input with the individual observations that have negative effect on the extracted model. As declared in equation (3), a difference function has been defined so that if the difference between observed and calculated values is smaller than the error range, the difference function will be set to zero; otherwise, it will be the difference between the difference function and error function: (3)in which , , *δt*_{obs,i} and *δt*_{calc,i} are *i*th observed and calculated fast polarization direction and observed and calculated delay times, respectively. and *E*_{δt,i} are errors on fast polarization direction and delay time for *i*th observation. and mean_{δt} are mean values of fast polarization direction and delay time for all observations, respectively.

There are also two additional options to reduce the effects of possible back‐azimuthal errors. Normally, there is no back‐azimuthal correction, but considering the possibility of back‐azimuthal error contained in observations, using tolerance on back azimuth can lead to better solutions. Also, if number of observations in one back‐azimuthal range compared with others is too high, then back‐azimuthal normalization can improve the final result.

In addition, the residual calculation algorithm proposed by Kaviani *et al.* (2013) is also included in the program. Unlike others, this algorithm is based on calculation of the penalty function and provides an alternative approach for misfit assessment. It is worthwhile to note that implementing such a variety into M‐Split gives us the opportunity to compare results from various algorithms, especially in the sensitivity analysis stage.

## PROGRAM MODULES

The M‐Split workflow is shown in Figure 1. It can be seen that this program is composed of three major steps: data preparation, grid search, and sensitivity analysis. In this program, there are a few points that should be considered. First, the main assumption is that input data have 90° of back‐azimuthal periodicity; hence, this program plots the data in mod 90 form. Therefore, all back azimuths are initially converted to mod 90. Second, applied grid search and misfit calculations for simple anisotropy with a one‐layered model are exactly the same as two‐layered models. Thus, quantitative comparison of best‐fitting model misfits for one‐ and two‐layered models aids users to decide about the number of anisotropic layers that should be used.

### Data Conversion

As shown in Figure 2a, the start window of the program consists of two options: with and without data conversion. If data conversion is needed, then the conversion window will open (Fig. 2b). Input of this program is compatible with the output of SplitLab, which is in tab‐delimited text format. M‐Split generates two sets of data formatted for different types of algorithms that the user may use within the program. The program is adapted for three different shear‐wave splitting analysis techniques, which are SC method (Silver and Chan, 1991), RC method (Bowman and Ando, 1987), and eigenvector method (Silver and Chan, 1991), and hence any of these methods can be used to generate the input files.

As mentioned before, the null observations can be detected using SplitLab and can be included in the analysis in M‐Split without weighting because the large error ranges given for null observations will lead to minimal weights in weighted and band‐fit approaches, causing null data to have no effect on the final model. Therefore, the input of weighted and band‐fit algorithms are arranged so that they cannot contain the null observations. It is also important to note that including null observations might lead to systematic error and deviate the final result from true model parameters. As an example, in the RC method, null observations are always 45° off from the correct direction (Wüstefeld and Bokelmann, 2006); thus, including null data will lead to wrong solutions. On the other hand, if SC (Silver and Chan, 1991) method is selected, null directions are expected to be parallel or perpendicular to fast polarization direction, so ones that are perpendicular will still contaminate the input data. In this respect, including null observations in the analyses is not recommended by the authors, but is not taken out entirely from the program for testing purposes. Therefore, including null data should be done with caution.

### Grid Search

Having data prepared in proper format, the next step is the comprehensive grid search to identify the best‐fitting model. As shown in Figure 3, the inversion window consists of four steps. The first step is the data quality selection where users will decide the quality of the observations to be used in the misfit calculation and model extraction. It must be noted that in most of the cases, the number of good quality data for a reliable estimation is not enough, therefore using the fair quality observations will be necessary.

The next step is back‐azimuthal correction options. Users can add a tolerance to the back azimuth so that program will calculate the misfit for the back‐azimuth range that can be defined as back azimuth plus–minus the tolerance with the increments of step size. Then, the minimum of the calculated misfits for the selected tolerance range will be considered as the misfit. It is necessary to mention that adding tolerance will drastically increase the process time because instead of a single point, misfit will be calculated for a vector. As another option in this part, back‐azimuthal normalization can also be selected. Using this option, the observations will be divided into nine clusters of 10° of back‐azimuthal range and each cluster will be normalized by the number of observations in that cluster. During adjustment of grid‐search parameters, the maximum and minimum values and their increments of the splitting parameters for both layers will be set. Like the tolerance, if increments are set to a small value, the processing time will increase because of the increase in number of the models that should be tested.

The final step is the misfit calculation selection, where users can choose from one of the formerly mentioned methods to be used as the misfit calculation formula. The advantage of this program is that all models with corresponding calculated parameters (i.e., misfit, *R*^{2}, etc.) will be written into a text file for future use. All models with positive *R*^{2} values will be written in another text file to use as an input file for the sensitivity analysis program. The final model will be displayed on the results displaying screen (Fig. 4). In some cases, two‐layered anisotropy parameters resulting from different algorithms can be very different but their resultant splitting can be very similar. Therefore, to visualize and compare the competing models, a model plotting console that plots splitting parameters of multiple models along with the observations is established (Fig. 5).

### Sensitivity Analysis

Nonuniqueness of the resultant models is common in such grid‐search applications involving multiple variables. In many cases, the quantified misfits between sets of models can be very close, outlining the inherited nonuniqueness in identification of complex anisotropy and vitality of sensitivity analysis. To handle this issue, a sensitivity analysis program has been developed (Fig. 6). For this purpose, the program uses all models with positive *R*^{2} values and generates contour plots to identify present local maxima (or minima) in the model parameters space for various misfit functions. The best‐fitting model parameters associated with each local maximum are also plotted with their standard deviations. Eventually, if misfits calculated for multiple local maxima are similar, the data can be represented equally well with multiple models and thus other independent data from geology, geodesy, tectonics, etc., is needed to choose among these models.

## SYNTHETIC TEST

To evaluate the program, formulas for apparent splitting parameters introduced by Silver and Savage (1994) are used to generate a discrete synthetic data corresponding to defined model parameters. The dominant frequency of 0.1 Hz is used during the data generation phase and also in real and synthetic data tests. Consequently, 40 points are randomly selected and used as synthetic observations to test the program. In the first step, the synthetic data are used without noise and all three algorithms recovered the parameters correctly. In the second step, Gaussian white noise is included to the synthetic data. At this stage, three signal‐to‐noise ratios are used (i.e., 5, 10, and 20). The error bonds for synthetic data are assigned by multiplying the differences between noiseless and noisy data with 0.5, 1, and 1.5, which produced nine different synthetic data sets to test the program.

During synthetic tests using the unweighted method with high noise rates or low signal‐to‐noise ratios (i.e., 5 and 10), the recovered model showed noticeable deviations from input model (Table 1; Fig. 7). In contrast, both band‐fit and weighted methods performed well even for noisy data (Table 1; Fig. 5b) and in most cases led to the correct model.

## ANALYSIS OF REAL DATA

During the final step of the testing process, splitting measurements made for station BKS by Bonnin *et al.* (2010) are used. This data set is ideal for this purpose because it has a good back‐azimuthal coverage and had already been studied for complex anisotropy. Grid searches conducted using weighted and band‐fit approaches led to results similar to Bonnin *et al.* (2010). The used data display highly variable error bonds, which are ignored in the unweighted approach, causing noticeable deviations from the published model. In this respect, unweighting is not optimal for high‐quality data sets, including well‐defined error bonds.

In this data set, penalty function algorithm also selected a best‐fitting model different than the minimum misfit, *R*^{2}, and *R*^{2} adjusted, which are intrinsically linked (Fig. 8). During sensitivity analysis, contour plots of both given in Figure 9 reveal that for all parameters there are two model clusters characterized by equally low misfits. The first cluster corresponds to model parameters given in Bonnin *et al.* (2010) and the second cluster is the one picked by penalty function. The misfit distribution within model parameter space is nearly identical, indicating that the data can be reproduced by both models. At this stage, nonuniqueness may be avoided only by incorporating data from independent sources such as assuming known plate motion to be parallel/subparallel to the upper layer fast‐axis direction.

## CONCLUSION

We presented a program that analyzes shear‐wave splitting observations to find the best‐fitting multilayered anisotropy. M‐Split can perform data preparation, model fitting, and sensitivity analysis for given data input. The program’s interface provides a user‐friendly environment that automatically performs the calculations and produces visual outputs. To reduce the model ambiguity and test robustness, multiple misfit calculation procedures are embedded to the program and each has been subjected to rigorous testing. After successful tests conducted on both synthetic and real data, the reliability of the designed software is established. M‐Split runs on MATLAB platform (see Data and Resources) independent of operating system and successfully tested on Windows‐ and Linux‐based computers.

## DATA AND RESOURCES

All data used in this article came from published sources listed in the References. The source codes of the program can be obtained from www.m-split-code.com (last accessed May 2017). All programs and functions in this software pack are developed using MATLAB software (https://www.mathworks.com/products/matlab.html?s_tid=hp_products_matlab, last accessed May 2017). To accomplish this, a commercial version of this program (R2014a) is used.

## ACKNOWLEDGMENTS

The authors thank Mickael Bonnin for providing the data that are used in the tests and validation process. The authors are also grateful to Eric Sandvol, Jonathan Robert Delph, Katrina Burch, and Colton Lynn for their helpful comments.