Earthquake magnitude estimation from early radiated energy

Click Here GEOPHYSICAL RESEARCH LETTERS, VOL. 35, L22307, doi:10.1029/2008GL035576, 2008 for Full Article Earthquake magnitude estimation from early radiated energy Gaetano Festa,1 Aldo Zollo,1 and Maria Lancieri2 Received 4 August 2008; revised 14 October 2008; accepted 16 October 2007; published 27 November 2008. [1] From inspection of a large set of Japanese events, we investigate the scaling of the early radiated energy, inferred from the squared velocity integral (IV2) with the final magnitude of the event. We found that the energy can only discriminate whether the event has a magnitude larger or smaller than 5.8, and in the latter case it can allow for realtime magnitude estimation. However, by normalizing IV2 for the rupture area, the initial slip scales with the magnitude between 4 < M < 7 following the expected scaling laws. We show that the ratio between the squared peak displacement and IV2 is a proxy for the slip following the same scaling but it can be directly derived from the data, without any assumption on the rupture area. The scaling relationship between initial slip and magnitude can be used for early warning applications, when integrated in a probabilistic, evolutionary approach. Citation: Festa, G., A. Zollo, and M. Lancieri (2008), Earthquake magnitude estimation from early radiated energy, Geophys. Res. Lett., 35, L22307, doi:10.1029/ 2008GL035576. 1. Introduction [2] Earthquake early warning systems are real-time monitoring infrastructures designed to provide a rapid notification of the potential effects of an impending earthquake, through the fast telemetry and the processing of data from dense instrument arrays deployed in the source region or surrounding the target site. [3] As a first order approximation, the amplitude and the characteristic frequency of the seismic records depend on the event location and magnitude and on the attenuation mechanisms that the waves undergo during the propagation from the source to the site of interest. The latter effect is almost independent of the event size and it can be modeled with sufficiently high accuracy. On the other hand, the earthquake location can be very rapidly determined from early signals recorded at a few stations close to the hypocenter [Horiuchi et al., 2005]. Therefore the lead time of an early warning system critically depends on the system capability to predict the final earthquake magnitude from measurements on early recorded signals. [4] Correlations between the final magnitude of an earthquake and several parameters measured on the early portions of P-waves and S-waves have been widely investigated for seismic early warning applications. The analyses of the seismic databases from southern California [Allen and Kanamori, 2003; Olson and Allen, 2005; Wu and Zhao, 2006], Taiwan [Wu and Kanamori, 2005], the EuroMediterranean region [Zollo et al., 2006] and Japan [Zollo et al., 2007; Lockman and Allen, 2005] show that the predominant period and the peak of ground displacement scale with the final magnitude over a wide range (4 < M < 8), arguing that the energy available to break new asperities may already differ at the initial stage of the rupture [Olson and Allen, 2005; Zollo et al., 2006]. [5] Although the observations indicate that the early peak and the dominant frequency of seismic signals do increase with the magnitude, the limit at which the prediction is reliable remains debatable. Rydelek and Horiuchi [2006] and Rydelek et al. [2007] have claimed that neither the initial predominant period nor the peak ground displacement significantly increase with a magnitude beyond M = 6. Moreover geological observations support the hypothesis that the arrest mechanisms of large earthquakes are mainly associated with structural features in the fault geometry, irrespective of the event size [Wesnousky, 2006]. [6] For early warning applications, the ground motion estimations associated with the more destructive S wave rely on the information carried out by the first seconds of the P wave at the same recording site. When the target site is far from the hypocenter (distances larger than 80 km; regional early warning), the S wave recorded close to the fault can help in constraining the ground motion estimates without significantly increasing the system lead-time. Nearby the epicenter, S data can even provide the only estimation of the magnitude, when the S minus P time is smaller than the P window (few seconds) required for the estimation of the magnitude. [7] In this work we investigate the scaling of the radiated energy, inferred from squared velocity integral measured on early P- and S-waves signals, with the final size of the event. From the analysis of moderate to large earthquakes, recorded by the Japanese strong motion networks Kik-net and K-net, we discuss the behavior of the ‘‘macroscopic’’ slip at the beginning of the rupture in the magnitude range 4 < M < 7. 2. Data Analysis [8] Let us consider the integral of the squared velocity (IV2) measured in the early portions of P-waves and S-waves IV 2c ¼ 1 Dipartimento di Scienze Fisiche, Università di Napoli ‘‘Federico II’’, Naples, Italy. 2 Osservatorio Vesuviano, Istituto Nazionale di Geofisica e Vulcanologia, Naples, Italy. Copyright 2008 by the American Geophysical Union. 0094-8276/08/2008GL035576$05.00 tcZ þDtc v2c ðt Þdt tc where the subscript c refers to the P or S phase, tc is the corresponding first arrival, and vc is the particle velocity measured on the seismograms. Finally Dtc is the length of the signal along which the analysis is performed. IV2 has L22307 1 of 6 L22307 FESTA ET AL.: EARTHQUAKE MAGNITUDE ESTIMATION L22307 Figure 1. (top) One horizontal component (EW) of the velocity recorded at about 45 km from the hypocenter, for two events of magnitude M = 5 and M = 7.1 (the key gives the recording station references). (bottom) The squared velocity modulus as a function of time for the same events as plotted in Figure 1 (top). Between the two events the level of the energy drops by more than three orders of magnitude, with a steep change of the trend in correspondence of the P and S wave arrivals. the advantage that it includes information about the energy radiated by the advancing rupture on the fault plane, and therefore it provides direct, although partial, insights into the physics of the fracture [Boatwright and Fletcher, 1984]. [9] We have processed about 2500 high-quality strongmotion records from the K-Net (http://www.k-net.bosai. go.jp/k-net/index_en.shtml) and Kik-net (http://www.kik. bosai.go.jp/kik/index_en.shtml) Japanese databases. The records correspond to events that occurred from 1996 to 2005 and with magnitudes ranging from 4.0 to 7.0. We limited the analysis to stations close to the fault, for which the hypocentral distance was less than 60 km. The P and S phases were manually picked from the records, and then the velocities were obtained by integration and band-pass filtering in the frequency band of 0.05 – 10 Hz. The IV2 values were grouped into magnitude bins of 0.3, with the bin width related to the errors in the magnitude estimates. For the analysis, we selected Dtp = 4s of signal after the first P arrival and Dts = 2s beyond the S arrival. [10] We stress here that 2s and 4s on the records do not correspond to 2s and 4s of rupture when the observer is at near source distances. Since the velocity at which a dynamic rupture advances on a fault cr is comparable with the velocities at which the information travels in an elastic medium (cp or cs for P and S waves respectively), the image that an observer has of an advancing rupture is deformed and stretched in the direction of the observer. Moreover, different observers look at the fault from different views, giving an overall image that relates to the envelope of the scanned areas [Festa and Zollo, 2006]. Since cs < cp, the fault region spanned by the P waves is smaller than the one imaged by the S waves, for the same time window on seismograms after the initial P- and S-arrivals. Since cs is close to cr, the S wave explores shallow regions in the observer direction, while the P wave images are almost isotropic around the hypocenter. These are the reasons for which we need a larger P time window to capture the same information as is carried by the S wave. [11] Figure 1 shows one component of the velocity, integrated from the accelerometric records for two events of magnitude M = 7.1 and M = 5.0 respectively, which were recorded at about 45 km from the hypocenter. It also shows a logarithmic plot of the squared velocity modulus as a function of time. The energy increases by more than three orders of magnitude between the two events, with a steep 2 of 6 L22307 FESTA ET AL.: EARTHQUAKE MAGNITUDE ESTIMATION L22307 Figure 2. Scaling of the 4s velocity integral for the P waves (squares) and 2s velocity integral for the S waves (circles) as a function of the magnitude. Superimposed on the picture we also add the velocity integral over the total S phase for large earthquakes (diamonds). The error-bars represent the standard deviation error around the mean value. The total velocity integral points are aligned along a straight path with a least-squares slope compatible with standard scaling laws. Straight lines with this slope fit the first P-wave and S-wave early energy points up to M = 5.8. variation in correspondence of the P and S arrivals and an almost constant level for few seconds beyond them. [12] To compare records from stations located at different distances from the hypocenter, we normalized the measurements to the reference distance R0 = 10 km, by analytically removing the geometrical spreading term log(R2/R20). The final value of the integral of squared velocity is therefore . referred to as IV210km c 3. Radiated Energy and Average Slip (squares) and [13] Figure 2 shows the plot for the IV210km P IV210km (circles) as functions of the final magnitude of the S event. To discuss the behavior of these quantities, we also add the velocity integral evaluated for the whole signal S in the case of large events (diamonds). When accounting for the whole duration, the points are aligned along a straight line, with a slope (a = 1.41 ± 0.04) that is compatible with the expected scaling factor of 1.5 [e.g., Scholz, 2002]. Straight lines with this slope fit both the P and S data up to a magnitude M = 5.8. Beyond this, the early energy increases less, or does not increase at all, with respect to the final magnitude. A rupture size having a magnitude M = 5.8 is comparable with the area imaged by the back-propagation of the selected P- and S-windows. [14] By interpreting the velocity integral representations in the light of the scaling laws, for a window time length of 4s for P-waves and 2s for S-waves, we can conclude that below M = 5.8 the apparent duration is smaller than the investigation time window and the increase in the emitted radiation is associated with both the increasing fracture area and the increasing average slip. Beyond a magnitude of 5.8, the data provide a partial image of the advancing rupture, coming from a fault portion which has almost the same area, despite the magnitude. Any increase in the velocity integral has to be ascribed to the slip. [15] Hence, the velocity integral can be used only to discriminate whether an event has a magnitude larger or smaller than M = 5.8 and in the latter case to evaluate the magnitude of the event. We compute regressions laws through the first P and S points up to magnitude M = 5.8. ) = 7.7(±0.3) + The resulting curves are log(IV210km p ) = 6.3(±0.4) + 1.4(±0.1)M for P waves and log(IV210km s 1.4(±0.1)M for S waves, where the velocity is measured in cm/s. We remark that below M = 5.8 any deterministic method (such as local Ml or moment Mw magnitude) can be used to estimate the magnitude since the selected time windows entirely contain the direct waves emitted by the source. Beyond M = 5.8 the prediction becomes ill-posed and if a scaling with the magnitude exists, it has to concern the slip. [16] For this, we transform the velocity integral into the radiated energy E [Kanamori et al., 1993]: Ec ¼ 4p R2 rcIV 2c F 2 <2c where we use the average values r = 2.7 g/cm3 for the density, F = p 2 ffiffifor ffi the free surface coefficient, cs = 3.3 km/s and cp/cs = 3. <2c is the ratio between the actual and the average squared radiation pattern that we fixed to 1 [Kanamori et al., 1993]. We also assume the expected value for the ratio Es/Ep = 16.7 [Boatwright and Fletcher, 3 of 6 L22307 FESTA ET AL.: EARTHQUAKE MAGNITUDE ESTIMATION L22307 Figure 3. The slip as a function of the magnitude from the P (squares) and S (circles) early records. The dashed and solid lines represent the best-fit curves through the points P and S respectively. The dotted lines represent the prediction bounds for a new observation with a 95% confidence level in correspondence of the S best-fit curve. 1984] to obtain the total radiated energy E from the P and S energy estimates. The radiated energy is related to the slip hDi as: estimates from the P and S data very consistent each with the other. Superimposed on these data, there is the plot of the best-fit lines going through the data points (solid line for S waves and dashed line for P waves). The two curves are: 1 E ¼ Dsh DiA 2 where Ds is the stress drop. Since the velocity integral computed along the complete signal scales with the magnitude following the theoretical scaling law, we see that the stress drop must be nearly constant in the investigated magnitude range. Therefore, we can estimate the average stress drop from the level of the velocity integral, obtaining Ds = 6.7 MPa. We note here that under the hypothesis of constant stress drop, the latter only influences the absolute value of the slip, and does not affect the scaling of the slip with the magnitude. [17] Finally, to obtain the slip, we have to normalize the radiated energy for the rupture area as seen from the early seconds on the records. For small earthquakes (M < 5.8), we use the moment versus area scaling relationship for a circular crack [Madariaga, 1976]: logðh DiÞP ¼ 0:30ð 0:05ÞM  2:6ð 0:3Þ logðh DiÞS ¼ 0:40ð 0:05ÞM  3:2ð 0:2Þ with the slopes comparable with a = 0.5 of the theoretically expected scaling laws [Scholz, 2002]. The high correlation of the mean values can also be quantified through the values of the correlation coefficients which are r2p = 0.86 and r2s = 0.92, for P and S wave respectively. [19] We note that, assuming a constant stress drop, the radiated energy scales as the product hDiA, while the peak displacement PD scales as hDi1/2. Hence, if the scaling law holds also in the early portion of the signal, the ratio PD2/ IV2 should scale as the slip. The strength of this parameter is that it represents a proxy for the slip directly measurable on seismic signals. Figure 4 shows the ratio for 4s of P-wave and 2s of S-wave. The two best-fit curves are: 16 M0 ¼ 3=2 DsA3=2 7p   log PD2 =IV 2 P ¼ 0:38ð 0:02ÞM  3:28ð 0:13Þ   log PD2 =IV 2 S ¼ 0:48ð 0:03ÞM  3:82ð 0:17Þ where the moment M0 is known from the velocity integral. Then, beyond magnitude M = 5.8, a constant value for the rupture area is used, corresponding to the one estimated for an earthquake of M = 5.8. [18] Figure 3 shows the log of the slip, as imaged from the first few seconds of P and S wave records on seismograms, as a function of the magnitude of the event. We see here that the slip from both the P-waves and S-waves scales with the magnitude from M = 4 to M = 7, with the slip Again, the two slopes are comparable with those of Figure 3, as well as with the expected theoretical scaling. The values of the correlation coefficients r2p = r2s = 0.97 are representative of the high correlation of the mean values. [20] Although the variation of the slip with magnitude is statistically significant, the data dispersion in Figures 3 and 4 indicates that the inferred scaling relationships obtained cannot be considered as deterministic laws, since the errors 4 of 6 L22307 FESTA ET AL.: EARTHQUAKE MAGNITUDE ESTIMATION L22307 Figure 4. Scaling of the ratio PD2/IV2 as a function of the magnitude in the early portion of the P-waves and S-waves. The scaling is comparable to that obtained in Figure 3. The dotted lines represent the prediction bounds for a new observation with a 95% confidence level in correspondence of the S best-fit curve. on the slip are large. Hence, a single measure of slip cannot significantly discriminate an event of magnitude M = 6 from an event of magnitude M 7. When drawing on Figure 4 a horizontal line in correspondence of the value of slip expected for a magnitude M = 6, the probability that a value smaller than this will belong to the distribution of the slip data for M = 7 is 10% and 19% from S and P estimates respectively. [21] However, since the data are normally distributed, the error on the mean value decreases when extracting information from several receivers which are representative of the slip distribution corresponding to a given value of magnitude. In such a case, the uncertainty associated with the mean becomes comparable with the error from the bestfit curve, which is about 0.5 in magnitude. 4. Conclusions [22] The squared velocity integral measured on the first few seconds of P and S waves as a function of the final magnitude of the event shows a two-slope curve with separation limit at about M = 5.8. For M < 5.8, both 4s of P-wave and 2s of S-wave image the whole rupture process, and in this range IV2 can be used as an estimator of the magnitude of the event. Beyond M = 5.8 the area imaged by the first few P and S records remains almost the same despite the magnitude and any increase in IV2 vs M plots has to be ascribed to the variation of the slip. Under the hypothesis of constant stress drop, the initial slip can be achieved by normalizing the radiated energy by the rupture area as seen by the first P- and S-wave segments. The rupture size up to magnitude M = 5.8 was estimated using the moment versus area relationship for circular crack while a constant value, corresponding to M = 5.8, was used for large earthquakes. After normalizing by the rupture area, a linear relationship between the slip and the magnitude yields in the whole magnitude range, whose slope is comparable with that from standard scaling laws. This result implies that, for 5.8 < M < 7, the average slip measured in the early stage of the rupture looks similar to the final distribution of the slip on the fault plane, according to the hypothesis that large slip zones are likely to be located within or close to the hypocenter [Mai et al., 2005]. Additionally, the slip in the initial part of the rupture has reached or is closed to its final value, indicating that the rise time should be associated with a time scale smaller than 4s of the P-wave and 2s of the S-wave. [23] Under the same assumptions, the trend of the slip with magnitude can be directly inferred by measurements on the data, through the ratio PD2/IV2. We prove that this quantity follows the same scaling as the theoretical model, with smaller errors in the estimations of the best-fit coefficients. Although the errors on the single measurement of PD2/IV2 are too large for a deterministic prediction of the magnitude, several measurements can be combined in a probabilistic evolutionary framework for early warning applications. If PD2/IV2 values are representative of the probability function at the corresponding magnitude, the error in the final estimation can be pushed down to about 0.5 units. [24] Acknowledgments. Authors acknowledge AMRA Scarl, ‘‘Analisi e Monitoraggio del Rischio Ambientale’’ for financial support and two anonymous reviewers for their valuable contribution in improving the manuscript content. References Allen, R. M., and H. Kanamori (2003), The potential for earthquake early warning in southern California, Science, 300, 5620, 786 – 789, doi:10.1126/science.1080912. 5 of 6 L22307 FESTA ET AL.: EARTHQUAKE MAGNITUDE ESTIMATION Boatwright, J., and J. B. Fletcher (1984), The partition of radiated energy between P and S waves, Bull. Seismil. Soc. Am., 74, 361 – 376. Festa, G., and A. Zollo (2006), Fault slip and rupture velocity inversion by isochrone backprojection, Geophys. J. Int., 166, 745 – 756. Horiuchi, S., H. Negishi, K. Abe, A. Kamimura, and Y. Fujinawa (2005), An automatic processing system for broadcasting earthquake alarms, Bull. Seismol. Soc. Am., 95, 708 – 718. Kanamori, H., E. Hauksson, L. K. Hutton, and L. M. Jones (1993), Determination of earthquake energy release and ML using TERRAscope, Bull. Seismol. Soc. Am., 83, 330 – 346. Lockman, A. B., and R. M. Allen (2005), Single-station earthquake characterization for early warning, Bull. Seismol. Soc. Am., 95, 2029 – 2039. Madariaga, R. (1976), Dynamics of an expanding circular fault, Bull. Seismol. Soc. Am., 66, 639 – 666. Mai, P. M., P. Spudich, and J. Boatwirght (2005), Hypocenter locations in finite-source rupture models, Bull. Seismol. Soc. Am., 95, 965 – 980. Olson, D., and R. M. Allen (2005), The deterministic nature of earthquake rupture, Nature, 438, 212 – 215, doi:10.1038/nature04214. Rydelek, P., and S. Horiuchi (2006), Is earthquake rupture deterministic?, Nature, 442, E5 – E6, doi:10.1038/nature04963. Rydelek, P., C. Wu, and S. Horiuchi (2007), Comment on ‘‘Earthquake magnitude estimation from peak amplitudes of very early seismic signals on strong motion records’’ by Aldo Zollo, Maria Lancieri, and Stefan Nielsen, Geophys. Res. Lett., 34, L20302, doi:10.1029/2007GL029387. L22307 Scholz, C. (2002), The Mechanics of Earthquakes and Faulting, 2nd ed., Cambridge Univ. Press, Cambridge, U. K. Wesnousky, S. (2006), Predicting the endpoints of earthquake ruptures, Nature, 444, 358 – 360, doi:10.1038/nature05275. Wu, Y. M., and H. Kanamori (2005), Rapid assessment of damage potential of earthquakes in Taiwan from the beginning of P waves, Bull. Seismol. Soc. Am., 95, 1181 – 1185. Wu, Y. M., and L. Zhao (2006), Magnitude estimation using the first three seconds P-wave amplitude in earthquake early warning, Geophys. Res. Lett., 33, L16312, doi:10.1029/2006GL026871. Zollo, A., M. Lancieri, and S. Nielsen (2006), Predicting the earthquake magnitude from peak amplitudes of very early seismic signals, Geophys. Res. Lett., 33, L23312, doi:10.1029/2006GL027795. Zollo, A., M. Lancieri, and S. Nielsen (2007), Reply to comment by P. Rydelek et al. on ‘‘Earthquake magnitude estimation from peak amplitudes of very early seismic signals on strong motion records’’, Geophys. Res. Lett., 34, L20303, doi:10.1029/2007GL030560.  G. Festa and A. Zollo, Dipartimento di Scienze Fisiche, Università di Napoli ‘‘Federico II’’, I-80124 Napoli, Italy. (festa@na.infn.it) M. Lancieri, Osservatorio Vesuviano, Istituto Nazionale di Geofisica e Vulcanologia, I-80124 Napoli, Italy. 6 of 6