Teleseismic inversion for the second degree moments of earthquake space–time distributions

Summary

Determining detailed descriptions of the space–time distribution of an earthquake's rupture using only teleseismic data is often a difficult and non unique process. However, a description of an earthquake that goes beyond a point source is often needed for many seismotectonic problems. We have developed a method that inverts measurements of the shifts in traveltime and amplitude of various seismic phases at global stations for the second degree polynomial moments of an event's space–time distribution. Our values for the second moments of the 1995 M_w 8.0 Jalisco, Mexico, thrust event correspond to a characteristic rupture length of 121±10 km along strike, a characteristic rupture width of 76±11 km downdip, a characteristic duration of 29±2 s, and a direction of rupture propagation along the strike of the subduction zone to the NW. For the 1995 M_w 7.2 Gulf of Aqaba earthquake, our estimates of the second moments correspond to a characteristic rupture length of 53±2 km, a characteristic duration of 12±1 s, and rupture propagation to the NE along the strike of the Dead Sea Transform. Our estimates of the first and second moments agree well with values from local seismic and geodetic networks. For both the Jalisco and Aqaba events, we are able to resolve the fault plane ambiguity associated with the events' moment tensors using our estimates of the second spatial moment.

Introduction

The spatial and temporal distributions of earthquake moment release have a large effect on the displacement field observed at locations near the fault. This strong dependence has enabled seismologists to image the space–time evolution of the rupture of earthquakes recorded by nearby instruments in remarkable detail (Wald & Heaton 1994; Ide et al. 1996). However, abundant near field data have generally been lacking except in regions of dense regional seismic networks, such as California and Japan. This has motivated a large amount of work on imaging earthquake ruptures from teleseismic data. Such far away observations have limited resolution; however, Ihmlé (1998) recently showed that P wave inversions for the M_w 8.2 Bolivia deep earthquake, which had good teleseismic coverage and relatively simple P wave propagation, were unable to differentiate between a slip distribution composed of compact subevents and one composed of widely distributed slip. This inherent non uniqueness is often associated with a choice of the spatial smoothing parameter that is used to regularize these inversions (Ihmlé 1998). Imaging the slip distribution of large earthquakes is further complicated because P waves provide only weak constraints on the total moment of an event (Ekström 1989). Thus, attempts to determine the details of slip distributions for earthquakes that are only recorded teleseismically are at best extremely non unique. There are several important properties of an earthquake's space–time behaviour that do not require this type of detailed finite fault inversion, such as resolution of the fault plane ambiguity, approximate extent of the rupture area, duration and direction of rupture propagation. In this paper, we develop a method for robustly estimating these large scale features of the slip distribution from teleseismic data.

General description of seismic sources

In the normal mode formulation, the spatial and temporal extent of an earthquake's rupture enters the equation for the displacement u at a position x through the stress glut rate tensor Γ˙(r, t) (Gilbert 1970; Dahlen & Tromp 1998):

where

and s_k is the eigenvector of the kth mode with eigenfrequency ω_k, attenuation parameter α_k and associated strain tensor E_k. Γ˙(r, t), introduced and described in detail by Backus & Mulcahy (1976a,b), is non zero only in the region of non elastic deformation, i.e. the source volume. In this paper, we use a dot to represent the time derivative and assume that Γ˙(r, t) has the form

where M̂ is a unit norm moment tensor and f˙(r, t) is the scalar source space–time function. For slip on a simple fault, this assumption is equivalent to requiring a constant orientation of the fault plane and slip vector, but it also allows slip on parallel faults. In contrast to near field inversions that seek to determine a detailed parametrization of f˙(r, t), we will invert for the polynomial moments of f˙(r, t), μ^(i,j), which are tensors of order i and have spatial degree i, temporal degree j and total degree i+j. The zeroth degree moment

is the standard seismic moment, M₀, while the normalized first degree moments μ^(1,0)/M₀ and μ^(0,1)/M₀ give the spatial and temporal centroids, respectively:

The degree 2 moments about a point r₀ and time t₀ are called central moments when r₀=μ^(1,0)/M₀ and t₀=μ^(0,1)/M₀:

The coordinate system axes (r, θ, φ) correspond to radius, colatitude and longitude respectively. If the degree 2 moments are calculated about a point that is slightly offset from the true centroid location, there is a small difference between the second moment and the second central moment:

