Cascading elastic perturbation in Japan due to the 2012 M w 8.6 Indian Ocean earthquake.

Since the discovery of extensive earthquake triggering occurring in response to the 1992 M w (moment magnitude) 7.3 Landers earthquake, it is now well established that seismic waves from earthquakes can trigger other earthquakes, tremor, slow slip, and pore pressure changes. Our contention is that earthquake triggering is one manifestation of a more widespread elastic disturbance that reveals information about Earth's stress state. Earth's stress state is central to our understanding of both natural and anthropogenic-induced crustal processes. We show that seismic waves from distant earthquakes may perturb stresses and frictional properties on faults and elastic moduli of the crust in cascading fashion. Transient dynamic stresses place crustal material into a metastable state during which the material recovers through a process termed slow dynamics. This observation of widespread, dynamically induced elastic perturbation, including systematic migration of offshore seismicity, strain transients, and velocity transients, presents a new characterization of Earth's elastic system that will advance our understanding of plate tectonics, seismicity, and seismic hazards.

INTRODUCTION

The loading of quasi-static stresses by the differential motions of tectonic plates is the primary driver of earthquakes and the deformation of the Earth’s crust at or near plate boundaries. Where and how the crust deforms largely depends on the elastic and frictional properties of the crust both in fault zones and in the bulk crust. Dynamic stresses from seismic waves can influence the nature and timing of crustal deformation by perturbing elastic properties. On the basis of field () laboratory studies () and simulation studies (), we posit that seismic waves induce instantaneous, or near-instantaneous, perturbations in elastic moduli, stress state, and frictional properties, including the behavior of fault gouge material, hereafter referred to as “elastic changes,” which persist during a recovery period and which may be observed with the proper instrumentation and metrics. In this instance, the mechanism for elastic change is a disturbance in the manner stress is transmitted across a fault or region (). In the case of granular material within a fault zone, force chains are destabilized and rearranged (–). Force chains are chains of connected particles within a layer of compressed granular material that support stresses across the layer. Within the bulk crust, there could be similar behavior at the scale of grains and cracks based on observations from large numbers of laboratory studies (, ).

We can probe the evolution of elastic properties via spatial and temporal changes in seismicity, seismic velocities, seismic attenuation, and crustal deformation. Here, we analyze Global Positioning System (GPS), strain meter, and Hi-net short-period seismic data recorded in northeast Japan in the days after the 2012 Mw (moment magnitude) 8.6 Indian Ocean earthquake (IOE). In response to perturbation by seismic waves from the IOE, we observe induced faulting, increased seismicity, crustal deformation, and velocity changes. Each observation is a manifestation of dynamically induced crustal metastability, with each having its own time constant for recovery. Such a widespread perturbation due to seismic waves reveals the interconnectivity of the Earth’s elastic system.

RESULTS

Seismicity

We begin our investigation by examining two Mw >5.5–triggered events identified by Pollitz et al. () that occurred off the east coast of Japan in the days after the 2012 Mw 8.6 IOE (Fig. 1). We find that these two events are part of several earthquake clusters that initiate along a spatial and temporal trend and migrate from northeast to southwest at a rate of ~70 km/day (Fig. 2A). Statistically, the observed increase in seismicity after the IOE is a 1 in 358 event in the Tohoku aftershock period and region when considering the geographic extent of the increase (see Materials and Methods). Most of the triggered earthquakes were normal faulting events in the shallow accretionary wedge, consistent with the state of stress in this region after the 2011 Mw 9.0 Tohoku earthquake (TOE) (, ). The two largest events, labeled “P1” (Mw 5.5) and “P2” (Mw 5.7), initiated ~30 and ~50 hours after the IOE, respectively (Table 1). The delay between the IOE and the initiation of the first cataloged event in this sequence prompted us to investigate noncataloged seismicity during this delay.

Fig. 1

Maps.

(A and B) These maps show (A) the IOE and (B) our study area in Japan. Symbols in (B) are as follows: blue and green triangles are Hi-net stations used to calculate seismic velocities; green stations are used to calculate seismic coherence. Red circles are the epicenters of the 2011 TOE and shallow earthquakes discussed in the text. Black triangles are F-net station HRO and extensometer KTA. Inset shows location of shallow earthquakes; hatch marks indicate region where seismic velocities are calculated.

Fig. 2

Observations.

