A new method for determining small earthquake source parameters using short-period P waves

We developed a new technique of inverting short-period (0.5–2 Hz) P waveforms for determining small earthquake (M <3.5) focal mechanisms and moments, where magnitude ~4 events with known source mechanisms are used to calibrate the "unmodeled" structural effect. The calibration is based on a waveform cluster analysis, where we show that clustered events of different sizes, for example, M ~4 versus M ~2, display similar signals in the short-period (SP, 0.5–2 Hz) frequency band, implying propagational stability. Since both M ~4 and M ~2 events have corner frequencies higher than 2 Hz, they can be treated as point sources, and the "unmodeled" structural effect on the SP P waves can be derived from the magnitude 4 events with known source mechanisms. Similarly, well-determined magnitude 2’s can provide calibration for studying even smaller events at higher frequencies, for example, 2–8 Hz. In particular, we find that the "unmodeled" structural effect on SP P waves is mainly an amplitude discrepancy between data and 1D synthetics. The simple function of "amplitude amplification factor" (AAF) defined as the amplitude ratio between data and synthetics provides useful calibration, in that the AAFs derived from different clustered events appear consistent, hence stable and mechanism independent. We take a grid-search approach to determine source mechanisms by minimizing the misfit error between corrected data and synthetics of SP P waves. The validation tests with calibration events demonstrate the importance and usefulness of the AAF corrections in recovering reliable results. We introduce the method with the 2003 Big Bear sequence. However, it applies equally well to other source regions in southern California, because we have shown that the mechanism independence and stability of the AAFs for source regions of 10 km by 10 km are typical. By definition, the AAFs contain the effects from the station site, the path, and crustal scattering. Although isolating their contributions proves difficult, the mechanism independence and stability of the AAFs suggest that they are mainly controlled by the near-receiver structure. Moreover, the ratios between the AAFs for the vertical and radial components from various events at different locations appear consistent, suggesting that these AAF(v)/AAF(r) ratios might be simple functions of site conditions. In this study, we obtained the focal mechanisms and moments for 92 Big Bear events with M_L down to 2.0. The focal planes correlate well with the seismicity patterns, while containing abundant finer-scale fault complexity. We find a linear relationship between log(M_0) and M_L, that is, log(M_0) = 1.12M_L + 17.29, which explains all the data points spanning three orders of magnitude (2.0 < M_L < 5.5)

Determination of earthquake focal depths and source time functions in central Asia using teleseismic P waveforms

We developed a new method to determine earthquake source time functions and focal depths. It uses theoretical Green's function and a time-domain deconvolution with positivity constraint to estimate the source time function from the teleseismic P waveforms. The earthquake focal depth is also determined in the process by using the time separations of the direct P and depth phases. We applied this method to 606 earthquakes between 1990 and 2005 in Central Asia. The results show that the Centroid Moment Tensor solutions, which are routinely computed for earthquake larger than M5.0 globally using very long period body and surface waves, systematically over-estimated the source depths and durations, especially for shallow events. Away from the subduction zone, most of the 606 earthquakes occurred within the top 20 km of crust. This shallow distribution of earthquakes suggests a high geotherm and a weak ductile lower crust in the region

Validating tomographic model with broad-band waveform modelling: an example from the LA RISTRA transect in the southwestern United States

Traveltime tomographic models of the LA RISTRA transect produce excellent waveform fits if we amplify the damped images. We observe systematic waveform distortions across the western edge of the Great Plains from South American events, starting about 300 km east of the centre of the Rio Grande Rift. The amplitude decreases by more than 50 per cent within array stations spanning less than 200 km while the pulse width increases by more than a factor of 2. This feature is not observed for the data arriving from the northwest. While the S-wave tomographic image shows a fast slab-like feature dipping to the southeast beneath the western edge of the Great Plains, synthetics generated from this model do not reproduce the waveform characteristics. However, once we modify the tomographic image by amplifying the velocity contrast between the slab and adjoining mantle by a factor of 2–3, the synthetics produce observed amplitude decay and pulse broadening. In addition to the traveltime delay, amplitude variation due to wave phenomena such as slab diffraction, focusing and defocusing provide much tighter constraints on the geometry of the fast anomaly and its amplitude and sharpness as demonstrated by a forward sensitivity test and snapshots of the seismic wavefield. Our preferred model locates the slab 200 km east of the Rio Grande Rift dipping 70°–75° to the southeast, extending to a depth near 600 km with a thickness of 120 km and a velocity of about 4 per cent fast. In short, adding waveform and amplitude components to regional tomographic studies can help validate and establish structural geometry, sharpness and velocity contrast

Juan de Fuca subduction zone from a mixture of tomography and waveform modeling

