Evaluating Centroid-Moment-Tensor Uncertainty in the New Version of ...

23 downloads 2019 Views 5MB Size Report
finite-source representation and complete elastodynamic wave field. The KIWI software (KInematic Waveform Inversion;. Cesca et al., 2010) allows the automatic ...
E ○

Evaluating Centroid-Moment-Tensor Uncertainty in the New Version of ISOLA Software by Efthimios Sokos and Jirˇ í Zahradník Online Material: Additional figures.

INTRODUCTION Focal mechanism and moment tensor determinations based on regional and local waveforms have become routine tasks in seismology. In recent years, considerable effort has been dedicated to extend software capabilities, in particular towards the inclusion of finite-extent sources and easy, user-friendly use. Two examples of these efforts are the FMNEAREG and KIWI software packages. Focal Mechanism from NEAr source to REGional distance records (FMNEAREG; Delouis et al., 2008; Maercklin et al., 2011) performs a grid search in order to find pure double-couple sources, based on a 1D (line) finite-source representation and complete elastodynamic wave field. The KIWI software (KInematic Waveform Inversion; Cesca et al., 2010) allows the automatic retrieval of pointsource parameters, and for magnitudes above a given threshold kinematic finite-fault rupture models can also be obtained. KIWI relies on a pre-calculated database of Green’s functions. KIWI performs the moment tensor inversion in two steps, the first of which is based on the fit of amplitude spectra. The second step is based on time-domain waveform fits. The ISOLA, ISOLated Asperities, software package (Sokos and Zahradník, 2008) also performs waveform inversions to find source parameters. ISOLA is based on FORTRAN codes and offers a user-friendly MATLAB graphical interface (GUI; www.mathworks.com/products/MATLAB). ISOLA allows for both single- and multiple-point-source iterative deconvolution (Kikuchi and Kanamori, 1991) inversion of complete regional and local waveforms. The moment tensor is retrieved through a least squares inversion, whereas the position and origin time of the point sources are grid searched. The computation options include inversion to retrieve the full moment tensor (MT), deviatoric MT, and pure doublecouple MT. Finite-extent source inversions may also be performed by prescribing a priori the double-couple mechanism (which will remain homogeneous over the fault plane). Green’s functions, including near-field terms, are calculated 656

Seismological Research Letters

Volume 84, Number 4

using the discrete wavenumber method (Bouchon, 1981; Coutant, 1989). ISOLA has been implemented since 2006 at the University of Patras, Seismological Laboratory (UPSL), Greece, to routinely compute moment tensors for M > 3:5 events in western Greece, and, since 2012 also at the National Observatory of Athens, to characterize events that occur over the whole country; the solutions are reported to the European Mediterranean Seismological Center. In 2012, the Iranian Seismological Center (Institute of Geophysics, University of Tehran) began using the code for M w > 5 earthquakes in Iran. ISOLA has also been employed in a number of research applications, for example, to study 30 earthquakes of M w 2.5–5.3 of the Efpalio 2010 earthquake sequence (Sokos et al., 2012), the M w 6.3 Movri Mountain 2008 earthquake (Gallovič et al., 2009), the M w 6.3 L’Aquila 2009 earthquake (Gallovič and Zahradník, 2012), the M w 7.2 Van 2010 earthquake (Zahradník and Sokos, 2011), and the M w 9 Tohoku earthquake (Zahradník et al., 2011). Inversion artifacts (spurious slip patches) have also been studied theoretically using ISOLA by Zahradník and Gallovič (2010) and Gallovič and Zahradník (2011). A unique analysis of two microearthquakes, M w 0.2 and 0.4, using ISOLA was carried out by Benetatos et al. (2013). Fojtíková et al. (2010) were able to successfully use ISOLA to study earthquakes with magnitudes down to M w 1.2. The moment tensor determinations still need improvements, in particular as regards the assessment of the solution reliability, or uncertainty (e.g., Valentine and Trampert, 2012). In this paper we will focus on this aspect. Many users assess the quality of their solutions by exclusively looking at the waveform fits. However, if only a few waveforms are used, the waveform fit can be excellent, whereas the focal mechanism can be completely incorrect. In this paper we describe a few ISOLA innovations towards the simple assessment of solution reliability. This assessment combines the signal-to-noise ratio in the frequency range of the inversion, waveform fit (hence variance reduction), condition number, and also the space–time variability of the solution. We illustrate these innovations with examples based on a real data set from Greece (Fig. 1). July/August 2013

doi: 10.1785/0220130002

