Large tsunamis have been responsible for great loss of life over time, including recent catastrophic events such as the 2004 Sumatra-Andaman earthquake (e.g., Shearer & Bürgmann, 2010) and the 2011 Tohoku earthquake (e.g., Lay, 2018), along with numerous other events (e.g., Lorito et al., 2016). Mitigation of loss of life requires rapid tsunami warning and evacuation of both local and remote regions based on quick quantification of the earthquake process and resulting tsunami excitation. Good progress has been made in developing seismological approaches (e.g., Duputel et al., 2011;Guilhem et al., 2013;Hirshorn et al., 2013;Kanamori & Rivera, 2008) and geodetic or seismogeodetic approaches (e.g., Blewitt et al., 2006;Colombelli et al., 2013;Hoechner et al., 2013;Melgar & Bock, 2013;Melgar et al., 2013Melgar et al., , 2015Melgar et al., , 2016;;Wright et al., 2012) that can provide rapid assessment of the seismic moment tensor (seismic moment and fault orientation) and/or constraints on the fault slip magnitude and along-coast extent. This information supplements direct tsunami observation by ocean bottom pressure sensors operated by the National Oceanic and Atmospheric Administration (https://www.tsunami.noaa.gov) as part of the Pacific tsunami warning system (e.g., Bernard & Titov, 2015) along with an increasing number of seafloor pressure sensors (e.g., Kanazawa et al., 2016;Tanioka, 2018) and seafloor geodetic systems (e.g., Hannemann et al., 2016;Sathiakumar et al., 2016) being operated by other countries.
In recent years, very successful warning of tsunami amplitudes at distant sites has been achieved, as there is sufficient time for analysis of seismic, geodetic, and nearby ocean bottom pressure sensor data to robustly characterize the tsunami wavefield as it expands into the ocean (e.g., Bernard & Titov, 2015;Yamazaki et al., 2012). It remains very challenging to quantify the tsunami excitation rapidly for regional warning, which requires rapid determination of the faulting geometry, seismic moment, slip distribution along the fault, and updip limit of faulting. Figure 1 compares regional maximum tsunami water height measurements (National Oceanic and Atmospheric Administration, https://www.ngdc.noaa.gov/hazard/tsu_db.shtml) with moment magnitude M W for 52 large interplate thrust events (Table 1) that are considered in this paper. The maximum water height measurements are a complex mixture of local runup, tide gauge, and deep water pressure sensor observations, so they are only a first-order characterization of the tsunami for each event. The data show systematic increase in tsunami amplitude with M W overall, with all events larger than M W 8.0 generating significant tsunami. However, there is large variation of maximum water height for magnitudes M W < 8.0. In part, this variation can be attributed to variation in rupture depth, with the symbols in Figure 1 being color coded to indicate whether there is evidence from finite-fault inversions (e.g., Hayes, 2017;Ye et al., 2016) for at least some shallow coseismic slip on the megathrust or not (Table 1). The influence of rupture depth is further demonstrated in Figure 2, which plots the same maximum water height measurements as a function of the centroid of the slip distribution for the finitefault inversions for the thrust events (Table 1). Earthquakes with slip distribution centroid depths less than 20 km tend to generate large tsunamis; these events span a large range of M W and include great ruptures like the 2011 Tohoku (M W 9.1) event and moderate magnitude (M W ~7.5-7.9) tsunami earthquakes (e.g., Kanamori, 1972) highlighted in magenta. Slip under deep water near the trench can be particularly tsunamigenic because large displacements can occur in low rigidity material near the toe of the upper plate. This holds true whether the subduction zone is filled with sediment or is an erosional zone (e.g., Ranero & von Huene, 2000;Scholl & von Huene, 2009). Thus, rapidly determining the depth extent of rupture is important for regional tsunami warning.
Rapid seismic and geodetic data processing can provide critical information about the geometry and strength of the earthquake if sufficient real-time data in the region are available, but these data have limited Figure 1. Log 10 of the maximum water height reported by the National Oceanic and Atmospheric Administration (Table 1) versus M W for the 52 interplate thrust events considered in this study. The peak tsunami amplitude or peak runup is used to define maximum water height. No tsunami was reported for the five events located at the bottom. The data are color coded by independent estimates of whether the rupture has at least some slip at shallow depth on the megathrust or not (Table 1). resolution of the far-offshore (updip) extent of thrust faulting, even for very well instrumented areas (e.g., Lay, 2018). Centroids of long-period moment tensors determined from regional data that lack good azimuthal distribution and ground displacements detected by one-sided onshore GPS generally do not constrain the updip limit of slip. Radiated near-field and teleseismic P wave energy allows shallowrupturing tsunami earthquakes to be recognized by deficiency of their high-frequency radiation (e.g., Newman & Okal, 1998;Sahakian et al., 2019;Ye et al., 2013), as does comparison of short-period P wave magnitude and long-period moment magnitude, M W (e.g., Kanamori & Ross, 2019). However, some ruptures, like the 2011 Tohoku earthquake, involve slip across the entire seismogenic plate boundary width (Lay, 2018), complicating such approaches. If the occurrence of shallow slip in part, or all, of the rupture can be determined reliably, regional tsunami warnings can be enhanced in terms of recognizing the potential for a larger-than-normal tsunami to have been excited by slip extending out to near the trench. This study presents a strategy by which that can be done in a time frame of <15 min, which can help on regional tsunami amplitude warning.
The approach presented here builds upon an initial study by Lay and Rhode (2019), which introduced the concept of using the strength of teleseismic P coda relative to P as a proxy for the occurrence of some slip or no slip at shallow depth under deep water for large megathrust ruptures. That work follows on early studies of P coda generation for submarine faulting (e.g., Ihmlé & Madariaga, 1996;Ward, 1979;Wiens, 1987Wiens, , 1989)). The underlying idea is simple and is illustrated in Figure 3. For thrust faulting on dipping plate boundary megathrusts in either island arc or continental arc subduction zones, earthquakes radiate both downgoing P and upgoing P and S that will reflect off the rock surface as depth phases pP and sP, respectively. For upgoing P paths that reach the rock surface under 1 for the 52 interplate thrust events considered in this study. No tsunami was reported for the five events located at the bottom. The data are color coded as in Figure 1. Both shallow and deep thrust events will produce direct P and surface reflected pP and sP phases (black lines and dots) at all teleseismic distances. Shallow slip, below the toe of the accretionary prism, will also generate leaky acoustic wave reverberations in the water, comprising the pwP series (blue), with the time between arrivals being larger for deeper water. Deep slip will not produce pwP for stations at landward azimuths and will produce relatively high-frequency pwP at seaward azimuths. The coda ringing produced by pwP can thus indicate whether an interplate thrust has significant slip at shallow depth under deep water or not. This information can guide rapid assessment of the likely strength of tsunami generated by a large megathrust rupture.
water, some energy transmits into the water as acoustic waves that will reflect off the water surface and travel back down to the seafloor where some energy transmits into the rock as the phase pwP, while some reflects back into the water producing a series of leaking reverberations that cause teleseismic pwP arrivals to "ring" in a sequence of arrivals with time interval controlled by the water depth. For water deeper than 6 km, the period of pwP oscillations is about 12 to 15 s, with a gradually decaying signal that can persist for many tens of seconds as P coda after the direct (P+pP+sP) wave packet has completed its ground motions at a teleseismic station. If the source involves at least some slip at shallow depth on the megathrust under deep water near the trench, reverberative pwP with a period of ~13 s will be present in the P coda at all azimuths from the source. If the source involves only slip deeper on the megathrust, close to the coast, arcward azimuths will receive little/no pwP arrivals, while seaward azimuths can have short period (<10 s) pwP oscillations. Thus, observation of high P coda amplitudes relative to P amplitudes is expected for shallow slip under deep water, with source radiation patterns for thrust events tending to make the ratios highest in the seaward direction because at those azimuths the downgoing P is more nodal than the upgoing pP that generates the pwP series (see synthetic examples in Lay & Rhode, 2019).
Conceptually, one can infer the occurrence or nonoccurrence of shallow slip by simple measures of the signal strength of P coda relative to P for a given event. Essentially, this involves remotely "hearing" prolonged booming echoes of the source sounds when slip occurs under deep water. Lay and Rhode (2019) confirm this expectation for a particular metric of teleseismic P coda /P at large epicentral distances (>80°), providing a new probe of the updip extent of rupture for megathrust events. Here we extend this approach for the practicalities of rapid tsunami warning, which requires use of data at <80°distance to reduce the time lag for decision making.
Lay & Rhode, 2019 focus on teleseismic P waves recorded at epicentral distances from 80+°to 120°where there is a several minute-long time interval between the first P arrival and the arrival of secondary seismic phases such as PP and PKIKP. For a large earthquake rupture that radiates direct P waves for a duration of less than 2 min there is about 2 min of immediately ensuing P coda duration in which no other arrival from radial Earth structure is expected. Measures of P coda amplitudes relative to P amplitudes can then reliably assess whether the P coda is augmented by persistent pwP near-source reverberations that indicate rupture under water or just has typical levels of scattered energy from propagation through the heterogeneous 3-D Earth. Summing measures from many stations suppresses effects of individual station/path complexity and enhances sensitivity to energy emerging from the source region. Band-pass filtering the data to emphasize the period range of pwP expected for rupture under deep water (Lay & Rhode, 2019; use a passband of 7 to 15 s) near the trench increases the sensitivity of the measures to the updip limit of slip during the earthquake as suggested by Figure 3.
Use of large distance P waves results in significant delay after the earthquake, given that the first P arrival at 80°will have a travel time of about 12 min, while that at 120°will take about 15 min to reach the station. The duration of the P wave generated during rupture (composed of all the P+pP+sP arrivals from the space-time distribution of slip on the megathrust) plus an interval of P coda sufficient for stable measurement must be allowed for, giving an effective delay of ~18+ min for event-averaged assessment of the relative amplitude of P coda to P. This can certainly contribute to regional and remote tsunami warning but is clearly of limited use for local tsunami warning. P waves at an epicentral distance of 30°have travel times of about 5 min, so use of signals at closer distances can cut minutes off the contribution to tsunami warning. Even closer distance might be considered, but the strong pwP ringing involves steep raypaths (small ray parameters) and is more pronounced at teleseismic distances where the P waves dive steeply into the mantle (P n phases have weaker pP n signals and hence will have lower pwPn amplitudes). Triplications from upper mantle discontinuities add complexity and total duration to direct P arrivals at distances from 15°to 30°. At close distances, strong S phases and surface waves may obscure the P coda signal as well. Generation of seismic coda observed at distances less than 30°still does warrant further investigation, but here we only consider observations of P coda at distances larger than 30°, where direct P waves are simpler. We seek to reduce the measurement delays to <15 min for contribution to local tsunami warning, which is somewhat longer than the lag time required for stable moment tensor inversions using long-period signals for well-instrumented regions. The 10.1029/2019JB018221 Journal of Geophysical Research: Solid Earth LAY ET AL. 10,589
added complication of the closer distance range is that secondary phases such as PP, PPP, and PcP (usually small) cut through the P coda time interval as shown in Figure S5 of Lay and Rhode (2019).
Figure 4 shows representative profiles of global seismic P wave observations, involving ground velocities filtered in the 7-to 15-s period range, plotted over the distance range 30°to 120°. The 2003 Hokkaido earthquake has been shown to not have coseismic slip extending to near the trench (e.g., Miyazaki et al., 2004), while the 2013 Santa Cruz Islands event does have substantial, tsunamigenic slip at shallow depth below deep water (e.g., Lay et al., 2013). The profiles, aligned on the P wave, have red lines bracketing the interval for the direct P (+pP+sP) ground motions (first interval) as defined by the source duration of the moment rate functions for the earthquakes inferred from finite-fault inversions (Ye et al., 2016), with 20 s added to account for the last depth phase arrivals. The subsequent interval is the P coda duration, which was defined in Lay and Rhode ( 2019) by the time from the end of the P window to the arrival time of PP at ~80°. For the 2003 event, clear PP, PPP, and S phase branches are apparent in the ground velocities, highlighted by dashed blue lines. PP, in particular, moves through the P coda window in the distance range 30°to 70°. For the 2013 event, note that the discrete travel time branches are less apparent at all ranges and the ground motions in the P coda window are comparable to those in the P window, in contrast to the 2003 profile. The azimuthal distributions of global network broadband stations used (GSN, FDSN, and Antarctica BB) are shown in the maps relative to the hypocenters. The azimuthal distributions are nonuniform, and there is an intrinsic tendency for stations at seaward directions to be sparse, while clusters of continental stations occur at landward directions; this motivates an azimuthal binning procedure when generating eventaveraged parameters.
The data beyond 80°in Figure 4 are particularly simple to evaluate because the P coda levels are not contaminated by expected secondary phases. Measurements for each station, in a given passband, of the ratios of root-mean-square (RMS) amplitudes of P coda (for the time window following the P window up to arrival of the PP phase at 80°) to RMS amplitudes of P (for the time window spanning P+pP +sP arrivals from source slip), are denoted as RMS_C/P. The mean or median values for the stations in 10°azimuth bins are determined, and the overall mean or median of the binned values is computed. This provides easily obtained event-specific measures sensitive to the presence or absence of strong pwP. Lay and Rhode (2019) found that event mean RMS_C/P varies by a factor of 4 among 39 major and great subduction zone thrust events. RMS_C/P >0.61 was found for all events (19) with independent finite-fault models indicating at least some slip located at shallow depth on the megathrust below deep ocean, while all events (16) except for two with no indication of shallow slip have RMS_C/P <0.61. One of the exceptions has a readily recognized early aftershock in the P coda window (rather than ringing pwP), while the other has a relatively uncertain finite-fault model that is being investigated further. Two events have shallow slip, but lack a nearby deep trench, so they should not have long-period pwP; indeed, they have low RMS_C/P values. Computing medians of the 10°azimuthal bins gives similar results with slightly greater separation of the bimodal populations (see Figure S8 in Lay & Rhode, 2019).
The data for the 39 interplate thrust events considered by Lay and Rhode (2019) are augmented here to include broadband observations for the distance range 30°to 120°, and corresponding data for 14 more large megathrust earthquakes from 1990 to 2018 are added (although the 2004 Sumatra M W 9.2 earthquake is not used for coda measurements due to its extremely long P duration). For the new events we again draw upon published finite-fault models (e.g., Hayes, 2017;Ye et al., 2016) to constrain the P wave duration and the bimodal classification as having shallow coseismic slip or not. In addition, corresponding data for 34 large shallow earthquakes in different tectonic environments (outer rise compressional, intraslab compressional, back-arc compressional, continental thrusts, outer rise extensional, intraslab extensional, oceanic strike slip, and continental strike slip) are analyzed. Among the latter, the thrust-and normal-faulting dip-slip events in oceanic environments are also of concern for tsunami warning. Table 1 lists the hypocentral parameters and M W of the 87 earthquakes, along with estimates of the centroid depth of the finite-fault slip distributions (H c ) and a priori classification of whether or not the event has shallow slip. Figure 5 shows the global centroid moment tensor focal mechanisms for all of the events, color coded to separate the subduction zone megathrust events (blue) from the other categories of events (red). The data set includes all useable shallow events larger than M W 7.5-10.1029/2019JB018221
LAY ET AL. 10,590
7.6 for all mechanism types from 1990 to 2018 along with a handful of smaller subduction zone megathrust events included to provide a good balance of deep-slip and shallow-slip events with good sampling of the regions that generate most tsunamis. Intrinsically, regions that have experienced few or no recent large earthquakes like Cascadia, Izu, Marianas, and Ryukyus are not well represented.
Measurements are made for all events following the basic procedure introduced by Lay and Rhode (2019); for each new event the data from 80°to 120°are used to compute the RMS_C/P ratios at all stations using the 7to 15-s band-pass filtered velocity signals, and these are binned in 10°azimuthal windows. For the binning, we compute the median of the ratios in each bin (rather than average) and then we compute the median of the bin medians to give the final event median RMS_C/P value. For the original 39 events, we use the same data from Lay and Rhode (2019); some of those events had original distance windows of 90°to 120°or 100°to 120°and included all broadband data from temporary deployments yielding many observations. We compute median RMS_C/P using medians of 10°azimuthal bins for the original distance windows to provide reference measurements with continuity to the earlier study (designated with ranges 80+°to 120°), given the near-perfect bimodal separation of the two populations of slip depth it provided. Very small differences are found using all data in the 80°to 120°range from just GSN, FDSN, and Antarctic stations for those 39 events (the respective median RMC_C/P measures have a correlation coefficient r = 0.970).
The median RMS_C/P measures for all of the interplate thrust events are plotted in ascending order in Figure 6. The corresponding values are listed in Table S1 in the supporting information. The range of the 40th to 60th percentile measurements for each event is included as an indication of stability of the median values (these bounds generously approximate a standard error of the median for the skewed distributions).
Color coding is used to indicate the independent inferences of whether there was at least some shallow slip during the rupture (blue) or not (red). The distribution is very similar to that for the 39 events in Lay and Mechanisms are color coded for interplate thrust events (blue) and intraplate or strike-slip transform fault events of various faulting geometry (red). Based on examination of published finite-fault models, each interplate thrust event is characterized as either suspected to have some shallow coseismic slip (or entirely shallow coseismic slip for tsunami earthquakes) or suspected to not have significant shallow slip (Figure 6). Specific event information is given in Table 1.
10.1029/2019JB018221 Journal of Geophysical Research: Solid Earth LAY ET AL. 10,591 Rhode (2019), with the same outliers, but there are 13 new events and our use of median of bin medians slightly reorders some events. One new event (2002 Papua New Guinea) with shallow slip occurred in a region with no deep trench, giving a total of three such cases (cyan). Two new events have moderately high RMS_C/P but uncertain slip models; 17 December 2016 Solomon Islands (e.g., Lay et al., 2017; https://earth- quake.usgs.gov/earthquakes/search/; Lee et al., 2018) and 5 May 2015 Papua New Guinea (e.g., Ye et al., 2016; https://earthquake.usgs.gov/earthquakes/search/), giving three such cases. The two new events involve initially deep ruptures (the 2016 Solomon Islands event likely begins as an intraplate rupture) that
propagate or trigger updip, and while the finite-fault models do suggest that some slip may extend to near the trench, the resolution is limited and the amount of shallow slip varies among the finite-fault models.
The 2016 Solomon Islands event produced outer rise normal faulting aftershocks along part of the rupture, a common indicator of shallow coseismic slip (e.g., Sladen & Trevisan, 2018;Wetzler et al., 2017), while the 2015 Papua had shallow thrust-faulting aftershocks all the way to the trench (Wetzler et al., 2017). Such events, with relatively deep moment tensor centroid depths, present particular difficulties for tsunami warning. One new event (5 March 2002 Mindanao) does not have evidence for shallow slip in available models (e.g., Ye et al., 2016; https://earthquake.usgs.gov/earthquakes/search/) but locates within the population that does, giving one clear outlier to the general separation of the populations that cannot be readily accounted for unless the published finite-fault models are found to be misleading, which is being examined separately.
If we only consider the red values (events with no shallow slip) in Figure 6, the mean value of median RMS_C/P for that population is 0.464 with a standard deviation of 0.087. The mean value of the blue RMS_C/P values (events with some/all shallow slip under a deep trench) is 0.859 with a standard deviation of 0.277. Observation of correspondingly measured median RMS_C/P larger than ~0.65 is strongly suggestive of pwP enhancement of the coda, whereas values less than ~0.55 are strongly suggestive of no shallow slip in the rupture. Intermediate values are ambiguous, and this includes some modest size tsunami earthquakes (names highlighted in magenta). There are significant changes in the availability of global seismic data over time; more recent events have substantially larger data sets with more complete azimuthal sampling, so earlier events (pre-2000) tend to have larger uncertainty. Rather than diving into a complex statistical classification effort (the subpopulations means are statistically distinct, but the actual measurements present a continuum rather than a bimodal distribution), we now consider the prospect of achieving comparable first-order population segregation using data at closer distances, as desirable for more rapid tsunami warning application. 10.1029/2019JB018221 Journal of Geophysical Research: Solid Earth LAY ET AL. 10,592
We used the same P and P coda windows to compute RMS_C/P for the data for each event for distance ranges of 30°to 120°and 30°to 80°. The expanded sampling provided by these data allows us to compare 10°azimuthal binning with 10°distance binning in the computation of medians of azimuthally binned or distance-binned medians. Referring to the examples in Figure 4, it is clear that 10°distance binning will lead to P coda measures at less than 60°being enhanced for events without high coda but only modestly affected for events with high coda. Variation of area covered by distance bins at different ranges is inconsequential given the highly nonuniform distribution of stations. Using 10°azimuth binning will tend to mitigate the effects of PP contamination at closer distances by computing the median for azimuths with a wide range of distances and prevents dense continental station distributions from being overrepresented and sparse data at seaward directions from being underrepresented. Figure 7 compares the resulting median RMS_C/P values with those for the 80+°to 120°values from 10.1029/2019JB018221 Journal of Geophysical Research: Solid Earth LAY ET AL. 10,593 Figure 6. The corresponding measurements are all listed in Table S1. The same color coding as used in Figure 6 is used for the comparisons in Figure 7.
There is very strong correlation between the results for the 30°to 120°estimates and the 80+°to 120°values for both azimuthal binning (r = 0.934) and for distance binning (r = 0.966). The data do overlap at large distances although we cull the data set spanning 30°to 120°to remove all non-FDSN networks and transportable station deployments from the original large distance data set of Lay and Rhode (2019). The separation of populations is still very good, but there is more overlap for values around 0.6 for the 30°to 120°measures than for the 80+°to 120°measures. For the 30°to 80°estimates, increased influence of PP and PPP biases in coda amplitudes is expected, with less efficacious suppression by binning. However, the correlations with the 80+°to 120°measures are still quite good for both azimuthal binning (r = 0.845) and distance binning (r = 0.834).
The median RMS_P/C measures for data from 30°to 80°with 10°azimuthal binning are plotted in sorted order in Figure 8, for comparison with Figure 6. There is more mixing of the populations for values around 0.70 ± 0.05 but still very good separation of the populations for larger and small values. The Peru 21 February 1996 tsunami earthquake is now an outlier, with low coda level. This early event has very limited data at the closer distances in the seaward direction, so it is dominated by arcward data that lack strong pwP. The mean value of the red (no shallow slip) median RMS_C/P values in Figure 8 is 0.600, with a standard deviation of 0.087. The mean value of the blue (some/all shallow slip with deep trench) values is 0.934, with a standard deviation of 0.222. Values >0.8 are strongly indicative of presence of shallow slip with potentially enhanced tsunami amplitudes. This figure provides a summary of the potential for rapid determination of RMS_P/C from data at distances of 30°to 80°, giving 3-4 min faster contribution to tsunami excitation assessment than provided by the measurements in Figure 6.
It may be possible to enhance the performance for closer distance RMS_P/C measures by "masking out" the signals of PP and PPP or by using distance-varying time windowing for the P coda duration, but we have kept this initial demonstration very simple, as is desirable for operational consideration. The P signal window can be quickly estimated from twice the centroid time for a rapid W phase moment tensor inversion, from time energy integrals of the broadband waveforms, or from simple inspection of the broadband data profiles, and a subsequent, comparable P coda duration window defined with very straightforward data processing (instrument correction and band-pass filtering) with no other decisions, just time lagged for enough stations to provide good data averaging to constrain the event summary parameter. 10.1029/2019JB018221
Journal of Geophysical Research: Solid Earth
The waveform profiles in Figure 4 suggest that frequency domain approaches should be promising for characterizing the P coda /P behavior as an alternative to the time domain RMS amplitude measures discussed so far. Each seismogram could be converted to a spectrogram, seeking persistent narrowband features that correspond to prolonged pwP ringing. As a first approach to such frequency domain measures, intended to keep procedures simple, we compute the spectra for 5-to 20-s band-pass filtered ground velocity waveforms in the Figure 9. Stacked spectra and spectral ratios of P (blue) and P coda (red) ground velocity waveform segments of equal length band-pass filtered between 0.05 and 0.2 Hz for events with no shallow slip (top row) and events with shallow slip (bottom row). Individual station P coda spectra were normalized to unity at 0.083 Hz, retaining correct relative amplitude of P spectra, and the median values at each spectral point for medians of 10°azimuthal bins for all P and P coda measurements in the 80°to 120°distance range are plotted in each spectral comparison. Median values of the spectral ratios are shown below the spectra. The horizontal light magenta lines indicate the 7-to 15-s passband used for root-mean-square (RMS) calculations in RMS_C/P. The short horizontal bold magenta lines indicate the 12-to 13-s spectral window used for narrowband spectral ratio comparisons shown in Figure 10. 10.1029/2019JB018221 Journal of Geophysical Research: Solid Earth LAY ET AL. 10,595
P and P coda signal windows, stacking the observations to obtain event-specific spectra. We use the same P durations as in the RMS analysis, but we now make the P coda durations exactly the same as the P durations for each event, to ensure uniform signal window lengths for spectral power estimates. The spectral amplitude at ~12-s period for each P coda window is set to unity, with the relative scaling of the corresponding P window spectra being retained, and the spectra are stacked separately (using median values at each frequency for medians in 10°azimuthal bins) along with corresponding stacking of the P coda /P spectral ratios for all stations in various distance ranges. Figure 9 shows examples of the median spectra for four events with no shallow slip (Figures 9a to 9d Thin magenta horizontal lines in the lower part of each spectral panel of Figure 9 indicate the band pass used for the RMS calculations (7 to 15 s or 0.067 to 0.143 Hz). The relative P and P coda spectra in this passband control the measures in Figure 6. For the events with no shallow slip (top row), the P spectra (blue) are higher than the P coda spectra (red) across the full bandwidth of the RMS measures, whereas for the events with some shallow slip (lower row), the P coda spectra are higher than or similar to the P spectra over the lowerfrequency part of the passband. The large peaks in P coda spectra for the 1995 Chile and 2003 Aleutian events correspond to the high RMS_C/P values for these stations in Figure 6. Note that the stacked spectral ratios (Figures 9e′ and 9f′) have clear peaks around 0.07 Hz for these two events. The 2010 Mentawai earthquake has a spectral ratio peak near 0.10 Hz (Figure 9g′), consistent with the shallower depth of the trench in the source region. The spectra for the 2011 Tohoku event (Figure 9h) have relatively high spectral ratios (Figure 9h′) over the whole bandwidth but no single large peak; this is likely due to the very broad extent of large slip from 30-to 7-km depth for this event (e.g., Yamazaki et al., 2018), which generated pwP with a broad spectrum of interfering resonances. As in Figure 6, the P coda /P behavior clearly suggests the presence of shallow slip for 2011 Tohoku but not as dramatically as one might have expected given the 50+-m estimates of slip near the trench (e.g., Fujiwara et al., 2011;Sun et al., 2017). The lack of seaward observations with stable signals is possibly a factor in this, but the wide depth distribution of slip is likely the dominant factor.
Noting the large, narrowband peaks in the P coda spectra for the 1995 Chile and 2003 Aleutian events in Figure 9, targeting the corresponding frequency band more selectively may hold promise of improving the recognition of deep water pwP. We take the average P coda /P spectral ratio at periods from 12 to 13 s (0.077 to 0.083 Hz; the range of the short bold magenta lines in the lower portion of each spectral plot in 10.1029/2019JB018221 Journal of Geophysical Research: Solid Earth LAY ET AL. 10,596
Figure 9), as a very specific search for deep water pwP. For the Tohoku event there is a small peak around 0.083 Hz in the spectral ratios as well, and much lower spectral ratios are observed for the events without shallow slip. Again, computing medians of medians in 10°azimuth bins for these 12-to 13-s period spectral ratios, we compare their behavior with the median RMS_C/P measurements with similar binning for 80+°to 120°and 30°to 80°data sets in Figure 10. The event values are listed in Table S2. There is larger scatter than for comparison of RMS measurements as in Figure 7, but low or high spectral ratios and RMS_C/P do separate the populations well. The correlation coefficient for the 80°to 120°measures is r = 0.684, while that for the 30°to 80°measures is r = 0.770. There is somewhat more population overlap in this narrowband spectral ratio, but this is likely an outcome of variable frequency deep water resonances due to trench depth variations among the different events, combined with the greater intrinsic variance of narrowband spectral measures. Consideration of full spectrograms that retain broader bandwidth may provide stabilization with respect to both of these aspects, but that is best suited to a neural network-based classification approach rather than the simple parametric approaches in this paper; that approach is being pursued separately.
Magnitude-based measures are also attractive for rapid applications in tsunami warning, largely due to the simplicity of the measurements and their widespread familiarity. While the RMS_C/P measures are very straightforward, they do not provide absolute source strength levels due to lack of geometric spreading corrections or direct calibration with M W . Tsunami warning approaches currently make use of the broadband body wave magnitude M wP (e.g., Hirshorn et al., 2013) and W phase M W , but there is also value in the classic body wave magnitude m B when computed for P wave signals with periods around 3.5 s (e.g., Kanamori & Ross, 2019). As a relatively short period measure, comparison of m B with M W can robustly isolate tsunami earthquakes due to their deficient radiation of shorter period signals. We apply the same simple processing as used for m B to our data and measure event median m B values for each event. Then we use the same magnitude formulation to develop a corresponding coda magnitude, m_coda, for the period range of pwP, and consider a differential magnitude m_codam B , as a proxy for pwP excitation, comparing it with RMS_C/P.
The procedure used for m B calculation exactly follows that of Kanamori and Ross (2019); broadband vertical component recordings are corrected to Weichert-type seismograph recordings; the instantaneous amplitude (A P ) and period (T P ) are measured and corrected to ground motion amplitude. The centroid time relationship t c = 1.2 × 10 -8 M 0 1/3 (Duputel et al., 2013; with M 0 being the independently estimated seismic moment in dyne cm) is used to define a P wave duration window of 3t c (cutting it at the arrival time of PP, if necessary), and the magnitude is computed using m B = log 10 (A P /T P ) + Q (Δ, h), where Q(Δ, h) is the Gutenberg and Richter distance (Δ), depth (h) correction (Richter, 1958). The distance range for which the latter corrections are reliable does not extend to diffracted distances, so we limit magnitude measurements to the 30°to 80°r ange. Dominant periods of the measurements used in the m B magnitudes are ~3 to 4 s.
The same magnitude formalism is used to estimate a P coda magnitude, m_coda, but the broadband ground displacement data are filtered in the passband 8 to 15 s (typical T P are in the range 11 to 14 s), and the P coda time window is from 5t c to 8t c (or until the S arrival). We also consider durations of 3t c to 6t c , but for several events with long tails on their moment rate functions the later window better isolates the coda. As is the case for the RMS_C/P measures, PP and PPP contaminate the m_coda measures at closer distances; hence, median values are used. The m_coda measures have the added attribute of being magnitudes that retain source strength information, but they are also single peak measures rather than RMS values so they are more sensitive to contamination from large late arrivals like PP, whereas RMS values intrinsically emphasize sustained large amplitudes throughout the coda window. The m B and m_coda magnitudes and their differentials are listed in Table S2.
The differential median m_codamedian m B measures (30°to 80°) are compared with the 80+°to 120°m edian RMS_C/P values in Figure 11. As a result of measuring the ~12-s period coda strength versus the ~3.5-s period P strength in the magnitude differences, the differentials are particularly effective at separating the moderate size tsunami earthquakes from other events of similar size that rupture to shallow depth. The magnitude differences have more overlap of the bimodal population, and the correlation coefficient with the RMS measures is 0.598, but at the end of the range there is very good separation for shallow-slip (small magnitude difference) and deep-slip (large magnitude difference) events. The corresponding comparison using
10.1029/2019JB018221 Journal of Geophysical Research: Solid Earth LAY ET AL.
10,597 the 30°to 80°median RMS_C/P values is shown in Figure 12, with similar general behavior and an overall correlation coefficient of 0.730. Using the 3t c to 6t c window for the coda duration for m_coda (see Figures S2 and S3) gives measures with ~0.16 lower correlation coefficients for the two distance ranges.
Given that all of the RMS and spectral and magnitude differential measures can be processed very rapidly, various measures may be combined in artificial intelligence frameworks to optimize population classification. Again, this is left for future efforts; the bimodal classification of updip extent of rupture that we have adopted is very simplified and may not represent the continuum of actual earthquake behaviors that contribute to fuzziness and overlap of the parametric measures that we explore here. That aspect of the problem must be considered further in more complex classification schemes. Practically speaking, very sophisticated approaches (such as finite-fault modeling or other procedures requiring many user decisions) are of limited use for rapid tsunami warning decision making. The simple parametric behaviors observed here give the primary contribution of coda-based approaches; observation of high relative P coda levels (large RMC_C/P, small m_codam B , large 12-to 13-s P coda /P spectral ratios) suggests shallow slip which is expected to enhance tsunami excitation, and in combination with a W phase determination of mechanism geometry and seismic moment, the added information about rupture depth can add confidence in the warning and provides emphasis that a larger than typical tsunami might have been generated.
Thus far, our focus has been on known interplate thrust-faulting ruptures in subduction zones. There are other types of large earthquakes of concern for tsunami generation, including interplate thrust-faulting events, normal-faulting events, and some near-coastal strike-slip events. Figure 5 shows the global distribution of large earthquakes from 1990 to 2018, highlighting the nonmegathrust events in red. For all 34 of the latter events, we computed RMS_C/P, m_codam B , and 12-to 13-s period P coda /P spectral amplitude ratios (Tables S1 and S2), using the same procedures applied to the megathrust events. Comparison of m_coda - 10.1029/2019JB018221 Journal of Geophysical Research: Solid Earth LAY ET AL. 10,598
m B and 80-120°RMS_C/P measures in Figure 13 establishes that there is a large range in relative coda levels, comparable to that for interplate thrust events, and there is good correlation (r = 0.734) between these parametric measures. The events are color coded according to mechanism type.
Strike-slip events systematically tend to have larger relative coda levels, particularly for the RMS measures, which is probably a manifestation of both close to nodal downgoing P radiation and strong near-source scattering due to the P and SH radiation patterns having radiation lobes in the horizontal directions. Note that the strike-slip events with the largest relative coda levels are interplate ruptures on transform faults, which rupture primarily within the crust. Intraplate strike-slip events have lower coda, notably for the large events in the Wharton basin, which have relatively large centroid depths (Lay, 2019). The two continental strike-slip events (2001 Kunlun and 2002 Denali) have moderate to low coda levels that must be generated by scattering rather than pwP. The same is true for the six continental thrust events, and it is notable that the coda levels for these events span a range similar to that for suboceanic thrusts in backarcs, outer rise oceanic lithosphere, and subducted slabs. The latter three populations have measures comparable to those for interplate thrusts that lack shallow slip. The back-arc thrust events are similar to the outliers in Figures 6 and 8 for thrust events with shallow slip but no deep trench in which to generate long-period pwP. The 1999 Vanuatu event did generate local 7-m tsunami runup (Regnier et al., 2003), so these events can be tsunamigenic, but their locations immediately identify them as distinct from megathrust events. The intraplate thrusts in the outer rise and in slabs tend to have greater depths than the megathrust ruptures, accounting for their weak excitation of pwP contributions to P coda . Tsunamigenic normal-faulting events in the outer rise (e.g., 2007 Kuril and 2009 Samoa) have higher coda levels than deeper intraslab normal fault ruptures such as the 2001 El Salvador and 2017 Chiapas.
Given that some of these events are tsunamigenic, we compare the RMS_C/P measures obtained for different distance intervals (30°to 120°and 30°to 80°) with the results for 80°to 120°(free of any PP contamination) in Figure 14. The correlations are high (r = 0.932 to 0.978), indicating that measures made at closer distances 10.1029/2019JB018221 Journal of Geophysical Research: Solid Earth LAY ET AL. 10,599
can help in tsunami warning assessment for these nonmegathrust events, informed by the moment tensor faulting geometry. Given an offshore faulting mechanism with significant dip-slip component, the P coda excitation levels can again contribute to assessment of likelihood that shallow faulting under deep water has enhanced tsunami excitation or not. Using P coda to recognize the shallow rupture of the thrust events under shallow water will require calibration of signals in the frequency band 5-7 s. For strike-slip events, the indication of shallow depth rupture for the transform faulting environments versus the intraplate environments also provides useful information, although tsunami excitation will intrinsically likely be modest with possible exception of events that rupture from offshore to onshore.
Simple parametric measures of seismic wavefields from large offshore earthquakes can be obtained in near real time due to open access to and telemetry of global broadband stations. This data availability underlies the ability to quickly determine broadband magnitudes, to invert long-period ground motions for the moment tensor, to evaluate radiated energy, and even to formulate rapid finite-fault inversions. However, resolving the updip extent of slip in thrust faulting events, particularly on megathrusts, remains challenging for rapid analyses. Onshore geodetic data have expanding global deployment and rapid processing but also intrinsically do not resolve the far offshore slip distribution for megathrust failures. Real-time data from offshore geodetic sites can resolve this problem but remain very limited in coverage globally. When rupture extends to near the trench, there is enhanced generation of tsunami signals for an event of a given size (Figure 2). This is the essence of shallow-rupturing tsunami earthquakes (Kanamori, 1972), which may be 10.1029/2019JB018221 Journal of Geophysical Research: Solid Earth LAY ET AL. 10,600
quickly recognizable based on their deficient short-period amplitude spectrum, but it also holds for ruptures spanning a wide portion of the megathrust such as the 2011 Tohoku and 2010 Chile earthquakes (e.g., Maksymowicz et al., 2017;Yue et al., 2014), which do not register as tsunami earthquakes. As shown here, assessment of long-period pwP seismic wave arrivals, which are expected to be preferentially generated by shallow slip under deep water, can help to resolve this issue.
When a large earthquake occurs, if regional seismic and seismogeodetic ground motions are available for rapid processing, the faulting mechanism and seismic moment or average fault slip and the along-strike faulting dimensions can be determined within ~3 to 15 min after the hypocentral time, depending on event location and station distribution. To supplement that information, the various parametric P coda /P measures described here can add the important nuance of whether the faulting extended into deep water (near the subduction zone trench), and if so, that tsunami excitation may be stronger than typical for the event's magnitude by up to the factor of 2 to 3 increase observed for tsunami earthquakes relative to comparable 10.1029/2019JB018221 Journal of Geophysical Research: Solid Earth LAY ET AL. 10,601
seismic moment events deeper on the megathrust (Figure 2). The RMS_C/P and other coda measures alone do not provide direct prediction of tsunami amplitudes; that can be inferred best in combination with the W phase moment tensor and/or geodetic observations. What is added by these measures is information about the updip extent of rupture, which is hard to rapidly constrain otherwise (and intrinsically, improved constraint on the overall rupture depth, which can be quite uncertain in rapid W phase inversions using sparse data). Figure 15 compares the median RMS_C/P measures for interplate thrusts using the 80+°to 120°data set with the H c values from finite-fault models (Table 1). Large coda levels are strong indicators of shallow slip under deep water, which is connected to shallow H c . Low RMS_C/P is strongly indicative of deeper H c . The coda information can be obtained within ~15 min after the event starts, and as more distant data become available, increasing confidence can be attached to the resulting measures and inferences. This is much sooner than the determination of reliable slip models from finite-fault inversions needed for direct estimation of H c .
Figure 16 directly compares median RMS_C/P measures with log 10 (maximum water height; logMWH). There is a general tendency for larger coda to be associated with larger tsunami amplitude, but it is important to note that absolute event size information is not represented by the coda ratio (the highly scattered variation of these RMS_C/P measures with M W is shown in Figure S4). In order to account for the contribution of event size, we consider the two populations in the interplate thrust event data set in Figure 1. Figure 17 shows separate regressions of logMWH on M W for the event subpopulations with at least some shallow slip (blue points) or no shallow slip (red points), omitting the ambiguous events with no tsunami reported, uncertain faulting models, or shallow slip but no deep water trench. The two populations have some scatter and overlap, but the influence of shallow rupture is strongly indicated, with convergence of the curves for very 1 for the 52 interplate thrust events considered in this study.
The data are color coded as in Figure 1.
Journal of Geophysical Research: Solid Earth large ruptures near magnitude 9. The regression for the events with no shallow slip has a slope close to 1, indicating that for deeper rupture, maximum water height roughly scales directly with M W for this subpopulation.
We use the behavior in Figure 17 to adjust the maximum water height measures for effects of varying M W , allowing us to compare directly with the coda measurements. To remove the magnitude effect, we assume the form shown in Figure 17 for the events lacking shallow slip:
Since MWH given in Table 1 is only a first-order characterization of the tsunami, here we use a slope a = 1.0 for simplicity, with corresponding b = -7.9. Since Abe's (1979) tsunami magnitude, M t , defined as a linear relationship, M t = logH+B (H is the tsunami amplitude), is generally compatible with M W , the choice of a = 1 is reasonable. Then, Δ(logMWH) ≡ logMWH -(aM w +b) can be interpreted as the enhanced water height caused by shallow slip.
Figures 18-20 show how Δ(logMWH) depends on the various measures of shallow slip. Figure 18 plots the resulting magnitude-corrected maximum water height values versus H c , demonstrating the clear influence of slip depth on enhancement of the water height. This figure can be compared to Figure 2 to see that the trend is slightly enhanced by the magnitude correction. The correlation coefficient for only the red and blue populations in Figure 2 is -0.529, whereas for the same subgroup in Figure 18 it is -0.575. Figure 19 Figure 16. Log 10 of the maximum water height reported by the National Oceanic and Atmospheric Administration versus median RMS_C/P for medians of 10°azimuthal bins of the 80+°to 120°data set for the 52 interplate thrust events considered in this study. No tsunami was reported for the five events located at the bottom. The data are color coded as in Figure 1.
Figure 18. Log 10 of the maximum water height reported by the National Oceanic and Atmospheric Administration, scaled for event size by subtracting (1.0 × M W -7.9), versus H c , the centroid of the slip distribution for the finite-faults reported in Table 1 for the 52 interplate thrust events considered in this study. No tsunami was reported for the five events located at the bottom. The data are color coded as in Figure 1. 10.1029/2019JB018221 Journal of Geophysical Research: Solid Earth LAY ET AL. 10,604
makes a similar comparison of magnitude-corrected maximum water height values versus median RMS_C/P for the 80+°to 120°data with 10°azimuthal bin medians. This can be compared with Figure 16, with the effects of M W scaling of the maximum tsunami height now being suppressed. The scatter is substantial, and the correlation does not change significantly (for the red and blue points it is 0.361 in Figure 16 and 0.346 in Figure 18). It is quite possible that large pwP can be generated by modest slip at shallow depth, so there need not be strong tsunami excitation for a large RMS_C/P. The 2003 Aleutian event is an example; it has a high coda level but only small reported tsunami. For this event slip does occur at shallow depth but it is only ~1 m (Ye et al., 2016). The strong coda produced indicates high efficiency of generation of pwP by relatively minor slip. Alternatively, tsunami can be enhanced by submarine slumping; this appears to hold for the 1998 Papua tsunami earthquake (e.g., Synolakis et al., 2002;Tappin et al., 2008). Figure 20 compares the magnitude-corrected maximum water height values with m_codam B measurements. This comparison has a relatively high correlation of 0.591 for the subset of red and blue points, improved from a 0.302 correlation with the uncorrected values. The enhanced separation of the tsunami earthquakes by the differential magnitude measurement is responsible for the stronger correlation.
We can perform linear regressions on how Δ(logMWH) depends on the various parameters, p:
to determine the coefficients c and d. The linear representation is entirely empirical, not physics based. Once c and d are determined, we can estimate the enhanced water height using (2) for a measured p. Alternatively, 10.1029/2019JB018221 Journal of Geophysical Research: Solid Earth LAY ET AL. 10,605
since the slope of equation ( 1) is 1, Δ(logMWH) can be also interpreted as an effective increment of magnitude, ΔM w = cp+d, if we are to estimate the enhanced water height in terms of an increased M W .
For example, the linear regression of ΔlogMWH on median RMS_C/P (80°+/120°) has slope c = 0.7702 and intercept d = -0.0273, and the same regression on m_coda-m B has c = 1.99 and d = 1.14. These regressions omit events with no tsunami detection, uncertain slip models, shallow slip under shallow water, or early aftershock contamination. Using these values of c and d, we get ΔM W = 0.54, 0.54, and 0.22 (for p = RMS_C/P) and 0.51, 1.02, and -0.06 (for p = m_coda-m B ) for the 2011 Tohoku (M W = 9.1), 1992 Nicaragua (M W = 7.7), and the 2003 Hokkaido (M W = 8.2) events, respectively. This illustrates how the occurrence of coda-generating shallow slip can be mapped into an effective magnitude increase that controls maximum water height. We tabulate ΔM W for all of the earthquakes used in the regressions in Table S3.
This regression-based approach can be formalized with propagation of errors, but the ambiguous nature of the maximum water level values and the scatter in the correlations with coda parameters may not warrant very precise calculations. As a conservative approach, tsunami warnings based primarily on rapid faulting mechanism, seismic moment, and coastal deformation observations can be given added emphasis on the tsunami likely being higher than normal if the coda measures are indicative of at least some slip on the shallow megathrust having occurred.
It is also important to recognize that large coda can be observed for some events with only modest shallow slip (and low tsunami excitation). The 2013 Aleutian and 1995 Chile events are clear examples of this. The generation of pwP is, of course, more complex than simple vertical reverberations above the source, and in detail it depends on 3-D geometry of the bathymetry and the presence of low-velocity sediment structures. Detailed 3-D modeling of subduction environments (e.g., Qian et al., 2019;Wu et al., 2018) may help to 10.1029/2019JB018221 Journal of Geophysical Research: Solid Earth LAY ET AL. 10,606
identify improved metrics sensitive to occurrence of shallow slip and pwP excitation along with recognition of structural scattering that can generate pwP laterally offset from the slip zone (e.g., Fan & Shearer, 2018;Yue et al., 2017).
Use of ~12-s P coda amplitude relative to P as a measure of pwP generation by slip at shallow depth under deep water during megathrust ruptures has potential to enhance tsunami warnings, given that such slip can be particularly tsunamigenic. The initial work of Lay and Rhode (2019), using observations at large distances (>80°) for which there is almost always a seismogram interval free of expected arrivals between the end of the direct P wave signals from the rupture and the ensuing PP and PKIKP phases, established that higher P coda levels relative to P correlate well with independent determinations of slip having occurred at shallow depths under deep water. This idea is extended here to closer distance data (30°to 80°) to reduce the lag time necessary to evaluate the relative level of P coda /P to less than 15 min, making it more useful for tsunami warning decision making. Simple parametric measures, such as 7-to 15-s period RMS P coda /RMS P (RMS_P/C), m_codam B , and 12-to 13-s P coda /P spectral ratios, provide strong indication of the presence of some shallow slip under deep water (high values) or lack of slip under deep water (low values), so that nuance can be added to tsunami warnings based on seismic moment and faulting geometry (shallow slip under deep water suggesting potential for stronger than typical tsunami for a given event size). Calibration of the coda measurements for shorter periods may allow shallow slip under shallower maximum water depths to be recognized. It would be straightforward to precalculate the optimal pwP passband to consider for regions with different trench depths. Information about the slip distribution on the megathrust gleaned from P coda levels has many other applications: for parameterization of finite-fault models, for evaluation of potential for ensuing tsunami earthquakes updip of megathrust failures that do not reach to the trench, and for assessment of frictional properties as a function of depth on subduction zone megathrusts.