Seismic tomography images of the upper mantle structures beneath the Pacific Northwestern United States display a maze of high-velocity anomalies, many of which produce distorted waveforms evident in the USArray observations indicative of the Juan de Fuca (JdF) slab. The inferred location of the slab agrees quite well with existing contour lines defining the slab's upper interface. Synthetic waveforms generated from a recent tomography image fit teleseismic travel times quite well and also some of the waveform distortions. Regional earthquake data, however, require substantial changes to the tomographic velocities. By modeling regional waveforms of the 2008 Nevada earthquake, we find that the uppermost mantle of the 1D reference model AK135, the reference velocity model used for most tomographic studies, is too fast for the western United States. Here, we replace AK135 with mT7, a modification of an older Basin-and-Range model T7. We present two hybrid velocity structures satisfying the waveform data based on modified tomographic images and conventional slab wisdom. We derive P and SH velocity structures down to 660 km along two cross sections through the JdF slab. Our results indicate that the JdF slab is subducted to a depth of 250 km beneath the Seattle region, and terminates at a shallower depth beneath Portland region of Oregon to the south. The slab is about 60 km thick and has a P velocity increase of 5% with respect to mT7. In order to fit waveform complexities of teleseismic Gulf of Mexico and South American events, a slab-like high-velocity anomaly with velocity increases of 3% for P and 7% for SH is inferred just above the 660 discontinuity beneath Nevada

Upper mantle P velocity structure beneath the Midwestern United States derived from triplicated waveforms

Upper mantle seismic velocity structures in both vertical and horizontal directions are key to understanding the structure and mechanics of tectonic plates. Recent deployment of the USArray Transportable Array (TA) in the Midwestern United States provides an extraordinary regional earthquake data set to investigate such velocity structure beneath the stable North American craton. In this paper, we choose an M_w5.1 Canadian earthquake in the Quebec area, which is recorded by about 400 TA stations, to examine the P wave structures between the depths of 150 km to 800 km. Three smaller Midwestern earthquakes at closer distance to the TA are used to investigate vertical and horizontal variations in P velocity between depths of 40 km to 150 km. We use a grid-search approach to find the best 1-D model, starting with the previously developed S25 regional model. The results support the existence of an 8° discontinuity in P arrivals caused by a negative velocity gradient in the lithosphere between depths of 40 km to 120 km followed by a small (∼1%) jump and then a positive gradient down to 165 km. The P velocity then decreases by 2% from 165 km to 200 km, and we define this zone as the regional lithosphere-asthenosphere boundary (LAB). Beneath northern profiles, waves reflected from the 410 discontinuity (410) are delayed by up to 1 s relative to those turning just below the 410, which we explain by an anomaly just above the discontinuity with P velocity reduced by ∼3%. The 660 discontinuity (660) appears to be composed of two smaller velocity steps with a separation of 16 km. The inferred low-velocity anomaly above 410 may indicate high water concentrations in the transition zone, and the complexity of the 660 may be related to Farallon slab segments that have yet to sink into the deep mantle

Lower mantle tomography and phase change mapping

A lower mantle S wave triplication (Scd) has been recognized for many years and appears to be explained by the recently discovered perovskite (PV) to postperovskite (PPV) phase change. Seismic observations of Scd display (1) rapid changes in strength and timing relative to S and ScS and (2) early arrivals beneath fast lower mantle regions. While the latter feature can be explained by a Clapeyron slope (λ) of 6 MPa/K and a velocity jump of 1.5% when corrected by tomographic predictions, it does not explain the first feature. Here, we expand on this mapping approach by attempting a new parameterization that requires a sample of D" near the ScS bounce point (δ VS) where the phase height (hph) and velocity jump (β) are functions of (δ VS). These parameters are determined by modeling dense record sections collected from USArray and PASSCAL data where Grand's tomographic model is the most detailed in D" structure beneath Central America. We also address the range of λ to generate new global models of the phase boundary and associated temperature variation. We conclude that a λ near 9 MPa/K is most satisfactory but requires β to be nonuniform with a range from about 1.0 to 4.0% with some slow region samples requiring the largest values. Moreover, the edges of the supposed buckled slabs delimitated by both P and S waves display very rapid changes in phase boundary heights producing Scd multipathing. These features can explain the unstable nature of the Scd phase with easy detection to no detection commonly observed. The fine structure at the base of the mantle beneath these edges contains particularly strong reflections indicative of local ultralow velocity zones, which are predicted in some dynamic models

A Comparative Study of the Sumatran Subduction-Zone Earthquakes of 1935 and 1984