MAIN INNOVATIONS Signal-to-Noise Ratio The signal-to-noise ratio (SNR) is evaluated based on amplitude spectra (Fig. 2). First, the user selects the signal time window (STW) on a waveform plot. The noise time window (NTW) is then chosen automatically as to have the same length as the STW and to end at the start of the STW. Thus, the noise window immediately precedes the signal window. Amplitude spectra are computed for both windows, smoothed and plotted individually for each component of ground motion. The SNR is computed as a function of frequency (SNR curve), averaged over the components, plotted and saved in a file (one file per station). Users may consult the SNR curve when choosing the frequency band that will be used to filter the waveforms during the inversion. A single SNR value (SNR average) is calculated by averaging the SNR curve over the entire frequency range chosen for the inversion. Whenever the frequency range chosen for the inversion is changed, the SNR-average value is updated. At the end, the SNR-average value is used to assess the quality of the solution. Caution is needed when records contain low-frequency disturbances (Zahradník and Plešinger, 2005, 2010), representing a kind of signal-generated noise. In such a case, high SNR values at low frequencies do not indicate the usable frequency range. In ISOLA, waveforms are band-pass filtered before the actual inversion. The frequency range for band-pass filtering is specified by four values: f 1, f 2, f 3, and f 4. The filter, which is applied both to real and synthetic waveforms, is flat ( 1) between f 2 and f 3, and cosine-tapered between f 1 and f 2, and

22°

22.5° 38.6°

PVO ANX

EFP

38.4°

SERG

DSF KALE

TRIZ

38.4°

38.6°

21.5°

DRO

21.5°



38°

38°

38.2°

38.2°

LAKA

GUR

22°

22.5°

Figure 1. Test event of 25 April 2012 (10:34 UTC; M w 4.3, Lat 38.4045° N, Lon 21.9877° E). Stations are marked by triangles. The beachball, plotted at the epicenter position, corresponds to the reference fault-plane solution, Case 1 of Figure 3. See also Ⓔ Figure S1 (in the electronic supplement).

again between f 3 and f 4. For simplicity, we discuss here mainly the f 1 and f 4 values. The SNR curves help mainly to define f 1, because the noise level (either natural or instrumental) limits the usable low-frequency range. The f 4 value is basically dictated by our limited knowledge of the Earth structure; indeed, with existing crustal models we can usually invert near stations (< ∼1 km) up to ∼1 Hz, whereas near-regional stations (∼100 km) can be inverted up to ∼0:1 Hz, and regional stations ∼1000 km up to ∼0:01 Hz. Example 1: SNR Case 1 of Figure 3 is the reference case (see Fig. 1). We invert all stations in the same frequency range: 0.04–0.05–0.15– 0.16 Hz. However, as shown by the SNR analysis, it is more appropriate to filter the waveforms recorded at SERG (Fig. 2) using a lower corner of 0.1 Hz due to the low-frequency noise of the accelerographs below this frequency, and similarly also at TRZ. For the broadband station KALE (also in Fig. 2) and for stations more distant than KALE, a proper lower corner is around 0.05 Hz. As regards the high-frequency limit at PVO and more distant stations we should decrease the maximum inverted frequency from 0.16 to, say, 0.09 Hz to ensure reasonable modeling. These specifications are used in Case 2 of Figure 3 in which we employ the new ISOLA feature: the station-dependent frequency range. The waveform match is shown in Ⓔ Figure S2 (available in the electronic supplement to this paper). In practice, we anticipate that the use of stationdependent frequency ranges will more remarkably affect specific cases, for example, when the station distribution is not favorable. Variance Reduction and Condition Number The variance reduction VR and condition number CN are complementary measures of solution quality. Definition: Variance reduction is defined as VR  1 − Σo − s2 =Σo2 , where o and s are the observed and synthetic waveforms, and the L2 norms in the nominator and denominator are evaluated as summation over all stations, components, and time samples. VR is commonly used to assess solution quality, but has two limitations: (i) Although VR is a global measure of waveform match, it may be biased by stations with the largest amplitudes, which usually are best fit. In other words, VR may have a large value in spite of poorly fit waveforms at some stations. Therefore, besides assessing global VR , ISOLA users also have the possibility to check VRs for individual station components. (ii) VR may attain very high values when only few stations are used in the inversion; then, naturally, a large VR value is not indicative of a reliable solution. VR should always be used together with other indicators to assess solution quality. The condition number (CN) is derived from the Green’s functions matrix G, which relates data u (observed waveforms) and model parameters m (coefficients of the linear combination of six elementary seismograms corresponding to six basic focal mechanisms), u  Gm. Definition: CN is the ratio of maximum to minimum singular values of G. In

Seismological Research Letters

Volume 84, Number 4 July/August 2013

657

(a)

(c)

Signal Noise

108

Amplitude (count*s)

104

106

SERG 102 104

102 100 NS,EW,Z Mean 100 10−2

10−1

100

101

102

10−2

10−1

Frequency (Hz)

100

101

102

Frequency (Hz)

(b)

(d)

Signal Noise

108

Amplitude (count*s)

104

106

