Effect of directivity on estimates of radiated seismic energy
Abstract
[1] To accurately estimate the seismic energy radiated in an earthquake, it is important to use appropriate corrections for path and source effects. Using slip models obtained by inversion of seismic data, we examine the effect of directivity on estimates of radiated energy and develop a method to correct for this source effect. From our calculations we suggest that the directivity correction for the earthquakes we studied is less than a factor of three; we also conclude that the directivity correction for dip-slip earthquakes with rupture propagating along strike is less than a factor of two at teleseismic distances. Thus teleseismic energy estimates of large subduction zone earthquakes will not be significantly affected by directivity.
1. Introduction
[2] To understand the physics of earthquakes, we have to relate seismologically observable parameters to the dynamics of faulting. One of the key seismological parameters that will help us achieve this objective is radiated energy. Radiated energy is defined as the wave energy that would be transmitted to infinity if an earthquake occurred in an infinite, lossless medium [Haskell, 1964]. When seismic waves that are generated by an earthquake travel through the earth, the earth structure acts like a filter and modifies these waves. To recover the radiated energy of the earthquake, the recorded seismic waves have to be corrected for propagation path effects such as attenuation, site effects and geometric spreading, and also for source effects such as radiation pattern and directivity. To achieve this, we require information of the earth structure over the entire frequency band, typically between 0.01 Hz to a few Hz for large earthquakes, since radiated energy is measured over a broad range of frequencies. Though significant progress has been made in the study of the structure of the earth, we still have insufficient information on the structural details that would be important at higher frequencies. This leads to inaccuracies in the corrections applied at higher frequencies which result in uncertainties in energy estimates.
[3] With the advent of broadband networks, estimates of radiated energy by direct integration of velocity records have improved considerably [Boatwright and Choy, 1986; Boatwright and Fletcher, 1984; Choy and Boatwright, 1995; Houston, 1990; Houston and Kanamori, 1990; Kanamori et al., 1993; Singh and Ordaz, 1994; Winslow and Ruff, 1999]. Despite the increased availability of broadband data, it has been observed that for the same earthquake the estimates of radiated energy from regional data differ from those obtained from teleseismic data by as much as a factor of 10 [Singh and Ordaz, 1994]. Venkataraman et al. [2002] used an empirical Green's function method to estimate radiated energy from regional data and showed that the energy estimates from regional data are in agreement with those obtained from teleseismic data for the 1999 Hector Mine earthquake. However, they did not correct their energy estimates for directivity.
[4] The primary focus of this paper is to examine the effects of directivity on energy estimates. First, we briefly discuss the method of estimating radiated energy from observed data followed by the theoretical formulae involved in calculating radiated energy from seismic waves. This will help us understand how directivity enters the formulation. Next, we use kinematic as well as earthquake slip models to illustrate the effect of directivity on estimates of radiated energy.
2. Estimating Radiated Energy From Observed Data
[5] Radiated energy can be estimated by integrating far field velocity records in a homogeneous medium. Seismic waves in the earth, however, are modified as they travel through the heterogeneous, attenuating structure of the earth. To recover the radiated energy of the earthquake, we have to correct the observed velocity records for these propagation path effects. If we have seismic stations at closely spaced points on the earth, we can determine the energy flux from these corrected velocity records, and then sum the flux over all these points to obtain the total radiated energy. Since we have limited station locations on the earth, the observed data cover only portions of the focal sphere. To overcome this limitation, we backproject the data to cover the entire focal sphere by removing the radiation pattern factor at each station and then using an average radiation pattern factor (details explained later). Thus we determine the total energy of the earthquake from the corrected data at an individual station: a method that we refer to as the single-station method.
[6] By backprojecting the corrected direct P and S wave data to cover the entire focal sphere, we account for energy leaving the source from all the rays taking off in all directions. P and S waves are sufficient to describe the earthquake source, as surface waves are only formed later as a result of the interaction between the P and S waves at the surface of the earth. We do not include the energy from surface waves since this energy is already measured when the energy is propagating as body waves. Thus by calculating the radiated energy from the direct P and S wave arrivals, we account for the total energy radiated by the earthquake.
(f)) from the displacement spectrum (
(f)) at a station:
(f) is the instrument response. We then use the moment rate spectrum, (f), at each station to compute radiated energy using the single-station method. We briefly discuss the corrections below, but more detailed discussions can be found in the work of Boatwright and Choy [1986], Newman and Okal [1998], and Venkataraman et al. [2002].
2.1. Attenuation Correction
[8] Anelastic processes in the earth attenuate seismic waves as they travel from the source to the receiver. When we determine radiated energy, we correct the observed data for this loss of energy due to attenuation. However, the attenuation structure of the earth, especially at higher frequencies, is not very well known. To compute energy from teleseismic data, we used the best available attenuation model and modified it to include the effects of lateral heterogeneities. The details of the approach are outlined by Venkataraman et al. [2002]. We also examined the effect of attenuation on energy estimates and concluded that uncertainties in attenuation corrections affect the energy estimates of large earthquakes (MW > 7.5) by less than a factor of 1.5.
2.2. Source Radiation Pattern Correction
[9] The radiation pattern factor determines the azimuthal variation in the amplitude of the seismic waves and is a function of the takeoff angle of the seismic ray and the azimuth of the station. To account for radiation in the entire focal sphere from data at limited seismic stations, we backproject the data from the seismic station to the source by dividing the observed data at a station by the radiation pattern factor and then multiplying it with the average radiation pattern (4/15 for the P wave and 2/5 for the S wave).
[10] Radiated energy from teleseismic data is computed using direct P waves. In shallow earthquakes recorded at teleseismic distances, the direct phases radiated by the earthquake cannot be separated from some of the reflected phases. Thus the teleseismic waveforms radiated by shallow events are usually modeled as a group of phases [Kanamori and Stewart, 1978; Boatwright and Choy, 1986]. For example, the P wave group comprises the direct P phase, and the depth phases, pP and sP. When the energy carried by these phases is small, the correction for the radiation pattern factor is large. This is especially the case for shallow strike-slip earthquakes, as has been pointed out by several investigators [e.g., Boore and Boatwright, 1984; Boatwright and Choy, 1986] and discussed in detail by Newman and Okal [1998]. In our study, we select a time window such that we minimize the inclusion of scattered arrivals in our energy estimates; also, we use a radiation pattern factor of 0.2 as a cutoff for the P wave group, so as to exclude nodal stations in our energy estimates.
2.3. Source Structure
[11] In computing radiated energy, we backproject the seismic waves to the source and determine the wave energy that radiates from the focal sphere surrounding the source. The underlying assumption in this formulation is the homogeneity of the medium surrounding the source; however, the material surrounding the earthquake source is far from homogenous. For example, in subduction zone earthquakes that occur at the plate interface, energy is radiated into the subducting oceanic plate as well as the overlying continental plate. For a reasonable contrast in density and velocity of the two layers, the energy estimates would change by a factor of two or less. Thus except in the case of sharp density and velocity contrasts, radiated energy estimates would not be significantly affected by heterogeneous source structure if the average values of density and velocity at the source region were used.
2.4. Scattering
[12] Scattering of seismic energy is caused when the seismic waves interact with small-scale heterogeneities. Owing to scattering, a part of the high-frequency energy arrives after the direct arrivals in waves that are called coda. Scattering can also decrease the amplitude of a seismic phase by shifting energy from the direct arrival into the coda [Newman and Okal, 1998]. Thus scattered energy can add or remove energy from the direct arrivals and affect the estimates of radiated energy. In most cases, by choosing a time window to include all the phases of interest (the P wave group for shallow earthquakes or the direct P for deep earthquakes), but excluding most of the scattered arrivals, we can limit the amount of scattered energy that is included in the energy estimates. However, at nodal stations, the signal-to-noise ratio is low and hence the scattered energy can be a significant part of the total energy; to avoid this problem, we exclude the nodal stations and only use stations that have radiation pattern factors 0.2 or larger in our energy estimates.
3. Theoretical Framework
is the density, α and β are the P and S wave velocities of the medium, 
α(t) and β(t) are the far-field P and S wave velocity records, respectively, and the integrations are over time and over a spherical surface surrounding the source.
(
, t) is the slip rate function on the fault plane, is the location of slip on the fault, r is the distance from to the station, the position of the observation point is given by θ and ϕ, and Rα(θ, ϕ) and Rβ(θ, ϕ) are P and S wave radiation pattern factors, respectively.
α(t, θ, ϕ) and β(t, θ, ϕ) are the moment rate functions for P and S waves observed in the direction given by (θ, ϕ). Thus equations (2a) and (2b) can be written as
3.1. Total Radiated Energy
e(θ, ϕ)dS.
3.2. Single-Station Method
[17] Ideally, we would like to determine the radiated energy flux (equations (7a) and (7b)) at closely spaced points on the earth, correct the flux for attenuation and geometric spreading, and then sum the flux over all these points to obtain the total radiated energy. However, since we have limited station locations on the earth, we use data at individual stations to determine the total energy of the earthquake: a method that we refer to as the single-station method.
3.2.1. Point Source
α(t, θ, ϕ) = (t) and β(t, θ, ϕ) = (t). Thus we can rewrite uα(t) and uβ(t) as
is the average radiation pattern coefficient for P waves and is given by
is the average radiation pattern coefficient for S waves and is given by
(t), is estimated from the displacement at a single station using equations (8a) and (8b).
(ω) is the moment rate spectrum and the hat is used to denote a quantity in the frequency domain (since ∣
(ω)∣ = ∣ω(ω)∣). We can also write the radiated energy as
3.2.2. Finite Source
(t, θ, ϕ) at a station from the observed displacement using equations (4a) and (4b); that is, the radiation pattern is removed, but the directivity effect is still included in (t, θ, ϕ). Then we substitute these estimates of (t, θ, ϕ) in equation (13) to estimate the total radiated energy from data at a single station and call it the “single-station estimate” of energy, 
R, given as
(t), is replaced by the P and S wave moment rate functions for a finite source, α(t, θ, ϕ) and β(t, θ, ϕ). Thus in equation (16), the dependence of the energy, R, on θ and ϕ is only due to directivity. Ignoring the contribution from P waves (which carry less than 5% of the radiated energy), R can be written as
R and the azimuthal variation of R is due to directivity. For a point source R should be exactly the same (= ER) at all stations. However, for a finite source, because of directivity, R varies with azimuth. An azimuthal average of R is not necessarily equal to the total radiated energy, ER. Thus to determine ER, we have to correct the average of R for directivity (by a method that will be described in a later section).
4. Observations of Directivity Effects
[21] To illustrate our method of evaluating the effect of directivity on energy estimates, we first use a theoretical Haskell source [Haskell, 1964]. We consider a longitudinal shear fault in which the rupture propagates along the slip direction, and the final displacement along the length of the fault is constant; we refer to this model as a unilateral strike-slip model (strike = 0°, dip = 90° and rake = 0°). For this unilateral strike-slip model we use the following fault parameters: fault length, L = 70 km, width W = 15 km, a ramp slip function with a rise time of 2.4 s, maximum slip = 3.2 m, rupture velocity = 0.9 times shear wave velocity (parameters similar to the 1992 Landers earthquake) with the rupture propagating along strike to the north (at an azimuth of 0°). For this model, the total radiated energy, ER, computed using Haskell [1964, equation (31)], is 3.8 × 1015 J.
R, using equation (16′) and the energy flux, eβ, using equation (7b), for different azimuths and takeoff angles. The estimates of energy flux are multiplied by 4πr2, i.e.,
R at a station to R at the same station. As mentioned earlier, R depends on both radiation pattern and directivity, whereas R depends only on directivity.
[23] To compute R and R as a function of angle (θ) and azimuth (ϕ), the focal sphere is divided into regions of equal area; we use about 50 values of θ where θ varies from 0° to 180°, and 50 values of ϕ at the equator, where ϕ varies from 0° to 360°. For the unilateral strike-slip model described above, Figure 1a shows R (dark lines) and 
R (gray lines) as a function of azimuth; for each azimuth we compute the energy estimates at a range of θ. Thus each line represents the variation in energy estimates with azimuth for a particular θ. All the computations are done using the analytical solution for radiated energy given by Haskell [1964]. Since the rupture is propagating to the north, the radiated energy is focused at azimuths close to 0°, and the radiated energy is minimal at stations away from the rupture direction (i.e., stations at an azimuth of 180°). For a point source, R is exactly the same at all stations; however, the finiteness of the source introduces strong directivity effects, and the resulting differences in R at different stations are quite pronounced (two orders of magnitude or more in this case).
R (equation (17)) calculated for different azimuths and takeoff angles; plotted in black are estimates of R obtained by using equation (16′) for the same azimuths and takeoff angles. The azimuthal distribution of R is due to the combined effect of the radiation pattern factor (the lobes) and the directivity factor (the curve with maximum in the direction of rupture propagation, i.e., maximum at 0°) (see equation (6b)), whereas the radiation pattern effects have been removed in R, and its azimuthal variation reflects the directivity effect only; the different lines in the plot show the variation in the energy estimates with azimuth for a particular takeoff angle. The takeoff angles for R are marked on the right. The total radiated energy ER is shown by the dark line. (b) Bilateral strike-slip model, with the rupture starting at the center and propagating to the north and south. (c) Transverse shear model (dip slip), with rupture propagating along strike (to the north).
[24] The azimuthal average of R is not exactly equal to the azimuthal average of R. Thus even if the station coverage is good, in the presence of strong directivity effects, we cannot recover the total radiated energy (ER) from R by simply taking the average. In the case considered in Figure 1a, the ratio (ER/average of R) = 1.53.
[25] To study the effects of fault length on directivity, we modified the length of the fault from L = 80 km to L = 110 km, keeping all other parameters the same. For the Haskell model, the total energy and the average energy from single-station estimates do not depend on fault length (as long as the smallest duration of the source time function is larger than the rise time, i.e., L(1/Vr − 1/β) > rise time). Thus directivity, a factor that is thought to be more important for longer faults, is independent of fault length in the Haskell model (within the limitations pointed out). This is because, in these models, energy is radiated only at the edges of the fault, i.e., at the beginning and end of rupture.
[26] We also developed an analytical solution for the energy radiated by a bilateral rupture model and used it to compute R at different stations. Shown in Figure 1b is an example for a bilateral strike-slip fault of total length L = 40 km (all other parameters are the same as those used for the unilateral strike-slip model); the rupture propagates from the center of the fault to the north (azimuth = 0°) and south (azimuth = 180°). From Figure 1b, we observe that the effect of directivity is less pronounced in the case of bilateral rupture. In the case considered in Figure 1b, the ratio: ER/average of R = 1.21, and the total radiated energy, ER, is ∼4.0 × 1015 J.
[27] The third case plotted in Figure 1c shows transverse shear (dip slip) faulting (strike = 0°, dip = 90°, and rake = 90°) on a fault of length = 70 km with rupture propagating along strike to the north and all other parameters the same as before. In this case, the average of R is larger than the total radiated energy, ER, and thus the ratio: ER/average of R = 0.51, and the total radiated energy, ER, computed by Haskell [1964, equation (66)] is 1.28 × 1015 J. Since R depends only on directivity and the directivity for both dip-slip and unilateral strike-slip cases is along the strike on a vertically dipping fault, both cases have the same R (black lines in Figures 1a and 1c). However, R depends on both radiation pattern and directivity and is thus different for the two cases (the gray curves in Figures 1a and 1c).
[28] For the three cases considered above, we compared estimates of the average R for regional (takeoff angles close to horizontal, i.e., ih ≈ 90°) and average R for teleseismic (takeoff angles close to vertical, ih varying between 18° and 25°) distribution of stations to the total radiated energy. Regional stations refer to stations that are within 500 km form the epicenter, whereas teleseismic stations are stations at distances between 30° and 90°. Directivity effects cause the average of the regional and teleseismic R to be different from the total radiated energy, ER. The problem is most severe for teleseismic estimates from unilateral vertical strike-slip faults with rupture propagating along strike, because directivity causes the energy to be focused in the direction of rupture propagation. Additionally, the S wave radiation pattern factor is also large in this direction, and the combined effect of these two factors results in strong focusing of radiated energy in the direction of rupture propagation; very little energy is received at stations away from the direction of rupture propagation. Thus R (shown in gray in Figure 1a) at stations in the direction of rupture propagation is very large. For the unilateral strike-slip model considered above, ER could be a factor of 10 larger than the average of teleseismic R. Hence directivity corrections for teleseismic estimates of energy from unilateral strike-slip faults are very important; moreover, it is essential to have a good azimuthal coverage of stations because if there are no stations at azimuths close to the direction of rupture, the energy estimates will be severely underestimated.
[29] On the other hand, in transverse shear faulting (vertical dip slip) with rupture propagating along strike, directivity focuses the energy in the rupture direction; however, the S wave radiation pattern factor at takeoff angles close to horizontal (i.e., regional distances) is small in this direction. Hence R (shown in gray in Figure 1c) at these stations is not as large as R at azimuths close to 0° (and 360°). For the case considered ER is about a factor of 5 smaller than the average of the regional R.
[30] Thus directivity corrections are very important for teleseismic single-station estimates for vertical strike-slip faults with the rupture propagating along strike, and for regional single-station estimates for vertical dip-slip faults with rupture propagating along strike.
[31] The models considered so far are simple kinematic models, and though they are useful in understanding the effects of directivity on radiated energy estimates, it is important to study the more complex models of real earthquakes. Hence we studied the effects of directivity on radiated energy estimates by using slip models that were determined from inversion of seismic data. We used the rupture model of Dreger [1994] and Wald and Heaton [1994] for the 1992 Landers earthquake, the rupture model of Ji et al. [2002] for the 1999 Hector Mine earthquake and the rupture model of Ji et al. [2003] for the Taiwan earthquake. For each of these earthquakes, we calculated the single-station estimates, R from the slip models. The computations were done using a simple algorithm in which we divided the focal sphere into surfaces of equal area and used equation (16′) to calculate R at a station for a given slip model.
[32] The Landers earthquake was a strike-slip earthquake; the rupture was primarily unidirectional and propagated along strike (azimuth ∼ 340°). We used the slip model of Dreger [1994] and Wald and Heaton [1994] and to calculate R over a range of azimuths and θ. For the Dreger [1994] model, shown in Figure 2a, the total radiated energy, ER, is computed to be 4.2 × 1015 J. From Figure 2a, we observe that most of the energy is focused along the fault strike, similar to the unilateral strike-slip model of Figure 1a, but the energy distribution is more complicated than the simple unidirectional Haskell model. The large R estimates at azimuths of about 80° and 180° are due to small radiation pattern factors at nodal stations. Later, when we compare the model energy estimates with the energy estimates from data, we do not use these nodal stations (stations at which the radiation pattern factor is smaller than 0.2).
R (equation (16′)). The total radiated energy ER is shown by the dark line. (b) Hector Mine earthquake: estimates of R from the slip model of Ji et al. [2002]. (c) Taiwan earthquake: estimates of R from the slip model of Ji et al. [2003].
[33] Shown in Figure 2b are the radiated energy estimates for the Hector Mine earthquake determined using the model of Ji et al. [2002]. The Hector Mine earthquake was essentially a bilateral rupture with a fault about 40 km long [Ji et al., 2002; Trieman et al., 2002]. However, the rupture broke three fault segments and has a complicated slip history [Ji et al., 2002] probably far more complex than the Landers earthquake. Thus as seen in Figure 2b, the variation of R at stations is quite complicated as compared to the bilateral Haskell model of Figure 1b. From the model, the total radiated energy, ER, for the Hector Mine earthquake, is computed to be 1.0 × 1015 J. A similar computation for the Chi-Chi, Taiwan, earthquake of 20 September 1999, using the model of Ji et al. [2003] is shown in Figure 2c. The Taiwan earthquake was a thrust earthquake with a complicated slip distribution. The earthquake ruptured from south to north along strike (∼20°) and also downdip. Using this model, the total radiated energy, ER, of the Chi-Chi, Taiwan, earthquake is computed to be 4.9 × 1015 J.
[34] Thus in the three earthquake models mentioned the slip distribution is more complicated than in the simple kinematic models considered earlier. Variations of fault strike, rake direction, amount of slip and direction of rupture as the rupture propagates along the fault cause additional complications that are captured in the earthquake models. Hence complex earthquake sources could result in complicated azimuthal variation of radiated energy. If slip models are available, we can correct for these effects using the procedure outlined in the next section.
5. Corrections for Directivity
[35] To determine a correction for directivity, we compute R from observed data and from slip models for the same set of stations. First, we compare these estimates for the Hector Mine and Taiwan earthquakes. Figures 3a and 3b show a comparison between the data and model estimates of R at several stations for the Hector Mine earthquake where we use the slip model of Ji et al. [2002] to determine the model estimates. From the figure, we observe that there are differences between the energy estimates predicted by the model and those determined from data. These differences arise due to the following reasons: the synthetics created by the model do not exactly predict the data at frequencies larger than 0.5 Hz; the model is nonunique and may not be completely representative of the actual slip distribution especially because the Green's functions at higher frequencies are not well determined; the estimation of energy from data may have inaccuracies. Though the two estimates of energy are significantly different at some stations, the average of R from the model and regional data are almost the same. Also, the average of R determined from the model and teleseismic data differ by less than a factor of two.
[36] Figure 3c shows a comparison between R determined from the model of Ji et al. [2003] and from teleseismic data for the Taiwan earthquake. In Figure 3c, we observe that R at teleseismic stations can vary by a factor of ten between stations. Hence it is important to have a good azimuthal distribution of stations, so that the average of the single-station estimates yields a realistic estimate of the energy of the earthquake. With the azimuthal coverage of stations used in the Taiwan earthquake, the average of the R determined from the slip model differs from that determined from teleseismic data by less than a factor of two.
[37] From the above discussion, it can be seen that the effect of directivity on radiated energy estimates depends on the slip model and station distribution. For a given distribution of stations, we can compute the average R from the slip model and the average R from data. We can also compute the total radiated energy (ER) of the earthquake for this slip model. The ratio of (ER)/(average R) determined from the slip model would be representative of the directivity effect of the earthquake for the given slip model and station distribution. Thus we could use this ratio to correct the average R determined from data for directivity. By using this method, the actual radiated energy would still be determined from the data and the correction is only a factor that is applied to this observed estimate; also, the actual details of the slip model will not significantly affect the estimate of total radiated energy.
R at each station from the observed records at the station and denote it as EiD. We use slip distribution models from literature and compute R from the model at the same stations and denote it as EiM. We also compute the total radiated energy, ER, for the model and denote it as EM. The correction factor is given by the ratio of the actual total energy for the model to the average of single-station energy estimates for the model, i.e.,
| Earthquake Identification | Energy, J | Directivity Correction | Source for Slip Model | |
|---|---|---|---|---|
| Data (R/T)a | Directivity Corrected | |||
| 920628, Landers | 8.6 × 1015 (T) | 2.6 × 1016 | 3.08 | Dreger [1994] |
| 2.9 × 1016 | 3.34 | Wald and Heaton [1994] | ||
| 930115, Kushiro-oki | 4.2 × 1016 (T) | 2.9 × 1016 | 0.68 | Takeo et al. [1993] |
| 941004, Shikotan | 1.5 × 1017 (T) | 1.4 × 1017 | 0.93 | Kikuchi and Kanamori [1995] |
| 940609, Bolivia | 1.3 × 1017 (T) | 1.3 × 1017 | 0.98 | Kikuchi and Kanamori [1994] |
| 950730, Chile | 2.6 × 1016 (T) | 2.1 × 1016 | 0.80 | Ruegg et al. [1996] |
| 990920, Taiwan | 8.8 × 1015 (T) | 6.6 × 1015 | 0.76 | Ji et al. [2003] |
| 991016, Hector | 2.0 × 1015 (T) | 1.0 × 1015 | 0.49 | Ji et al. [2002] |
| 3.0 × 1015 (R) | 0.9 × 1015 | 0.31 | Ji et al. [2002] | |
- a R, regional data used to compute energy; T, teleseismic data used to compute energy.
[39] Since the slip models used for the Landers earthquake, the Hector Mine earthquake and Taiwan earthquake were determined by inverting both regional and teleseismic data (and GPS data in some cases), we have better constraints on the slip distribution of these earthquakes. For the other large earthquakes (Shikotan, Kushiro-oki, and Bolivia), slip distributions are not as well constrained, so the directivity corrections estimated are less reliable. However, these earthquakes are dip-slip earthquakes with rupture propagating along strike and the directivity corrections for these earthquakes at teleseismic distances are less than a factor of 2; moreover, as mentioned earlier, the actual details of the slip model will not significantly affect the estimate of radiated energy. Thus we can be confident that radiated energy estimates of large subduction zone earthquakes in which the rupture propagates along strike are not significantly affected by directivity.
6. Ratio of S to P Wave Energy
(t) is the same for both the S and P wave, from equation (13) we can write
[41] To determine the appropriate value of the ratio of S to P wave energy to use in our calculations, where the single-station estimates are used to determine the average, we define q1 = (average of Rβ)/(average of Rα), where Rα and Rβ are the P and S wave energy estimates as given by equation (16). Figure 4 shows a plot of the ratio q1, for a range of θ (where Rα and Rβ are averaged over azimuth at each takeoff angle). At teleseismic takeoff angles, i.e., ih between 18° and 25°, the ratio, q1 ≈ 25–30 for the three Haskell models studied earlier (unilateral strike-slip, bilateral strike-slip, and dip-slip models). Since the rupture propagates along strike in both the unilateral strike-slip and the dip-slip cases, and since the radiation pattern factor is averaged in the single-station method (equation (16)), we observe that q1 is the same for both the unilateral strike-slip and dip-slip cases.
Rβ)/(average of Rα), where Rα and Rβ are the P and S wave energy estimates as given by equation (16), is plotted for a range of takeoff angles (where Rα and Rβ are averaged over azimuth at each takeoff angle). At teleseismic takeoff angles, i.e., ih between 18° and 25°, the ratio q1 ≈ 25–30 for the three Haskell models studied earlier (unilateral strike-slip, bilateral strike-slip, and dip-slip models). The values obtained for the unilateral strike-slip and dip-slip cases are the same because R depends only on directivity.
[42] For the Landers earthquake, using the slip model of Dreger [1994], we obtained q1 ≈ 21; for the Hector Mine earthquake, we used the slip model of Ji et al. [2002] to obtain q1 ≈ 25 and for the Taiwan earthquake we used the model of Ji et al. [2003], and calculated the value of q1 ≈ 28. From the above calculations, we observe that the value of q1 for teleseismic data is close to the value obtained for a point source, and hence in this study we use the point source value of q1 ≈ 23.2. When a larger number of reliable slip models are available, the approach outlined above could be used systematically to understand the variation of q1 as a function of directivity.
7. Discussions and Conclusions
is the uncertainty of the energy estimate given by a factor, and 
s are the uncertainty factors of each of the corrections. As mentioned earlier, for most earthquakes, our estimates of the uncertainty in attenuation correction, radiation pattern factor plus directivity and ratio of S to P energy are 1.5 each, and the uncertainty due to source structure is a factor of 2. Then,
= 2.7. Thus we can state that the uncertainty in the resultant energy estimate is a factor of 2.7. Admittedly, this is a subjective statement for an average event. For some earthquakes (e.g., the events for which the source structure and the detailed rupture pattern are known), the uncertainty is smaller, and for earthquakes for which the source structure is laterally heterogeneous (e.g., subduction zone boundary), the uncertainty is larger.
[44] We use simple kinematic models as well as slip models obtained by inversion of seismic data to examine the effects of directivity on estimates of radiated energy and observe that directivity could cause a difference of more than three orders of magnitude in energy estimates between different stations. We observe that source complexities could result in complicated azimuthal variation of radiated energy. For example, the Hector Mine earthquake was mostly bilateral, but the complex rupture on pattern causes a complex azimuthal distribution of energy and could potentially bias the energy estimates. Thus to obtain accurate estimates of radiated energy, it is essential to have a good azimuthal coverage of stations and also apply a correction for directivity. Using the slip models available, we develop a correction for directivity. For a particular slip model of an earthquake, we can calculate the total radiated energy and the average single-station energy estimates for a particular station distribution; the ratio of these two estimates would be representative of the directivity effect of the earthquake for the given slip model and station distribution. We can use this ratio to correct the average of the single-station energy estimates determined from data for directivity and rupture complexity. By using this method, the radiated energy would still be determined from the data and the correction is only a factor that is applied to this observed estimate; hence, the actual details of the slip model will not significantly affect the estimate of radiated energy. Using the above method and slip models from literature, we calculated the directivity corrections for a few large earthquakes. From our calculations, we suggest that the directivity correction for the earthquakes we studied is less than a factor of three; we also conclude that the directivity correction for dip-slip earthquakes with rupture propagating along strike is less than a factor of two at teleseismic distances. Thus teleseismic energy estimates of large subduction zone earthquakes will not be significantly affected by directivity.
Acknowledgments
[45] We wish to thank Luis Rivera for helpful discussions on directivity. We thank Jim Mori and Ralph Archuleta for comments and suggestions that considerably improved our manuscript. This research was partially supported by NSF Cooperative Agreements EAR-9909371 and EAR-0125182 and also USGS-HQGR0035. This is contribution 8947 of the Caltech Division of Geological and Planetary Sciences.
