Source rupture process of the 2011 Fukushima-ken Hamadori earthquake: how did the two subparallel faults rupture?

Abstract

The 2011 Fukushima-ken Hamadori earthquake (M_W 6.6) occurred about a month after the 2011 Great Tohoku earthquake (M_W 9.0), and it is thought to have been induced by the 2011 Tohoku earthquake. After the 2011 Hamadori earthquake, two subparallel faults (the Itozawa and Yunodake faults) were identified by field surveys. The hypocenter was located nearby the Itozawa fault, and it is probable that the Itozawa fault ruptured before the Yunodake fault rupture. Here, we estimated the source rupture process of the 2011 Hamadori earthquake using a model with two subparallel faults based on strong motion data. The rupture starting point and rupture delay time of the Yunodake fault were determined based on Akaike’s Bayesian Information Criterion (ABIC). The results show that the Yunodake fault started to rupture from the northern deep point 4.5 s after the Itozawa fault started to rupture. The estimated slip distribution in the shallow part is consistent with the surface slip distribution identified by field surveys. Time-dependent Coulomb failure function changes (ΔCFF) were calculated using the stress change from the Itozawa fault rupture in order to evaluate the effect of the rupture on the Yunodake fault. The ΔCFF is positive at the rupture starting point of the Yunodake fault 4.5 s after the Itozawa fault started to rupture; therefore, it is concluded that during the 2011 Hamadori earthquake, the Yunodake fault rupture was triggered by the Itozawa fault rupture.

Findings

Introduction

The 2011 Fukushima-ken Hamadori earthquake (M_W 6.6) occurred in the southeastern part of Fukushima Prefecture, northeastern Japan (Figure 1a) at 17:16 Japan Standard Time (JST; 08:16 UTC) on 11 April 2011. In northeastern Japan, particularly in northeast Ibaraki and southeast Fukushima prefectures, seismicity associated with swarm-like, shallow normal faulting has increased as a result of the 2011 Great Tohoku earthquake (M_W 9.0) (Kato et al. 2013). The 2011 Hamadori earthquake was the largest event in this period. During a field survey, Toda and Tsutsumi (2013) identified two subparallel surface ruptures along the previously mapped Itozawa and Yunodake faults that exhibit a predominantly normal sense of slip. The hypocenter determined by the Japan Meteorological Agency (JMA) has been shown to be close to the Itozawa fault (Figure 1b), which implies that the Itozawa fault ruptured before the Yunodake fault. The difference between strikes of the Itozawa and Yunodake faults is approximately 25°, and they are not conjugated. Interferometric synthetic aperture radar (InSAR)-derived crustal deformation also shows that the Itozawa and Yunodake source faults are not conjugated (e.g., Fukushima et al. 2013; Kobayashi et al. 2013). There are some cases where parallel faults have ruptured during an earthquake (e.g., the 2008 Wenchuan earthquake, China (Tong et al. 2010)). Moreover, there are also some cases where continuous faults have ruptured during an earthquake (e.g., the 1992 Landers earthquake, California (Wald and Heaton 1994) and the 2002 Denali earthquake, Alaska (Asano et al. 2005)). However, no cases have been previously observed where subparallel and unconjugated faults have ruptured during an earthquake, such as the case during the 2011 Fukushima-ken Hamadori earthquake. Therefore, it is important to estimate the spatiotemporal slip distributions on the source faults with their complicated fault geometry in order to understand the relationship of the two fault ruptures.

Figure 1

Map showing the location of the study area. (a) The stars indicate the epicenters of the 2011 Tohoku earthquake (green star) and the 2011 Hamadori earthquake (yellow star). The slip distribution of the 2011 Tohoku earthquake by Kubo and Kakehi (2013) is also shown on the map. (b) The stars indicate the epicenters of the 2011 Hamadori earthquake (yellow star) and the aftershock (13:39 JST on 18 June 2011) (light blue) with the F-net moment tensor solution and the active faults (AIST 2007) (red and orange lines). The aftershock is used for testing velocity structure models. The white inverted triangles indicate the strong motion stations located within 50 km of the epicenter. The black inverted triangles are the stations used for the kinematic waveform inversion.

In previous studies (Hikima 2012; JMA 2012; Shiba and Noguchi 2012), the source rupture process of this event has been estimated using a two-fault model (the Itozawa and Yunodake faults) based on strong motion data. However, different slip distributions have been identified, particularly in relation to the Yunodake fault. It is considered likely that such differences were caused by variations in the authors’ assumptions used for the fault model, particularly in relation to the rupture starting point and rupture delay time of the Yunodake fault.