KALE 102 104

102 100 NS,EW,Z Mean 100 10−2

10−1

100

101

102

10−2

Frequency (Hz)

10−1

100

101

102

Frequency (Hz)

▴ Figure 2. Signal and noise spectra for stations (a,c) SERG and (b,d) KALE. (a) and (b) The signal and noise Fourier amplitude spectra. (c) and (d) The signal-to-noise ratios per component and mean values. Dashed vertical line is the lowest usable frequency.

ISOLA, the singular-value decomposition of G is not actually computed; instead, CN  is obtained indirectly as CN  p maxeigenval = mineigenval , where max_eigenval and min_ eigenval are eigenvalues of matrix GT G, and GT denotes transposition of G. CN measures the reliability of the inversion from the viewpoint of the source-station configuration, frequency range, and crustal model used. CN does not depend on the particular data used in each inversion; that is, CN can be calculated without using actual waveforms (Zahradník and Custódio, 2012). In this paper we evaluate CN for the bestfitting source position. The most informative cases are those in which comparable applications provide contrastingly very large or very small values of CN, indicating ill- and wellconditioned problems, respectively. What is called large and small depends on the application, and we will give examples later. Ill-conditioned problems are always alarming, whereas well-conditioned problems may still yield solutions with considerable uncertainty (e.g., none of the singular values are much smaller than the others, but all are small). 658

Seismological Research Letters

Volume 84, Number 4

Focal-Mechanism Variability Index, FMVAR A useful indicator of solution quality can be obtained by quantifying the variability of the focal-mechanism solution in the vicinity of the best-fitting source position and time. The so-called correlation plot, in which the correlation between observed and synthetic waveforms is analyzed as a function of the trial source position and time, was already available in the previous version of ISOLA. The new version uses the correlation plot to better focus on the focal-mechanism variability around the optimum solution. Let corropt denote the largest correlation value (the one corresponding to the best-fit, or optimal, solution). We set up a correlation threshold in this paper equal to 0:9 × corropt , and define the acceptable solutions as those having their correlation between 0:9 × corropt and corropt . Each acceptable solution is compared with the optimal solution using the so-called Kagan angle (Kagan, 1991). This angle, hereafter K-angle, expresses the minimum rotation between two focal mechanisms. We can now introduce the variability measure FMVAR as follows: Definition: FMVAR is the July/August 2013

Case

Stations

1 reference

As Case 2

Frequency Band (Hz) 0.04-0.05-0.15-0.16

Strike/Dip/Rake and K- angle from reference (°)

Focal STVAR Mechanism

Mw

Depth (km)

DC %

SNR

VR

CN

FMVAR (°)

327/32/-45 0

4.3

8

87

97

0.69

4.9

10

0.11

322/30/-45 5

4.3

7

96

79

0.78

4.3

17

0.11

SERG

0.10-0.10-0.15-0.16

EFP none TRIZ KALE LAKA ANX only NS PVO only EW,Z DSF DRO only NS,Z GUR

none 0.10-0.10-0.15-0.16 0.05-0.06-0.15-0.16 0.05-0.06-0.15-0.16 0.04-0.05-0.15-0.16 0.05-0.06-0.08-0.09 0.05-0.06-0.08-0.09 0.05-0.06-0.08-0.09 0.05-0.06-0.08-0.09

3

All except SERG, TRIZ, KALE

As Case 2

311/26/-70 14

4.4

8

60

60

0.52

4.8

12

0.13

4

SERG TRIZ

As Case 2

217/23/-144 45

4.6

12

74

7

0.88

13.7

23

0.20

5

TRIZ GUR

As Case 2

284/20/-93 23

4.5

8

72

6

0.98

14

38

0.25

6

DSF GUR

As Case 2

316/33/-53 6

4.4

11

99

9

0.89

2.1

16

0.42

7

DRO GUR

As Case 2

321/37/-46 7

4.4

11

75

7

0.82

2.4

22

0.49

8

GUR

As Case 2

279/84/79 44

4.6

4

20

2

0.98

4.8

63

0.81

9

TRIZ

As Case 2

256/23/-129 37

4.6

9

67

4

0.99

21

45

0.45

10

LAKA

As Case 2

301/31/-102 37

4.4

9

95

10

0.98

9.7

35

0.21

2



Figure 3. Summary of the experiments performed with the real-data test event. The individual cases simulate scenarios in which only some stations of Figure 1 are used, the reference solution (Case 1) is not known and user wants to judge the quality of the solution by means of the indices printed in bold (SNR, VR, CN, FMVAR, and STVAR), using also plots in Figures 4 and 5. The underlined values indicate an unreliable solution.

