A Decade of Short‐Period Earthquake Rupture Histories From Multi‐Array Back‐Projection

Teleseismic back‐projection imaging has emerged as a powerful tool for understanding the rupture propagation of large earthquakes. However, its application often suffers from artifacts related to the receiver array geometry. We developed a teleseismic back‐projection technique that can accommodate data from multiple arrays. Combined processing of P and pP waveforms may further improve the resolution. The method is suitable for defining arrays ad‐hoc to achieve a good azimuthal distribution for most earthquakes. We present a catalog of short‐period rupture histories (0.5–2.0 Hz) for all earthquakes from 2010 to 2022 with MW ≥ 7.5 and depth less than 200 km (56 events). The method provides automatic estimates of rupture length, directivity, speed, and aspect ratio, a proxy for rupture complexity. We obtained short‐period rupture length scaling relations that are in good agreement with previously published relations based on estimates of total slip. Rupture speeds were consistently in the sub‐Rayleigh regime for thrust and normal earthquakes, whereas a tenth of strike‐slip events propagated at supershear speeds. Many rupture histories exhibited complex behaviors, for example, rupture on conjugate faults, bilateral propagation, and dynamic triggering by a P wave. For megathrust earthquakes, ruptures encircling asperities were frequently observed, with downdip, updip, and balanced patterns. Although there is a preference for short‐period emissions to emanate from central and downdip parts of the megathrust, emissions updip of the main asperity are more frequent than suggested by earlier results.