In this study, in order to obtain slip distributions on the two faults, we objectively determine not only a hyperparameter describing the strength of smoothing constraint but also a set of additional parameters such as rupture starting point, rupture delay time, and first-time window triggering velocity based on Akaike’s Bayesian information criterion (ABIC) (Akaike 1980), whereas previous source inversions determined only hyperparameters (e.g., Ide and Takeo 1997; Sekiguchi et al. 2000). We then calculate the stress field from the estimated moment release history of the Itozawa fault and get temporal changes in the Coulomb failure function (time-dependent ΔCFF) in order to discuss whether the Yunodake fault rupture was triggered by that of the Itozawa fault.

Estimation of the source rupture process using kinematic waveform inversion

Data

The observed ground motions are represented by the convolution of the source, the propagation path, and site effects. In order to use appropriate Green’s functions in kinematic waveform inversion analysis, it is necessary to check whether observed waveforms can be reproduced using a given velocity structure model. In this study, we selected strong motion stations based on a comparison between the theoretical and observed waveforms of an aftershock (13:39 JST on 18 June 2011, M_W 4.4).

The observed data from 26 strong motion stations of K-NET, KiK-net, F-net, and JMA, which are located within 50 km of the epicenter of the 2011 Hamadori earthquake (Figure 1b), were used. The original acceleration data from these stations (except F-net stations, which only provide original velocity data) were integrated into velocity; velocity seismograms were then band-pass filtered from 0.1 to 1.0 Hz. Assuming a one-dimensional velocity structure model for each station (extracted from the three-dimensional Japan Integrated Velocity Structure Model (Koketsu et al. 2012)), the theoretical waveforms of the aftershock were then calculated using the discrete wave number method (Bouchon 1981) and the reflection and transmission coefficient matrix method (Kennett and Kerry 1979). Theoretical waveforms were then compared with the observed waveforms (Figure 2), and we selected 18 stations (Figure 1b) where the observed waveforms were found to be well reproduced.For the kinematic waveform inversion, three components of the velocity seismograms recorded from the 18 selected strong motion stations were used (Figure 1b), and these were band-pass filtered from 0.1 to 1.0 Hz and resampled at 10 Hz. The dataset comprises 35 s of the time series from 1.0 s before the S-wave arrival.

Figure 2

Examples of a comparison between observed waveforms (black lines) and theoretical waveforms (red lines). R, T, and U are the radial, transverse, and up-down components, respectively. In the first row, the observed waveforms are well represented; therefore, this station was selected for the kinematic waveform inversion. In the second row, the observed waveforms are not well represented; therefore, this station was not selected.

Method

Green’s functions were calculated in the same way as the theoretical waveforms of the aftershock. We assumed a model with two planar faults (using the Itozawa and Yunodake faults) based on the source fault model inferred from the InSAR information analyzed by Fukushima et al. (2013). The setting of the assumed fault model is shown in Table 1. The Itozawa fault model is divided into 77 subfaults of 2 × 2 km² (11 along the strike and 7 along the dip). The Yunodake fault model is also divided into 63 subfaults of 2 × 2 km² (9 along the strike and 7 along the dip).

Table 1 Setting of the assumed fault model

We used the multiple time window linear source inversion method (Hartzell and Heaton 1983) with a spatiotemporal smoothing constraint on slips (Sekiguchi et al. 2000). The rake angle variations were limited within −90° ± 45° by the non-negative least squares problem (Lawson and Hanson 1974). The temporal moment release history at each subfault is expressed by a series of six smoothed ramp functions: each has a duration of 1.0 s and shifts from the preceding one by 0.5 s. We used 63 possible rupture starting points in the Yunodake fault (the center of all the subfaults) and eight cases of the rupture delay time for the Yunodake fault at 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, and 6.5 s. In addition, we used six first-time window triggering velocities of 2.04, 2.21, 2.38, 2.55, 2.72, and 2.89 km/s, which correspond to between 60% and 85% of the S-wave velocities in the velocity structure model at the rupture starting point of the Itozawa fault. By using the minimum ABIC condition, we selected the strength of smoothing and the best set of model parameters such as the rupture starting point, the rupture delay time of the Yunodake fault, and the first-time window triggering velocity. The observational equation and smoothing constraint are combined as

$$\left[\begin{array}{c}\hfill \mathbf{G}\hfill \\ \hfill \lambda \mathbf{S}\hfill \end{array}\right]\mathit{m}=\left[\begin{array}{c}\hfill \mathit{d}\hfill \\ \hfill 0\hfill \end{array}\right]\text{,}$$

where G is the M_G × N matrix of theoretical Green’s functions, λ is the strength of smoothing, S is the M_S × N matrix of the elements of the smoothing constraint equations, m is the model vector, and d is the vector of the velocity data obtained by integrating the original acceleration data. In this study, ABIC is defined as

