# Seismological Research Letters

- © 2014 by the Seismological Society of America

*Online Material: *Figures of complete stress inversion results.

## INTRODUCTION

The estimation of the stress‐field orientation from focal mechanisms of earthquakes is a relevant tool to understand crustal mechanics and the physics of earthquakes. In global seismology, Formal Stress Inversion (FSI) is a well‐established technique to study tectonic processes associated with the occurrence of large earthquakes (e.g., Hardebeck and Michael, 2004; Yoshida *et al.*, 2012). Information on the stress‐field orientation is also relevant for the exploitation of hydrocarbon and geothermal reservoirs. Knowledge of the orientation of the maximum horizontal stress (*σ*_{Hmax}) is crucial for reservoir development, such as drilling and leakage risk assessment (Terakawa *et al.*, 2012; Martínez‐Garzón *et al.*, 2013). Additionally, the FSI technique may be useful for understanding the physics of rupture processes at a microscale in the frame of rock deformation experiments in the laboratory.

In seismological practice, the stress tensor is obtained either from inverting earthquake focal mechanisms or directly from first‐motion polarities. Most of the developed FSI methods share two main assumptions:

The stress field is homogeneous within the considered rock volume.

The slip of the fault is parallel to the direction of the tangential traction (Wallace, 1951; Bott, 1959).

Estimation of stress‐field orientation is a nonlinear inverse problem that can be solved either directly or linearized. When solving the nonlinear inverse problem, the best‐fitting stress tensor to the group of focal mechanisms is obtained using grid‐search algorithms (Gephart and Forsyth, 1984; Gephart, 1990; Arnold and Townend, 2007) or Monte Carlo sampling‐based optimization methods (Angelier, 1984; Xu, 2004). The linearized inversion scheme is solved by a generally less computationally extensive least‐squares technique (Michael, 1984; Hardebeck and Michael, 2006).

Because of the specific characteristics of the inversion problem, the uncertainty assessment is generally a necessary part of FSI. In particular, the following topics are essential:

Because a number of faults with different orientation can be reactivated under the same state of stress (McKenzie, 1969), a single focal mechanism only constrains the stress‐field orientation inside the corresponding dilatational quadrant. Therefore, the input dataset must contain certain variety in the focal mechanisms to constrain the solution.

When inverting for the stress tensor, the knowledge of the fault plane on which the slip has occurred is explicitly assumed. However, according to the concept of the seismic double‐couple, there are two possible fault planes: the true fault plane and an auxiliary fault plane (Aki and Richards, 2002). When lacking additional information, both fault planes are equally possible.