The second central moments describe the spatial and temporal extent of the source. Following Backus (1977a) and Silver & Jordan (1983), we can define the characteristic dimension, x_c(n̂), of the source region in a direction n̂, the characteristic duration τ_c, the characteristic (or apparent) rupture velocity v_c and the average velocity of the instantaneous spatial centroid, v₀:

where L_c is the maximum value of x_c, that is, the characteristic rupture length along the direction of the eigenvector associated with the largest eigenvalue of μ̂^(2,0).

A number of theoretical papers on second moments have examined the non uniqueness in source inversions and the relationships between second moments and standard source models (Backus 1977b; Okal 1982; Doornbos 1982a,b; Pavlov 1994; Das & Kostrov 1997), but prior observational work on determining the polynomial moments of earthquakes has been limited to the inversion of body waveform data (Gusev & Pavlov 1988) or surface wave amplitude spectra assuming a fault plane a priori (Bukchin 1995). In this paper, we develop a formalism for inverting frequency dependent traveltime and amplitude measurements for the source moments of polynomial degrees 0, 1 and 2. Data are obtained for arbitrary wave groups, such as surface waves and body waves, as differential measurements relative to synthetics computed for a point source in a 3 D earth model.

2010 Frequency-dependent phase and amplitude data

We measure the phase and amplitude of a target wave group relative to a synthetic seismogram using the GSDF technique of Gee & Jordan (1992). The synthetic seismogram, s˜, is represented as a sum over travelling wave groups, indexed by mode number n. In the frequency domain, it is assumed that the data seismogram, s, can be related to the synthetic by small differences in phase delay, δτ_pⁿ, and amplitude reduction time, δτ_qⁿ:

GSDF measurements of the generalized delay time, δτ_pⁿ(ω), have been utilized in a variety of structural inverse problems (Gaherty et al. 1996; Katzman et al. 1998) because the relationships between δτ_pⁿ(ω) and structural model parameters are easily linearized. In those studies, synthetic seismograms were calculated for short duration earthquakes (M_w∼6) using a 1 D earth model, and the phase delays of the data were attributed to structural variations along the path between the source and station. Here, we study large earthquakes (M_w≥7.0) that cannot be represented adequately as point sources, even at low frequencies. We make approximate corrections for Earth structure using a 3 D earth model, so that our measurements of δτ_pⁿ and δτ_qⁿ can be attributed primarily to properties of the seismic source.

The synthetic seismograms in this study were calculated for a point source in both space and time through normal mode summation for the PREM earth model (Dziewonski & Anderson 1981), corrected for 3 D elastic structure using the degree 12 aspherical model of Su et al. (1994) and the asymptotic approximations of Woodhouse & Dziewonski (1984). We also corrected fundamental modes above 7 mHz for smaller scale heterogeneity using the degree 40 phase velocity maps of Ekström et al. (1997). The 3 D corrections only offset the phase of the surface wave arrivals; no corrections are made to the amplitudes for the effects of 3 D anelastic structure or focusing and defocusing. The source is specified by a centroid location, centroid time and moment tensor from the Harvard CMT catalogue (Dziewonski et al. 1997). To limit any contributions to our measurements by lateral heterogeneity unaccounted for by the 3 D models, we restricted our measurements of surface waves to low frequencies (≤20 mHz), and we omitted any stations for which our synthetic seismograms did not closely resemble the observed seismograms.

The first step in the measurement process was the creation of an isolation filter by a weighted summation of normal modes. The isolation filter was designed to approximate the complete synthetic in the time interval containing a target wave group and decay to zero outside that interval (Fig. 1). The isolation filter was then cross correlated with both the complete synthetic and the data seismograms. The resulting cross correlagrams were windowed to localize the signals in the time domain and then put through a series of narrow band filters (with centre frequencies ω_i) to localize them in the frequency domain. The windowed and filtered cross correlagrams were then fitted with two five parameter Gaussian wavelets from which the phase and amplitude parameters, δτ_p(ω_i) and δτ_q(ω_i), were calculated (Gee & Jordan 1992). This procedure corrects for the interference between the phase targeted by the isolation filter and other arrivals on the seismogram, as well as the biases introduced by the windowing and filtering operations. For shallow earthquakes, we typically made measurements on vertical component seismograms for fundamental mode Rayleigh waves at frequencies between 3 and 20 mHz and for P waves between 10 and 50 mHz.

Effects of an extended source on GSDF measurements