(A) Earthquakes (origin 31.32°N, 134.10°E). Red triangles represent P1 and P2 (top to bottom). Blue triangle represents P0. Black vertical line extending through all subfigures is IOE. (B) Stacked interstation seismic coherence. (C) Cumulative aftershocks for P0 and P2 (black stars) with map view inset. Blue and red curves are cumulative number of earthquakes in the P0 and P2 clusters, respectively. The P0 main shock and northern aftershocks (blue dots) generally precede the P2 main shock and aftershocks (red dots), and the P0 southern aftershocks (red dots). Vertical dashed lines mark P0 (blue, left) and P2 (red), and initiation time of earthquake cluster following P0 (blue, right). (D) Stacked fractional change in interstation travel times with 1 SD; red shaded area indicates time period where travel time decrease may have initiated. (E) Volumetric strain at extensometer KTA.

Table 1

Hypocenter information for earthquakes P0, P1, and P2 (Fig. 1).

JMA, Japan Meteorological Agency; NEIC, National Earthquake Information Center; UTC, Coordinated Universal Time; AQUA, Japan Focal Mechanism Catalog System.

	JMAtime(UTC)	JMAlatitude(degrees)	JMAlongitude(degrees)	JMA(NEIC)magnitude	JMA(AQUA)depth (km)
P0	2012-04-1211:19:57	36.829	141.343	5.6 (5.4)	28.8 (5.0)
P1	2012-04-1214:50:50	37.452	141.734	5.9 (5.5)	26.8 (13)
P2	2012-04-1310:10:00	36.947	141.422	6.0 (5.7)	31.5 (8.0)

Because many small noncataloged events can be missed, a far more sensitive measure of localized seismicity can be obtained by measuring the high-frequency (5 to 10 Hz) coherence of seismic energy across nearby Hi-net () borehole stations (Fig. 1B). Most energy in this band is necessarily produced locally (within tens of kilometers) because seismic energy from more distant sources is highly attenuated. Coherence is a robust and elegant correlation metric in that it detects all seismic emission regardless of its waveform, including traditional earthquakes, tremor, microseismicity, and clusters of small events that overlap in space and time (movies S1 and S2).

We find that seismic coherence begins increasing in the geographic region surrounding P1 and P2 immediately after the passing of seismic waves from the IOE, preceding both event sequences, and continues at an elevated level for about 10 days (Fig. 2B). Coherence increases in the north and then migrates south (fig. S1), in agreement with cataloged seismicity. Thus, seismic coherence along with the earthquake catalog establishes a spatial and temporal correlation between the IOE and local seismicity in the 10 days after the IOE.

The third largest event in the seismicity after the IOE, labeled “P0” (Mw 5.4), initiated before the observed spatial-temporal trend predicts seismicity, indicating that it was triggered independently of the observed northeast-to-southwest trend in seismicity. However, it has an unusual aftershock sequence that connects it to our other observations. If we plot the cumulative earthquake count for the aftershock sequences of P0 and P2, we see that they both exhibit a typically sharp increase in seismicity immediately after their respective main shocks (Fig. 2C). However, P0’s aftershock sequence abruptly terminates, and a new cluster of events begins to emerge to the south at the time the spatial-temporal trend predicts that seismicity should increase (Fig. 2C, inset). This new cluster of events is not itself an aftershock sequence of a new event because its largest event is not at the sequence beginning. Migrating seismicity and seismicity near P0 can be explained by Coulomb stresses near the leading edge of a slip front migrating from the northeast to southwest, consistent with our seismicity and coherence observations.

Dynamic strain

Nonlinear elastic effects induced by dynamic strains exhibit amplitude and wave duration dependence (, –). Peak dynamic strains during Love waves generated by the IOE are 50 μstrain at a period of 35 s measured at F-net station HRO (Fig. 1B) (). P2 generates peak strains of 1.5 μstrain at a period of 12 s. Surface wave amplitudes remain above 10 μstrain during the IOE and its Mw 8.2 aftershock for ~190 s. In contrast, there are three magnitude 5+ events in the 10 days before the IOE in Japan; the closest spatially, a Mw 5.7 deeper event (52 km) on the plate interface 63 km from P1 that occurred on 1 April, produced dynamic strains similar to P2. The other two smaller events occurred on 2 and 3 April and are 127 and 322 km away from P1, respectively. Strain amplitudes are at least 30 times larger from the IOE than any local event, and thus, nonlinear effects that should depend nearly linearly on the dynamic strain amplitude () from the IOE should far exceed those from local events. Wave durations are also much longer from the IOE and may “condition” the crust, changing the elasticity further during wave excitation, in the manner commonly observed in laboratory studies (, ) and suggested in a simple model for earthquake triggering ().