mean K-angle of all acceptable solutions (as specified by the user-defined correlation threshold) with respect to the best-fit solution. Because Kagan angles address only the double-couple (DC) variability of the acceptable solutions, ISOLA also constructs histograms of their DC-percentage (DC%). Deviation of the DC% from 100 reflects the amount of the non-DC components (i.e., compensated linear vector dipole [CLVD]

and isotropic, or just CLVD, in case of the full and deviatoric MT inversion, respectively). The non-DC components have been often observed, especially in geothermal fields (e.g., Cannata et al., 2009; Foulger et al., 2004). It is important to study the DC-percentage in relation with the correlation plots because the amount of DC% trades off with the optimum space–time position of the source (e.g., Křížová et al., 2013).

Seismological Research Letters

Volume 84, Number 4 July/August 2013

659

1

2

3

4

5

6

7

8

9

10

▴ Figure 4. Nodal lines and P, T axes corresponding to the space–time grid search in Cases 1 to 10 of Figure 3. The same earthquake is studied in various scenarios, simulating the availability of a different number of stations. Acceptable solutions, whose correlation is between 0:9 × corropt and corropt , are plotted (black nodal lines) along with the best-fit solution (white nodal lines). Space–Time Variability Index, STVAR This index measures the size of the space–time region corresponding to the given correlation threshold, independently of the focal-mechanism variation. Definition: STVAR is the surface of the area in the space–time plot occupied by the solutions within a given correlation threshold, normalized by the total area of the investigated space–time region. STVAR is complementary to FMVAR . The most favorable applications are those in which both FMVAR and STVAR indexes are small. Example 2: VR, CN, FMVAR, and STVAR This example again refers to the same earthquake as Example 1. We shall create various scenarios, that is, simulate situations for which only a few stations are available, a solution is obtained, and we need to decide whether this solution is reliable. We will be able to check whether our assessment was correct a posteriori, by comparing each solution with the reference solution. Let us consider all cases in Figure 3. The nodal lines corresponding to the acceptable solutions specified by the correlation threshold are shown in Figure 4. Two groups of solutions can be distinguished and identified by FMVAR in Figure 3: Group A solutions (Cases 2, 3, 4, 6, and 7) are relatively concentrated nodal lines, not deviating much from the bestfit solution, characterized by FMVAR < 30°, and the remaining Group B solutions (Cases 5, 8, 9, and 10). Our goal is to guess which of these inversions presents a reliable solution. Besides FMVAR , other parameters such as SNR , CN, or STVAR may also be used to assess solution quality. One possibility in a blind situation, that is, when the reference solution is not known, is to nominate as reliable all solutions in Group A, provided they have a relatively low CN (e.g., < 10). Solutions 2, 3, 6, and 7 satisfy this requirement. 660

Seismological Research Letters

Volume 84, Number 4

We may also require a good space–time resolution, for example, STVAR < 0:30. We would then remove solutions 6 and 7, and the only remaining reliable solutions would be 2 and 3. If we further request a high VR (e.g., > 0:6) we would remove solution 3. We could also request a good SNR (e.g., > 5), which in this example would eliminate solutions 8 and 9, both of which had already been rejected by the previous criteria. To check whether our blind guess of reliable focal mechanisms is good we adopt a criterion that the deviation between the best-fit solution and the reference solution is characterized by a K-angle < 15° (in the fourth column of Fig. 3). We find that the first guess, which resulted in keeping the solutions 2, 3, 6, and 7 was correct, the solutions indeed are close to the reference solution, and, vice-versa, no solution was lost. In the second and third guess we eliminated solutions 6, 7, and even 3, which, however, were close to the reference solution, so we lost them because the adopted criteria were too strict. We should include here a remark about solution 4, rejected due to its large CN. Such a rejection was very important because, although belonging to Group A and having acceptable STVAR , the best-fit solution is very far from the reference one, as seen in Figure 4 and in the fourth column of Figure 3 (K-angle  45°). The large CN in this case (and analogously in solutions 5 and 9) is due to the small epicentral distance of the involved stations. The MT solution is vulnerable to small errors of any kind, for example, the epicenter position. This example shows that in unfavorable situations even the solutions with acceptable nodal lines concentrated around the best-fit solution may still be incorrect. This example demonstrated the usefulness of quality measures, but it has also shown that a more extended practical July/August 2013