To relate the differences between the point source synthetic and data seismograms to our source parametrization, we determined the partial derivatives of δτ_pⁿ and δτ_qⁿ with respect to the polynomial moments of f˙(r, t). A Taylor series expansion truncated to second order gives the expression relating a point source seismogram, s˜, to the seismogram resulting from an extended source, s (e.g. Dahlen & Tromp 1998):

where r_R and r_S are the receiver and point source locations respectively, t_S is the point source time, and μ^(1,0) and μ^(0,1) represent any deviation of the centroid (r₀, t₀) from the point source location. The operator ∇_S acts on the source spatial coordinates (r_S, θ_S, φ_S) only. In the GSDF formalism, s˜ and s are represented as sums over travelling modes:

where we have dropped the argument r_R for convenience. To obtain the explicit expression for s˜, we converted the expression for the synthetic seismogram from a sum over normal modes to a sum over travelling modes using the Poisson sum formula (Appendix A):

where Δ is the epicentral distance and A_n(ω) and B_n(ω) depend on the eigenfunctions of the earth model, the assumed centroid moment tensor and the epicentral distance. The perturbations to the continuous wavenumber, δλ_n, account for the corrections made to the synthetics for 3 D Earth structure (Appendix A). Eqs (12) and (14) yield expressions for τ_pⁿ(ω) and τ_qⁿ(ω):

Substituting eqs (12) and (13) into eq. (11) and making the low frequency approximation,

gives the following formulae for δτ_pⁿ(ω) and δτ_qⁿ(ω):

Expressions for the terms involving τ_pⁿ and τ_qⁿ are calculated from eqs (15) and (16). The partial derivatives of δτ_pⁿ and δτ_qⁿ with respect to the 14 independent elements of μ^(0,1), μ^(1,0), μ^(0,2), μ^(1,1) and μ^(2,0) can be simply read off eqs (18) and (19). We have also calculated the partial derivatives of the GSDF measurements with respect to changes in the moment tensor elements, M_ij, using the expressions

In practice, isolation filters for most phases involve contributions from a large number of normal modes ranging over multiple values of the overtone number n. Thus the measured values of δτ_p and δτ_q will be a weighted sum over the δτ_pⁿ and δτ_qⁿ for the individual branches that contribute to the isolation filter. The partial derivative of a GSDF measurement can also be calculated by a weighted sum of the partial derivatives of the individual branches (Gee & Jordan 1992, eqs 94–95). Fig. 1 shows examples of the isolation filters for the P wave and the Rayleigh wave from a shallow source along with the normal mode weights used to calculate the filters. Only the fundamental mode branch is needed to match the Rayleigh wave, while contributions from several branches are needed to match the P wave. While all of the GSDF measurements in this paper were made on vertical component records of either fundamental mode Rayleigh waves or P waves, the technique is completely generalizable to an arbitrary waveform. Thus, intermediate and deep earthquakes can be studied by the inclusion of higher mode surface waves. Additionally, second orbit surface waves can be used to increase azimuthal coverage of the source (Appendix A). In the future, we will incorporate Love and other S_H waves to improve our resolution of the higher moments.

Examples of partial derivatives

Several intuitive properties can be seen in the partial derivatives (Fig. 2) for the Rayleigh wave isolation filter shown in Fig. 1. The effect of a shift in centroid time, μ^(0,1), on the frequency dependent arrival times δτ_p(ω) is independent of azimuth, and this type of shift has no effect on the frequency dependent amplitude parameter δτ_q(ω). ∂δτ_p/∂μ^(1,0) shows a cosine dependence on azimuth with a maximum at 270°. This pattern corresponds to late arrivals at stations to the west and early arrivals at stations to the east owing to an eastward shift in the centroid location. A similar pattern, shifted by 90°, is found for ∂δτ_p/∂μ^(1,0). These traveltime shifts due to a mislocated r_S are independent of frequency. The centroid depth, μ_r^(1,0), is the only component of the first degree moment that significantly perturbs the amplitude reduction time, δτ_q. In contrast, the degree 2 moments primarily perturb the δτ_q measurements, having almost no influence on δτ_p. ∂δτ_q/∂μ̂^(1,1) approximates a cosine in azimuth and increases in amplitude with frequency. This behaviour is just the directivity pattern of Ben Menahem (1961), which has been used in many studies to infer directivity from surface waves. The φφ component of μ̂^(2,0) is non zero for sources with finite extent in the east–west direction. ∂δτ_q/∂μ̂^(2,0) has positive values (i.e. decreased data amplitudes) at stations to both the east and west due to destructive interference (Fig. 3). ∂δτ_q/∂μ̂^(0,2) is independent of azimuth, but increases in amplitude with frequency, corresponding to the spectral roll off expected for a source of finite duration.