Introduction
Following the extraordinary 26 December 2004 Sumatra-Andaman (M W 9.3) earthquake, back-projection rupture imaging has become an important technique to complement finite-fault source inversions for determining the kinematics of very large ruptures (e.g., Ishii et al., 2005;Krüger & Ohrnberger, 2005;Walker et al., 2005).Backprojection exploits the coherence of (usually) teleseismic P waveforms with limited prior assumptions on the fault geometry.Applications have targeted, for example, megathrust subduction earthquakes (e.g., Koper et al., 2011;Lay et al., 2010;Meng et al., 2011;Palo et al., 2014), aftershock detection (e.g., Feng et al., 2020;Kiser & Ishii, 2013;Tilmann et al., 2016), complex ruptures along multiple faults (e.g., Lay et al., 2018;Meng, Ampuero, Stock, et al., 2012;Ruppert et al., 2018), intermediate-depth earthquakes (e.g., Kiser et al., 2011;Ye et al., 2020), moderate size earthquakes (e.g., D'Amico et al., 2010;Taymaz et al., 2021), and near-field tsunami prediction (e. g., An & Meng, 2016;Xie & Meng, 2020), among others.Back-projection is generally applied to earthquakes that are large enough to reveal source complexity in teleseismic recordings.The technique makes use of the coherency of short-period filtered waveforms, which are able to resolve details of the earthquake rupture better than long-period waveforms (e.g., Kiser & Ishii, 2017).The frequency band used must subsequently be taken into account in the interpretation of the inferred rupture history.Short-period radiation is related to variations of slip and rupture velocity (e.g., Madariaga, 1977;Madariaga, 1983;Marty et al., 2019).This implies that the back-projected rupture history can be related to fault heterogeneities.For example, subduction zone megathrust earthquakes were found to radiate short-period energy predominantly in deeper parts of the megathrust, downdip of the maximum slip regions derived from finite-fault inversions, which are sensitive to the low-frequency emissions of the rupture (e.g., Lay et al., 2012;Yao et al., 2011Yao et al., , 2013)).The apparently dominant downdip short-period radiation is proposed to occur at the transition between brittle and ductile regions (e.g., Lay et al., 2012;Simons et al., 2011).Recently, K. Wang et al. (2020) correlated the 27 February 2010 Maule (M W 8.8) earthquake short-period rupture with downdip segmentation at the base of the overriding mantle wedge.In contrast to this view of the along-dip segmentation of seismic emission during large earthquakes, Meng et al. (2018) provided the first evidence for an encircling rupture, that is, high-frequency rupture emissions both updip and downdip of the slip patch for the 16 September 2015 Illapel (M W 8.3) earthquake.Similarly, Vera et al. (2023) presented rupture emissions outlining the main slip asperities (updip and downdip) for the 14 November 2007 Tocopilla (M W 7.7) earthquake.Moreover, the observations of Meng et al. (2018) for the 2015 Illapel event are in contrast to earlier back-projection studies, which had shown the expected, predominantly downdip rupture pattern (e.g., Melgar et al., 2016;Tilmann et al., 2016;Yin et al., 2016).It raises the question if this type of complex behavior is actually more common and just requires higher resolution rupture images.
Fracture mechanics predicts that instabilities such as earthquakes either propagate at sub-Rayleigh speed or exceed the shear wave speed (e.g., Andrews, 1976;Burridge, 1973;Das & Aki, 1977), where the latter is only expected for mode II cracks.In mode II cracks, the rupture propagates in the direction of displacement, which is also the direction of the initial shear stress resolved onto the fault.Such supershear earthquakes are associated with simple fault geometry and homogeneous stress-strength conditions (e.g., Bouchon et al., 2010), but also to damage zones under relatively low stress leading to unstable supershear ruptures, that is, between the shear wave speed and
A disadvantage of the back-projection method is that the array configuration can cause notable artifacts in the recovered rupture (e.g., Meng, Ampuero, Luo, et al., 2012).Back-projection often leads to a persistent time-space trade-off of the earthquake rupture toward the seismic array, often called "swimming" artifact.This effect arises from the low curvature of the time-distance travel time curve.Sources closer to the receivers but activated later will have the same arrival time as earlier sources farther away, resulting in a point source in space and time appearing as an extended source drifting toward the array.In animations showing the evolution of the backprojected energy with time, this slightly irregular drifting looks like a swimming motion, giving the artifact its name.Because of the dependence on the array azimuth, it is easy to understand why the swimming artifact can be reduced by combining multiple array images.
Depth phases can cause additional artifacts in the form of "ghost" emitters corresponding approximately to the bounce points of the surface-reflected phase, but they also can carry additional information.For large intermediate-depth earthquakes, the time-delay between P and depth phases (e.g., pP and sP) allowed to improve the resolution in depth by combining both P and depth phase backprojections (e.g., Kiser et al., 2011).For more shallow earthquakes (40-100 km), depth phases can contribute significantly to uncertainties (e.g., Zeng et al., 2020).To our knowledge, however, a systematic imaging method that integrates depth phases and multiple arrays for the shallow depth range has not been reported yet.
In this study, we extended the multi-array approach of Rössler et al. (2010) by including depth phases depending on the radiation pattern (for earthquakes deeper than 40 km) and weighted seismic array images.We provide an algorithm for automatically estimating rupture length, directivity, speed, and aspect ratio from backprojections.Short-period rupture lengths are additionally used to calculate scaling relations.This modified method is applied to derive a catalog of rupture histories of recent large earthquakes (01/2010-12/ 2022) in the 0.5-2.0Hz frequency range, which is complete for M W ≥ 7.5 and depths less than 200 km.The analysis focused on complex ruptures, depth-varying short-period radiation for large subduction earthquakes, and the detection of supershear ruptures for strike-slip events.We also showed that short-period ruptures encircling asperities are frequent in subduction megathrust earthquakes.The results suggest that encircling rupture patterns around slip patches are primarily related to the high-stress gradient around the asperity (and seismogenic barriers).

Multi-Array Multi-Phase Back-Projection
The back-projection method is similar to beamforming in maximizing the coherency of time-shifted waveforms at an array.Unlike in beamforming, there is no assumption of a planar wavefield.Instead, the time shifts are calculated from the predicted travel times for a grid spaced around the hypocenter.In practice, the grid is usually two-dimensional, chosen to be either a horizontal plane at the hypocentral depth, or a plane aligned with one of the nodal planes of the focal mechanism, or an a priori known fault surface, for example, the slab interface in subduction zones.The waveforms are back-projected onto this grid.The theoretical arrival of a target seismic wave (e.g., P wave) based on a reference velocity model controls the beamforming delays.For each grid point, the resulting array beam is: where N is the number of stations, b i (t) is the beam for the ith grid point, u k the vertical component waveform recorded at station k, t ik the travel time between the grid-point i and station k in a reference velocity model, and Δt k the station-specific static correction term accounting for differences between the reference and true velocity model.
The station correction terms can, to a large extent, absorb the effect of 3D Earth heterogeneities on arrival times.They are usually determined by cross-correlating the first few seconds of the rupture recorded by each receiver.The resulting time-shifted arrivals are then compared to those predicted for the catalog hypocenter, and the differences correspond to the necessary correction terms.Thus, the back-projection image retrieves the rupture nucleation at the catalog hypocenter by definition.For very large earthquakes, aftershocks can alternatively be used to correct source-receiver paths away from the hypocenter (e.g., Ishii et al., 2007;Meng et al., 2016;Palo et al., 2014).The advantage over the hypocenter-based calibration is that the location errors of several events are averaged, such that a possible bias from mislocation of the mainshock hypocenter is reduced, and furthermore, spatially varying station terms due to 3D structure effects can be accommodated.The calibration with aftershocks is particularly important when the rupture pattern is compared with aftershocks.In contrast, the mainshock calibration offers advantages in near-real-time applications since only the hypocenter is required.Alternatively, travel time calculations could be performed in a 3D Earth model (e.g., Liu et al., 2017).For simplicity, in this work, we adopt the hypocenter calibration method.
In order to carry out the actual rupture tracking, the maxima of beamformed energy (E i ) and semblance (S i ), defined in Equations 2 and 3 below, are used to locate the most intense emission at each time step.The energy represents the amount of radiation emitted, and semblance provides a measure of the coherence of waveforms which is not affected by the amplitudes of individual traces and is, therefore, more effective for tracking the location of the earthquake rupture (e.g., Neidell & Taner, 1971;Palo et al., 2014;Rössler et al., 2010).For both measures, a time window of length W needs to be defined, which should accommodate at least two periods of the longest period analyzed.
Both measures provide an image of the earthquake rupture.Given a grid of sources, tracking the local semblance maxima provides a proxy of the rupture propagation.The energy peak in each time window provides a relative measure of the source time function of short-period seismic energy, which is related to but not necessarily proportional to the moment rate, for example, as derived from finite fault solutions.Thus, the time-integrated energy maps provide a summary view of the high-frequency energy radiation.
For the processing of data from multiple arrays, we initially follow the approach of Rössler et al. (2010) where semblance and energy maps from N a arrays are multiplied to provide the rupture and energy radiated, respectively: Here, we introduce two modifications: (a) we use exponents (equivalent to weights in log-space) to balance the contributions of individual arrays (based on their azimuthal distribution to avoid artifacts due to clustering of arrays in certain azimuthal ranges), and (b) for earthquakes deeper than 40 km, we combine P and pP backprojections to reduce artifacts related to the depth phase (see Figure 1a).The modified expressions are: We assign p = 1 to P and p = 2 to pP waveforms, and N P is the number of phases used.Here N P = 1 is used for shallow earthquakes and 2 for deep earthquakes, but the method is open to experimenting with alternate seismic wave arrivals (e.g.,: sP, PP, and PKIKP; e.g., Kiser et al., 2011;Koper et al., 2012).The weighting exponent γ is set proportionally to the sum of the two half-angles between the azimuths of target and neighboring arrays (see Figure 1b), where the median of the array receiver coordinates is assumed to be the reference location for the weighting, and the normalization is chosen such that ∑ N a i γ i = 1.Therefore, for a single array γ = 1, and for two, we always have γ = 0.5, but for irregularly shaped arrays, the weighting depends on the distribution.This weighting ensures internal consistency in the sense that for a perfectly coherent signal, a large array could be split into two sub-arrays, and the resulting peak energy values would be the same as would have resulted from a single array.The same principle can be extended to combine backprojections from arrays symmetrically well-located around the source.However, the resulting source image might be biased due to a non-symmetric configuration, that is, clustering of arrays in a specific azimuthal range.Here, a weighting scheme based on the azimuthal coverage allows for balancing the contribution of each array, preventing potential artifacts from affecting the solution.The term f(t t 0 ), with t 0 the origin time of the earthquake, is a taper function adapted from Kiser et al. (2011) to mute the waveform before the first arrival when multiple seismic phases are included:

Journal of Geophysical Research: Solid Earth
If τ = 0, the rupture has initiated, and the results are fully incorporated to image the source; otherwise, they are suppressed.T is the period of the cosine taper function controlling the transition between admitted and suppressed intervals.

Processing Details
Fifty-six large earthquakes in the time range 01/2010-12/2022 (all events with M W ≥ 7.5, depth ≤200 km according to Global CMT catalog, https://www.globalcmt.org)were back-projected.Figure 2 shows the distribution of the earthquakes color-coded by faulting type.Thrust earthquakes dominate the catalog (57%), followed by strike-slip (27%) and normal (16%) faulting mechanisms.Tables S1-S3 in Supporting Information S1 present the earthquake source parameters for strike-slip, thrust, and normal faulting earthquakes, respectively.
The number of (ad-hoc) arrays weighted and combined in the analysis depended on their availability at epicentral distances between 30 and 100°, where the P waveforms are not subject to complications like triplications as observed at smaller distances.The array selection and weighting prioritize azimuthal distribution while incorporating as many arrays as possible.Dense arrays were formed depending on the configuration and geometry of permanent broadband networks, for example, North America (e.g., US: United States National Seismic Network; AK: Alaska Regional Network), Japan (NIED Hi-net Network), Europe (many national and global networks distributed by ORFEUS-EIDA; Strollo et al., 2021) and Africa (e.g., AF: AfricaArray).Additionally, local and regional networks (often temporary deployments) were evaluated (e.g., data availability, waveform coherence) to form smaller aperture arrays and maximize the azimuthal coverage.Therefore, the extent to which a good azimuthal distribution can be achieved differs between regions and also depends on the time and time gap after the event (as temporary network data are often only openly available a few years after the experiment).Earthquakes in Indonesia were favorably located to be within range of several arrays, for example, networks in Asia, Europe, Africa, Australia, and Antarctica.In contrast, earthquakes in northern and central Chile suffered from limited coverage, although a somewhat reasonable azimuthal distribution can be achieved by combining networks in North America, Africa, and Antarctica.For events in southernmost Chile, Australian networks improved the coverage.
The data pre-processing involves mainly the removal of the instrument response from the waveforms and basic quality assessment.The mainshock location, origin time, and theoretical P and pP wave arrival times controlled the rupture imaging, allowing a reference frame to merge multiple array backprojections.P and pP wave arrival times were predicted based on the IASP91 velocity model (Kennett & Engdahl, 1991).For P waves, static corrections were determined by measuring the relative time shifts of first arrivals with the adaptive stacking method of Rawlinson and Kennett (2004).We considered the first 15 s after the P-wave onset on bandpass-filtered (0.4-3.0 Hz) vertical velocity waveforms.This frequency band is a little wider than the band used for the back-projection (0.5-2.0 Hz) and optimizes coherence while enough high frequencies are retained for a precise alignment.For pP waves, theoretical arrival times were calibrated with the P wave corrections.We also removed anomalous traces that could impact the waveform stack during the adaptive stacking as part of the input quality control.
We combined P and pP backprojections for earthquakes deeper than 40 km and only if the theoretical pP wave amplitude reached a minimum of 40% of the direct P wave amplitude; see Section 3.2 for synthetic tests and extended discussion on the radiation pattern.This provides a minimum separation time between P and pP wave arrivals of ∼11 s (depth ≥40 km; arrays at epicentral distances of 30-100°) in the IASP91 Model (Kennett & Engdahl, 1991).The ≥40 km depth threshold is additionally aimed to mitigate the impact of depth phases in the back-projection, as observed for intermediate-depth (40-100 km) earthquakes (e.g., Zeng et al., 2020).To calculate the pP/P amplitude ratio, we considered the radiation pattern for P and pP waves (based on the respective Global CMT focal mechanism), including the P-P wave reflection coefficient of the free surface derived from the Aki and Richards approximation (Aki & Richards, 2002).The P wave slowness was calculated with the TauP Toolkit (Crotwell et al., 1999).To determine the take-off angles at the source and calculate the surface reflection Earthquakes back-projected with the combined P and pP approach are highlighted in magenta color.Labels designate the earthquake ID, with details given in Figure 4 and Tables S1-S3 in Supporting Information S1.Red segments show major plate boundaries (Bird, 2003).
coefficient, we used averages from upper, middle, and lower crustal layers of the CRUST2.0model (Vp = 6.4 km/ s, vs. = 3.5 km/s, and Density = 2.9 g/cm 3 ; Laske et al., 2013).We considered the peak amplitudes evaluated over the array to calculate the pP/P amplitude ratio.Figure S9 in Supporting Information S1 shows the P-P wave reflection coefficient as a function of the incidence angle.Table S7 in Supporting Information S1 presents the pP/ P amplitude ratios for the events and arrays included in the combined P and pP back-projection.
For strike-slip, normal, and thrust (non-megathrust) earthquakes, the target grid was placed on a horizontal plane at a reasonable depth to minimize artifacts, that is, the hypocentral depth, with grid points every 5 km.For megathrust earthquakes, the grid followed depth variations of the SLAB1.0model (Hayes et al., 2012).We remind that depth is not resolved explicitly, but the horizontal position is translated to depth to obtain accurate theoretical travel time predictions.For bilateral ruptures of a large extent (>400 km), we additionally mapped semblance maxima over a pre-defined sub-region of the grid to probe secondary rupture patterns, that is, the 2010 Maule (M W 8.8), 2011 Tohoku-Oki (M W 9.1), and 2017 Komandorsky Islands (M W 7.8) earthquakes.We also added the 2015 Illapel (M W 8.3) earthquake and defined sub-regions along-dip to simultaneously track updip and downdip rupture emissions as reported by previous back-projection studies, that is, Meng et al. (2018).
In our back-projection analysis, we considered a moving window of length W = 6 s, chosen to be three times the dominant period (back-projection frequency band: 0.5-2.0Hz) and moving forward in 1 s increments.Longer values for W would have over-smoothed the rupture image, and shorter ones result in instability.The time reference is chosen such that the time axis refers to the beginning of the sliding window used in the source region with zero sets to the nominal origin time of the event.This choice of reference means a coherent arrival will become visible at nominally negative times.Prior to back-projection, we applied a 0.5-2.0Hz bandpass filter, corresponding to the highest frequency range for which consistently sufficient waveform coherency has been obtained in previous rupture imaging studies (e.g., Meng et al., 2015Meng et al., , 2016Meng et al., , 2018;;Palo et al., 2014).
Finally, the end of the rupture was visually determined, but an automatic estimation was included for comparison.Reactivation of earlier fault segments or scattered semblance maxima were considered to indicate the end of the rupture.The apparent source time function of radiated energy was also considered in the visual inspection, as small values compared to its peak indicate the end of the rupture.For the automatic rupture end determination, we followed the approach of Xie and Meng (2020), which also relies on detecting scattered maxima.The earthquake rupture is declared to have ended when the standard deviation of the rupture emission points within a 30 s window exceeds a threshold of 50 km.Additionally, we defined three criteria to avoid under-and overestimated rupture end times.First, we included lower and upper boundaries to restrict the search window after the earthquake initiation relative to magnitude and faulting type.Specifically, we assume a linear rupture propagation and estimate the rupture length following the scaling relations of Blaser et al. (2010).We, therefore, convert the upper and lower bounds of rupture speed (6.0 and 1.5 km/s, respectively) into upper and lower bounds for rupture durations.Second, if the threshold in the variance can not be passed, it will be reduced by 1 km until a peak value can be reached.This could be necessary for events with compact rupture lengths that would not allow passing a significant threshold in the spatial variance.Third, complex events, for example, bilateral ruptures or complex fault systems, might present high variance in the location of rupture emissions, which would wrongly indicate the end of the rupture.After inspecting the processed events, we ignore early variance peaks up to a specified time delay, specific to each event.

Estimation of Basic Source Parameters
We provide automated estimates of rupture length, directivity, and speed based on the rupture propagation derived from semblance maxima.Furthermore, we calculate aspect ratios to provide a proxy of rupture complexity.We emphasize that these parameters are only reflective of the high-frequency behavior of the earthquake and cannot always be assumed to coincide with the long-period properties based on total slip.Nevertheless, at least for length and directivity estimates, we expect a reasonable correlation (see discussion in Section 3.5.3).Clearly, if spurious scattered semblance maxima are included, errors in these estimates can become quite large, and it is, therefore, important to carefully determine the end of the rupture as accurately as possible.In our analysis, we used the manually determined rupture end times.Results for automatically determined rupture terminations are presented in Supporting Information S1.
To quantify rupture length (see Figure 1c), candidate lengths L are estimated by the projection of the semblance maxima (blue-red circles in Figure 1c) on lines passing through the epicenter (white star) for all azimuths 0-180°J (magenta line and open black circles show the realization for a 70°azimuth).The maximum value of L over all azimuths is chosen as the rupture length (red line in the bottom left subplot in Figure 1c).
The aspect ratio measures the sparsity of the short-period rupture emissions and is defined by the quotient between the minimum and maximum length estimates obtained from the rupture length grid search (aspect ratio close to 1: significant sparsity; aspect ratio close to 0: very low sparsity).We emphasize that the aspect ratio does not attempt to describe the earthquake dimension or rupture area, which is not well represented by high-frequency emissions, but the complexity, especially for larger earthquakes (M ≥ ∼8).Specifically, a large aspect ratio indicates a high degree of sparsity, for example, a conjugate fault rupture as observed for some oceanic intra-plate earthquakes (e.g., 2012 Wharton Basin and 2018 Gulf of Alaska) or rupture emissions outlining a simple megathrust earthquake asperity, that is, encircling rupture patterns.Subduction earthquakes tend to have shallowly dipping rupture planes, such that their surface projection represents approximately the fault plane.Here, the aspect ratio represents the actual aspect ratio of the high-frequency rupture and is sensitive to whether emissions occur, for example, downdip, updip, or encircling the main slip patch.Due to their much shorter rupture lengths and thus higher relative sensitivity to scattered emissions, the aspect ratio of events smaller than M ∼8 can not usually be reliably estimated, but for consistency, it is reported for all analyzed earthquakes.
Directivity corresponds to the average direction in which the rupture fronts propagate.It is measured similarly to rupture length, with lines of all azimuths 0-180°pivoting through the epicenter, but the quantity minimized here is the sum of the squares of the perpendicular distances of all semblance maxima to the line (i.e., average squared lengths of blue arrows in Figure 1c, which are shown symbolically for a few semblance peaks only).For simple ruptures, the azimuths returned by the length and directivity measurements will be very similar, but the length estimate is controlled by the end points, whereas the directivity estimate is controlled by all points simultaneously.The ambiguity in the actual directivity is resolved by considering an imaginary line perpendicular to the rupture direction and centered on the epicenter; we then choose the directivity based on which side of this line more rupture emission points are found.
To determine the rupture speed, we first calculate estimates of instantaneous rupture velocities by dividing the distance between the epicenter and a subsequent semblance peak by the time elapsed since the origin time, which really represents the average rupture velocity from nucleation to the current time step.The resulting time series is strongly affected by location uncertainties for small distances and times but stabilizes quickly (see Figure 1c).The time series of rupture velocities is smoothed within a 4 s wide time window, and its peak value is then taken as the event rupture speed; maxima during the first unstable 10 s are ignored.Although strictly speaking, this estimate physically represents the maximum average velocity (where the average is taken over all preceding times) and represents a reasonable estimate for most ruptures.Based on theoretical considerations, for supershear ruptures on simple faults, supershear velocities are thought to be attained quickly and maintained until the end of the rupture (Das & Aki, 1977), such that they can be detected straightforwardly by this estimation procedure.Nevertheless, the estimated rupture velocities can be significantly lower than the peak rupture speeds, particularly for complex ruptures with directional changes or those with a slow start.To mitigate the effect of a slow start caused by a small triggering nucleation event, for example, the 17 December 2016 Solomon Islands (M W 7.9) earthquake, instead of referring back to the epicenter for the whole rupture duration, we reset the reference point to an emission point at the beginning of the main phase of the rupture.
We also calculated average rupture velocities for comparison with our preferred maximum values.While average velocities can serve as a reference for assessing the typical estimated velocity in each earthquake rupture, the use of maximum velocity estimates enables the evaluation of rupture velocity transitions, such as those from subshear to supershear rupture.These transitions may be suppressed or underestimated when relying on averages.

Case Study: The 2020 Kuril Islands Earthquake
We introduce the presentation of our results based on the 25 March 2020 East of Kuril Islands (M W 7.5) earthquake as an example (Figure 3).Results for the complete catalog are presented in Figures S12-S67 of Supporting Information S1, including a summary of derived parameters in Figure 4 and Table S4 in Supporting Information S1.We selected the Kuril Islands event because of the relatively simple rupture and because, at the time of writing, to our knowledge, no other back-projection analysis had been published yet for this event.This event is an intraplate thrust earthquake, which probably was triggered by compressional bending stresses in the deep interior of the subducting Pacific plate (Ye et al., 2021).
The left inset of Figure 3a shows the distribution of arrays, where the arrays in North America and Australia are weighted more strongly as they cover a larger backazimuthal range.In Figure 3b, the waveform coherency near the rupture initiation is visualized.It is generally good except for the pP phase at the China array.The main plot compares the back-projected rupture with the USGS-NEIC finite fault slip solution (https://earthquake.usgs.gov/earthquakes/eventpage/us70008fi4/finite-fault); for other events, frequently finite slip models from the literature are shown instead.Because this event is not a subduction megathrust earthquake, we back-projected onto a horizontal plane at 58 km depth, the hypocentral depth of the event.
For the East of Kuril Islands earthquake, the finite fault solution of Ye et al. (2021) favored rupture on the northwest-dipping fault plane.A further advantage of a horizontal plane is that this assumption will not fail in case of more than one planar fault being activated (e.g., the 2018 Gulf of Alaska earthquake; Figure S52 in Supporting Information S1).For the Kuril event, the main slip patch appeared close to the epicenter, and at the beginning of the rupture, both finite slip and short-period rupture spatially agreed.After 15 s, the high-frequency rupture propagated to the east of the secondary slip patch.This time period revealed the most energetic short-period emission.Close to this area, the most intense aftershock activity was observed (see Figure 1   The rupture propagated unilaterally along a nearly straight line, making the length and directivity estimation straightforward (Figure 3c).The implied strike direction of the rupture track is clockwise rotated by some 10-20°w ith respect to the strike direction indicated by the moment tensor.At face value, this would imply that the shortperiod emission points did not strictly propagate along strike but moved additionally updip (or downdip, depending on which fault plane is the correct one).The automatically determined rupture duration was 2 s longer than the visually determined end time (see Figure S1 in Supporting Information S1 for automatic end times).Event ID (for cross-referencing to other figures, tables, and Supporting Information S1), moment magnitude, rupture length (colored by rupture time), aspect ratio, directivity, and rupture speed (colored by hypocentral depth).The events are categorized into normal (upper), thrust (middle), and strike-slip (bottom) earthquakes, with the order in each category determined by magnitude.Events in bold text are assumed to be subduction megathrust earthquakes.Events labeled with a black bullet indicate that a secondary rupture pattern was obtained over a predefined grid and included in the rupture parameter estimates.The vertical red bars in the rupture speed panel show the expected shear wave speed (V S ) at the hypocentral depth in the IASP91 velocity model (Kennett & Engdahl, 1991).Rupture speeds in the range 95-105%V S and >105%V S define "likely' and "strongly" supershear earthquakes, respectively.
Finally, we can compare the time history of emitted short-period energy with the moment rate obtained from the SCARDEC and USGS databases (highlighted in Figure 5).The functions show two peaks, and the timing of the first peak and subsequent trough agree.However, at short periods, the second energy peak exceeds the first one, unlike for the USGS and SCARDEC solutions, and the rupture appeared to continue for longer.Interestingly, in the source time function estimated by Ye et al. (2021) (their Figure 1), the moment rate of the second peak also exceeded the first one.

Synthetic Tests
To evaluate the lateral resolution, we considered the array configuration used in the Kuril Islands earthquake back-projection (see Figure 3b) and simulated a unilateral rupture scenario.Synthetic seismograms were generated using a Ricker wavelet with a central frequency of 1 Hz and 50% of Gaussian random noise (e.g., Kiser et al., 2011;Ricker, 1953).The maximum amplitude of the signature defined the standard deviation of the Gaussian noise.Five sources placed every 15 km at the hypocentral depth of the Kuril event were used to simulate  the rupture propagation.Although using only five separated sub-sources does not lead to realistic-looking seismograms, this scenario allows us to visualize the accuracy of the recovered rupture easily.
To observe the effect of merging arrays, we assumed P and pP waves with the same amplitude and simulated a rupture propagating to the northeast at 2.5 km/s.Figure 6 presents tests using only P (upper subplots) and combined P and pP (lower subplots) backprojections with a variable number of contributing arrays.The resolution improves notably with the number and wide azimuthal distribution of combined arrays.Furthermore, the combined P and pP back-projection can better resolve the rupture compared to the P-only back-projection, also resulting in more accurate rupture length and directivity estimates.Rupture speeds are often mildly overestimated with variations up to 22% with respect to the true speed (pP⋅P back-projection; Figure 6).We remind that the rupture speed corresponds to the maximum average velocity evaluated over the rupture time.Overestimated rupture speeds, that is, the maximum average speed, may arise or be affected by artifacts or smearing effects between the activation of the source points and the initial unstable first few seconds of the rupture tracking.We preferred maximum over average speeds to quantify the uppermost speed limit from short-period back-projection.This also allows us to identify supershear ruptures that could be underestimated with average values while providing a standard frame to evaluate the earthquake database.It should also be noted that in the tests of Figure 6 and Figures S4-S6 in Supporting Information S1 (described below), P and pP Ricker wavelets have the same amplitude.The variations in the amplitude of the waveforms (Figure 6b) correspond to interference between P and pP wavelets, which arises from the set of emission points used for simulating the rupture propagation.
We tested other scenarios with rupture propagating to the southwest (Figure S4 in Supporting Information S1) or bilaterally (Figure S5 in Supporting Information S1) and observed similar improvements.For the bilaterally propagating rupture, the addition of the AU and CN arrays continues to improve the estimate, while for the northeast propagating rupture, only marginal improvements were observed from adding the third and fourth arrays.For the northeast-propagating rupture scenario, we also tested the pP back-projection alone to assess the distorting effect of P waves from the later part of the rupture on the pP back-projection image (Figure S6 in Supporting Information S1).Tracking the rupture initiation is challenging for the pP back-projection but overall the result is comparable with the P back-projection by merging multiple arrays.
We also conducted synthetic experiments with different pP/P amplitude ratios, which in reality depend on the radiation pattern.We tested the northeast-propagating rupture scenario shown in Figure 6.If the pP wave amplitude is at least as large as that of the P wave, the combined P and pP approach performs better than only the P back-projection with regard to the accuracy of rupture length and directivity estimates and the tracking of emission points (see Figure 7).If the pP amplitude drops to 40% of the P wave amplitude, the resolution decreases, but consistency is still comparable to the input.When the pP amplitude ratio drops to 30% of the P wave amplitude, the combined P and pP back-projection fails to separate emission points properly, and the results of the P-only back-projection are preferable.
We also tested slower and faster rupture speed scenarios of 1.5 and 4.0 km/s.If the rupture propagates at 4.0 km/s, the advantages over the P back-projection remained, but both P and combined pP and P backprojections failed to recover the beginning and end of the rupture, thus underestimating its length (Figure S7 in Supporting Information S1; the input rupture length is underestimated by up to 26%).The variations in the recovered speed are 17% of the input speed.The back-projection for a slow rupture scenario of 1.5 km/s is more prone to artifacts (Figure S8 in Supporting Information S1).Scattered emissions points are frequent, and rupture speeds are severely overestimated (up to 71%).
Following the tests, as previously mentioned in the description of the processing workflow (Section 2.2), we exclusively applied the combined P and pP approach for events deeper than 40 km depth and only for those arrays where the predicted pP wave amplitude is at least 40% of the P wave; otherwise, the P-only back-projection is used.Additionally, the tests suggest resolving emission points ∼15 km apart (horizontally) should be possible in most cases.Furthermore, we tested the back-projection with a selection of aftershocks to provide a more realistic rupture simulation in terms of noise and complexity of the waveforms.We considered a sequence of eight aftershocks that approximately covered the 2011 Tohoku-Oki (M W 9.1) earthquake extent.Please refer to Table S8 in Supporting Information S1 for details on the aftershock source parameters and Figures S68-S69 in Supporting Information S1 for their visualization (aftershocks marked with stars, outlined in cyan color).Before stacking, we pre-filtered the data within the 0.1-5.0Hz frequency range for display purposes, which is a wider band than the one later implemented in the back-projection (0.5-2.0 Hz).We also normalized the waveforms to prevent an image dominated by the largest aftershocks, which could be larger than realistically expected for a single rupture.Thus, normalization allows for a balanced contribution of individual aftershocks and ensures that the back-projection tests focus on rupture tracking.Due to the significant distance between the emission points, ranging from 38 to 85 km, we conducted tests with propagation speeds of 2.5 and 4.0 km/s, that is, unilateral and bilateral rupture scenarios.This provides gaps up to ∼90 s in the rupture propagation (2.5 km/s: gaps between 21 and 93 s; 4.0 km/ s: 13-58 s), allowing to assess the performance of the back-projection when dealing with interruptions in rupture emissions, that is, potential smearing artifacts.Following the processing of rupture synthetics, we performed P and combined P and pP backprojections in the 0.5-2.0Hz frequency range.We conducted P⋅pP backprojections without considering the threshold depth of ≥40 km and the radiation pattern (aftershock depth ranging from 19 to 53 km; Table S8 in Supporting Information S1).This evaluation aimed to assess the performance of the method in challenging rupture tracking scenarios, complementing the previous Ricker wavelet analyses while keeping the same aftershocks and station coverage used for the P back-projection tests.Because of the temporal span of the  S8 in Supporting Information S1) and the teleseismic waveform coherence for moderate magnitude events (M6.1-7.1;Table S8 in Supporting Information S1), we obtained only three arrays, well distributed in azimuth, that recorded all of the aftershocks, that is, EU (58-86 stations), USA (41-89 stations) and IND (4 stations).See Figures S70-S78 in Supporting Information S1 for the waveforms and array distribution used in the tests.Additionally, Figures S71-S73 in Supporting Information S1 presents an example of the stepwise inclusion of aftershocks to generate the waveforms shown in Figure S70 of Supporting Information S1 and back-projected in the tests of Figures S68a-S69a in Supporting Information S1.This illustration highlights the impact of the coda waves in generating complex rupture synthetics through a sequence of aftershocks.The limited number of stations and the complexity of the waveforms provided an additional challenging input to be back-projected and recovered.
Figure S68 in Supporting Information S1 shows the P wave backprojections for the selection of aftershocks, that is, southwest (left column), northeast (central column), and bilateral (right column) propagating ruptures at 2.5 km/s (upper subplots: a, b, and c) and 4.0 km/s (lower subplots: d, e, and f).The array coverage and waveform coherence for each scenario are presented in Figures S70-S75 of Supporting Information S1.The difference in the recovered rupture length and directivity relative to the input was relatively small, with an absolute variation of up to 18% for rupture length and 4°for directivity (rupture parameter estimates inset in Figure S68 of Supporting Information S1).The recovered rupture speed is overestimated by up to 31% compared to the input speed (as seen in the time series inset in Figure S68 of Supporting Information S1).However, it is important to consider the effect of artifacts when interpreting peaks in rupture speed, specifically, the smearing of emissions when the rupture stops propagating between the source points.Overall, the rupture propagation can be recovered in both space and time, even in more complex bilateral rupture scenarios (Figures S68c and S68f in Supporting Information S1; see the comparison between the back-projected rupture and the input sources, color-coded relative to activation time).The results for combined P and pP backprojections are presented in Figure S69 of Supporting Information S1.Compared to the P back-projection, the main difference in the rupture parameters is an overestimation of speeds (up to 52% of the input speed).However, it is worth noting that most of the rupture speed time series are below the input speed (see magenta dashed line in the speed time series; Figures S68 and S69 in Supporting Information S1).Potential smearing effects due to gaps in the emission points and the observed tracking instability during the first few seconds of the rupture should be considered when interpreting rupture speeds and overall rupture parameter estimates.
The Ricker wavelet and aftershock sequence tests showed the reliability of combining different array and P and pP wave backprojections (depth ≥40 km; pP wave amplitude ≥40% of the P wave) for tracking the earthquake rupture.This is observed even in complex bilateral rupture scenarios and when arrays of limited aperture complement the processing.

Short-Period Source Time Functions
The short-period energy time functions for all events are compared with USGS and SCARDEC (Vallée & Doeut, 2016) moment rates in Figure 5.The time history of radiated high-frequency seismic energy frequently agrees (e.g., 2021 Chignik, ID 43;2020 Caribbean, ID 8;2016 Kaikoura, ID 32;2013 Scotia Sea, ID 11;2012 Wharton Basin, ID 15;2010 Mentawai, ID 33) with the seismic moment rate functions, but it is sometimes apparently time-shifted, following (e.g., 2016 Melinka, ID 23;2013 Craig, ID 1;2011Tohoku-Oki, ID 47) or, less often, preceding (e.g., 2016 Solomon, ID 35;2012 Philippine Islands, ID 25) the seismic moment rate.A few earthquakes presented quite different shapes (e.g., 2020 Kuril Island, ID 16).However, it has to be noted that the USGS and SCARDEC moment rates also do not always agree.One striking example is the Iquique-Pisagua (M W 8.1) earthquake (ID 41; Figure 5), where the SCARDEC estimate suggests a much shorter source time function.For the Chiapas (M W 8.2) earthquake (ID 56; Figure 5), SCARDEC (and the short-period energy function) completely lacks the secondary peak seen in the USGS estimate.Because of this variability, it is difficult to judge whether the differences between moment rate and short-period energy function are physical (i.e., related to the time-varying ratio of slip to high energy emissions), methodological artifacts, or a combination of both.A possible physical explanation is that the ratio of short-period to long-period spectral energy is known to vary significantly between earthquakes (or equivalently, the seismic energy to moment relation, which is related to the variability in stress drop).It seems likely that this variability can also extend to the case of different asperities within a given rupture, which would then cause different time histories of moment rate and short-period energy function.2020), (j) Heidarzadeh et al. (2017), and (l) Moreno et al. (2018).Bottom inset: High-frequency rupture pattern classification.The dominant type is noted in the upper of each plot.

Megathrust Rupture Patterns
We now consider the rupture patterns of subduction megathrust earthquakes, specifically the spatial relationship between the short-period emissions and the finite slip models characterizing the long-period behavior.Previous observations indicated a strong preference for short-period emissions downdip of the main slip asperities (e.g., Lay et al., 2012).Figure 8 presents the rupture patterns for the largest subduction megathrust events in our catalog.
We can distinguish a variety of earthquake rupture patterns referred to as downdip, updip, and balanced (segmented or encircling) determined by the relation between short-period radiation and the asperity (see Figure 8 for a graphical definition; the same events are presented in Figure S10 of Supporting Information S1 color-coded by rupture time).We note that "encircling rupture patterns" are overall a purely kinematic description to highlight high-frequency emissions around megathrust asperities and does not necessarily mean a bifurcation of the rupture but could be, for example, the result of an extending rupture crack.Here, we define, somewhat arbitrarily, the asperity as the area enclosed by the contour corresponding to 75% of the maximum slip and classify each emission point whether it is downdip, updip or within the main asperity.For the downdip and updip classification, the emissions have to be outside the asperity contour and downdip or updip of a transect through the peak slip in the direction of the fault strike (see schematic overview in the lower-right part of Figure 8; the yellow-red color represents the slip patch, the encircling black contour define the asperity, and the black dashed line provides the reference for the updip and downdip emissions relative to the strike).If more than 75% of the short-period emission points are located downdip or updip from the asperity, then the rupture is classified as downdip or updip; otherwise, it is defined to be balanced.For balanced rupture patterns and based on visual inspection, we classify a rupture as balanced-encircling for emission patterns clearly outlining the asperity with a semicontinuous ring, and as balanced-segmented otherwise, that is, where there are noticeable gaps in the pattern of emission points (see visual representation in Figure 8).
The 2010 Maule earthquake represents a classical downdip rupture with a bilateral propagation (∼29% of the emission points are southward of the epicenter; Figure 8b and Figure S12 in Supporting Information S1).The 2011 Tohoku-Oki earthquake ruptured bilaterally and mainly downdip (Figure 8a and Figure S16 in Supporting Information S1).However, emissions were also observed near the trench where a large shallow slip occurred (up to 50 m).Toward the end of the rupture (135-180 s; Figure S16a in Supporting Information S1), a subsidiary pattern with updip and downdip emissions encircled the second larger slip area off the coast of Fukushima Prefecture (see the southern 10 m contour in Figure 8a).The 2010 Mentawai earthquake propagated predominantly along the downdip edge of the main asperity and then ruptured updip near the trench region (Figure 8d and Figure S15 in Supporting Information S1).
Notably, the observations include several events with updip and balanced-encircling patterns (see graphical definition inset in Figure 8).For the 2016 Pedernales earthquake, the reverse of the conventional pattern occurred, first with short-period emissions updip of the main asperity, followed by a downdip propagation near the southern tip of the rupture (Figure 8j and Figure S41 in Supporting Information S1).The 2021 Kermadec Islands earthquake initially propagated to the northeast along with the main asperity, followed by a rupture moving updip in the final phase (Figure 8h and Figure S63 in Supporting Information S1).The 2014 Iquique-Pisagua earthquake presented a balanced-encircling rupture around the larger slip area (Figure 8f and Figure S31 in Supporting Information S1); a similar sequence defined the 2016 Melinka earthquake (Figure 8l and Figure S47 in Supporting Information S1).The 2021 Chignik earthquake (Figure 8g and Figure S64 in Supporting Information S1) ruptured eastward and encircled the slip area.More or less parallel updip and downdip emissions outlined the main slip asperity for the 2015 Illapel earthquake (Figure 8c and Figure S39 in Supporting Information S1), consistent with the back-projection study of Meng et al. (2018) for this event.The main rupture propagated northward along the downdip.A secondary front in the shallow part of the megathrust (less than ∼15 km depth) followed the region parallel to the trench for over ∼150 km.
Other events exhibited balanced-segmented ruptures.The 2014 Iquique earthquake (main M W 7.7 aftershock) showed a unilateral northeast short-period rupture partitioned around the asperity but still outlining the maximum slip region (Figure 8k and Figure S32 in Supporting Information S1).Offshore of the Alaska Peninsula, the 2020 Simeonof Island earthquake first propagated downdip to the northwest and west (Figure 8i and Figure S60 in Supporting Information S1), but the last short-period emissions are observed east of the epicenter.
To evaluate the depth distribution of short-period emissions more systematically, we present the along-dip horizontal distance distribution of radiation for the megathrust earthquakes shown in Figure 8; see Figure 9. Emissions located inside the asperity represented a median of only ∼5% of the total emissions, with some exceptions, for example, the 2014 Iquique-Pisagua earthquake with ∼48% of the radiation emitted within the asperity.For most events, short-period radiation is observed dominantly downdip, but the large variation in rupture patterns also demonstrated that the classical pattern of dominant downdip high-frequency radiation has exceptions.Notable examples are dominant updip short-period ruptures, for example, the 2016 Pedernales and 2021 Kermadec Islands earthquakes, and balanced ruptures, for example, the 2015 Illapel, 2016 Melinka, and 2021 Chignik earthquakes.There was no apparent relation between the distribution of short-period radiation and the centroid depth.
In addition, we compared multi-array backprojections for major megathrust earthquakes with IRIS backprojections.Specifically, we considered IRIS solutions from the Global Seismographic Network (GSN; http:// ds.iris.edu/spud/backprojection),which provides the tracking of the earthquake rupture from a wide azimuthal coverage.Figure S79 in Supporting Information S1 presents the IRIS backprojections (blue-red squares) compared to multi-array back-projections (blue-red circles, color-coded by rupture time) for the events shown in Figure 8.In general, the rupture tracking is comparable in both solutions regarding timing and directivity.Blue and red dots present updip and downdip short-period emissions relative to the earthquake asperity, which is defined by the area where the slip exceeds 75% of the maximum slip.Gray circles show emissions inside the asperity.White stars mark the horizontal distance from the epicenters.Earthquake magnitude and Global CMT centroid depths are shown on the top; earthquake rupture pattern classification on the bottom.The events are ordered by centroid depth.
However, the IRIS solutions have shown a higher susceptibility to artifacts, which is evident in the broader extent and sparsity of the rupture emissions.Thus, the combination of multiple arrays, weighted according to their azimuthal distribution, presented greater control over artifacts and, as a result, over the rupture tracking and the image of potential complexities.In the analysis, it is also important to consider the lower frequency band used in the IRIS-GSN solutions (0.05-0.25 Hz) relative to the multi-array back-projection (0.5-2.0 Hz).This lower frequency band could imply more significant smearing effects or artifacts due to the lower frequency range and, therefore, lower resolvability.The impact of the lower resolution, particularly the presence of artifacts, is more evident when estimating rupture parameters from IRIS backprojections.In effect, the results show similar directivity relative to this study but significantly larger rupture lengths and highly overestimated speeds (for a detailed comparison, refer to the rupture parameter estimates in Table S9 of Supporting Information S1).

Complex Megathrust Ruptures and Depth-Varying Short-Period Radiation
Seismicity prior to a great earthquake is often demonstrated to be active on the rupture zone edge while a seismic quiescence dominates the region later experiencing the maximum slip according to the "Mogi donut" hypothesis (e.g., Kelleher & Savino, 1975;Mogi, 1969).According to this hypothesis, the spots at the asperity edge become seismically active as the stress on the asperity increases, outlining the potential rupture area (e.g., Kanamori, 1981).For example, Schurr et al. (2020) observed a Mogi donut in the years prior to the 2014 Iquique-Pisagua earthquake.Sippl et al. (2021) showed microseismicity forming three half-ellipses along the central Chile megathrust, with one of them encircling the 2015 Illapel earthquake area.In back-projection, short-period rupture emissions are related to rupture velocity variations, in turn related to stress variations (e.g., Marty et al., 2019).They are often seen to concentrate around the main slip asperity and thus reflect presumably the high-stress gradient around the asperity, equivalently to the Mogi donut pattern of seismicity.
Many megathrust earthquakes exhibited balanced rupture patterns, where the back-projection highlights propagation to shallower and deeper regions relative to the asperity (peak slip) (Figure 9).This pattern could also be related to the presence of seismogenic barriers under the assumption of a strong rupture deceleration or stopping phase (e.g., Madariaga, 1983).For example, the 2015 Illapel event showed a balanced-encircling rupture (Figure 8c and Figure S39 in Supporting Information S1); it was first observed by Meng et al. (2018) from backprojection and interpreted as a splitting of rupture fronts surrounding a large asperity or barrier.Recently, a geometrical kink in the slab interface and the continental Moho depth were considered to explain rupture emissions both updip and downdip of the slip asperities for the 2007 Tocopilla earthquake; see Vera et al. (2023).The 2010 Mentawai (Figure 8d and Figure S15 in Supporting Information S1) and 2011 Tohoku-Oki (Figure 8a and Figure S16 in Supporting Information S1) earthquakes were tsunamigenic with a rupture that propagated up to the trench (e.g., Hill et al., 2012;Iinuma et al., 2012;Lay et al., 2011).This is compatible with the backprojections, where rupture emissions were partly located in shallow regions below accretionary structures.The transition from the seismogenic to a shallow aseismic sliding region might explain a strong rupture velocity variation and the shallow source of short-period radiation.For the 2016 Pedernales event, Agurto-Detzel et al. ( 2019) proposed a barrier mechanically controlled by the subduction of a rough oceanic relief in the updip limit of the earthquake rupture.Similarly, the back-projection showed a rupture updip of the main slip patch (at 15-20 km depth; Figure S41 in Supporting Information S1) that agrees with high residual (>∼1 km) bathymetry data (see Figure 9 of Agurto-Detzel et al. ( 2019)), indicating that the strong short-period earthquake radiation might be related to this barrier structure.
In summary, short-period emissions from the shallow part of the megathrust interface, updip of the main slip asperities, are more frequent than suggested by earlier results.

Impact of Automatic and Manual Rupture Duration Estimates
Automatic and visually determined rupture end times show some agreement (see inset in Figure S1 of Supporting Information S1), but there is substantial scatter among the shorter rupture times, for example, for manually picked rupture times less than 50 s, the average difference relative to automatic times is 20 s.Automatically determined rupture durations are, on average, longer than manually determined durations, but underestimations also occur.Thus, manually picked end times are sufficiently superior to prefer them even though they introduce subjectivity.Figure S2 and Table S5 in Supporting Information S1 summarize the rupture parameters derived from the automatic end times.Figure S3 in Supporting Information S1 compares the rupture parameters from the manually picked and automated results.

Rupture Speed
Figure 4 shows the rupture speed estimations derived from back-projection, additionally expressed as a percentage of the shear velocity at the hypocentral depth in the IASP91 model.The underlying rupture speed time series are summarized in Figure S11 of Supporting Information S1.Because the rupture speed estimates are based on the horizontal projection, rupture velocity might be underestimated if rupturing updip or downdip dipping faults, for example, for a dip angle of 30°, this effect would amount to ∼15%.However, this effect will only affect a small subset of intra-plate thrust or normal faults.
Thrust earthquakes averaged a rupture speed of 2.1 km/s and 56% of the shear wave speed, with a range between 25% and 86% of the shear wave speed, placing propagation speeds firmly in the expected sub-Rayleigh regime.The 2010 Mentawai, 2011 Tohoku-Oki, and 2021 South of Sandwich Islands earthquakes were slower than the average.The Mentawai and Tohoku-Oki events are characterized by large displacement at very shallow depths, implying a low normal stress environment and large tsunamis.The Sandwich Islands event also presented a shallower rupture for at least part of the rupture propagation.This observation agrees with the multiple independent technique results presented in Metz et al. (2022), that is, a shallow and slow rupture (1.5-2.1 km/s) that follows the varying strike direction along the trench, as part of a complex rupture also involving a faster propagation on deeper parts of the plate interface.
The 17 December 2016 Solomon Islands (M W 7.9) earthquake is a special case (Figure 8e and Figure S45 in Supporting Information S1).A megathrust rupture released the bulk of the moment, but its hypocenter is at 105 km depth, and the first motion mechanism indicated downdip intraslab faulting, which is thought to have triggered the main megathrust event (e.g., Lay et al., 2017;S.-J. Lee et al., 2018).In the back-projection, this sequence was visible as a ∼55 km gap between the first emissions and the dominant megathrust rupture pattern.The short time interval needed to cross the gap would imply an extremely fast rupture speed of 5.7 km/s, equivalent to ∼71% of the P wave speed at the hypocentral depth.This indicates that, most likely, P waves from the intraslab sub-event triggered the megathrust rupture (note that the relative P wave speed is calculated in the reference model, and P wave speeds in the down-going oceanic crust would be slower; thus, higher percentages are possible).To explore rupture propagation in the main phase, we placed the reference point at the megathrust rupture initiation (Figure S46 in Supporting Information S1); then, a more typical rupture propagation speed of 2.3 km/s was obtained.
For normal-faulting earthquakes, the average was 2.5 km/s and 60% of the shear wave speed, ranging between 38% and 82% of the shear wave speed.The 2019 Northern Peru (106 km depth; Figure S58 in Supporting Information S1) earthquake showed the fastest rupture (3.7 km/s and 82%V s ), which is at the upper limit of the sub-Rayleigh speed.Although this intermediate-depth event was fast, the depth does not appear to affect rupture velocity systematically.For example, the 2014 Rat Islands event at nearly the same depth (109 km; Figure S35 in Supporting Information S1) and the 2019 Ecuador-Peru Border Region earthquake (at 146 km depth; Figure S56 in Supporting Information S1) ruptured slowly at 2.0 and 1.7 km/s.For the Rat Island event, previous backprojections and finite fault modeling showed a similar rupture velocity in the range of ∼1.5-2.5 km/s (e.g., Twardzik & Ji, 2015;Ye et al., 2014).However, the lack of solid evidence for a shallow or steeply-dipping fault plane preference hampers the evaluation of the projection effects in the rupture velocity estimation.
Strike-slip earthquakes showed the greatest variability.The average rupture speed was 2.7 km/s and 78% of the shear wave speed, but the range spans from 37% to 137%.This is comparable with the observations of D. Wang et al. (2016), where backprojections showed strike-slip earthquakes propagating at average speeds greater than 3.0 km/s, faster than dip-slip events.Two events were clearly within the unstable supershear range (between V s and ̅̅ ̅ 2

√
V s ): the 2013 Craig (4.6 km/s and 137%V s ; Figure S26 in Supporting Information S1) and 2018 Palu (4.2 km/s and 125%V s ; Figure S54 in Supporting Information S1) earthquakes.The Palu event was supershear from early on as reported by Bao et al. (2019) using the high-resolution MUSIC back-projection, while the oceanic interplate Craig initiated at subshear speed and transitioned to supershear after ∼20 s (see time series and supershear reference in Figure S11 of Supporting Information S1).The average rupture speed for the Craig earthquake inferred by us agrees with the finite-fault modeling results (4-5 km/s) of Aderhold and Abercrombie (2015).Previously, Yue et al. (2013) provided faster peak velocity estimates from Sg and Sn arrivals (>4.5 km/s) and finite-fault inversion (5.5-6 km/s), but as we do not consider time-variable propagation rates, these observations are not in conflict.The 2013 Balochistan event propagated at a speed of 3.3 km/s (Figure S29 in Supporting Information S1), representing 98% of the shear wave speed, that is, around the shear velocity but (very likely) faster than Rayleigh waves.Additionally, the 2018 North of Honduras (95%V s ; Figure S51 in Supporting Information S1) and 2020 Caribbean (104%V s ; Figure S59 in Supporting Information S1) earthquakes propagated at "likely" supershear (95-105%V s ; by assuming expected rupture velocity errors in the estimates).For the Caribbean earthquake, Tadapansawut et al. (2021) inferred supershear rupture fronts with peak velocities larger than 5 km/s from finite-fault modeling, that is, significantly faster than our average-based estimates.
For average velocity estimates based on manually and automatically determined rupture end times, please consult Tables S4-S5 in Supporting Information S1 and Data Sets S2-S3.Additionally, refer to Figure S80 in Supporting Information S1 for an illustrative example of the estimation of average velocities based on manually determined rupture end times.Overall, the average rupture speeds are approximately 70% and 65% of their respective maximum values.This is in relation to the and automatic rupture end times, respectively.The decrease in estimated rupture speed based on the averages could potentially lead to an underestimation or complicate the evaluation of supershear ruptures.This effect may be particularly pronounced when a rupture transitions from subshear to supershear, as observed in the 2013 Craig earthquake (shifting from 4.6 km/s or 137%V s to 2.9 km/s or 86%V s ; Table S4 in Supporting Information S1).Therefore, we preferred maximum velocity values over averages.

Rupture Length
The rupture length estimates from backprojections were used to derive magnitude-length scaling relations for each faulting slip type (normal, thrust, and strike-slip), based on the logarithm of rupture length and moment magnitude (Figure 10 and Table S6 in Supporting Information S1).Due to the composite nature of the earthquake and variability in moment estimates, the 2021 Sandwich Islands event was not included in the regression for thrust magnitude-length scaling relation.Thrust and normal earthquakes showed similar magnitude-length dependencies, while strike-slip earthquakes had a large data dispersion on average.We additionally compared the newly short-period scaling relations to established scaling relations (e.g., Blaser et al., 2010;Wells & Coppersmith, 1994) derived from aftershocks, geodetic modeling, and empirical  S6 in Supporting Information S1).For normal faulting earthquakes (Figure 10a), our scaling relation is very similar to those derived from both Wells and Coppersmith (1994) and Blaser et al. (2010).For thrust earthquakes (Figure 10b), our scaling relation predicts shorter rupture lengths than Wells and Coppersmith (1994) but close to the predictions of Blaser et al. (2010), which is based on a much larger database.Lengths for strike-slip earthquakes presented a large dispersion, with the regression pointing to slower growth of rupture length with magnitude than either Blaser et al. (2010) or Wells and Coppersmith (1994) but still within the expected, albeit large, uncertainty range (Figure 10c).The first limitation here is the sparsity of strikeslip events with M W > ∼8.Second, several of the large strike-slip earthquakes were complex and included multiple conjugate faults, leading to large ambiguity in the interpretation of the rupture length as determined by our algorithm (as the cumulative fault length ruptured will be much larger than this estimate).This pattern is often found for intraplate earthquakes in the oceanic lithosphere, for example, the great 2012 Wharton Basin (e.g., Duputel et al., 2012;Hill et al., 2015;Meng, Ampuero, Stock, et al., 2012;Satriano et al., 2012) and 2018 Gulf of Alaska (e.g., Krabbenhoeft et al., 2018;Lay et al., 2018;Ruppert et al., 2018) earthquakes, where backprojections were broadly compatible with the more comprehensive assessments for those events, for example, aftershocks, seafloor topography and/or finite fault modeling.

Aspect Ratio
We remind that manually determined rupture end times were preferred over automatic values.Therefore, the following analysis benefits from excluding potential outliers after the rupture ends.We also considered the full back-projected rupture, that is, we have not removed potential outliers from the middle part of the rupture.
We explored the relationship between aspect ratios and source complexity.Simple elongated ruptures were characterized by small aspect ratios (aspect ratio close to 0: very low sparsity in the rupture emissions), for example, the 2017 Komandorsky Islands (ratio ∼0.19; Figure S49 in Supporting Information S1), 2016 Wharton Basin (ratio ∼0.25; Figure S40 in Supporting Information S1), and 2011 Kermadec Islands (ratio ∼0.17; Figure S19 in Supporting Information S1) earthquakes (all of them bilateral ruptures).In contrast, complex events with no straightforward interpretation of "rupture length" tended to have large values (aspect ratio close to 1: high sparsity in the rupture emissions).For strike-slip earthquakes, aspect ratios close to 1 often indicate the activation of multiple faults.Prominent examples are the 2012 Wharton Basin (ratio ∼0.81; Figure S21 in Supporting Information S1) and Gulf of Alaska (ratio ∼0.70; Figure S52 in Supporting Information S1) earthquakes that involved rupture of conjugate fault systems (e.g., Krabbenhoeft et al., 2018;Satriano et al., 2012).
It should be noted that some strike-slip earthquakes exhibit fluctuations in high-frequency rupture emissions that are not reflected in their small aspect ratio estimates due to their total rupture length.For instance, the 2018 North of Honduras (M W 7.5) strike-slip earthquake presented a low ratio (ratio ∼0.29; Figure S51 in Supporting Information S1).This low ratio indicates an elongated rupture, as observed in the back-projection (Figure S51 in Supporting Information S1).However, it is also possible to observe NS fluctuations in the rupture emissions that are not reflected in the overall ratio due to the large EW extent of the rupture (rupture length: 143 km; Figure S51c in Supporting Information S1).Aspect ratios for sub-segments of the North of Honduras rupture would indicate higher values relative to the total evaluation of the earthquake extent, likely attributed to fluctuations in the highfrequency radiation, rather than having a physical meaning.Supershear earthquakes were found to be associated with the most mature parts of the rupture zone (Perrin et al., 2016).Since fault length can be considered a proxy for the fault maturity (e.g., Manighetti et al., 2007;Perrin et al., 2016), back-projected speeds in the supershear range combined with linear faults (small aspect ratios < ∼0.3), indicate a large degree in fault maturity, for example, Palu (ratio ∼0.21; Figure S54 in Supporting Information S1) and Craig (ratio ∼0.31; Figure S26 in Supporting Information S1) supershear earthquakes.Not all strike-slip earthquakes with a small aspect ratio were supershear.For example, the 2016 Wharton Basin is a sub-shear event with a small aspect ratio (ratio ∼0.25; Figure S40 in Supporting Information S1).For thrust earthquakes, the aspect ratio also quantified complexity.However, the complexity usually does not arise from a rupture on different faults but from the distribution of high-frequency emissions on the megathrust, for example, the 2015 Illapel earthquake and its balanced rupture pattern (ratio ∼0.84; Figure S39 in Supporting Information S1).
Regardless of the faulting slip type, a compact rupture can also result in aspect ratios close to 1 due to the greater importance of variance in the estimated locations of emission points for smaller ruptures, for example, the 2021 Loyalty Islands (ratio ∼0.84; Figure S62 in Supporting Information S1) and 2014 Solomon Islands (ratio ∼0.61; Figure S34 in Supporting Information S1) earthquakes.Therefore, rupture complexity is linked to aspect ratios > ∼0.5, but such a large value does not necessarily imply complexity but can also indicate compactness, particularly for earthquakes with M W < ∼8.

Conclusions
We present a catalog of short-period rupture histories for 56 large earthquakes (0.5-2.0 Hz; M W ≥ 7.5; 01/2010-12/2022) based on a new implementation of the teleseismic back-projection method.The approach includes multiple arrays and combines P and pP waveforms for earthquakes deeper than 40 km, which also considers the radiation pattern, that is, in selecting only those arrays where the pP depth phase is expected to have a significant amplitude.Based on the backprojections, rupture length, directivity, speed, and aspect ratio were estimated algorithmically.The main findings are as follows: 1. We found differences in rupture patterns between finite fault slip models and short-period emissions for subduction megathrust earthquakes.We confirmed the preference for short-period seismic energy to be emitted downdip of the main slip asperity (as previously reported in the literature), but additionally identified many examples of balanced rupture patterns around the main asperity, including updip short-period ruptures.2. Short-period ruptures for megathrust earthquakes, for example, balanced patterns around slip patches, might result from the high-stress gradient around asperities and seismogenic barriers (strong gradients in the slip are needed to generate high-frequency emissions).3. Earthquake rupture speeds were consistently in the sub-Rayleigh regime for thrust and normal faulting earthquakes, with a median of 56% and 60% of the shear wave speed, respectively.Strike-slip earthquakes showed the greatest variability with a median of 78% and a range between 37% and 137% of the shear wave speed.The 2013 Craig (137%V s ) and 2018 Palu (125%V s ) events propagated in the unstable supershear range.
The 2013 Balochistan (98%V s ), 2018 North of Honduras (95%V s ), and 2020 Caribbean (104%V s ) earthquakes were "likely" supershear (95-105%V s ; by assuming expected rupture velocity errors in the estimates) and in any case faster than Rayleigh waves.4. We presented new scaling relations from short-period backprojections comparable to finite difference methods.Thrust and normal earthquakes showed a similar magnitude-length relationship compared to established "long-period" relations.Strike-slip events presented a large dispersion but still within the expected uncertainty of previously published scaling relations.Limitations for strike-slip earthquakes were the lack of events with a large moment magnitude (M W > ∼8) and the difficulty of interpreting rupture lengths for ruptures on complex conjugate fault systems.
Whereas overall, the method is suitable for near-real-time processing, the array selection and the definition of secondary search grids need to be done manually.Furthermore, the automatic determination of the end of the rupture is not always robust and needs to be confirmed manually.This paper has analyzed the fault rupture histories for 56 earthquakes, but we intend to continue to analyze future earthquakes with M W ≥ 7.5 and depths less than 200 km.

Data Availability Statement
The back-projection catalog is available in Supporting Information S1 as map views and in machine-readable format (Data Set S1), including the earthquake parameters derived from back-projection (Data Set S2) and parameter estimates based on automatic rupture end times (Data Set S3).The data files used in this paper are available at Vera et al. (2024): GFZ Data Services (https://doi.org/10.5880/GFZ.2.4.2024.001).Earthquake source parameters obtained from the CSN, GEOFON, ISC, SSN, USGS-NEIC catalogs, and Sippl et al. (2018), Lay et al. (2013), Nocquet et al. (2017), Bilek et al. (2011), Zhao et al. (2011), and Adhikari et al. (2015); see Tables S1-S3 in Supporting Information S1 for additional details.Fault geometry for the 2013 Balochistan (Figure S29 in Supporting Information S1) and 2016 Kaikoura (Figure S43 in Supporting Information S1) earthquakes extracted from Avouac et al. (2014) and Hamling et al. (2017).Major plate boundaries from Styron and Pagani (2020).Finite fault slip solutions for the 2015 Gorkha (Figure S37 in Supporting Information S1) and 2017 Chiapas (Figure S50 in Supporting Information S1) earthquakes from Yagi and Okuwaki (2015) and Okuwaki and Yagi (2017).Hi-net data was provided by the National Research Institute for Earth Science and Disaster Prevention (NIED).NIED/Hi-net data was obtained with the HinetPy Python package (Tian, 2020).The

Figure 1 .
Figure 1.Multi-array multi-phase back-projection.(a) Workflow for rupture imaging using P and pP waveform arrivals and multiple seismic arrays.pP backprojections are only included for earthquakes with depth ≥40 km and are weighted according to radiation pattern (see text).The last step involves the extraction of a few rupture parameters based on the timing and locations of rupture maxima.(b) Array weighting.The sum of azimuthal half-angles between the target and its two neighboring arrays is proportional to the weights γ.(c) Automatic estimation of the rupture aspect ratio, speed, length, and directivity; see text for details.

Figure 2 .
Figure 2. Global distribution of back-projected earthquakes colored by focal mechanism.Earthquakes back-projected with the combined P and pP approach are highlighted in magenta color.Labels designate the earthquake ID, with details given in Figure4and TablesS1-S3in Supporting Information S1.Red segments show major plate boundaries(Bird, 2003).

Figure 3 .
Figure 3.The 25 March 2020 East of Kuril Islands (M W 7.5) earthquake back-projection (0.5-2.0 Hz).ID: 16-Origin Time: 2020-03-25 02:49:21.2UTC-Lon: 157.70°E-Lat: 48.96°N-Depth: 58 km-Type: thrust.(a) Rupture pattern.Blue-red dots show semblance maxima tracking the earthquake rupture color-coded by time and scaled by energy radiated.Black dotted contours outline the short-period energy radiated (normalized to 1).The yellow-red polygons present the USGS-NEIC finite fault slip solution for comparison.Trench line derived from SLAB1.0 model (Hayes et al., 2012).Focal mechanism from Global CMT catalog.Inset: Array weights and energy radiated source time function.(b) Multi-array configuration and time-shifted P and pP waveforms.(c) Automatic rupture parameter determination; see Figure 1c for further information on the format.

Figure 4 .
Figure 4. Earthquake rupture parameters from back-projection (0.5-2.0 Hz).Event ID (for cross-referencing to other figures, tables, and Supporting Information S1), moment magnitude, rupture length (colored by rupture time), aspect ratio, directivity, and rupture speed (colored by hypocentral depth).The events are categorized into normal (upper), thrust (middle), and strike-slip (bottom) earthquakes, with the order in each category determined by magnitude.Events in bold text are assumed to be subduction megathrust earthquakes.Events labeled with a black bullet indicate that a secondary rupture pattern was obtained over a predefined grid and included in the rupture parameter estimates.The vertical red bars in the rupture speed panel show the expected shear wave speed (V S ) at the hypocentral depth in the IASP91 velocity model(Kennett & Engdahl, 1991).Rupture speeds in the range 95-105%V S and >105%V S define "likely' and "strongly" supershear earthquakes, respectively.

Figure 5 .
Figure 5. Short-period (0.5-2.0 Hz) energy radiated source time functions derived from back-projection.The energy radiated maxima evaluated over the rupture time provides the source time function (in gray).Additionally, SCARDEC (in blue) and USGS-NEIC (in red) moment rate functions are shown for comparison.The event highlighted with a dashed box is used as an example and discussed in detail in the text.

Figure 6 .
Figure 6.Back-projection of northeast-propagating rupture synthetics with variations in the number of combined arrays.(a) Synthetic rupture scenario.(b) Synthetic waveforms for the European (EU), North American (USA), Australian (AU), and China (CN) array.A pP/P amplitude ratio of 1 is assumed.(c) P wave (top panel) and combined P and pP (bottom panel) backprojections.The number of combined arrays increases from left to right.The white stars mark the location of the emission points.Inset: Array weights and energy radiated source time function.Rupture parameters estimates; rupture speed, length, directivity, and aspect ratio.

Figure 7 .
Figure 7. Back-projection of northeast-propagating rupture synthetics with variations in the pP/P amplitude ratio.P (top panel) and combined P and pP (bottom panel) backprojections.The pP/P amplitude ratio increases from left to right, expressed in percentages.Input rupture speed equivalent to 2.5 km/s.The white stars mark the true emission points.Inset: Array weights and energy radiated source time function.Rupture parameters estimates; rupture speed, length, directivity, and aspect ratio.

Figure 8 .
Figure 8. Back-projected earthquake rupture classification for major megathrust earthquakes (0.5-2.0 Hz).Colored dots show semblance maxima, scaled by energy radiated: blue for maxima downdip of the main asperity, red for those updip, and gray for emissions inside the asperity (see text for details of the definition of asperity).The dashed black line indicates the transect along with the fault strike and maximum slip (magenta cross mark) used to classify the rupture pattern relative to the asperity (black contour).The subplots are sorted by rupture duration.The dark-magenta arrows (lower-right corner) indicate the multi-array distribution relative to the epicenter (white star).Trench location from SLAB1.0 model (Hayes et al., 2012); focal mechanisms from Global CMT catalog.The yellow-red background shows the slip distribution of (a) Iinuma et al. (2012), (b) Moreno et al. (2012), (c) Tilmann et al. (2016), (d), (e), (g) and (h) USGS-NEIC finite fault solution, (f) and (k) Schurr et al. (2014), (i) Crowell and Melgar (2020), (j) Heidarzadeh et al. (2017), and (l) Moreno et al. (2018).Bottom inset: High-frequency rupture pattern classification.The dominant type is noted in the upper of each plot.

Figure 9 .
Figure9.Along-dip horizontal distance distribution of short-period radiation (0.5-2.0 Hz) for large megathrust earthquakes.Blue and red dots present updip and downdip short-period emissions relative to the earthquake asperity, which is defined by the area where the slip exceeds 75% of the maximum slip.Gray circles show emissions inside the asperity.White stars mark the horizontal distance from the epicenters.Earthquake magnitude and Global CMT centroid depths are shown on the top; earthquake rupture pattern classification on the bottom.The events are ordered by centroid depth.

Figure 10 .
Figure 10.Rupture length scaling relations for short-period back-projection based estimates.(a) Normal, (b) thrust, and (c) strike-slip earthquakes.The black line shows the back-projection relationships (solid line) and the 95% of confidence interval predicted (dashed line).The squares, circles, and triangles present the data distribution for each faulting type.The gray circle shows the 2021 South of Sandwich Islands earthquake not included in the regression.Blue and red lines show the well-known scaling relationships of Blaser et al. (2010) andWells and Coppersmith (1994) for comparison.

ournal of Geophysical Research: Solid Earth
VERA ET AL.