Strain

Coincident in time with the arrival of waves from the IOE and the increase in rate of microseismicity observed with waveform coherence is a change in slope in the volumetric strain time series measured at station KTA (Figs. 1B and 2E). After a period of persistent extension after the TOE due to afterslip, this region of Honshu Island is in a period of transition from extension back to compression during 2012 (Fig. 3), the predominant behavior in the interseismic period. During this transition, strain fluctuates between volume increase and decrease. There are no other strain meters in our study area. Strain meters to the north do not show a transition from positive to negative volume change coincident with the arrival of waves from the IOE, indicating that this signal is not simply a response of the instrument to the seismic waves, but instead indicates a regional strain signal. The purpose of Fig. 3 is not to show a signal associated with the IOE because the data are highly smoothed in both time and space compared to strain measured at station KTA (Fig. 2E), but rather to show that the regional strain rate in the direction of plate convergence is neutral (zero slope) at the time of the IOE. Therefore, the strain field produced by local phenomena, such as the observed cluster of earthquakes, can force the strain rate to go positive or negative without being overprinted or dominated by the direction of the regional strain field.

Fig. 3

Geodetic.

(A) Two sets of GPS stations (blue and green). The mean position of the blue stations is subtracted from mean position of the green stations along the direction of plate convergence indicated by the black arrow. (B) Times series of the difference. The red line indicates the IOE. The two sets of stations are rapidly diverging after the TOE, but long period strain rate is near zero at the time of the IOE and slightly negative by the end of 2012. During 2012, the strain alternates between contraction and dilation.

Seismic velocity

We next characterize seismic velocities to determine whether there is a bulk crustal change in elasticity at the time of the IOE (, ). Our measurements cover the onshore crust where the Hi-net seismic array is located (Fig. 1B). We use ambient noise in the passband 0.1 to 1 Hz recorded at 113 three-component stations over 58 days, and use 10-day stacks and all nine component pairs to calculate temporal changes in seismic travel times. Each point on the curve represents a stack of the previous 10 days in the figure (Fig. 2D). We find that normalized travel times δt/t are reduced by up to 10−4 after the IOE, corresponding to increased velocities. The reduction in travel times only becomes apparent 4 days after the IOE; however, because of the moving-window averaging procedure, it is not possible to precisely discern when the velocity change begins within the 4-day interval after the IOE. Regardless, the minimum travel time is first observed ~10 days after the IOE, indicating that all of the velocity changes took place at the time of, and/or soon after, the IOE. The perturbation persists for about 3 weeks. No simultaneous, statistically significant change in seismic coda-wave attenuation was observed.

These seismicity, strain, and seismic velocity observations cover two adjacent geological regions: the shallow accretionary wedge offshore and the shallow crust onshore, suggesting a widespread elastic effect (Fig. 1B). All observations initiate coincident with or shortly after the passing of waves from the IOE but have different slow dynamical recovery times (2) (note that the term “slow dynamics” refers to the recovery process independent of the sign of the elastic change). The increase in seismicity and seismic coherence recovers over ~14 days, the travel times recover over ~21 days, and the recovery time of the strain is unknown because of the unavailability of data. Thus, it appears that the seismic velocities, seismicity, and seismic coherence reflect recovery processes that are different and short-lived.

DISCUSSION

The most plausible explanation for the timing of the increase in seismicity and seismic coherence is a seismic wave–induced weakening of the normal faults in the shallow accretionary wedge. Because most of the triggered seismicity is delayed, the mechanism is not simply instantaneous Coulomb failure. Other explanations include random coincidence, local interactions between the magnitude 5+ events and nearby clusters, changes in loading rate due to afterslip on the plate interface, and some other process that is affecting seismicity locally such as increased pore pressures. With the exception of random coincidence, which is not supported by our statistical analysis (see Materials and Methods), none of the explanations exclude a triggering role for the dynamic stresses from the IOE.

All nine component combinations for these three-component instruments (ZZ, ZT, ZR, TT, TZ. TR, RR, RZ, and RT) show a reduction in travel times, which is likely dominated by surface waves traveling between stations. We conclude that the velocity increase is widespread throughout the uppermost crust. We justify this statement based on the following evidence. The result shown in Fig. 2D represents the full frequency range of 0.1 to 1 Hz. As a proxy for the depth of the velocity change, we tested two frequency bands: 0.1 to 0.2 Hz (mid-crust) and 0.2 to 1 Hz (shallow crust). The reduction in travel time is present but weaker for the frequency range of 0.1 to 0.2 Hz. On the basis of Rayleigh dispersion characteristics over these frequency intervals, we conclude that the travel time reduction is strong in the upper 5 km and weak but present below 5 km. Further evidence suggests the observations cannot be solely a result of near surface influences because all Hi-net stations are in boreholes at least 100 m in depth. Additionally, all previous observations of near surface effects exhibit a velocity decrease (, ).