Fig. 3 shows the GSDF measurements that result from two simple models of an extended earthquake source, unilateral and bilateral rupture. In these calculations, the centroid and moment tensor are known exactly, and f˙(r, t) is symmetric about the centroid, so the dominant contribution comes from the second moments. In both cases, almost no deviation from the point source synthetic is observed at low frequencies, but a characteristic pattern emerges in the δτ_q measurements as frequency increases. For unilateral, westward rupture (Fig. 3a), the effects of μ̂^(1,1) and μ̂^(2,0) approximately cancel in the direction of rupture propagation, while they add constructively in the opposite direction. Thus a ‘half cosine’ pattern emerges in the δτ_q(ω) measurements, with the smallest data amplitudes occurring 180° away from the direction of rupture propagation. We observe this azimuthal pattern, which is diagnostic of a predominately unilateral rupture, for both the 1995 Jalisco and Gulf of Aqaba earthquakes, which are discussed in the Results section. In the bilateral case (Fig. 3b), μ̂^(1,1) is zero, so the μ̂^(2,0) term dominates, producing small amplitudes in both directions along the strike of the fault. For both of these examples, the finite duration of the source shifts the δτ_q measurements to positive values with increasing frequency, owing to spectral roll off.

Directivity effects are often assumed to average to zero in source time function and source spectrum estimation techniques if sufficient azimuthal coverage is available. The amplitude perturbations owing to finite source effects shown in Fig. 3 have only positive values of δτ_q, indicating a reduction in the data amplitude relative to the point source synthetic. While the perturbations due to the propagation of the centroid, ∂δτ_q/∂μ^(1,1) (Figs 2f and g), average to zero as a function of azimuth, the perturbations owing to μ̂^(2,0) do not (Figs 2f and g). The amplitude increases in the direction of rupture propagation resulting from μ̂^(1,1) are cancelled in the unilateral case by the decreases resulting from μ̂^(2,0) (Fig. 3a). In the perfectly bilateral case (Fig. 3b), μ̂^(1,1)=0, so only the amplitude reductions resulting from μ̂^(2,0) occur. Thus techniques that assume that amplitude effects resulting from directivity average to zero will underestimate the earthquake's seismic moment.

Inversion method

Optimization technique

We can formulate a linear inverse problem to minimize ||Am−b|| for a vector m, which contains the 14 independent elements of μ^(0,1), μ^(1,0), μ^(0,2), μ^(1,1) and μ^(2,0) as well as the six independent elements of M. The data vector, b, contains the GSDF measurements and A contains the partial derivatives calculated using eqs (18)–(21). For sources that obey the assumption in eq. (3), the moments of f˙(r, t) are uniquely determined by the displacement field (Bukchin 1995), making the linear inverse problem theoretically solvable. In practice, however, owing to errors in the assumed moment tensor, the effects of unmodelled lateral heterogeneity and the small relative size of the minimum eigenvalue of μ^(2,0), a linear least squares inversion does not result in a robust solution. Das & Kostrov (1997) incorporated linear inequality constraints among the moments to stabilize their inversion; however, their constraints, which are based on the Hausdorff inequalities, required a priori information about the maximum size of the source region. The power of these constraints depended strongly on the a priori limits on the source region size (Das & Kostrov 1997). To avoid setting a priori limits on the extent of the source region, we stabilized our inversion by incorporating the physical constraint that the 4 D source region have non negative volume. This constraint, which is an inherent property of all earthquakes, is equivalent to requiring the 4×4 matrix

to be positive semi definite. This constraint is non linear in the elements of m and encompasses the constraints placed on the second degree moments of f˙(r, t) by other authors through the use of Bessel inequalities (Bukchin 1995, see also Appendix B). Methods for solving linear problems subject to matrix inequality constraints such as (22) have recently been developed in the optimization literature. We use the semi definite programming approach of Vandenberghe & Boyd (1996). To put our problem into this form, it is necessary to restate the non linear (least squares) objective function as a linear objective function subject to a non linear (quadratic) constraint. We introduce a dummy variable, c, such that ||Am−b||≤c. Then our problem can be restated as

