- © Seismological Society of America
We present a graphical user interface in the MATLAB environment for teleseismic shear‐wave splitting analysis. In view of increasingly large seismological experiments, SplitRacer is designed to process large data sets quickly for which it integrates automatic and user‐active steps. This enables a straightforward analysis based on the user’s own quality criteria. SplitRacer offers a number of features designed to ensure precise measurements as well as user‐friendliness. Data can be downloaded directly from all data centers that are part of the International Federation of Digital Seismograph Networks Web Services. Existing data can easily be used when cut into 1‐hr three‐component seismograms. A first automatic data check is followed by visual assessment of the phases in question. We account for a possible misalignment of the sensor by analyzing the particle motions and comparison with the theoretical back azimuths of the events recorded at a station. The shear‐wave splitting analysis is based on the minimum (transverse) energy method (Silver and Chan, 1991), which constrains the fast‐polarization direction and the delay time. Measurements can be repeated for a number of time windows to enable a statistical error analysis. SplitRacer provides an overview of single‐phase splitting measurements per station that can be used to detect azimuthal dependencies. To infer the anisotropic properties beneath a station, SplitRacer offers the possibility to fit models of one or two anisotropic layers or a model of smoothly varying anisotropy to the observed data. Furthermore, SplitRacer also offers a multiphase joint‐splitting analysis for one and two hypothetical layers in which all events at a given station are used simultaneously to constrain the anisotropic model. This significantly reduces the influence of noise and makes the measurement more robust. As an example, we apply SplitRacer to constrain seismic anisotropy beneath the permanent station BKS of the Berkeley Digital Seismic Network.
The analysis of shear‐wave splitting has been used to constrain the anisotropic properties of the Earth’s interior for more than three decades now. As such, it belongs to the standard analysis tools for many seismic experiments and is applied worldwide, continuously improving our understanding of the dynamic processes of the Earth’s interior (e.g., Silver and Chan, 1991; Savage, 1999; Long and Silver, 2009).
Anisotropy is thought to be controlled by two different mechanisms: crystallographic‐preferred orientation (CPO) of minerals such as olivine, which forms by, for example, simple shear through dislocation creep induced by mantle flow (e.g., Savage, 1999; Long and Silver, 2009) and shape‐preferred orientation caused by aligned cracks or melt‐filled lenses (e.g., Crampin et al., 1984; Holtzman and Kendall, 2010). If a shear wave travels through an anisotropic medium, it is split into two orthogonally polarized fast‐ and slow‐wave components. The polarization of the fast‐wave component is used to infer the fast‐axis direction of the anisotropic layer; its thickness and anisotropic strength depend on the delay time between the two wave components.
Previous shear‐wave splitting studies in conjunction with laboratory experiments have been used to derive anisotropic properties at crustal, midlithospheric, asthenospheric depths, as well as the core–mantle boundary region (e.g., Savage 1999; Wookey et al., 2005; Long and Silver, 2009; Ko and Jung, 2015). In ocean settings, shear‐wave splitting is mostly interpreted in terms of CPO in relation to flow approximately aligned with absolute plate motion directions and/or as frozen fossil spreading directions preserved in the lithosphere (e.g., Savage, 1999; Long and Silver, 2009). In continental settings, shear‐wave splitting is often more complicated to interpret since anisotropy might be caused by active fault activity, frozen‐in fossil anisotropy at lithospheric depth characteristic of past and present orogenic events, as well as mantle flow at asthenospheric depth (e.g., Silver and Chan, 1991; Fouch and Rondenay, 2006; Flesch and Bendick, 2012).
To investigate the anisotropic properties of the different regions beneath a station, the shear‐wave splitting analysis is usually based on teleseismic phases. This includes all core shear phases from events between 85° and 180° distance to the receiver. Most often, SKS, SKKS, and PKS phases are used, and we refer to them as XKS phases. As these phases travel through the liquid outer core, their conversion from P to S at the core–mantle boundary leads to a polarization in the direction of back azimuth. A single teleseismic shear‐wave splitting measurement is always a path‐integrated measurement (Long and Silver, 2009), characteristic of all anisotropy along the receiver side of its path. Thus, it has a poor depth but good lateral resolution, in which the latter can be obtained from Fresnel zone estimates (e.g., Rümpker and Ryberg, 2000).
A number of software codes have been published to perform shear‐wave splitting measurements on selected XKS phases. These mostly rely on choosing a certain time window for the analysis and often use the technique suggested by Silver and Chan (1991), which is a grid search for the splitting parameters, the fast polarization , and delay time dt. Although some codes focus on individual measurements (single event at a single station; Silver and Chan, 1991; Wüstefeld et al., 2008), others favor an automatic algorithm (e.g., Teanby et al., 2004; Savage et al., 2010; Wüstefeld et al., 2010; Liu and Gao, 2013) for the selection of a suitable time window for the analysis. This decreases the manual input and has the potential to objectify the categorization of measurements (Savage et al., 2010). Automatic processing routines have been advanced further by automatic selection of suitable events (e.g., Walther et al., 2014), which lead, due to the rigorous quality criteria, often to the discharging of many events (Evans et al., 2006; Wüstefeld et al., 2010).
In this study, we present a semiautomatic processing approach to select usable events and time windows that seems most suitable in view of the increasing number of permanent seismic stations as well as large temporary seismic experiments. SplitRacer aims to automatize repetitive steps but requires the user’s manual input to ascertain that the maximum amount of data is used. To ensure user‐friendliness, as in the perhaps most widely used code in teleseismic shear‐wave splitting analysis SplitLab (Wüstefeld et al., 2008), SplitRacer is a graphical user interface (GUI) in the commonly used MATLAB environment (see Data and Resources). This may be especially beneficial for users unfamiliar with the method and processing. Advanced users will benefit from the possibilities in tuning the pre‐processing and analysis criteria to their needs. All processing steps can be repeated, different settings applied, and compared. Data can be added, and processing be halted or continued. The interpretation of single‐splitting results is strengthened by the theoretical layer fitting of single‐splitting measurements and multiphase joint‐splitting analysis for one and two layers. To illustrate SplitRacer’s workflow, we apply it to the permanent station BKS of the Berkeley Digital Seismic Network (BDSN) and compare our findings with those of previous publications.
PROGRAM CONTENT AND STRUCTURE
SplitRacer provides a stand‐alone GUI and analysis system in the MATLAB environment. It has been tested for MATLAB 2015b‐2016b. For full functionality, it requires a Linux or Windows operating system and the GNU software Wget (see Data and Resources). The program is available from the authors’ institutional website (see Data and Resources), in which we also provide a data example as discussed in this article. SplitRacer is an open source software that can be modified by the user for different purposes such as teleseismic or local S shear‐wave splitting. Please consult the user guide for instructions. SplitRacer consists of five subfolders that store the subroutines for its three main functions: download, pre‐processing, and analyzing teleseismic data for shear‐wave splitting. The defaults folder comprises all necessary nonscript files such as the text files for data centers and event lists, travel time tables, and the user guide. The folder main contains the main GUI functions.
SplitRacer builds comprehensive data libraries for each processing step while the user works with the program. This ensures the possibility of repeating processing steps, adding data at a later point, or processing the same data set with different settings and comparing results. Figure 1 shows SplitRacer’s workflow. First, data must be downloaded or imported. Then, the pre‐processing routines help to select usable phases suitable for further analysis. Finally, the shear‐wave splitting analysis is applied. The reader is referred to the user guide for specific instructions and complete descriptions of available program options. Statistics, results, and waveforms can be exported at intermediate steps. To start the program, the user must simply run the MATLAB function start_SplitRacer.m and follow the instructions in the menu. Figure 2 shows the initial interface and main layout.
Given recent development within the seismological community, the program package offers a data download per International Federation of Digital Seismograph Networks Web Services (FDSNWS), which offers data from worldwide networks. The user can choose between different data centers offering FDSNWS, their networks and stations, as well as start and end time, distance range, and threshold magnitude. Data centers can be added in the appropriate text file, data_centers.dat. Once a data center is selected, a request is sent and available networks are shown in the GUI. This also applies for choosing a channel name for broad band (BH) and high broad band (HH) streams, which promptly lists the available stations in its own menu. SplitRacer automatically downloads 1‐hr three‐component miniSEED (mseed) files based on an event list containing earthquakes above a given magnitude worldwide (e.g., default is magnitude 6). The event list is stored in defaults, and the user can alter it or add a new one. The reader is referred to the user guide for how to use existing data.
The pre‐processing of the data aims to distinguish between phases that can and cannot be used for the shear‐wave splitting analysis. We designed several subroutines in which the user controls the criteria to select suitable phases efficiently. SplitRacer’s pre‐processing steps consist of three interdependent subroutines: initial pre‐processing, visual quality check, and check station misalignment (see Fig. 1). Criteria for the automatic pre‐processing routines that are applied during the initial pre‐processing are set by the user, and their impact can be viewed beforehand via the view waveforms option. Events and phases for the analysis are selected based on a chosen band‐pass filter to enhance the signal‐to‐noise ratio (SNR), as well as an SNR‐cutoff value below which events (phases) are discarded. During the visual quality check, only events (phases) that fulfilled the previous criteria are displayed to be checked upon by the user. Finally, a possible misorientation of the sensor can be analyzed using the check station misalignment option.
When analyzing a new data set, it is important to learn what kind of quality the data actually has, before choosing the parameters for the automatic processing, which might otherwise lead to an unfavorable diminishing of the data set. Using the view waveforms option, SplitRacer offers a seismogram and XKS phase viewer as illustrated in Figure 3 using an example from station BKS. The user chooses an event, upon which the vertical, north, east (ZNE) components of the seismogram are plotted in full length. Next to each trace, spectral amplitudes are also shown. To the bottom, all XKS phases are displayed separated by tabs. Phases are marked using a travel‐time table look‐up originally generated with TauP for the iaspei91 model (Crotwell et al., 1999). Each tab shows a 100 s zoom‐in of the ZNE seismogram around the phase, as well as the particle motion of the northeast (NE) components for the current filter settings and a longer period filter of 15–50 s. In the center, a button panel composed of the same options as in the pre‐processing GUI is shown. The information box to the bottom right displays which criteria are currently applied and if the phase is selected or discarded using these settings. This enables the user to tune the pre‐processing to the data set by automatically selecting only those events (phases) that will be used for further analysis. As default, we suggest applying a broadband filter between 4 and 50 s, which is usually appropriate for the analysis of XKS phases.
Once automatic pre‐processing criteria are chosen, the entire data set is processed. During initial pre‐processing, SplitRacer applies a 1 s low pass and resamples to 20 Hz upon reading the mseed files. Only phases that are not preceded by another XKS phase within a 10‐s time window are allowed for further analysis to avoid interference of neighboring phases. For each phase, SplitRacer applies the chosen broadband filter and calculates an SNR based on a 20‐s‐long noise time window ending 5 s prior to the theoretical XKS arrival of the phase. The signal time window of 25 s starts 5 s before the theoretical XKS arrival to account for variations of origin time or focal location estimation as suggested by Liu and Gao (2013). We use horizontal components to calculate the SNR. Only phases with an SNR value above the chosen threshold are kept for further analysis. The remaining phases are saved for the visual quality check.
Once the initial pre‐processing is completed, the user needs to view the remaining phases individually to decide which waveforms should be kept for the splitting analysis. To simplify the decision of which phases should be kept, SplitRacer shows the phase/event in several panels. Figure 4 shows the respective GUI with an example for an event on 16 September 2006. The filtered 1‐hr ZNE seismogram (top left) is cut to incorporate a time window of 100 s around the phase (bottom left). The amplitudes of usable phases should be seen clearly in both plots, which can serve as a first indication of the quality of the entire events. The horizontal components are then rotated by the initial back azimuth that leads to radial and transverse components (top right). As XKS phases are polarized in the direction of back azimuth, there should be little energy around the phase arrival time on the transverse component unless the phase traveled through an anisotropic medium. The user should then carefully check if a clear phase is visible just after the ray‐theoretical arrival time. If the phase is affected by anisotropy, its particle motion becomes elliptical (bottom right). In case of absence of anisotropy or if the back azimuth and fast axis of the anisotropic medium are parallel or perpendicular to each other, the particle motion is linear. We show the particle motion for the entire 100 s NE components, for the XKS signal time window and for the XKS signal time window filtered with a bandpass of 15–50 s (from left to right). This allows SplitRacer to check for a possible misalignment of the sensor by automatically calculating the difference between the theoretical back azimuth and the long axis of the ellipse of the long‐period particle motion (Rümpker and Silver, 1998). We use these values to correct for a possible misalignment of the sensor.
The user decides whether the phase should be kept for further analysis using the button panel placed below the particle motions taking the above and following criteria into account. The automatically chosen signal time window can be altered by the user if it does not capture the phase accurately. This is of significant importance in the presence of noise or other phases within the automatically chosen time window. Established criteria for choosing appropriate phases are the phase’s SNR, a clear phase arrival on both radial and transverse components, a nearly elliptical particle motion for the signal window, and the nonexistence of energy from other phases in the signal time window (e.g., Savage, 1999; Liu and Gao, 2013). Additionally, we also suggest to check whether a reasonable correction value for a possible misalignment of the sensor is calculated for the displayed phase.
Once all phases have been checked for their usability for the shear‐wave splitting analysis, the function check station misalignment assesses the distribution of the sensor’s possible misorientation (see Fig. 5 for an example) using the phases saved in the visual quality check step. The top left histogram in Figure 5 shows the distribution of misalignment values for station BKS, using ∼25 years of data starting in August 1991. The top right plot shows the distribution of misalignment values over back azimuth. This also shows the main axes of incoming events. The bottom plot displays the distribution of misalignment over time. If values change significantly with time, a possible change in sensor orientation can be detected. The general usability of a station might be called into question if misalignment values show scattering over the entire data range. A mean misalignment value is saved for further analysis; however, the user decides if it is appropriate to use a mean value in the splitting analysis itself. If not, a misalignment value for each phase is available for further analysis.
SPLITTING ANALYSIS (SINGLE PHASE)
Once the pre‐processing is completed, the actual shear‐wave splitting measurements are performed. SplitRacer’s Splitting Analysis offers three different shear‐wave splitting measurement functions. The first is a standard single‐phase shear‐wave splitting analysis tool, and tool for simultaneous inversion of all phases at a station for one and two layers. In SplitRacer, shear‐wave splitting analysis is based on the energy minimization method of Silver and Chan (1991) that has been shown to be the most robust technique for SKS‐type waves in which the back azimuth is readily known (e.g., Long and Silver, 2009). We choose to perform the splitting analysis in the frequency domain. First, NE components are rotated into radial and transverse components (see Fig. 6, upper left plot for the implementation in SplitRacer) using the back azimuth including a possible misalignment correction. This is followed by a grid search to find the splitting parameters (fast polarization) and dt (delay time) that best remove the energy on the transverse component. This corresponds to linearizing the particle motion. Here, we use a resolution of 0.1 s for dt with a maximum value of 4 s and 1° for between 0° and 180° for the fast polarization.
The analysis is performed using the settings chosen in the GUI. Available settings are the filter range, which can be imported from pre‐processing routines or put in manually, the number of time windows, and the range used in the analysis. For consistency, we suggest using the same settings as in pre‐processing, in which we opt for a wide broadband filter. However, as shear‐wave splitting measurements may vary with frequency (e.g., Rümpker et al., 2003), the data set can be processed using different filter settings. To test for frequency dependencies, we always show the particle motion for the entire trace and the XKS time window for the currently used filter, as well as for the XKS window filtered with 10–50 s and 15–50 s (Fig. 6, center top).
As default, we suggest using up to 50 slightly shifted time windows for the analysis (see, e.g., Walker et al., 2004), because this allows for a thorough statistical analysis of the results. Starting with the original time window as chosen in the automatic or manual pre‐processing steps, time windows are randomly generated around the original. As several studies have demonstrated the influence of the time window on the measurement (e.g., Teanby et al., 2004; Savage et al., 2010; Wüstefeld et al., 2010; Liu and Gao, 2013), this feature offers the possibility of checking upon consistency between measurements, average results over many time windows, and possibly discarding events. The splitting parameters of each measurement are shown in the histogram (see Fig. 6, upper right plot).
SplitRacer calculates a mean energy grid over all time windows (Fig. 6, bottom left plot) and the 95% confidence level (Fig. 6, bottom right), for which we adopted the corresponding routines from Wüstefeld et al. (2008). These are based on the original error calculation by Silver and Chan (1991). It assumes that the signal is combined with a Gaussian white noise process, in which the energy on the transverse component for a test fast polarization and delay time is a random χ2 variable with n degrees of freedom that can be estimated from the seismogram. Confidence levels are then calculated using the F‐test for two model parameters. We altered these by the findings of Walsh et al. (2013), who demonstrated errors made in the derivations for the error calculation by Silver and Chan (1991). These overestimated the degrees of freedom and led to smaller standard errors. Here, we derive our final pair of splitting parameters from the minimum value of the 95% confidence level. We apply these splitting parameters to formulate the inverse‐splitting operator (e.g., Rümpker and Silver, 1998) that effectively removes the energy on the transverse component and linearizes the particle motion (Fig. 6, bottom center).
SplitRacer first calculates shear‐wave splitting measurements using the chosen settings for the entire station or data set before the user views each phase with the different plot windows as described above and decides how to categorize each measurement. For this purpose, we plot the results in tabs for all phases in one window. The user can compare the different measurements by clicking through all tabs, which provides a straightforward overview of final data and measurement quality. If a time window does not capture the phase accurately, it can be altered again. Instead of categorizing measurements automatically, the user categorizes each measurement (see Fig. 6, bottom left). This offers more flexibility and control for the user but also requires accurate handling. Automatic categorization algorithms as suggested by, for example, Liu and Gao (2013) or Savage et al. (2010) are independent of the user and thus more objective, though some cases require the user’s interaction. We believe that users should be familiar with their own data and be able to adjust the categorization depending on the data set. Hence, we suggest keeping in mind several factors when deciding a category, which can be good, average, poor, or null for null measurement. A good quality measurement has a clear phase onset on both radial and transverse, near elliptical particle motion, none or only minor scattering of splitting parameters shown in the histogram, small 95% confidence level, as well as good removal of the energy on the transverse component as seen in the linearization of the particle motion. Usually, the energy of the XKS phase in question dominates the particle motion plot even for the entire 100 s. To judge the effectiveness of the application of the inverse‐splitting operator, we calculate the percentage of energy reduction on the corrected transverse component relative to the original component. From extensive testing, we found that a good measurement should have above 85% energy reduction. Figure A1 shows an example of a good measurement.
We realize, however, that there are cases of null measurements or complex anisotropy that might not fulfill these criteria. Null measurements can occur when there is no anisotropy along the ray path, or when the initial polarization is parallel or perpendicular to the fast axis of the anisotropic medium (see Fig. A2). To categorize as a null measurements, a clear phase arrival must be visible on the radial component, but not on the transverse component (e.g., Liu and Gao, 2013). Hence, the particle motion is linear from the start. Often, the splitting analysis will indicate an unusually high delay time as well as an energy grid showing minimum energy on two different fast polarization axes and a large confidence level. The application of the inverse‐splitting operator reproduces a similar particle motion to the uncorrected one, and the energy reduction is rather low or slightly negative. We suggest that measurements are not null if they scatter over the entire delay time and/or fast polarization range and if the application of the inverse‐splitting operator leads to an energy increase on the transverse component.
We suggest the label average measurement if the waveform appears slightly noisier than that of a good measurement, but still displays noticeable energy on the transverse component (Fig. A3). Particle motion may stray from the typical elliptical form, but its corrected counterpart should visibly become more linear. In cases of complex (e.g., 3D) anisotropy beneath a station, a near‐perfect removal of the energy on the transverse component may not be possible. The same is true for phases contaminated with noise, which cannot effectively be removed. We suggest that energy reduction should still be above 50%, whereas the histogram should show none or only minor scattering. The quality poor should be designated to all remaining cases, for example, nonelliptical particle motion only or in combination with high scattering of splitting parameters in the histogram and/or energy removal of less than 50%, or an increase in transverse‐component energy (see Fig. A4).
MODELING OF SPLITTING PARAMETERS
To visualize the results of the single phase analysis, SplitRacer offers an overview (see Fig. 7). To the right, the pie charts show which phases were used and how they were categorized. As shear‐wave splitting measurements are often checked upon their dependency on back azimuth, we plot the measured fast orientation and delay time as function of back azimuth, separately below. In cases of single‐layer anisotropy, the splitting parameters should be consistent within their error bars except for null measurements that (in theory) should be characterized by a back azimuth that is parallel or perpendicular to the measured fast polarization from other events with different back azimuths. Here, we offer to calculate a simple mean value.
If two (or more) anisotropic layers exist along the ray path, splitting parameters display a distinct π/2 periodicity that are usually denoted apparent splitting parameters (e.g., Savage and Silver, 1993; Silver and Savage, 1994; Rümpker and Silver, 1998). The option to fit two layers calculates the resulting apparent splitting parameters as a function of back azimuth for a period of 10 s for all possible two‐layer models using steps of 10° and 0.2 s in both layers (based on Silver and Savage, 1994). The 10 models whose calculated apparent splitting parameters yield the best fit to the observed/measured data are plotted together with the observed data in the lower graphics. Here, the parameters for the best‐fitting model are printed in black (Fig. 7). In the upper left plots, we show the distribution of layer parameters, in which each model is denoted by a distinct symbol in both panels. The fit between model and data highly depends on the number of measurements used and their back azimuthal distribution. Users should check carefully if the distribution of splitting parameters warrants a two‐layer case. From experience, we note that if a two‐layer case is well constrained, the 10 best models should be clustered within a narrow range of parameters for each layer. If the different models do not cluster closely together, as they do, for example, in Figure 8 for the continuous model, the observed apparent splitting parameters can be fit equally well by multiple two‐layer models. Hence, the fit is nonunique.
As Rümpker and Silver (1998) have shown, apparent splitting parameters as a function of back azimuth resulting from two anisotropic layers can often be fitted by a three‐parameter model with continuously changing fast‐axis directions between the bottom and top of the anisotropic region and an integrated delay time. Using the option fit continuous model, SplitRacer calculates apparent splitting parameters as a function of back azimuth for a period of 10 s resulting from models constructed by 10 layers, in which the fast‐axis direction varies by equal steps between fast axes chosen for the uppermost and lowermost layers and in which the delay times of the individual layers are equal and add up to give the integrated delay time. We show an example for station BKS in Figure 8.
As an alternative to the fitting of (apparent) splitting parameters resulting from theoretical models and observed data, SplitRacer also offers a simultaneous inversion for all waveforms. We apply this multiphase joint‐splitting approach for one layer if the splitting measurements do not show significant dependence as function of back azimuth, and for two layers if measurements show significant dependence as function of back azimuth indicative of a π/2 periodicity. Joint splitting relies purely on the waveforms themselves and inverts all waveforms, even null measurements, at a given station simultaneously. For this analysis, we use the waveforms previously categorized as good, average, or null, which were automatically stored for further analysis during the categorization process.
For one layer, we use the same approach as for the conventional (single waveform) splitting analysis that is a grid search over the splitting parameters (fast polarization) and dt (delay time). In the inversion process, we stack the energy grid of each phase. The smallest energy value of the stacked grids yields one pair of splitting parameters that should minimize the transverse energy component on all phases simultaneously. This is similar to the approach by Wolfe and Silver (1998) and was also used in a similar form by Homuth et al. (2014) and Wölbern et al. (2014). Joint splitting reduces the influence of noise significantly and has the potential to yield more robust measurement compared with a single‐phase splitting result. We calculate the total energy reduction from the amount of each phase’s individual energy reduction and visually check the linearization of all particle motions after application of the same inverse‐splitting operator to the waveforms. If the total energy reduction is high and most particle motions show significant linearization, we conclude that the anisotropy beneath a station may be best characterized by a single anisotropic layer. For stations BKS, this seems invalid as we will discuss later.
For two layers, we apply a grid search with four independent parameters, the fast polarizations and delay times for both layers. In the inversion, we apply a two‐layer inverse apparent splitting operator in the frequency domain (see Rümpker and Silver, 1998, their equation 10). The apparent splitting operator yielding the best energy minimization on all transverse components simultaneously is then used to correct the particle motions of the signal time window. Similarly, Özalaybey and Savage (1994) also used a two‐layer grid search on single waveforms and later stacked the transverse energy misfit spaces. Here, we calculate the errors using the same expressions as for the single and joint‐splitting one‐layer approach, but change the model parameters in the F‐test to four. Again, we calculate the total and individual energy reduction of the transverse component and check the linearization of the particle motions visually (see Fig. 9, right panel, for a full discussion see the Application to station BKS section). We also calculated synthetic waveforms based on a two‐layer anisotropic model for station BKS after Özalaybey and Savage (1994) to test this approach. In Reiss et al. (2016; see their appendix), we show that the joint splitting for two layers, based on synthetic waveforms, recovers the model exactly.
The user decides whether a one or two‐layer approach is more appropriate based on the outcome of the single‐phase splitting analysis. The joint‐splitting analysis can only be applied after calculating and categorizing the single‐phase splitting measurements (see Fig. 1). The analysis can either be calculated by including all successful measurements (good, average, and null) or by omitting the average phases. Hence, it is important to carefully categorize the results of the measurements. We suggest using the same settings for the joint‐splitting analysis as for the single‐splitting measurements. However, for two layers, the calculations are very time intensive. We thus suggest doing the analysis for one or a few time windows only. Given the simultaneous inversion for four parameters, this technique is slightly less robust if waveforms are noisy and different time windows introduce slightly different waveforms. Please see the Application to station BKS section for further discussion.
APPLICATION TO STATION BKS
We use SplitRacer to determine the anisotropic properties beneath the station BKS of the BDSN and to test the proposed workflow. The station is located in California, U.S.A., in close vicinity of the San Andreas fault (SAF), a transpressional dextral strike‐slip plate boundary that separates the Pacific from the North American plate (e.g., Wallace, 1990). In the last three decades, shear‐wave splitting in this region has been analyzed by a number of studies (Savage and Silver, 1993; Özalaybey and Savage, 1994; Silver and Savage, 1994; Hartog and Schwartz, 2001; Polet and Kanamori, 2002; Bonnin et al., 2010) that all observe variations of splitting parameters with back azimuth for stations located close to the SAF. BKS has been used in all these studies and is thus ideal for testing purposes, especially in the presence of complex anisotropy.
First, Silver and Savage (1994) observed azimuthal variations of splitting parameters at several stations close to SAF, one of which was BKS. They described these as apparent splitting parameters resulting from a two‐layer anisotropic model: the lower layer exhibits a fast‐axis direction of 90° and 0.9 s delay time, whereas the anisotropy of the upper layer is characterized by a fast‐axis direction of 140° (−40°) and 1.1 s delay time. The analysis of Özalaybey and Savage (1994) yielded a similar result with splitting parameters of 90°, 1.4 s and −45°, 1 s for the lower and upper layers, respectively. They used a four‐parameter grid search on the original waveforms, whereas Silver and Savage (1994) applied the grid search to the measured splitting parameters based on two‐layer models and compared their fit. Both found the upper layer to be parallel to the strike of the SAF. Özalaybey and Savage (1994) attributed the lower layer anisotropy to a remnant slab or to local effects due to the northward motion of the Mendocino triple junction. Hartog and Schwartz (2001) used both splitting parameters and waveforms in their inversion. However, this did not yield a unique solution for all parameters. They thus constrained the upper layer at −35°, assuming that anisotropy in this layer is caused by the relative motion between the Pacific and North American plates and found an approximate east/west direction of the lower layer. This might either be caused by past or present active convection or by motion of the North American lithosphere.
In this study, we use ∼25 years of data recorded at station BKS starting in August 1991. Using SplitRacer, we requested data above magnitude 6 and used a band‐pass filter of 4–50 s as well as a cutoff SNR of 2.5 during pre‐processing. Of 2003 downloaded seismograms, we obtained 301 XKS phases from 281 events that yielded an SNR above 2.5 and were used in the visual quality check. Of these, 86 events and 89 phases were kept to calculate a mean misalignment that resulted in a value of 0.74°. We applied this correction and performed the shear‐wave splitting analysis for all phases and for 50 (randomly chosen) time windows. The results of the analysis are summarized in Table 1. Using a very strict categorizing approach, we obtained 20 good, 39 average, 28 poor, and 2 null measurements (see Fig. 7). Our mean errors are 23° and 1 s for and dt, respectively, which are considerably larger than those obtained by Hartog and Schwartz (2001), which are 10° and 0.3 s for and dt, respectively. This is due to the more appropriate error formulation of Walsh et al. (2013), which is used in SplitRacer. The single‐splitting measurements exhibit a distinct π/2 periodicity, which we try to fit by calculating apparent splitting parameters resulting from two anisotropic layers, as described previously. Using a least‐squares approach, the best match between observed data and apparent splitting parameters is obtained by a model with lower layer anisotropy of ∼80° and 0.8 s and an upper layer of approximately −50° and 1.2 s. However, the distribution of the 10 best‐fitting two‐layer models shows that the observed splitting results can just as well be matched by various two‐layer models with slightly variable characteristics. However, the best‐fitting model agrees well with the one proposed by Silver and Savage (1994). We also try to fit our data using a continuous modeling approach that leads to a best‐fitting model with a lowermost fast axis of ∼80°, an uppermost fast axis of −40°, and an integrated delay time of 2 s (see Fig. A1). Here, the distribution of the 10 best‐fitting models indicates same variations along the delay time axis, but only very minor variations along the fast‐polarization axes.
As discussed previously, the fit between observed data and modeling results is nonunique. Hence, we also attempt to further constrain the anisotropy by applying the two‐layer joint‐splitting approach described in the previous section. We apply this analysis to all waveforms except for those that were categorized as poor. We also repeat the analysis by including only good and null measurements. Given the number of 61 phases, we use one time window only. For both cases, the analyses yield the same results with a lower anisotropic layer of 75°, 1 s, and an upper layer of −55°, 1.4 s, except that the error bars are slightly larger when using only phases categorized as good and null. This likely results from using only one‐third of the data in the inversion. When using all phases, confidence levels are very small in the order of one grid step for both parameters. With five time windows, the results remain unchanged (see Table 2). We check the linearization of particle motions visually to verify that the transverse‐component energy has been reduced. Overall, the energy reduction is ∼72% using all phases and ∼67% for using only good and nulls. In comparison with Özalaybey and Savage (1994), we find the delay times of the two layers to be reversed though they are within their error bars.
Overall, the best‐fitting two‐layer model and the results of the joint‐splitting analysis are in close agreement and corroborate the findings of Özalaybey and Savage (1994) and Silver and Savage (1994). However, we find that the inversion of waveform data leads to more stable results than the fitting of observed and calculated (apparent) splitting parameters. The joint‐splitting waveform analysis seems better suited to constrain two‐layer anisotropy, especially if the data coverage and its back‐azimuthal distribution is pronounced. Our results suggest that the anisotropy of the upper layer is slightly stronger than assumed previously. However, a tectonic interpretation of our results would require a more detailed analysis of data from neighboring stations, and this is beyond the scope of the current article.
We present MATLAB code and a GUI for the semiautomatic analysis and interpretation of teleseismic shear‐wave splitting, which is designed for the rapid processing of large amounts of data. The code is tested by analyzing 25 years of data recorded at station BKS of the BDSN. Using SplitRacer, seismogram data can be downloaded directly from all data centers offering FDSNWS. The mseed files are automatically pre‐processed based on the user’s settings, such as SNR and frequency range, which can be adjusted depending on the quality of the data set. Visual inspection of phases fulfilling the previously applied criteria is used to check or alter the time window for which the splitting analysis is done. Furthermore, SplitRacer allows to constrain a possible sensor misorientation and to apply an appropriate correction angle. The shear‐wave splitting analysis is performed per phase and can be repeated for, for example, 50 slightly alternating (randomly chosen) time windows that enable statistical error analysis. We also suggest guidelines on how to categorize the results of the splitting measurements. SplitRacer offers the possibility to fit modeling results based on one or more anisotropic layers to the observed splitting parameters. This, in conjunction with a joint inversion of all waveforms at one station in terms of one or two anisotropic layers, can be used to interpret results in cases of complex anisotropy. We show that for station BKS, the splitting measurements are best explained by a two‐layer anisotropic model with a lower layer of 80°, 1 s, and an upper layer of −55° and 1.4 s.
DATA AND RESOURCES
Waveform data for station BKS were accessed through the Northern California Earthquake Data Center (NCEDC), doi: 10.7932/NCEDC. SplitRacer is available at http://www.geophysik.uni-frankfurt.de/64002762/Software (last accessed date October 2016). MATLAB is available at http://www.mathworks.com/products/matlab/ (last accessed September 2016). Wget is available at https://www.gnu.org/software/wget/ (last accessed September 2016).
We thank all persons who tested earlier versions of SplitRacer and/or partly contributed to the development of the MATLAB code, especially Ingo Wölbern, Benjamin Homuth, Ulrike Löbl, and Manvendra Singh. We are indebted to Paul Silver for providing the impetus for this work during G.R.’s postdoctoral fellowship at the Department of Terrestrial Magnetism, Carnegie Institution of Washington (1996–1998). Some parts of the code are based on Fortran subroutines provided by Paul. We thank the two anonymous reviewers for constructive comments on the submitted article. This study has been supported by Deutsche Forschungsgemeinschaft.
Here, we show examples for different categories of shear‐wave splitting measurements as mentioned in the Splitting Analysis (Single Phase) section. This appendix is meant to serve as a guide on how to categorize shear‐wave splitting measurements. Figure A1 shows a splitting measurement we suggest to characterize as good. The NE as well as RT components of the event show a favorable signal‐to‐noise ratio (SNR), and the particle motions are elliptical. The histogram shows no scattering and the 95% confidence level is small. The energy reduction is very high (95%) and the particle motion is nicely linearized when applying the found pair of splitting parameters. Figure A2 shows a null‐measurement. The SNR is good but there is very little energy on the transverse component. Accordingly, the particle motion is nearly linear.
Figure A3 shows a phase we suggest to characterize as average. Though the SNR is good and the histogram shows no scattering, the particle motion is not as nicely linearized. Compared to a phase characterized as good, the energy reduction is lower (74%) and the confidence level larger.
Figure A4 shows a phase we would categorize as poor. The energy on the transverse component at the time of the phase arrival is not above the surrounding noise level and the particle motion is neither null nor elliptical. The histogram shows string scattering of results over the used time windows. Thus, errors become infinite.