The sense of velocity change (positive) is in contrast to most laboratory (, ) and field observations (negative) (). We posit that the observed velocity increase is due to a closing of cracks and fractures in the upper crust as a result of negative strain (contraction) conditions induced by IOE-triggered faulting. The increase in normal faulting events occurring offshore produces a negative strain perturbation onshore, which can be observed on strain meter KTA (Figs. 1B and 2E). The contraction after the IOE precedes the larger cataloged offshore events, but its initiation temporally corresponds with the burst of coherence shown in Fig. 2B, suggesting widespread microseismic normal faulting offshore during this time. Stress is transferred from the plate interface to the shallow, offshore fore arc quasi-continuously during the afterslip of the TOE. Stress is periodically transferred to the interior of Honshu Island by normal faulting offshore and, in this case, is due to a triggered burst of seismicity (Fig. 4).

Fig. 4

Interpretation.

(A to C) Cycle of stress transfer from (A) plate interface to shallow fore arc (B) resulting in normal faulting earthquakes and (C) stress transfer to interior Honshu Island. Blue triangles represent Hi-net stations. Red stars are normal faulting earthquakes.

The volumetric strain observed at station KTA continues to trend negative after velocities are recovering and may exhibit a different recovery time (data do not exist to test this concept). We note that the observations at station KTA are a point measurement and the velocity observations represent an average of the area beneath the stations used to calculate them (Fig. 1B). We expect that station KTA, being close to normal faulting earthquakes offshore, represents an upper bound in magnitude and duration for strains occurring in the area where increased velocities are observed. Nevertheless, we can use the observed strains to estimate the sensitivity of velocities to stress using the strain perturbation observed at the time of the highest observed velocities. The velocity sensitivity is calculated by dividing the fractional change in velocity (8.0 × 10−5) (the negative of the fractional change in travel time) by the stress, obtained from the volumetric strain (3.6 × 10−8), multiplied by the bulk modulus (5.0 × 104 MPa).

These calculations highly depend on how sensitively velocity variations are measured.

The velocity sensitivity is 4.4 × 10−2 MPa, which is higher than that observed during a slow slip event in Mexico (7 × 10−3 MPa) () but lower than that observed at SAFOD (San Andreas Fault Observatory at Depth) (2.4 × 10−1 MPa) (), and due to tidal deformation in California (0.5 MPa) ().

Over the past 20 years and most notably in the last 10 years, seismologists have focused on earthquake interaction manifest only by dynamic earthquake (or tremor/slow slip) triggering (, –). The overarching hypothesis supported by this study is that broad regions of Earth’s crust may be perturbed by dynamic stresses, leading to cascading elastic effects where crustal material is forced into a metastable state, followed by a slow dynamical (, ) recovery to either the original or a new equilibrium state. Strain focusing may occur in fault zones where effects are expected to be larger than in the surrounding bulk rock mass as observed in laboratory studies (). Thus, the influence of a perturbation may be highly heterogeneous and favor highly damaged regions associated with faults.

To summarize, we posit that IOE surface waves trigger a migrating elastic disturbance from north to south. The disturbance is manifest by large numbers of small, triggered events as well as several large earthquakes. The triggering of dominantly normal faults transfers stress, leading to contraction beneath Honshu, causing a widespread increase in velocity. The recovery rates of all of these phenomena after the disturbance vary, but none are observed to last more than several weeks.

Many of the phenomena we observe bear strong resemblance to laboratory observation in bulk and fractured materials (, , , ). We also note that similar effects are induced from quasi-static (and potentially dynamic) induced stress changes in aftershock regions of earthquakes (). To our knowledge, a broad, dynamically induced effect at a large distance from an earthquake source has not been observed before.

Our study region is in the aftershock zone of the TOE and is subject to ongoing aftershocks, afterslip, viscous relaxation, and other postseismic processes; thus, it is logical to expect that many components of this elastic system are highly damaged or in a critical state near failure, where small perturbations can activate observable changes. Characterizing the full elastic response to wave disturbance will ultimately affect our understanding of earthquake nucleation, triggered earthquakes, earthquake forecasting, induced seismicity associated with hydraulic fracturing, and seismic hazards, in addition to improving our understanding of the state and behavior of Earth’s elastic systems at all scales.