where I^N is the identity matrix with dimension equal to the number of GSDF measurements, N, and ≥0 indicates that the matrix is positive semi definite. The equivalence between the linear least squares problem and (22) can be seen by calculating the eigenvalues of the N+1×N+1 matrix in (23), which are non negative when the matrix is positive semi definite. This restatement of the problem is known as using Schur complements to represent a non linear constraint as a linear matrix inequality (Vandenberghe & Boyd 1996). The semi definite programming algorithm also allows us to incorporate a variety of linear inequality constraints, such as requiring the source to lie on one of the fault planes of the moment tensor or constraining the source to occur below the free surface, when appropriate for a particular event. The physical constraint (22) greatly reduces the non uniqueness in the minimization problem without requiring any a priori information about the extent of the source.

Any difference between the assumed point source location (r_S, t_S) and the true centroid location(r₀, t₀) will result in a slight difference between the estimated value of the second moments and the values of the second central moments (see eq. 8). Correcting the estimates of the second moments for a change in centroid location can violate the positive definite constraint (22). Thus we prefer to iterate the inversion (by updating the CMT parameters) until the correction term is negligible. Usually only one additional iteration is necessary for the elements of m corresponding to the changes in the CMT parameters to be negligible.

For events of moderate size, M_w≈7, we usually need to include high frequency (20–50 mHz) P wave measurements to resolve the second moments accurately. These inversions also require the inclusion of low frequency Rayleigh wave measurements to constrain the moment tensor. We have noticed that combining these two types of measurements often leads to smaller than expected values of μ^(0,2). This bias results from a small discrepancy (≤20 per cent) in the relative amplitudes of P and R1 waves in our synthetics as compared to the data. Iterating the centroid depth often removes most of this anomaly, but there appears to be a residual bias even after this process. To remove this bias from our inversions, we include a denuisance parameter that allows for a small, frequency independent excitation anomaly in the P wave measurements. The estimated value of this parameter usually corresponds to δτ_q values of 1–2 s, and hence is a smaller signal than the observed directivity.

Error analysis

We have also incorporated the non negative volume constraint (22) into our estimates of the errors in our solution. This is important because the optimal solution to a semi definite programming problem will often lie on the boundary of the feasible region defined by the matrix inequality constraints (Vandenberghe & Boyd 1996). Standard linear estimates do not account for the existence of an infeasible zone, but one method that can is the grouped jackknife (Tukey 1984). The jackknife technique relies on the assumption that different elements of a data set are independent. In our data sets of GSDF measurements, errors in the measurements at a given station are likely to be correlated from one frequency to the next. For instance, all the arrivals at one station may be early owing to fast velocity anomalies along the path. Thus, to provide a more conservative (larger) estimate of the variance of our model, we created M subsets of data by deleting all the measurements in each of M azimuthal bins (usually 12 bands of 30°) to create the ith subset. The jackknife method then gives an estimate of the covariance matrix, C_mm, of our model parameters:

where m̂_i is the estimate produced by the semi definite program technique using the ith subset of the data, and m̂_(.) is the average of the M values of m̂_i. The diagonal elements of C_mm give the variances of the individual model parameters.

The jackknife estimate of C_mm can also be used to determine the variances of derived quantities such as the characteristic rupture dimension, x_c, in a direction n̂. The equation for our estimate of x_c, eq. (9), can be rewritten in vector form as

where z is a function of n̂, and m̂ contains the elements of our estimate of μ̂^(2,0). Using standard error propagation equations (Bevington & Robinson 1992, p. 41), we can also estimate the variance in our estimate of x_c:

Results

1995 Kobe earthquake

The 1995 M_w 6.9 Kobe earthquake was recorded extremely well by the Japanese strong motion network, allowing for a good comparison of our teleseismic results with more accurate estimates determined from the local data. In particular, we use it to check for biases in the zeroth and first degree moments of f˙(r, t). For this event we initially calculated synthetics using the Harvard CMT (Dziewonski et al. 1996) (20:46:59.4, 34.78°N, 134.99°E, 20 km depth). We then made measurements, shown in Fig. 4, of the fundamental mode first orbit Rayleigh wave at 46 global stations for which our synthetics closely resembled the observed seismograms. The δτ_p(ω) measurements showed a pattern of late arrivals on the data seismograms at stations to the north and early data arrivals at stations to the south, while the δτ_q measurements were essentially a flat function of azimuth, indicating that the data amplitudes were slightly larger than the synthetics. The lack of a clear directivity signal in the amplitude measurements is because the measurements were made at wavelengths for which an event of this size is well represented by a point source.