This paper discusses a software package, MSATSI^{1}, that allows the FSI to be used in the MATLAB (http://www.mathworks.com/products/matlab/; last accessed May 2014) environment. The new algorithm has been redesigned from the SATSI^{2} algorithm (Hardebeck and Michael, 2006). In our package, we include a visualization tool to represent 0D–3D FSI results using a variety of plots. We test MSATSI using datasets with different characteristics to test the validity of the FSI and the visualization tools for every range of 0D–3D input datasets. The package is licensed by GNU General Public License (GPL) and freely available to all users.

## SATSI STRESS INVERSION

In nature, the transition between two different states of stress must be continuous. However, during the investigation of a distribution of stress‐field orientations the obtained stress‐tensor solutions might depend on the way that the input data was bound (e.g., Hardebeck and Michael, 2004; Townend and Zoback, 2004). The SATSI (Hardebeck and Michael, 2006) algorithm was created to handle this problem. Here, the input focal mechanisms can be grouped prior to the inversion into a number of subareas (we use the term “grid points”) distributed over a number of dimensions, from 0D (single FSI in one grid point is performed) and 1D (e.g., temporal changes of stress field) up to 4D (e.g., spatiotemporal distribution; Fig. 1), with time treated as an additional coordinate. Then, a stress tensor for each grid point is inverted simultaneously using a damped least‐squares inversion scheme in order to obtain a smoothed solution. We have implemented SATSI due to the robustness of the FSI method combined with the damped least‐squares optimization (Menke, 1989) and because of its computational performance. For a complete description of SATSI, we refer to Hardebeck and Michael (2006) and Michael (1984, 1987).

The inversion provides the orientations of the three principal stress axes and a quantity () that reflects the relative stress magnitude *R*: (1)in which *σ*_{1}–*σ*_{3} are the magnitudes of the three principal stress axes obtained from the deviatoric stress tensor. The relative stress magnitude quantifies whether the magnitude of the intermediate principal stress *σ*_{2} is closer to the magnitude of the most compressive (*σ*_{1}) or the least compressive principal stress (*σ*_{3}). The bootstrap resampling method applied to the input focal mechanisms results in the uncertainties of the stress‐axes orientations and the relative stress magnitude.

Figure 2 shows the SATSI flowchart that is common regardless of the number of dimensions in which the grid points are distributed. However, the actual programs called during the FSI are different for the 0D/1D/2D and 3D/4D cases. The main steps in the calculation of stress‐field orientation using SATSI are as follows:

Define the optimum damping parameter for the inversion. This is performed by multiple calls to the subroutine

*satsi_tradeoff.exe*, resulting in a function depicting the dependence between data misfit and model length. The combination of data misfit and model length that minimizes these two values is then chosen.The damped FSI is performed simultaneously for all grid points using the subroutine

*satsi.exe*. The best stress‐tensor orientation and relative stress magnitude for each grid point is then determined.The uncertainty assessment is performed by bootstrap resampling of the original input data for each grid point. This is done by the

*bootmech.exe*subroutine, which creates new subsets by randomly sampling the original dataset with replacement and then performing new FSIs.Finally, the

*bootuncert.exe*subroutine selects the subset of bootstrap data within the specified confidence intervals.

The source code of SATSI package is written in C programming language. As there is a certain sequence of programs that must be called (compare with Fig. 2), Hardebeck and Michael (2006) provided additional Perl scripts to semiautomatize the calculation process. However, significant improvements can be made in terms of program workflow homogenization and optimization for different input data types, more efficient and flexible user input–output data handling, and features extensions.

## THE MSATSI PACKAGE

MSATSI is composed of two routines. The first (*msatsi.m*) wraps the whole FSI workflow (compare with Fig. 2). It handles data exchange between consecutive programs of the original SATSI package and provides a well‐structured output. The second (*msatsi_plot.m*) is a routine to graphically represent the FSI results. This routine is capable of handling FSI solutions of up to three dimensions and generating a variety of plots.

Setup of initial parameters for FSI, as well as the properties of graphical output, are handled using MATLAB properties (PropertyName–PropertyValue pairs). Here, we only discuss the most important features of both routines; the complete list of the properties is available in the package documentation.

### msatsi.m

The *msatsi.m* program calculates the stress orientation from focal mechanisms and provides the bootstrap resampling uncertainty assessment (compare with Fig. 2). A few aspects have been added and/or modified with respect to SATSI, the most relevant ones are as follows:

*msatsi.m*is capable of performing the FSI for the 0D case, that is, to calculate stress‐field orientation for a single grid point (the static stress field is calculated from provided data, with no spatial or temporal variations).The routine handles input/output data in the same way, regardless of their dimensionality. Therefore, the same programs are used for any dimension.

*msatsi.m*generates the distribution of P and T axes from the given set of focal mechanisms. This is useful for a visual inspection of input data to elaborate if the diversity in focal mechanisms is adequate to perform the FSI.The routine automatically assesses the appropriate damping parameter based on input data. The corresponding trade‐off plot is also generated (Fig. 3). Note that small variations in the value of the damping parameter do not have a significant impact on the results.

*msatsi.m*handles the FSI and input/output errors generated by SATSI executable programs.All relevant FSI results are gathered in one output for user convenience. Additional information regarding the bootstrap resampling data contained in the selected confidence interval is also included.

Table 1 lists the most important parameters controlling the FSI procedure, together with suggested intervals and default values. A minimum number of seismic events per grid point is necessary to constrain the FSI (see Introduction). If the minimum number is not reached, the inversion will not be performed for this particular grid point. The fraction of correctly picked fault planes *F* deals with the ambiguity of the fault‐plane solution from seismic data. With MSATSI, the user can explicitly specify the probability that each input fault plane is believed to be the real one. This is useful, for example, if the input fault planes come from direct geologic measurements (in this case *F*=1). The number of iterations for the bootstrap resampling is important for the assessment of the uncertainty. A rule‐of‐thumb would be to perform approximately 20 times the number of input data (Efron and Tibshirani, 1986). Last, the desirable confidence interval of the uncertainties should be specified. Generally, the default value equal to 95% is used.

### msatsi_plot.m

The purpose of the *msatsi_plot.m* routine is to provide a variety of graphical representations of the stress‐tensor inversion results calculated with *msatsi.m*. The program is capable of displaying all varieties of 0D–3D stress‐inversion output data. The routine also gives the opportunity to represent a reduced number of dimensions of a performed FSI set by so‐called “slices.” For example, a 3D distribution can be sliced into a set of 2D maps along an arbitrary coordinate. This feature substantially simplifies the analysis of FSI results.

The routine includes five basic types of plots. Generally, they display the orientation of the principal stress axes and the relative stress magnitude (*R*), together with the uncertainties. The uncertainties can be represented as confidence intervals or scattered clouds of bootstrap‐resampled solutions. The five available plots are as follows:

The

*Stereonet*plot is suitable to represent the 0D–1D FSI results (e.g., only one FSI/temporal changes). The direction of the principal stress axes is plotted using stereonet and lower hemisphere equal‐area/angle projection (Fig. 4a).*R*values for each grid point are plotted in a separate figure.*Profile*is appropriate to represent 1D FSI results. This plot is composed of three figures showing changes in plunge and trend of the principal stress axes and*R*for each grid point (Fig. 4b).*Stereomap*is appropriate for visualization of 2D FSI results (e.g., surface distribution). Orientations of principal stress axes for each grid point forming a map are plotted using the stereonet (Fig 5b). The 2D distribution of*R*values is shown as a separate figure.The

*WSM*(World Stress Map) option creates a 2D plot compatible with those created by World Stress Map project (Heidbach*et al.*, 2010). This map plot shows the orientation of the maximum horizontal stress according to the respective stress faulting regime classification criteria (Zoback, 1992) applied to the best stress‐inversion solution and corrected by the estimation of the maximum horizontal stress published in Lund and Townend (2007; Fig. 5a). The estimation is performed by the subroutine*SH.m*(available for public use) attached inside the main code of*msatsi_plot.m*.*Stereovolume*is appropriate for 3D FSI data. The results are presented as a 3D plot with cardinal stress‐axis directions pointing according to their trends and plunges. The idea is analogous to*Stereomap*but extended to 3D (Fig. 6).

The generated plots can be modified using over 30 properties to adjust the plots to the user preference. The complete list of properties is provided in the MSATSI documentation. Multiple code examples are provided, together with detailed descriptions of how the graphical representation of the results was obtained.

### Content of Software Package MSATSI

The MSATSI package is composed of the following elements:

The root directory contains

*msatsi.m*and*msatsi_plot.m*, together with modified SATSI executable files compiled under the Microsoft Windows platform.The /bin folder contains modified SATSI executable files compiled for Windows, Linux, and Mac platforms.

The /examples directory contains illustrations of the usage of the MSATSI package. Examples 1D–3D cover the cases presented in this paper. Example 0D presents another case study not described here.

Last, /src contains modified C codes from the original SATSI package. Package compilation instructions under different platforms are also provided.

## APPLICATION OF MSATSI TO DIFFERENT CASE STUDIES

### Temporal Stress‐Field Orientation Changes at the North Anatolian Fault Zone (1D)

The North Anatolian fault zone (NAFZ) represents the plate boundary between the Anatolian and Eurasian plates and is one of the best‐studied transform fault zones in the world. It extends along more than 1200 km from eastern Anatolia toward the North Aegean Sea.

Our study was designed to better understand whether local rotations of the stress field can be used as an indicator to characterize the seismotectonic setting along individual fault segments during the seismic cycle. For this purpose, we applied MSATSI to investigate temporal stress‐field variations associated with the occurrence of the 1999 *M*_{w} 7.4 İzmit earthquake in the NAFZ in northwest Turkey (Bohnhoff *et al.*, 2006; Ickrath *et al.*, 2014).

The dataset used here covers the preseismic, aftershock, and postseismic phases of the 1999 *M*_{w} 7.4 İzmit earthquake, which ruptured a 140 km long segment reflecting a right‐lateral strike‐slip faulting mechanism. We made a compilation of 989 fault‐plane solutions of seismic events distributed along the İzmit rupture segment of the NAFZ with magnitudes between 0.8 and 7.4, determined from various local and regional seismic networks (Bulut *et al.*, 2007; Görgün *et al.*, 2010; Örgülü, 2011). The seismicity covers the time interval between 1943 and 2007.

The input dataset was subdivided into a period before the İzmit earthquake, two months of data of İzmit aftershocks, and the postseismic period. The aftershock period, which contains a higher number of focal mechanisms, was also further subdivided temporally, following the segmentation of the İzmit event using a moving time window of 70 events with overlap of 20 events. The subsets were inverted for the stress tensor using MSATSI. We used a damped FSI and generated 2000 bootstrap resamplings of the original dataset. The uncertainties were determined assuming 95% confidence intervals.

Figure 4 shows the results from the FSIs using Stereonet and Profile plots. We found distinct temporal variations of the principal stress‐axes orientations. The regional stress field prior to İzmit and following the three‐month aftershock sequence reflects a stable strike‐slip regime with an ∼N130°E‐trending subhorizontal maximum principal stress (*σ*_{1}), in accordance with the dextral east–west‐trending NAFZ. In contrast, the early coseismic period is dominated by a northeast–southwest extensional normal‐faulting regime. For a detailed description and interpretation of these results, refer to Ickrath *et al.*, (2014).

### Surface Stress‐Field Orientation at The Geysers Geothermal Field (2D)

The Geysers (TG) geothermal field in northern California is the world’s largest producing geothermal field. It covers an area of 280 km^{2} and is composed of approximately 70 reinjection and 350 steam‐production wells. The existing microseismicity, with more than 16800 events (*M*≥1) recorded between 2007 and 2012, is induced by water injection and steam production (e.g., Martínez‐Garzón *et al.*, 2013).

We selected 4859 focal mechanisms (*M*_{D} 1–4.5) distributed over the entire geothermal field, available through the North California Earthquake Data Centre, to analyze the surface distribution of the stress field at reservoir depth (located on average between 2 and 2.5 km). TG (20 km×14 km) was divided into grid points covering an area of 2 km×2 km. We then performed the 2D damped FSI for grid points containing more than 30 microearthquakes. The uncertainty assessment was performed by bootstrap resampling of the input data 2000 times, and the confidence interval was set to 95%.

Figure 5a shows the orientation of the maximum horizontal stress for each of the FSIs performed using the *WSM* plot. For most of the grid points, the direction of *σ*_{Hmax} is in good agreement with the regional stress field in northern California (Provost and Houston, 2003). Although the local stress regime within the reservoir is generally normal faulting/transtensional, the grid point with coordinates easting = 515 km and northing = 4295 km is already within the strike‐slip criteria.

Figure 5b presents the distribution of principal stress axes using the *Stereomap* plot. For the cases in which the bootstrap resampled points are scattered over large areas (e.g., *σ*_{1} and *σ*_{2} stress axes with northing = 4295 km), the representation of the confidence intervals with the *bootstrap* option might provide a more realistic representation. For the grid point showing strike slip in Figure 5a, here we see that the uncertainties for the plunges of *σ*_{1} and *σ*_{2} are very large. This fact, together with the low *R* value, shows that the stress magnitudes of these axes are similar. Finally, in the southern part of the field (inversions with easting = 525 km), the bootstrap resampling for *σ*_{2} and *σ*_{3} axes are distributed over a larger range of azimuths, which suggest that here the magnitude of these axes are small compared to *σ*_{1} and, thus, the uncertainties in the trend of *σ*_{Hmax} are larger. These observations are supported by the distribution of *R* values for each grid point.

At a global scale, the normal stress regime at TG may be related to the volumetric contraction of the reservoir and consequent reduction of horizontal stresses within the reservoir. However, the normal stress regime observed at some particular areas of the field where geothermal exploitation is relatively new may also be related to the tensile opening of new fresh fractures to increase the permeability of the reservoir Enhanced Geothermal Systems.

### Synthetic Test (3D)

We have created a synthetic catalog of focal mechanisms with the objective of investigating how the fault geometry affects the final FSI result and to check the correct performance of the 3D procedure. We generated a 3D grid network (3×5×5) by modification of the three fault‐plane parameters: dip direction , dip angle *δ*={30°,45°,60°,75°,90°}, and rake *λ*={−90°,−45°,0°,45°,90°}. Therefore, the different combinations of these parameters cover the whole range of faulting styles from pure strike slip, through a variety of oblique faulting types to a pure thrust or normal faulting. For each grid point, 100 focal mechanisms were generated. The *m*‐th randomly generated fault plane for a grid point with subindexes *i*, *j*, *k* was described by the triplet of fault‐plane parameters in which *N*(·) are random numbers drawn from a normal distribution with means equal to the given values of dip direction, dip angle, and rake, and , *σ*_{δ}, and *σ*_{λ} are the standard deviations of the sampled focal mechanisms.

In this test, all standard deviations were settled to . This number was selected in order to ensure a correct variability in the data. Other performed synthetic tests using a uniform distribution found that, to properly constrain the stress tensor, a minimum variability of 30° should be included in the fault‐plane solutions (Hardebeck and Hauksson, 2001). As a result, the tension T and pressure P axes of the fault‐plane solutions oscillate around some preferred orientation. Therefore, the orientation of the most compressive (*σ*_{1}) and least compressive (*σ*_{3}) principal stress axes should correspond to the average orientation of the P and T axes, respectively. This could be seen as a case in which most of the earthquakes occur on optimally oriented faults with respect to the stress field.

In this analysis, we did not use the damped scheme because we did not want to obtain results influenced by the neighboring grid points. The uncertainty assessment was calculated by 2000 bootstrap resampling of the original input data.

Although we know *a priori* that the fault plane specified here is the appropriate one to be inverted, we tried two values for the fraction of correctly picked fault planes. In the first case, we assumed that all specified fault planes were correct, which would be analogous, for example, to inverting geologic slickensides (*F*=1); whereas, in the second case, we randomly selected one of the two possible options (*F*=0.5). This is the most common case in seismology when no additional data is available. Therefore, this allows us to test how different the solutions are by randomly selecting one fault plane.

The visualization of all the inverted cases for each of the two parameters plotted with the option *Stereovolume* is available in the electronic supplement to this paper. To facilitate the visualization, Figure 6 shows the stress‐inversion results for the two *F* values tested and the rake *λ*=0°, which includes the case of pure strike slip. The results obtained in both cases are remarkably similar, which indicates that randomly selecting the fault plane also provides a valid solution. However, certain obliquities are introduced in the horizontal axes for the *F*=0.5 case, given that the inversion is capable of distinguishing better between certain pairs of axes; therefore, for this case, the uncertainties appear slightly smaller than for the case of *F*=1.

As expected, the least crucial parameter to drive the stress regime is the dip direction, because this parameter systematically rotates the best solution of the three principal stress axes along the azimuth.

Once variations with the fraction of correct fault planes and with the dip direction were checked, we reduced the initial 3D grid to 2D by removing the dip direction to more readily visualize the other two parameters. In Figure 7, we show the data slice for and *F*=1. The three main faulting regimes according to Anderson (1951) are visible together with intermediate states. As it was expected, purely normal faulting (*σ*_{1} vertical) is found at *δ*=45°, *λ*=90°. Pure strike‐slip (*σ*_{2} vertical) regime appears at *δ*=90°, *λ*=0°. Finally, pure thrust faulting (*σ*_{3} vertical) is found at *δ*=45°, *λ*=−90°. For the same values of the rake, variations in the dip angle will result in introducing certain obliquity in the orientation of the stress axes that will be increased for the largest dip angles. However, the rake is the main influence in determining the stress regime.

The distribution of *R* values corresponding to the stress inversions shown in this slice is plotted in Figure 8. We noticed that for the stress regimes closer to pure normal faulting, the *R* value tends toward 1; in the regimes closer to pure inverse, the *R* value tends toward 0. This is explained by the fact that in the pure cases (less oblique) the magnitude of the horizontal stresses is very similar compared to the magnitude of the vertical stress. This difference is even more prominent for the case of *F*=1 (see electronic supplement).

## CONCLUSIONS

The MSATSI software package conducts a thorough stress‐tensor inversion under the MATLAB environment based on the SATSI algorithm (Hardebeck and Michael, 2006). Within MSATSI, the entire FSI procedure has been automatized, improved, and made more user‐friendly. With the second program, a whole suite of graphical representations of the FSI output data has been provided to appropriately visualize the stress‐tensor inversion results and their confidence intervals. The MSATSI package is GPL‐licensed and is freely available.

MSATSI has been applied to different datasets to test the performance of the developed routines over a wide variety of circumstances. For analyses of the natural seismicity along the NAFZ in northwest Turkey and the induced seismicity from The Geysers geothermal field in California, the results are consistent with other studies and/or field observations, indicating that the performance of MSATSI is correct. Additionally, we have also validated the 3D FSI using synthetic subsets of focal mechanisms.

Among the numerous FSI packages already existing, MSATSI has the advantage of unifying most of the potentially desirable circumstances. Although it can be used to analyze variations of the stress‐field orientation by using different grid points, it is also possible to perform only one stress inversion for the input data. In addition, the procedure can be damped (smoothed) to obtain more realistic results from natural observations, but this can be also deactivated for desirable cases (e.g., our synthetic example). Finally, we have demonstrated that the performance of the method is correct regardless of the size or nature of the input earthquakes.

## ACKNOWLEDGMENTS

We acknowledge funding within the Helmholtz‐Alberta Initiative (HAI) and from the Helmholtz Association in the framework of the Helmholtz Young Investigator Group ”From microseismicity to large earthquakes.” We thank Managing Editor Mary George and Associate Editor John N. Louie for useful remarks. We thank Jeanne Hardebeck for her valuable comments that helped to improve the manuscript. We acknowledge Oliver Heidbach and John Townend for providing reviews to the initial version of the manuscript.

## Footnotes

↵1 MSATSI is available to download from http://www.gfz-potsdam.de/en/research/organizational-units/departments/department-3/geomechanics-and-rheology/projects/software/msatsi/ and http://www.induced.pl/msatsi (last accessed May 2014).

↵2 SATSI is available to download from http://earthquake.usgs.gov/research/software/#SATSI/ (last accessed May 2014).