testing is needed to set up criteria for distinguishing between potentially reliable and unreliable solutions. Example 3: Correlation Plots and Histograms of the K-angle To better analyze FMVAR and STVAR , Figure 5a–f shows (a, c, e) complete space–time correlation plots and (b, d, f ) histograms of K-angle and DC-percentage. Three cases of Figure 3 are selected, Cases (a, b) 2, (c, d) 4, and (e, f ) 6. For comparison, the optimal correlation is normalized to unity, and the scale is adjusted to easily identify the acceptable correlation range (0.9, 1.0). True maxima of the correlation can be obtained from the VR values in Figure 3 (VR equals the correlation squared). (a, b) Case 2 is a favorable situation with a perfect space–time resolution and a small variability of the mechanism. (c, d) Case 4 has a remarkably worse space–time resolution and shifts the formally best-fitting solution to the depth of 12 km, that is, 4 km below the reference solution. The space–time elongation of the correlation strip is unrelated to CN, as CN just refers to the moment tensor, not to the space and time coordinates. The CN value is large, none of the acceptable mechanisms (with the correlation between 0.9 and 1.0) are correct, but their variability is only modest. Note that, besides the main strip in Case 4, we also find a secondary strip situated in the right top corner of (c) the correlation plot. However, solutions on the secondary strip have a polarity opposed to those on the main strip. The solutions on the secondary strip yield a few large values in the K-angle histogram, in Figure 5d, seen also as a mixture of the P and T axes for Case 4 in Figure 4. This effect is more common in relatively narrow frequency ranges and high-frequency inversions. In such cases, the user should check the MT solutions against the first-motion polarities in at least a few stations (this can also be easily performed in ISOLA). (e, f ) Case 6 contains the correct mechanism among the accepted solutions, but the space–time variation is even larger than in Case 4. As the source inversion is also dependent on the assumed crustal structure, calculations should be repeated in ISOLA for a suite of models available for the investigated region. Effects of the velocity model upon the correlation plots, FMVAR and STVAR are illustrated for Case 2 in Ⓔ Figure S4 (see supplement). In the studied frequency range (below 0.16 Hz) the structural effect is marginal. Theoretical Uncertainty Without Waveforms In the preceding section we analyzed the moment tensor uncertainty by means of variations in solution detected by grid searching in space and time. In the most favorable cases, the focal mechanisms have a simple statistical distribution for which the mean is close to the best-fit solution. In such cases (and only in such cases), the solution uncertainty might also be quantified by means of the correlation matrix (numerically calculated from the acceptable solutions of a given correlation threshold). The error ellipsoid can then be constructed, and the confidence intervals can be calculated in 8D space (6 components of the moment tensor, plus depth, and time). ISOLA

provides a simple tool of a similar physical meaning; that is, an uncertainty analysis of the focal mechanism by means of a 6D theoretical error ellipsoid, calculated for a fixed source position and time. When inverting for the moment tensor with a given (assumed) space and time position of the source, the inverse problem is linear and it is easy to calculate a theoretical uncertainty (Press, 1997). The MT uncertainties reported by some agencies rely on this concept. In practice, however, users are mostly interested in the uncertainty of strike, dip, and rake, which we address below. Similarly to the well-known problem of the hypocenter location, for which a 3D error ellipsoid is constructed by means of the eigenvectors and eigenvalues of the involved design matrix, here we compute a 6D error ellipsoid (because the full MT inversion has 6 model parameters), which reduces to 5D in the case of deviatoric MT inversions. The ellipsoid is constructed numerically by regularly sampling the parameter space around an assumed moment tensor. The ellipsoid shape and orientation are fully determined by the source-station configuration, frequency range, and the crustal model used. The absolute size of the error ellipsoid is codetermined by the assumed variance of the data (waveforms) σ 2 , for which σ is the standard deviation of data. The data (waveform) variance is rarely known, but this problem is not critical if we make only relative assessments of the uncertainty, for example, when comparing solutions of the same earthquake in two different networks. Then we may assume that data error σ is proportional to a simple order-of-magnitude estimate of the data (using equation 13 of Zahradník and Custódio, 2012), σ  Ce GM 0 , where C is a constant, e G is the peak displacement of the medium in response to a unit-moment source dislocation (in the studied source-station configuration and frequency band), and M 0 is the seismic moment. Finally, having numerically sampled the interior of the error ellipsoid, that is, obtained a statistical set of the moment tensors, we map them into strike, dip, and rake angles, and construct their histograms. The new version of ISOLA offers a special tool for this purpose. After calculating Green’s functions and performing the MT solution, it is possible to compute an estimate of σ 2 , controlled by the user’s choice of constant C. Reasonable results are often found using C  1. Alternatively, the user can prescribe his or her preferred value of σ 2 . For example, σ 2 may be formally assumed as the posterior misfit Σo − s2 =ndf (equation 5.30 of Shearer, 2009). Here ndf is the number of degrees of freedom, which equals the number of independent data minus the number of the inverted parameters. This ISOLA tool automatically makes all computations associated with eigenvectors and error ellipsoid, and at the end plots histograms of the source angles and nodal lines. The K-angle is used once again to indicate discrepancies with respect to the best-fit solution. Example 4: Theoretical Uncertainty Without Waveforms This method is applied to the best-fit solution (position, time, and moment tensor) of Case 2, without using any seismograms. ISOLA takes the source parameters and the station

Seismological Research Letters

Volume 84, Number 4 July/August 2013