Owing to the lack of a clear directivity signal in the amplitude measurements, we only inverted for changes to the first moments and moment tensor. The results of the inversion are given in Table 1 and Fig. 5, and the fit to the data is displayed in Fig. 4. The centroid location moved 24 km to the south, which fits the δτ_p pattern, and 5 km shallower. The changes in the moment tensor slightly increased in the seismic moment while retaining approximately the same focal mechanism. However, the changes to the M_rθ and M_rφ components are not well constrained. The new centroid location is about 10 km from the value determined by (Ide & Takeo 1997) (34.61°N, 135.075°E, 9.2 km) using local data and provides a better estimate than the starting Harvard CMT values (Fig. 5). We have calculated synthetic seismograms for the revised CMT and remeasured the data (Fig. 6). The results show essentially flat functions of azimuth at a value of zero for both δτ_p and δτ_q. Inversion of these remeasured data produces extremely small (≤∼1 km) shifts in the centroid location and correspondingly small shifts in the moment tensor, which are below our resolution limit. The residual scatter in Fig. 6, which has an average rms amplitude of 7 s depending on frequency band, is indicative of the amount of error in our measurements due to unmodelled propagation effects. This level of scatter is significantly smaller than the signal produced by directivity in large earthquakes, as will be shown in the next section.

1995 Jalisco earthquake

The 1995 M_w 8.0 Jalisco, Mexico, earthquake occurred on a shallow thrust plane from the subduction of the Rivera Plate under the North American Plate (Fig. 7). A good test of our method comes from comparing our estimate of μ̂^(2,0) with the slip distribution determined from a local GPS study by Melbourne et al. (1997). We performed an initial inversion for changes in the centroid and moment tensor. The new centroid location (15:36:28.80, 19.43°N, 104.92°W, 8.0 km) corrected the early R1 arrivals at stations to the southeast found for the Harvard CMT (Dziewonski et al. 1997). For most shallow sources, the M_rθ and M_rφ components were not well constrained, so we settled on a moment tensor (M_rr=4.6, M=−3.2, M=−1.3, M_rθ=5.0, M_rφ=−3.0, M=1.9×10²⁰ N m) through a process of forward modelling. Measurements for synthetics calculated from this CMT are shown in Fig. 8. The primary signals visible in the data are low R1 amplitudes at stations to the SE. The magnitudes of these signals increase with frequency, as expected for a finite source. The baseline level of the amplitude reduction times, δτ_q, also shifts upwards from about 8 s at 5 mHz to about 12 s at 10 mHz. This overall decrease in amplitude with frequency is indicative of the spectral roll off caused by the finite duration of the source time function. We inverted the data in Fig. 8 under the non negative volume constraint (22) and under the constraint that the characteristic length in the radial direction, x_c(r̂), be ≤10 km. This added constraint ensured that the source volume did not extend above the free surface.

Our inversion (Tables 1, 2 and 3) yielded a characteristic duration τ_c of 29±2 s, a centroid velocity v₀ of 3.7±0.2 km s⁻¹ at an azimuth of 307°, a characteristic rupture length L_c of 121±10 km along the strike of the subduction zone, a characteristic width w_c of 76±11 km downdip, and a characteristic rupture velocity v_c of 4.1±0.1 km s⁻¹. The three eigenvalues of μ̂^(2,0) specify the semi axis lengths of an ellipsoid whose axes are oriented in the direction of the associated eigenvectors. This characteristic source volume, shown for the Jalisco event in Fig. 9, is an estimate of the region that contributed significantly to the moment release. The source ellipsoid is nearly planar, with its longest dimension oriented NW–SE along the strike of the subduction zone. The ellipsoid dips slightly to the NE as expected for a shallow thrust plane. Fig. 10 compares our estimate of the characteristic source dimensions (solid ellipse) to the GPS slip inversion of Melbourne et al. (1997). Our centroid location agrees well with the value calculated from the slip distribution and the Harvard CMT. The NEIC epicentre is displaced about 80–100 km to the SE. The position of the centroid relative to the epicentre agrees with the direction of the characteristic velocity. We calculated the characteristic dimensions for Melbourne et al.'s distribution (dashed ellipse in Fig. 10), which was confined to a single plane (strike 302°, dip 16°). Its larger dimension is along strike, L_c=91 km, and its smaller dimension is oriented downdip, w_c=50 km. Thus, our estimate of μ̂^(2,0) is consistent with the local inversion.