A M_s 7.7 earthquake struck the western, equatorial coast of Sumatra in December 1935. It was the largest event in the region since the two devastating giant earthquakes of 1833 and 1861. Historical seismograms of this event from several observatories around the world provide precious information that constrains the source parameters of the earthquake. To more precisely quantify the location, geometry, and mechanism of the 1935 event and to estimate the coseismic deformation, we analyze the best of the available teleseismic historical seismograms by comparing systematically the records of the 1935 earthquake with those of a smaller event that occurred in the same region in 1984. First we constrain the source parameters of the 1984 event using teleseismic records. Then, we compare the records of the 1935 event with those of 1984 from the same sites and instruments. To do this, we choose several time windows in the corresponding seismograms that contain clearly identifiable phases and deconvolve the modern event from the older one. The deconvolutions result in very narrow pulses with similar sizes, thus confirming similar locations and mechanisms for the events. The initiation of the 1984 event was on the subduction interface at a depth of 27 ± 2 km; its M_0 is 6.5 x 10^(19) N m (M_w is 7.2). The sense of slip was nearly pure thrust, on a plane dipping 12°. The 1935 event also involved rupture of the shallow subduction interface, but was about five times larger (M_0 3.3 x 10;^(20) N m, M_w 7.7) and initiated a few kilometers to the southeast, along strike. The 1935 rupture propagated unilaterally toward the southeast. The along-strike rupture length was about 65 km. From these source parameters, we calculate the surface deformations, assuming an elastic multilayered medium. These deformations compare favorably with those actually recovered from paleoseismic data in the form of coral microatolls

Rupture Directivity Characteristics of the 2003 Big Bear Sequence

We have developed a forward modeling technique to retrieve rupture characteristics of small earthquakes (3 < M < 5), including rupture propagation direction, fault dimension, and rupture speed. The technique is based on an empirical Green’s function (EGF) approach, where we use data from collocated smaller events as Green’s functions to study the bigger events. We tend to choose smaller events with similar focal mechanisms for EGFs; however, we show that the events with different focal mechanisms can work equally well when corrected for radiation pattern effect. Compared to deconvolution, this forward modeling approach allows full use of both the shape and amplitude information produced by rupture propagation. Assuming a simple 1D source model, we parameterize the source time function of a target event as the convolution of two boxcars, featuring the rise time τ_r and the rupture time τ_c; we solve for τ_r and τ_c in a grid search manner by minimizing the waveform misfit between the three-component data and the synthetics constructed from the EGFs. The rupture propagation direction, fault length, and rupture speed can then be estimated by fitting the observed azimuthal pattern of τ_c from P and S waves. We apply the approach to the 12 largest events (M_w ≥ 3:3) of the 2003 Big Bear sequence (excluding the mainshock) in southern California. Among them, seven events are found to exhibit robust rupture directivity. The fact that the ruptures of these events propagate in all directions reveals complicated fault geometry at depth. We compute the stress drop Δσ ~ 2 M_0/π L^3 for each event using the resolved fault length. The results show large variations ranging from ~1 to 90 Mpa, with no dependence on moment. However, Δσ appears inversely correlated with rupture speed V_r; in particular, events with larger Δσ tend to propagate at smaller V_r, whereas events with smaller Δσ propagate faster

Upper-mantle structures beneath USArray derived from waveform complexity

Tomographic imaging of the crust and upper mantle beneath the western United States has greatly improved with the addition of USArray data. These models display many detailed images of both fast and slow blobs penetrating into the transition zone. To study such features, we apply a newly developed technique, called MultiPath Detector analysis, to the SH waveform data. The method simulates each observed body waveform by performing a decomposition; by [S(t)+C×S(t–Δ_9LR))]/2, where S(t) is the synthetics for a reference model. Time separation Δ_9LR) and amplitude ratio C are needed to obtain a high cross-correlation between a simulated waveform and data. The travel time of the composite waveform relative to the reference model synthetics is defined as Δ_T. A simulated annealing algorithm is used to determine the parameters Δ_(LR) and C. We also record the amplitude ratio (Amp) between the synthetics for the reference model relative to the data. Generally, large Δ_(LR) values are associated with low Amp's. Whereas the conventional tomography yields a travel time correction (Δ_T), our analysis yields an extra parameter (Δ_(LR)), which describes the waveform complexity. With the array, we can construct a mapping of the gradient of Δ_(LR) with complexity patterns. A horizontal structure introduces waveform complexity along the distance profile (in-plane multipathing). An azimuthally orientation Δ_(LR) pattern indicates a vertical structure with out-of-plane multipathing. Using such maps generated from artificial data, we can easily recognize features produced by dipping fast structures and slow structures (DSS). Many of these features display organized waveform complexity that are distinctly directional indicative of dipping sharp-edges. Here, we process the array data for events arriving from various azimuths and construct maps of multipathing patterns. The similarity between tomographic features and complexity maps is striking. When features are dipping such as the slab structures beneath the Cascade Range and Nevada, strong complexity is observed from Southeastern events arriving along these ray paths with split pulses separated up to 6 s for both. This requires extended slab segments to at least 600/300 km with a 4/8 per cent velocity jump along the edges. One of the most dramatic set of DSS observations is associated with a slow northwest dipping conduit beneath Yellowstone that extends into the transition zone. A number of forward modelling experiments are included for the strongest patterns formed by sharpening present tomographic images