661

▴ Figure 5. (a–f) Variability of the focal mechanisms (beachballs) given different space and time position for Cases 2, 4, and 6 (a, b, c, d, and e, f, respectively). (a, c, e) Space–time correlation plots, with maximum correlation normalized to unity. The acceptable solutions are those with normalized correlation between 0.9 and 1.0. Vertical axis refers to the trial source depth. Horizontal axis refers to the temporal grid search between −3 and 3 seconds with respect to origin time. The largest beachball is the best-fit solution. (b, d, f) Histograms of the K-angle and the DC-percentage expressing variability of the acceptable solutions in (a, c, e; with respect to the best-fitting mechanism). Color version of this plot is available as Ⓔ Figure S3 (see supplement).

662

Seismological Research Letters

Volume 84, Number 4

July/August 2013

▴ Figure 6. (a–f) Uncertainty of source angles for Case 2 (theoretical assessment for a fixed depth). (a–c) Histograms of strike, dip, and rake angles, respectively. The optimal solution is shown by thick line, whereas the dotted-line vertical strips mark the histogram width, interactively determined by the user. (d,f) K-angles and nodal lines associated with the source-angle histograms. The optimal solution is plotted using white nodal lines. (e) The uncertainty of the DC-percentage is also shown.

distribution, evaluates the Green functions and provides e G  2:949 × 10−20 N −1 . Choosing C  1, we get σ 2  1:2 × 10−8 m2 . The code then provides graphic output shown in Figure 6a–f. Jackknifing In the theoretical uncertainty assessment, we studied an effect of fictitious data errors upon model parameters. In a real inversion the data contains noise and the Green’s functions are imperfect. These are true errors and we can partly analyze their effect upon the source parameters if we invert an artificial, systematically changed, data set. This is the idea behind jackknifing. Jackknifing consists of performing MT inversions of real waveforms in which we repeatedly remove the same amount of data, for example, one station component or one complete station. Here we demonstrate systematic removal of a single station component. For each removal we perform a complete space–time grid search and least-squares MT calculations. The new version of ISOLA computes the repeated inversions fully automatically and reports the solution family by means of nodal lines and histograms. Histograms are constructed for the strike, dip, and rake angles, and also for the source position and time. K-angle is calculated to characterize deviation of each individual solution with respect to the best-fitting allcomponent solution.

Example 5: Jackknifing Figure 7a–h shows the jackknifing result for Case 2 of Figure 3. Although the uncertainty assessment of the position and time is an advantage over the theoretical analysis (fixed source depth), the source depth is constrained to fall within a narrow range (7 or 8 km). In particular, this range is much narrower than that found with using the correlation threshold in the grid search (2–15 km, see Fig. 5a). The same holds true for the source time and focal mechanism. In this sense, jackknifing has a clear tendency to reduce the potential uncertainty. This happens because jackknifing uncertainties are based exclusively on the parameter variation that results from that part of the data and modeling error which differs from one reduced data subset to the other, which are generally small. In practice, jackknifing may sometimes display a larger parameter variation than shown in the presented example. This may happen in situations for which we use a 1D crustal model, but the stations are situated in two or more distinctly different geological provinces, for example, inside and outside a sedimentary basin. Then, if the 1D model does not account for the basin structure, model errors for stations inside the basin are large, hence keeping or removing the basin stations when varying the data set may yield a large variation of the model parameters. However, in such situation the user would probably prefer to use station-specific crustal models, which is a possibility offered by ISOLA. Jackknifing may also yield large solution uncertainties if gross data problems exist at some stations, for

Seismological Research Letters

Volume 84, Number 4 July/August 2013

663

(b)

(c)

20

15

15

15

15

10

10

50 100 150 200 250 300 350

0 10 20 30 40 50 60 70 80 90

Dip

Strike

(f)

25

0 0

60

120 180

2

4

6

Rake

(g)

15

10

5

0 −180 −120 −60

0 0

10

5

5

0

Count

20

Count

20

5

(e)

(d)

20

Count

Count

(a)

8 10 12 14 16 18 20

Source position

(h)

10

20

10

Count

Count

Count

10 15

5

5 5 0

0 −3

−2

−1

0

1

2

3

Source time (s)

0 0 10 20 30 40 50 60 70 80 90 100

0

2

4

6

8

10

Kagan angle

DC%



Figure 7. (a–h) Uncertainty assessment of the source angles, source position, time, and DC-percentage for Case 2, using jackknifing (removing single-station components). The nodal lines and K-angles with respect to the solution with all stations and components used (thick nodal lines) are also shown.

example, presence of a long-period disturbance due to tilt or clip (e.g., Zahradník and Plešinger, 2005, 2010). Implementation of the automatic detection and removal of the disturbances is under way.