To resolve the fault plane ambiguity for the Jalisco focal mechanism (Fig. 7), we calculated the characteristic rupture length in the direction normal to each of the candidate fault planes. For the shallowly dipping plane (strike 302°, dip 19°), x_c(n̂)=21±2 km and for the steeply dipping plane x_c(n̂) = 76±10 km. The shallowly dipping plane has a much smaller value of x_c(n̂), and we therefore infer that it is the true rupture plane. This agrees with Melbourne et al.'s study using local data. The small value of x_c(n̂) is due to the dip of the characteristic source volume to the northeast. This dip is a robust feature of the inversion. Fig. 11 shows histograms for each element of the first and second moments for the 12 different inversions carried out for the jackknife calculation (corresponding to deleting each of the 12 30° azimuth bins from the full data set). In general, the distributions are peaked near the value for the full data set with a relatively small range of values. The positive values of μ_rθ^(2,0) (most trials) and the negative values of μ_rφ^(2,0) (11 of 12 trials) both correspond to a source volume that dips to the NE. The characteristic volume has a dip of about 7° to the NE, while the moment tensor (19°) and the fault plane determined from the local network (17°) have slightly steeper values. Thus, even for this shallow thrust event, it is possible to recover the approximate dip direction from the second moments.

1995 Gulf of Aqaba earthquake

The 1995 M_w7.2 Gulf of Aqaba earthquake ruptured the southern extension of the Dead Sea Transform (Fig. 12a). First orbit Rayleigh waves at low frequencies (≤10 mHz), measured using synthetics computed for the Harvard CMT, showed no systematic signal in either δτ_p or δτ_q. However, higher frequency measurements of teleseismic P waves did display systematically low amplitudes at stations to the south of the event (Fig. 13). We inverted measurements of δτ_p and δτ_q for R1 at 5 and 10 mHz and P waves at 25, 35 and 45 mHz for the first and second moments and the moment tensor elements. The characteristic dimension in the radial direction was constrained to be less than the Harvard CMT centroid depth of 18 km. The results of our inversion (Tables 1, 2 and 3) correspond to a characteristic duration, τ_c, of 12±1 s, a centroid velocity, v₀, of 4.3±0.3 km s⁻¹ in the direction N05°E and a characteristic rupture length, L_c, of 53±2 km in the direction N06°E. The centroid location (Fig. 12) is shifted slightly to the west of the gulf of Aqaba, presumably due to incorrectly modelled 3 D elastic structure.

The characteristic rupture volume corresponding to our values of μ̂^(2,0) is shown in map view in Fig. 12(b). One fault plane of the moment tensor strikes NNE–SSW while the other strikes WNW–ESE. The NNE striking plane has a value of x_c(n̂)=12±1 km while the ESE striking plane has a value of x_c(n̂)=50±2 km. Thus, we prefer the NNE striking plane as the true rupture plane. This fault plane is consistent with the orientation of the Dead Sea Transform in this region. The small difference in the azimuth of the characteristic rupture length, 006°, from the strike of the fault plane of the moment tensor, 018°, may be related to this event rupturing more than one fault. Pinar (1997) and Klinger et al. (1999) showed that this event's rupture jumped north across a series of faults in the Dead Sea Transform system. This jump would cause a difference between the azimuth of the fault plane in the focal mechanism and the azimuth of the characteristic rupture volume of the sense we determine, but the difference may also be a result of the limited station coverage in southern azimuths.

Discussion

We have developed an inversion procedure that can estimate the first and second degree moments of earthquake source distributions without any a priori constraints on the orientation or extent of the source. The second moments are a desirable quantity to determine because they are independent of the types of parametrizations and smoothing operations typically employed in teleseismic finite fault inversions, and they can be used as integral constraints in determining the space–time history of rupture. The technique presented in this paper can be applied to any seismic phase for which an isolation filter can be constructed. Future studies will take advantage of the different samplings of the source provided by using measurements from a variety of phases/components. Our estimates of the first and second moments agree well with values derived from slip inversions using local data for both the Kobe and Jalisco events. Additionally, the Jalisco and Aqaba events demonstrate the ability of the technique to resolve the fault plane ambiguity.