MATERIALS AND METHODS

Dynamic strain method

We calculated the dynamic strain from the IOE and P2 by dividing the transverse component of the velocity seismogram from F-net station HRO by the Love wave phase velocity and then multiplying it by the normalized Love wave displacement eigenfunction for a depth of 10 km. We used the velocity model of Nakamura et al. ().

Seismic coherence method

We obtained three-component data from Hi-net borehole stations between 9 March and 31 May 2012 (purple stations, Fig. 1B). We preprocessed the data by applying an anti-alias filter at 10 Hz, down-sampling the time series to 20 samples per second, and applying a 5-Hz high-pass filter. Then, we calculated the envelope function using a Hilbert transform.

Each station in our subnet that has at least three other stations within 100 km is processed as a “central” station. For each central station, we took 10-s templates on each component and padded them with 10 s of zeros on both sides. We then selected 30-s intervals on each component from nearby stations centered in time on the template from the central station. Then, we calculated the correlation coefficient with a maximum lag of ±10 s. We stacked the correlations for all station pairs and all components for the central station. To produce Fig. 2B, we stacked correlations for all central stations and then used a slot filter to remove frequencies around 1/day to remove a strong diurnal signal. Because of higher levels of anthropogenic, incoherent noise during the day, we systematically measured higher levels of seismic coherence during the nighttime.

Earthquake statistics

The observed increase in seismicity in the TOE aftershock region after the IOE was geographically widespread. To capture the geographic component of this increase in seismicity, we developed the following statistical method. We classified all earthquakes by their distance from latitude 31.32°N and longitude 134.10°E because a great circle path from this point to earthquake P1 roughly follows the plate interface in the vicinity of P1. After selecting a region that encompasses the entire TOE aftershock zone (500 km long), we binned all earthquakes by distance using a bin size of 25 km. Rather than having 20 bins in fixed positions, we put 300 bins in random locations along the plate interface to avoid biasing our results by a subjective selection of bin locations. At a given time horizon, we counted the number of bins where the rate for 6 days after the time horizon doubles over the rate for 20 days before the time horizon.

First, we used a time horizon equal to the hypocentral time of the IOE. We calculated 1 million realizations of 300 randomly placed bins and determined that, on average, 48.06 of the 300 bins have a doubling of their seismicity rate. Then, we selected 1 million time horizons between 21 days after the TOE to the end of 2013 and determined that, on average, 7.38 of 300 bins have a doubling of seismicity rate. In this distribution, only 1 in 358 of the realizations had 48.06 or more bins with a doubling of seismicity rate.

Strain method

Strain data from three-component extensometer KTA before 27 March 2012 and after 3 May 2012 were not of sufficient quality to use in this study. We manually fixed nontectonic offsets, outliers, and glitches from three strain components of extensometer KTA and accompanying barometric recordings. The tidal and barometric signals recorded on KTA were removed using the program cleanstrain+. To determine areal strain, we calculated the principal components by treating the extensometer as a rectangular rosette. A Poisson’s ratio of 0.3 was used to determine the volumetric strain from the areal strain.

Seismic velocity method

To calculate temporal variations in seismic velocities, we used the software package MSNoise (www.msnoise.org). A full description of the method and software is described by Lecocq et al. (), which we summarize here. There are three steps to the calculation of δv/v: (i) computing cross-correlation functions (CCFs) of ambient seismic noise time series at different dates for individual pairs of seismic sensors, (ii) measuring travel time delays of different arrivals (direct or coda waves) between these individual CCFs and a defined reference CCF, and (iii) averaging these travel time delays for different correlation lag times over different sensor pairs and interpreting these travel time delays using a simple model of uniform relative velocity change within the studied area (δv/v = constant). MSNoise uses moving-window cross-spectrum analysis () to estimate time delays in the CCFs. Uncertainties are estimated by determining the weighted mean and SD of time delays across all component pairs.

We stacked CCFs and time delays across 113 three-component stations using all nine component combinations in the frequency band of 0.1 to 1 Hz. Each point in Fig. 2D is the travel time delay estimated from a 10-day stack of 56,952 component pairs relative to the initial 10-day recording period. The initial 10-day recording period is shown with zero time delay in Fig. 2D. Our robust observation of increased velocities is revealed only through this massive stacking procedure.