$$\begin{array}{l}\mathrm{ABIC}=\left(M_{\mathrm{G}}+M_{\mathrm{S}}–N\right)log\left({\left|\mathit{d}-\mathbf{G}\mathit{m}\right|}^2+{\lambda }^2{\left|\mathbf{S}\mathit{m}\right|}^2\right)\\ -M_{\mathrm{S}}log{\lambda }^2+log\left|{\mathbf{G}}^{\mathbf{T}}\mathbf{G}+{\lambda }^2{\mathbf{S}}^{\mathbf{T}}\mathbf{S}\right|+C\text{,}\end{array}$$

where C is a constant term.

Results

The best location of the rupture starting point on the Yunodake fault is indicated by the blue star in Figure 3a, and the rupture delay time is selected as 4.5 s. The best speed of the first-time window triggering is 2.04 km/s (60% of the S-wave velocity at the rupture starting point of the Itozawa fault). In Figure 4, the synthetic waveforms at all stations used in the inversion are compared with the observed waveforms, and these closely reproduced the observed waveforms.

Figure 3

The spatial relationship of the source model. (a) The spatial relationship between the rupture starting points (red and blue stars) and the fault planes. (b, c) Final slip distributions on the map. The stars indicate the rupture starting point on each fault. (c) The Yunodake fault is superimposed on the Itozawa fault to enable easy imaging of the slip distribution on the Yunodake fault.

Figure 4

Comparison between observed velocity waveforms (black lines) and synthetic velocity waveforms (red lines). The numbers next to the waveforms are the maximum amplitude (cm/s) of each component for each station. EW, east-west; NS, north-south; UD, up-down.

The estimated slip model shows that a large slip area is located in the northwestern shallow part of the Itozawa fault (Figure 3b). The estimated large slip areas on the Yunodake fault are located in the shallow part of the fault plane and in the vicinity of the rupture starting point (Figure 3c). The average amount of slip in the northwestern part of the Itozawa fault that lies within a depth of 5 km from the earth’s surface is approximately 1.3 m, and the average amount of slip of the Yunodake fault that lies within a depth of 5 km from the earth’s surface is approximately 1.0 m. The total seismic moment of this event is estimated to be 1.45 × 10¹⁹ N · m (M_W 6.7), and the seismic moment released from the Itozawa fault is almost the same as that released from the Yunodake fault.

Toda and Tsutsumi (2013) found a maximum vertical surface slip of 2.1 m at Shionohira, which is situated northwest of the epicenter along the Itozawa fault. They also found that the shape of the slip distribution of the Itozawa fault was an asymmetric chevron tapering linearly toward the edge. In addition, they found a slip range of 0.1 to 0.9 m all along the Yunodake fault and a plateau-like slip distribution. Our estimated slip distributions of the two faults in the shallow part are consistent with these observed surface slip distributions, although the observed surface slip along the Yunodake fault (0.1 to 0.9 m) is slightly smaller than the estimated average amount of slip (1.0 m). However, this difference could be explained by the fact that the estimated large slip area is not located in the shallowest part, which would cause the observed surface slip to be slightly smaller than the estimated average amount of slip.The temporal progression of the slip is shown in Figure 5. The rupture can be seen to propagate mainly to the north-western up-dip direction on the Itozawa fault. On the Yunodake fault, it can be seen to first propagate mainly to the up-dip direction and then to the southeast in the shallow part. In addition, the results show that the Yunodake fault begins to rupture at the time when the Itozawa fault is rupturing.

Figure 5

Temporal slip distribution on the two fault planes. The stars indicate the rupture starting points of the two faults, and the figures at the bottom right represent the final slip distributions. The times noted on each figure represent the time elapsed from the start of the rupture of the Itozawa fault.

The effects of the Itozawa fault rupture on the Yunodake fault plane

Calculation of the time-dependent ΔCFF

In a previous study, Hikima (2012) calculated the static ΔCFF using his final slip model on the Itozawa fault plane and found an area with a positive ΔCFF that included his rupture starting point when he assumed the apparent friction coefficient to be 0.8. In our results, which were obtained from the kinematic waveform inversion, the Itozawa fault was found to be in the process of rupturing when the Yunodake fault started to rupture. Therefore, we consider it important to calculate the time-dependent ΔCFF using the obtained moment release history of the Itozawa fault in order to evaluate the effect of the Itozawa fault rupture on the Yunodake fault plane.

Method