One potential source of error in our estimates comes from the affects of unmodelled lateral heterogeneity. We have limited our measurements to combinations of frequency bands and phases that minimize these affects. Also, our data selection process removes any waveforms with large, unmodelled anomalies. The clearest way to evaluate the amplitude of unmodelled heterogeneity is to make GSDF measurements for relatively small earthquakes that should show little to no directivity affects in the frequency bands we work in. Fig. 6 is an example of this using the Kobe earthquake. There are no systematic azimuthal anomalies for this event's Rayleigh waves at frequencies below 15 mHz, and the amplitude of the scatter, which presumably results from unmodelled 3 D heterogeneity, is much smaller than the directivity signal observed for the Jalisco event. We have made measurements for other relatively small earthquakes and find that the δτ_q measurements for P waves in the band between 20 and 50 mHz typically have rms scatter of 1–2 s, much smaller than the directivity signal observed for the Aqaba event. Future improvements in the accuracy of our synthetic seismograms, such as the inclusion of 3 D attenuation models, will expand the range of frequencies in which we can confidently measure the waveform anomalies resulting from the finite spatial extent of seismic sources.

While regions such as Mexico and Japan have extensive local networks that allow detailed study of their earthquakes, there are many interesting events that do not have sufficient local data sets. For instance, the 1998 M_w 8.2 Antarctic earthquake occurred over 1000 km from the nearest seismometer. We applied the technique described in this paper to measurements of the fundamental mode Rayleigh wave for this event to determine its first and second moments (McGuire et al. 2000). In addition to being able to resolve the fault plane ambiguity for this event, the estimates of the centroid and rupture length, l_c, were used by Toda & Stein (2000) to calculate the change in Coulomb failure stress at the locations of this event's aftershocks. In the future, we hope to utilize second moments to help constrain the details of slip distributions from teleseismic data. The moments can be used as integral constraints to help regularize teleseismic finite fault inversions without requiring arbitrary choices of damping parameters. We expect future studies of second moments using the technique developed in this paper to be useful in a wide range of earthquake source, seismotectonic and plate motion studies.

2010 Appendix A:Travelling wave representation of normal-mode summation synthetic seismograms

Here we present an outline of the derivation of eq. (14), the travelling wave representation of our synthetic seismograms. The expression for a synthetic seismogram due to a point source computed by the summation of normal modes can be written as

where n is the overtone number, l is the angular degree, P_lm is the associated Legendre polynomial and _nI_l^R and _nE_l^S are differential operators related to the excitation of the kth mode [k=(n, l)] at the receiver and by the source, respectively (Zhao & Jordan 1998, Appendix A). Using the Poisson sum formula (Dahlen & Tromp 1998, Chapter 11), the sum over angular degree l can be transformed into an integral over continuous wavenumber λ:

Substituting in the asymptotic expression for P(−cos Δ),

and performing the contour integral in the upper half plane gives

where λ_n is the wavenumber at which ω_n(λ)=ω, v_n(λ_n) is the group velocity at λ_n, and the index j relates to the orbit number. The corrections made to our synthetics for lateral heterogeneity lead to a slight shift in the eigenfrequencies of the modes such that λ_n→λ_n+δλ_n.

For first orbit waves, j=0 and only the first exponential term (corresponding to waves travelling in the direction from the source to the receiver along the minor arc) is retained, leading to eq. (14):

where the expressions for A_n and B_n are

with the source scalars

with

and the receiver scalars

where v=(v_r, v_θ, v_φ) is the unit directional vector of the instrument and a dot over an eigenfunction (U, V, W) represents a radial derivative. The source and receiver locations are r_S=(r_S, θ_R, φ_R) and r_R=(r_R, θ_R, φ_R). Results for other odd orbit waves, such as R3, are obtained in the same way with appropriate values of j. The expressions for even orbit waves, such as R2, are found from eq. (A4) by taking the second exponential term corresponding to waves travelling along the major arc with the appropriate values of j.

2010 Appendix B:Non-linear constraints on the moments

Here we show that the matrix inequality constraint we incorporate in our inversion (22) encompasses the constraint derived from Bessel inequalities by Bukchin (1995). For the planar fault case treated by Bukchin, eq. (22) is equivalent to

where x and y are the Cartesian spatial coordinates and t is the time coordinate. One property of a positive semi definite matrix is that its upper left determinates are non negative, i.e.

and

The determinant in eq. (B2) is

Solving for the second temporal moment assuming eq. (B3) holds gives

which can be written in terms of the second moments as

where Adj denotes the adjoint of a matrix. Eq. (B6) is equivalent to the bessel inequality, eq. (24), of (Bukchin 1995).

Acknowledgments

We thank Satoshi Ide for providing us with the centroid of his slip distribution for the Kobe event, and Tim Melbourne for providing his slip distribution for the Jalisco event. Figures in this paper were generated using the GMT software freely distributed by Wessel & Smith (1991). This research was sponsored by the National Science Foundation under grant EAR 9805202.

References