CONCLUSION This paper is a follow-up of Sokos and Zahradník (2008), in which the basic version of the ISOLA software was described. The present paper describes main recent innovations that allow the assessment of solution quality. Different stations may be used in the moment tensor (MT) inversion with a different frequency range. The low-frequency range suitable for the MT inversion is determined by analyzing the signal and noise spectra. The solution quality is measured based on the signal-tonoise ratio, variance reduction, condition number, and two new indices, FMVAR and STVAR . FMVAR (the focal-mechanism variability index) quantifies the stability of the focal mechanism in a space-time grid search within a prescribed correlation threshold. STVAR (space–time variability index) measures the stability of the source position and time in a grid search within the same correlation threshold. The uncertainty analysis is also provided by means of a 6D theoretical error ellipsoid, suitable mainly when designing a seismic network, when real waveforms are not available. Automated jackknifing is a possible complement of the standard MT inversion, in which real data and modeling errors are partly taken into account. The main innovations are illustrated by five examples. All of them concern the same (M w 4.3) earthquake, processed in 664

Seismological Research Letters

Volume 84, Number 4

ten different ways (Cases 1 to 10 in Fig. 3). Various scenarios have been simulated, in which just a few stations were available and the goal was to assess the reliability of the solution, while keeping the reference solution blind. Many of the solutions were successfully classified as unreliable because of their high condition number, large FMVAR or STVAR index. Successfully rejected in the blind simulation was also the most dangerous Case 4 of Figure 4, in which the mean of the accepted family of nodal lines coincided with the best-fit solution, but both were far from the correct solution. The latter was the most alarming case, showing that sometimes the true solution may be well beyond the estimated uncertainty limits. Although ISOLA still always gives some focal mechanism, we believe that with the current built-in uncertainty (quality) measures, the software provides enough warning against physically meaningless results. The software is under continued development and can be downloaded at http://seismo.geology .upatras.gr/isola/. Although ISOLA was developed under Windows, Linux and MacOS versions are now also available. Finite-extent source applications using ISOLA will be described elsewhere. Future work is planned towards speedup of the FORTRAN codes, automation of the procedures, and inclusion of complete Green’s function for 3D crustal models from (external) 3D codes.

ACKNOWLEDGMENTS Financial support was partly obtained from the grant agencies in the Czech Republic: GAČR 210/11/0854 and MSM July/August 2013

0021620860. The authors were also supported by the ERASMUS Greek–Czech bilateral agreement on the Staff Training Mobility project. We thank Ronnie Quintero who organized the ISOLA training course in Costa Rica for 15 participants from Central and South America, co-sponsored by IASPEI (3– 9 September 2011). The GMT software of Wessel and Smith (1998) was used in all figures of this paper. Susana Custódio helped with MacOS version, corrected English, and provided useful comments. Waveform data provided by the Hellenic Unified Seismic Network and Corinth Rift Laboratory Network (station TRIZ) are greatly appreciated.

REFERENCES Benetatos, C., J. Málek, and F. Verga (2013). Moment tensor inversion for two micro-earthquakes occurring inside the Háje gas storage facilities, Czech Republic, J. Seismol. 17, 557–577, doi: 10.1007/ s10950-012-9337-0. Bouchon, M. (1981). A simple method to calculate Green’s functions for elastic layered media, Bull. Seismol. Soc. Am. 71, 959–971. Cannata, A., P. Montalto, E. Privitera, and G. Russo (2009). Characterization and location of infrasonic sources in active volcanoes: Mount Etna, September–November 2007, J. Geophys. Res. 114, B08308, doi: 10.1029/2008jb006007. Cesca, S., S. Heimann, K. Stammler, and T. Dahm (2010). Automated procedure for point and kinematic source inversion at regional distances, J. Geophys. Res. 115, B06304, doi: 10.1029/2009JB006450. Coutant, O. (1989). Program of numerical simulation AXITRA, Tech. rep., LGIT, Grenoble, France (in French). Delouis, B., C. Charlety, and M. Vallée (2008). Fast determination of earthquake source parameters from strong motion records: Mw, focal mechanism, and slip distribution, in EGU General Assembly, Geophysical Research Abstracts 10, abstract 04939. Fojtíková, L., V. Vavryčuk, A. Cipciar, and J. Madarás (2010). Focal mechanisms of micro-earthquakes in the Dobrá Voda seismoactive area in the Malé Karpaty Mts. (Little Carpathians), Slovakia, Tectonophysics 492, 213–229. Foulger, G. R., B. R. Julian, D. P. Hill, A. M. Pitt, P. E. Malin, and E. Shalev (2004). Non-double-couple microearthquakes at Long Valley caldera, California, provide evidence for hydraulic fracturing, J. Volcanol. Geoth. Res. 132, 45–71, doi: 10.1016/s0377-0273(03) 00420-7. Gallovič, F., and J. Zahradník (2011). Toward understanding slipinversion uncertainty and artifacts II: singular value analysis, J. Geophys. Res. 116, B02309, doi: 10.1029/2010JB007814. Gallovič, F., and J. Zahradník (2012). Complexity of the Mw 6.3 2009 L’Aquila (Central Italy) earthquake: 1. Multiple finite-extent source inversion, J. Geophys. Res. 117, B04307, doi: 10.1029/ 2011JB008709. Gallovič, F., J. Zahradník, D. Křížová, V. Plicka, E. Sokos, A. Serpetsidaki, and G.-A. Tselentis (2009). From earthquake centroid to spatialtemporal rupture evolution: Mw 6.3 Movri Mountain earthquake, June 8, 2008, Greece, Geophys. Res. Lett. 36, L21310, doi: 10.1029/ 2009GL040283. Kagan, Y. Y. (1991). 3-D rotation of double-couple earthquake sources, Geophys. J. Int. 106, no. 3, 709–716, doi: 10.1111/j.1365246X.1991.tb06343.x. Kikuchi, M., and H. Kanamori (1991). Inversion of complex body waves, III, Bull. Seismol. Soc. Am. 81, 2335–2350. Křížová, D., J. Zahradník, and A. Kiratzi (2013). Resolvability of isotropic component in regional seismic moment tensor inversion, Bull. Seismol. Soc. Am. 103, no. 4 (in press).