The stress field was calculated by solving the elastodynamic wave equations with the finite difference method (FDM) using discontinuous grids from Ground Motion Simulator (GMS) (Aoi and Fujiwara 1999). The FDM model space was 60 × 60 km² in the horizontal direction and 30.1 km in the depth direction. The underground structure model for the stress field calculation was assumed to be a homogeneous half space, with V_P, V_S, density, and Q of 5.0 km/s, 2.9 km/s, 2.6 × 10³ kg/m³, and 680, respectively. The discontinuous grid system consists of two regions with different grid spacings. The grid spacing of the shallower region is small, whereas the grid spacing of the deeper region is three times coarser (Aoi and Fujiwara 1999). We used discontinuous grids that consisted of a grid spacing of 0.05 km in the shallower region up to the depth of 18.55 km (including the Itozawa and Yunodake faults) and of 0.15 km in the deeper region. The moment release history of the Itozawa fault plane is given by the spatially interpolated model of the obtained slip model using the kinematic waveform inversion with intervals of 2.0 to 0.4 km. The time-dependent ΔCFF is calculated according to ΔCFF = Δτ + μ′ Δσ (e.g., Toda et al. 2011), where τ is the shear stress, σ is the normal stress on the fault plane, and μ′ is the apparent friction coefficient. In this study, two cases of μ′ = 0.4 and 0.8 were tried. For shear stress τ, the positive direction is set to the slip direction of the first-time window at the rupture starting point of the Yunodake fault.

Results

For μ′ = 0.8, the time-dependent ΔCFF has a positive value of approximately 0.8 MPa at the rupture starting point of the Yunodake fault 4.5 s after the Itozawa fault begins to rupture (Figure 6a,b). Furthermore, for μ′ = 0.4, the time-dependent ΔCFF also has a positive value of approximately 0.3 MPa at the rupture starting point of the Yunodake fault, 4.5 s after the Itozawa fault begins to rupture (Figure 6c,d). Therefore, in consideration of such results, it is concluded that the Itozawa fault rupture triggered the Yunodake fault rupture.

Figure 6

The time-dependent ΔCFF. (a, c) The distribution of ΔCFF on the Yunodake fault plane at 4.5 s after the Itozawa fault begins to rupture. The stars indicate the rupture starting point of the Yunodake fault. (b, d) Time-dependent ΔCFF at the rupture stating point of the Yunodake fault. ‘0 s’ is the time at which the Itozawa fault starts to rupture. The top panels (a, b) are for μ′ = 0.8, and the bottom panels (c, d) are for μ′ = 0.4.

Discussion

In the previous section, it was demonstrated that the rupture of the Yunodake fault was caused by the dynamic effect of the rupture of the Itozawa fault. However, it is possible that the effects of other earthquakes such as the 2011 Tohoku earthquake could have impacted the Yunodake fault; accordingly, this possibility is also investigated.

Toda et al. (2011) calculated the static ΔCFF on known major faults and megathrusts using the source model of the 2011 Tohoku earthquake and the M_W 7.9 aftershock, and the results showed that the static ΔCFF on the fault near the Yunodake fault has a positive value of approximately 0.1 MPa when an apparent friction coefficient of 0.4 is adopted. However, the static ΔCFF on the Yunodake fault was not calculated in that study. Therefore, in order to discover whether the Itozawa fault rupture triggered the Yunodake fault rupture, we calculated the effect of the static ΔCFF of the 2011 Tohoku earthquake on the Yunodake fault; the static ΔCFF was calculated based on Okada (1992) using the final slip distribution of the 2011 Tohoku earthquake (Kubo and Kakehi 2013). The value of the static ΔCFF at the rupture starting point of the Yunodake fault was found to have a positive value of approximately 0.15 MPa, which is considerably smaller than the time-dependent ΔCFF of +0.8 MPa from the rupture of the Itozawa fault. Therefore, we consider this to indicate that the rupture of the Yunodake fault was triggered primarily by the rupture of the Itozawa fault.

Conclusions

The source rupture process of the 2011 Hamadori earthquake was estimated using fault model parameters such as the rupture starting point and the rupture delay time of the Yunodake fault. The results show that the Yunodake fault began to rupture from the northern deep point 4.5 s after the Itozawa fault began to rupture. The estimated final slip distributions are consistent with the surface slip distributions found by the field surveys.

The time-dependent ΔCFF was then calculated using the estimated moment release history on the Itozawa fault plane, and it was found to be positive at the rupture starting point of the Yunodake fault, 4.5 s after the Itozawa fault began to rupture. Furthermore, the obtained time-dependent ΔCFF was found to be larger than the static ΔCFF caused by the 2011 Tohoku earthquake, which implies that the effect of the Itozawa fault rupture on the Yunodake fault plane was larger than that of the 2011 Tohoku earthquake. Therefore, we conclude that during the Hamadori earthquake, the rupture of the Yunodake fault was triggered by the rupture of the Itozawa fault. Considering our results, it may be said that the rupturing of these unconjugated faults was due to the dynamic effect.