Maercklin, N., A. Zollo, A. Orefice, G. Festa, A. Emolo, R. De Matteis, B. Delouis, and A. Bobbio (2011). The effectiveness of a distant accelerometer array to compute seismic source parameters: The April 2009 L’Aquila earthquake case history, Bull. Seismol. Soc. Am. 101, no. 1, 354–365, doi: 10.1785/0120100124. Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery (1997). Numerical Recipes in Fortran 77: The Art of Scientific Computing, Second Ed., Cambridge University Press, 992 pp. Shearer, P. (2009). Introduction to Seismology, Second Ed., Cambridge University Press, 396 pp. Sokos, E. N., and J. Zahradník (2008). ISOLA a FORTRAN code and a MATLAB GUI to perform multiple-point source inversion of seismic data, Comput. Geosci. 34, 967–977. Sokos, E. N., J. Zahradník, A. Kiratzi, J. Janský, F. Gallovič, O. Novotný, J. Kostelecký, A. Serpetsidaki, and G.-A. Tselentis (2012). The January 2010 Efpalio earthquake sequence in the western Corinth Gulf (Greece), Tectonophysics 530–531, 299–309. Valentine, A. P., and J. Trampert (2012). Assessing the uncertainties on seismic source parameters: Towards realistic error estimates for centroid-moment-tensor determinations, Phys. Earth Planet. In. 210–211, 36–49. Wessel, P., and W. H. F. Smith (1998). New, improved version of generic mapping tools released, Eos Trans. AGU 79, 579. Zahradník, J., and S. Custódio (2012). Moment tensor resolvability: Application to southwest Iberia, Bull. Seismol. Soc. Am. 102, no. 3, 1235–1254, doi: 10.1785/0120110216. Zahradník, J., and F. Gallovič (2010). Toward understanding slip inversion uncertainty and artifacts, J. Geophys. Res. 115, no. B9, B09310, doi: 10.1029/2010JB007414. Zahradník, J., and A. Plešinger (2005). Long-period pulses in broadband records of near earthquakes, Bull. Seismol. Soc. Am. 95, no. 5, 1928– 1939, doi: 10.1785/0120040210. Zahradník, J., and A. Plešinger (2010). Toward understanding subtle instrumentation effects associated with weak seismic events in the near field, Bull. Seismol. Soc. Am. 100, no. 1, 59–73, doi: 10.1785/ 0120090087. Zahradník, J., and E. Sokos (2011). Multiple-point source solution of the M w 7.2 Van earthquake, October 23, 2011, Eastern Turkey, Report to European Mediterranean Seismological Center, http://www .emsc‑csem.org/Files/event/239856/Van_ISOLA.pdf (last accessed March 2013). Zahradník, J., F. Gallovič, E. Sokos, and G.-A. Tselentis (2011). Preliminary slip model of M 9 Tohoku earthquake from strong-motion stations in Japan: An extreme application of ISOLA code, Report to European Mediterranean Seismological Center, http://www .emsc‑csem.org/Files/event/211414/ISOLA_Report_tohoku.pdf (last accessed March 2013).

Seismological Research Letters

Efthimios Sokos University of Patras Department of Geology Seismological Laboratory 26504 Patras, Greece [email protected]

Jiří Zahradník Charles University in Prague Faculty of Mathematics and Physics 18000 Prague, Czech Republic [email protected]

Volume 84, Number 4 July/August 2013

665