Skip to main content

Ground motions characterized by a multi-scale heterogeneous earthquake model

Abstract

We have carried out numerical simulations of seismic ground motion radiating from a mega-earthquake whose rupture process is governed by a multi-scale heterogeneous distribution of fracture energy. The observed complexity of the Mw 9.0 2011 Tohoku-Oki earthquake can be explained by such heterogeneities with fractal patches (size and number), even without introducing any heterogeneity in the stress state. In our model, scale dependency in fracture energy (i.e., the slip-weakening distance Dc) on patch size is essential. Our results indicate that wave radiation is generally governed by the largest patch at each moment and that the contribution from small patches is minor. We then conducted parametric studies on the frictional parameters of peak (τp) and residual (τr) friction to produce the case where the effect of the small patches is evident during the progress of the main rupture. We found that heterogeneity in τr has a greater influence on the ground motions than does heterogeneity in τp. As such, local heterogeneity in the static stress drop (Δτ) influences the rupture process more than that in the stress excess (Δτexcess). The effect of small patches is particularly evident when these are almost geometrically isolated and not simultaneously involved in the rupture of larger patches. In other cases, the wave radiation from small patches is probably hidden by the major contributions from large patches. Small patches may play a role in strong motion generation areas with low τr (high Δτ), particularly during slow average rupture propagation. This effect can be identified from the differences in the spatial distributions of peak ground velocities for different frequency ranges.

Background

Near-field ground motion is strongly affected by the heterogeneity of earthquake rupture processes, such as multiple ruptures, rupture directivity, and the acceleration and deceleration of the rupture front. To quantitatively estimate the ground motion for seismic hazard assessment, characterization of this heterogeneity is essential. In kinematic earthquake models, an assumed target fault plane is usually decomposed by many rectangular sub-faults on which the spatio-temporal slip distribution is determined by a waveform inversion of seismic data (e.g., Hartzell and Heaton 1983). Numerous slip models have been presented for the Mw 9.0 2011 Tohoku-Oki earthquake (e.g., Ide et al. 2011; Koketsu et al. 2011; Simons et al. 2011; Suzuki et al. 2011). In these slip models, model parameters (e.g., total slip, rupture time, and rise time) have been attributed to each sub-fault.

The seismic waves inverted for slip models are usually displacements or velocities filtered in a low-frequency band, which is typically <1 Hz. There has been debate over whether the source of high-frequency waves is identical to the large slip areas in slip models (Hartzell et al. 1996; Kakehi et al. 1996; Nakahara et al. 1999). Although a general consensus has not yet been reached (e.g., Nakahara 2008), for many subduction earthquakes, high-frequency sources tend to be located deeper than the large slip areas (e.g., Lay et al. 2012). High-frequency sources are typically identified as several small isolated regions referred to as 'strong motion generation areas’ (SMGAs) (e.g., Kamae et al. 1998; Irikura and Miyake 2011). Some SMGA models have also been presented for the Tohoku-Oki earthquake (Kurahashi and Irikura 2011; Asano and Iwata 2012; Satoh 2012). Given that the number of parameters is much smaller when characterizing the kinematic earthquake rupture process in terms of SMGAs than when using a finely discretized finite-fault rupture model, the approach using SMGAs may be a more robust approach for ground motion simulations.

From a dynamic point of view, Ide and Aochi (2005) proposed that the heterogeneity in earthquake faults can be expressed as a superposition of different scales of heterogeneity. This is termed 'multi-scale heterogeneity’ and a seismic source can be represented as cascading ruptures of many circular patches of different sizes on a crack plane. In terms of kinematic finite source description, such heterogeneity may correspond to a class of composite models developed by Frankel (1991) and Zeng et al. (1994), and similar to those of Bernard et al. (1996). Aochi and Ide (2011) and Ide and Aochi (2013) attempted to characterize the Mw 9.0 2011 Tohoku-Oki earthquake based on this source expression and were able to dynamically simulate the growth of the rupture process from a small rupture through to the Mw 9 event. In these studies, a large elliptical patch corresponding to the large slip area of a Mw 9 event, which was unknown prior to this earthquake, was key in explaining the gross seismic moment release, whereas the distributed small circular patches governed the rupture growth before triggering the large elliptical patch. These small patches were introduced to reflect moderate earthquakes that had occurred previously in this region.

The apparent similarity in both the timing and locations of SMGAs (Kurahashi and Irikura 2011; Asano and Iwata 2012) and the circular patches of Ide and Aochi (2013) suggest that such small-scale heterogeneity corresponds to SMGAs. The dynamic models explored by Aochi and Ide (2011) and Ide and Aochi (2013) were neither constrained nor calibrated directly using observed data. Moreover, the radiation of ground motion based on this multi-scale heterogeneity has not yet been discussed. The purposes of this study are to show the characteristics of ground motions from the simulated fractal patch model for the 2011 Tohoku-Oki earthquake and to discuss the similarities (or differences) with respect to other synthetic heterogeneity models of the earthquake.

Methods

Model and method - reference model

The scaling issue in fracture energy or the slip-weakening distance of a slip-weakening relationship is considered to be an intrinsic feature of the fractal interface (Matsu’ura et al. 1992; Ohnaka 2003; Ide and Aochi 2005). However, there are other interpretations, such as an apparent effect from marked heterogeneities in fault strength (Campillo et al. 2001) or a posteriori obtained features caused by an increase in fault damage (Andrews 2005) or thermal pressurization (Wibberley and Shimamoto 2005). The multi-scale heterogeneity model on which this study is based has been proposed based on the scaling relationship of the earthquake source (i.e., fracture energy is proportional to earthquake size) and on fractal features (Ide and Aochi 2005). It has been proposed that the 2011 Tohoku-Oki earthquake was governed by a large patch of fracture energy, and to rupture on such a big patch, the rupture must have initiated on small, connected patches, including the foreshock area (Aochi and Ide 2011; Ide and Aochi 2013). Most of the small and medium patches are assumed to be located on deeper parts of the subduction interface, according to the pattern of seismicity over the last century. We first used this model (Figure six of Ide and Aochi 2013), which is shown in Figure 1, as a reference. This simulation was based on the 3D boundary integral equation method (BIEM) in a homogeneous medium with the free surface approximation.

Figure 1
figure 1

Snapshots of the dynamic rupture simulation originally shown in Figure six of Ide and Aochi (2013). The earthquake starts gradually at depth to the west and ruptures an asperity offshore of Miyagi. The rupture then turns to the east, which was the area of the foreshock of 9 March 2011, and accelerated towards the trench. Finally, the rupture propagates laterally to both the north and south. At that stage, moderate asperities located in the south (i.e., off Fukushima and Ibaraki) are ruptured.

We calculated seismic wave propagation using a finite difference method (Aochi and Madariaga 2003; Dupros et al. 2008) and introduced dynamic rupture scenarios. This is a sequential simulation procedure (Aochi and Madariaga 2003; Aochi and Douglas 2006). We focus our discussion on relatively low frequencies and hence adopted a 1D layered structure (Table 1) as used in the routine moment tensor inversion of F-net (National Research Institute for Earth Science and Disaster Prevention, Japan). The dynamically simulated earthquake source was implemented as a series of point sources through the tensors of seismic moment rate. The grid size (Δs) of the finite difference scheme was 500 m, and the time step of the calculations was 0.02 s. The maximum reliable frequency from the simulation was estimated as fmax = Vmin/(5Δs) = 1.3 Hz, where Vmin is the minimum wave velocity in the medium. A snapshot of the ground motion on Earth’s surface is shown in Figure 2, and seismograms at five selected stations are shown in Figure 3. The E-W component is the most dominant of the three components in the given geometric configuration, and these five stations (one from each prefecture) adequately represent the characteristics of the wave field along strike. At the nearest station (MYG003), we observed two principal wave packets at ca. 50 and 120 s, which are confirmed in the snapshots at 60 and 120 s. The first radiation at 60 s in Figure 2 is related to the initial phase of the rupture process observed at 30 s (Figure 1). It is also evident from Figure 2 that the later wave phases (120 s in Figure 2) appear to originate from the shallower part of the fault plane near the trench axis. The wave fronts that are visually identified on Figure 2 are limited. The first wave packet is due to the initial rupture propagating to the west (30 s in Figure 1), whereas the second wave packet is associated with the rupture break at the shallowest part of the fault (60 s in Figure 1). The bilateral propagation of the rupture (>70 s in Figure 1) then radiates a wave packet. In addition, some locally radiated wave fronts are observed from off the Fukushima coast (100 and 120 s in Figure 2), but the phases of such waves are difficult to identify in the waveforms. In the spectra of the waveforms (Figure 3), there is an observed lack of radiation energy at ca. 1 Hz. This implies that the rupture process in this situation is governed mainly at frequencies lower than 1 Hz.

Table 1 Crustal structure used in the ground motion simulations
Figure 2
figure 2

Ground motion on Earth’s surface. Snapshots of wave propagation on the ground surface simulated in a 1D structure from the source model of Ide and Aochi (2013). The three components are shown every 20 s. The epicenter is shown by a star, and the fault area (Figure 1) is the rectangle. The waveforms at the five stations (triangle symbols) are shown in Figure 3.

Figure 3
figure 3

Seismograms at five selected stations. Comparison of observed and simulated ground motions (Figure 2) for the stations that are shown in Figures 1 and 2. The E-W velocity components are compared. The observation times (top panels) are aligned from 14:46:00 (JST) for a duration of 250 s that begin at 16 s. The origin time in the simulation was set to zero. The velocity spectra are shown for a frequency range between 0.01 and 2 Hz.

However, according to the SMGAs reported for the 2011 Tohoku-Oki earthquake (Kurahashi and Irikura 2011; Asano and Iwata 2012; Satoh 2012), it is expected that it should be possible to identify certain phases with different frequency contexts. One possible explanation for this is that the difference is due to heterogeneity in stress conditions as our dynamic concept varies fracture energy (keeping stress conditions uniform), whereas the SMGA includes stress heterogeneity. From the point of view of the slip-weakening relationship for fault friction, both approaches imply that small patches may 'locally’ radiate more seismic energy compared with large patches. This study further explores this through parametric studies.

Parametric studies

We simplified our heterogeneous description to carry out parametric studies. We assumed that one big and a few small patches are aligned at the same depth (Figure 4a). The difference in patch dimensions varies fourfold, including a large elliptical patch equivalent to a circle of 100-km radius (Table 2). The rupture criterion and friction are described by a linear slip-weakening relationship used by Aochi and Ide (2011) and Ide and Aochi (2013), which includes three parameters, namely, slip-weakening distance Dc, peak strength τp, and residual stress τr:

τ Δu = τ r + Δσ 1 - Δu / D c for Δu < D c τ r otherwise ,
(1)
Figure 4
figure 4

Geometry of the parametric study. (a) Fault heterogeneity was defined by one large ellipse with high fracture energy and four small circles with low fracture energy. The hypocenter was set at (0, 0) on the figure. Around the hypocenter and in order to avoid the rapid onset of the initial crack, it was assumed that fracture energy increases linearly with hypocentral distance. The difference in peak strength is three times in such a case (case 2 in Figure 5) and was then varied in the parametric study. (b) Wave propagation was calculated for the same geometry as in Figure 2, with the same mechanism (strike N 205°E, dip 15°, and rake 90°) and hypocenter. However, the patch distribution is for synthetic studies, and so we renamed the five stations as R1 to R5 and focus on the E-W (X) component as it is the most significant. The panel shows a snapshot of the simulation corresponding to case (a).

Table 2 Model parameters of the parametric studies

where the dynamic strength drop Δσ is defined as Δσ ≡ τp - τr. The fracture energy Gc is defined as Gc ≡ Δσ × Dc/2. We assumed that Dc is proportional to the patch dimension (e.g., Ide and Aochi 2005). In this case, Dc is 3.2 and 0.8 m for a large and small patch, respectively. To initiate rupture on the large ellipse having a large Dc, we assumed that Dc also scales with the distance from the hypocenter around the initiation point and produces a small (but large enough with respect to the simulation element size of 2 km) finite initial crack (radius of 10 km) located at the coordinate origin in Figure 4. This helps to minimize an artificial initial phase of the dynamic rupture process (e.g., Aochi and Ide 2004; Ide and Aochi 2005; Aochi and Douglas 2006).

Given that we did not aim to discuss the effect of the ground surface on the rupture process, all the simulations were carried out in a 3D infinite, homogeneous elastic medium (P and S wave velocities of 6,000 and 3,464 m/s, respectively, and material density = 2,500 kg/m3) using a BIEM (Fukuyama and Madariaga 1995, 1998). In fact, as observed in Figures 2 and 3, the seismic wave radiation originating from the rupture arrival on the ground surface is obvious since the dynamic rupture significantly interacts with the ground surface. This masks the other radiation effects and may be exaggerated, as no shallow, slow layer was included (e.g., Goto et al. 2012; Ide and Aochi 2013). To better identify the radiation from the deep part of the fault as proposed in SMGA models, we ignored the interaction of the dynamic rupture with the ground surface.

We described the derivative stress with respect to the initial stress level (assumed to be zero, i.e., τ0 = 0) rather than to the absolute stress level. Therefore, instead of the two parameters τp and τr we refer to the stress excess Δτexcess(≡τp - τ0) and the static stress drop Δτ(≡τ0 - τr). A large Δτexcess may lead to a delay in the rupture onset. A large Δτ may produce a large fault slip. We firstly kept Δτ constant and changed Δτexcess (case 1 in Table 2). We then varied Δτexcess (case 2), and finally, we changed both parameters simultaneously (case 3). The ground motions were compared at several stations along the fault strike (the same locations as in Figure 1), assuming the same hypocenter location and fault geometry as for the 2011 Tohoku-Oki earthquake. As the focus of this study is the wave radiation from the causal source, we used the simple 1D structure given in Table 1.

Results

Figure 5 shows snapshots of the dynamic rupture propagation for the three cases (Table 2). In all cases, the dynamic strength drop Δσ on the small patches was set to 20 MPa and the fracture energy Gc was defined from the constitutive relationship of Equation 1. The rupture process on the small patches is significantly different in case 2, in which the southernmost isolated patch is triggered. For the other three small patches superimposed on the large patch, it is observed that the slip rate and fault slip become larger with increasing stress drop Δτ (case 2). No significant anomaly (rupture delay or large slip) on the small patches is observable in case 1, in which we varied only Δτexcess. Case 3 shows intermediate features between the other cases.

Figure 5
figure 5

Snapshots of dynamic rupture propagation for the three cases listed in Table1. The parameters used in the figure (case 1) are written in italics in Table 2. The frictional parameters vary for the small patches to maintain σ = (τp - τr) = 20 MPa. Case 1: Δτexcess for the small patches is higher. Case 2: Δτ for the small patches is higher. Case 3: both Δτexcess and Δτ are higher. The final fault slip is shown in the right column.

We now consider in detail how the value of Δτexcess of the small patches affects the rupture process. Given that little difference is observed for Δτexcess = 15 MPa (case 1 of Figure 5), we increased this value. When we assume Δτexcess = 17.5 MPa (3.5 times that of the large patch), the small patches do not rupture. This can lead to a different scenario, although this is not relevant to this study in which we expect the small patches ruptured. Studying more strictly the parameter range, we find that the small patches rupture when Δτexcess = 15.5 MPa, but not when Δτexcess = 16.0 MPa. Figure 6 shows the rupture times on the fault plane. In Figure 6b, the rupture onset on the small patches should have been most influenced. In fact, there is some delay in the isochrones when the rupture arrives (e.g., at 20 and 50 s). However, once the rupture starts, the rupture velocity is as rapid as in the other parts, and there is little impact on the subsequent isochrones.

Figure 6
figure 6

Rupture times on the fault plane. Rupture time contoured at every 5 s on the fault plane in case 1 for various values of Δτexcess for the small patches. (a) 15 MPa, (b) 15.5 MPa, (c) 16 MPa, and (d) 17.5 MPa. The rupture duration at each point is shown in color.

We next compared seismograms (Figure 7) for the stations defined in Figure 4b. The coordinates (X, Y, and Z) were taken along the east, north, and upward directions. We used the same station locations for the same hypocenter and fault orientation (strike N 201°E, dip 15°, and rake 90°) shown in Figure 2. However, the small patch distribution was no longer calibrated to the 2011 Tohoku-Oki earthquake, and we renamed the five receivers as R1 to R5. The same structural model (Table 1) was adopted so that the reliable maximum frequency fmax was 1.3 Hz. Of the scenarios simulated in Figure 5, the difference between cases 1 and 3 is the least significant. We observed that the changes made for case 2 have an influence at both high and low frequencies, particularly for R3 to R5. The amplification at high frequencies is significant in the early stages (ca. 70 s at R3 and ca. 100 s for R4), and this can also be recognized at low frequencies. This is because a small patch of high stress drop is ruptured simultaneously along with the ongoing rupture of the large patch, as inferred from the rupture time shown in Figure 6. However, at receiver R5, significant differences in high frequencies also appear in the later phase (140 s), as the furthest small patch slightly displaced from the large patch is ruptured, whereas no significant differences in low frequencies are found in this later phase. As in Figure 3, there is a lack of radiation energy at ca. 1 Hz (Figure 7). This reflects the limit of the dynamic rupture setting, as the smaller patches are represented by a dimension of 25 km and a Dc of 80 cm. This change appears at 0.2 to 0.5 Hz (case 2) and is due to the difference in radiation from the small patches. This influence is also confirmed by the distribution of peak ground velocities (PGVs) from case 2 (Figure 8). In Figure 8, we sampled the synthetic three-component seismograms at every 5 km on the ground surface and applied two different bandpass filters (0.02 to 0.1 Hz and 0.1 to 0.5 Hz). The PGV was then calculated for the X-component (the largest component), for the vector of the three components (X, Y, and Z), and for the geometric mean of the horizontal components (X and Y). The behavior of this isolated patch is similar to that of SMGAs. High PGV values are found at higher frequencies that correspond to both small patches and SMGAs.

Figure 7
figure 7

Ground motion comparison for the three scenarios presented in Figure5. The receiver locations are shown in Figure 4b. The E-W (X) component is shown for a bandpass between 0.1 and 0.5 Hz and for a band between 0.02 and 0.1 Hz. At the bottom of the figure, the Fourier spectrum is also shown for a frequency band of 0.01 and 2 Hz. A frequency of up to 1.3 Hz is numerically reliable, but there is no wave radiation around 1 Hz under the given parameters, and as such, we restrict our discussion to the results obtained at <0.5 Hz.

Figure 8
figure 8

PGVs from the simulation shown in case 2 of Figure7at two frequency ranges. (a) 0.02 to 0.1 Hz and (b) 0.1 to 0.5 Hz. From left to right, each panel shows a PGV value for the X-component, three-component vector, and geometric mean of horizontal components.

Therefore, we gleaned from these simulations that heterogeneity in the static stress drop Δτ (i.e., residual stress level τr) leads to more significant differences in ground motion than those produced by heterogeneity in the stress excess Δτexcess (or peak strength τp). Small patches are visible when they are isolated from a large patch, such that the small patches are in some way separated from the nearby ongoing rupture. Therefore, the locally identified SMGAs may be different from the main rupture front or characterized by a larger stress drop. A more complicated 'inhomogeneous’ fault model (e.g., Ripperger et al. 2008) might be more appropriate in such cases. Stochastically heterogeneous models (e.g., Ide and Aochi 2005) should be used in future ground motion estimation studies.

Discussion

The main effect of the multi-scale heterogeneity in Dc is on the rupture velocity, as well as the rupture directivity (Ide and Aochi 2013). Such small patches exert a role when the rupture starts and expands. It might be expected that small patches influence high-frequency generation sources even after the rupture has progressed. However, the dimension of the small patches is not significant compared with the rupture progressing on a large patch. Thus, small patches are difficult to recognize in the waveforms. By varying τp and τr for small patches, we found that it is difficult to create a delayed rupture on small patches by changing these values. We examined how a delayed rupture is possible on the smaller patches. Through increasing the Δτexcess of the large patch, we attempted to make the rupture slower so that stress accumulation on the small patches was also slower. Figure 9 shows snapshots of the simulations for Δτexcess = 5.5 MPa for the large patch and then assuming (a) Δτexcess = 16.280 MPa (2.96 times larger) for the small patches and (b) Δτexcess = 16.335 MPa (2.97 times larger). The three small patches are ruptured in the former case (a), whereas no small patches are ruptured in the latter case (b). In example (a), the first onset of the rupture of the small patches is delayed, but it propagates very rapidly once it has started. The waveforms are compared in Figure 10.

Figure 9
figure 9

Snapshots of dynamic rupture propagation for two other cases. Δτexcess for the large patch is 5.5 MPa, and Δτexcess for the small patches are (a) 16.280 MPa and (b) 16.335 MPa. These examples show the difficulty of controlling the rupture process for the small patches. In case (a), the rupture onset is slightly delayed for the small patches, but once rupture starts, it propagates very rapidly. In case (b), no rupture occurs.

Figure 10
figure 10

Ground motion comparison for scenarios (a) and (b) shown in Figure9. The E-W (X) component is shown for a bandpass between 0.1 and 0.5 Hz and for a band between 0.02 and 0.1 Hz. At the bottom of the figure, the Fourier spectrum is also shown for a frequency band of 0.01 and 2 Hz. For reference, case 1 from Figure 7 is also shown.

Our parametric study confirms the observations of Ide and Aochi (2005). We have determined that the role of small patches during rupture propagation makes the rupture front heterogeneous, not by delaying it, but through advancing it according to the low fracture energy given that the rupture propagation is governed by the energy balance. In our model, the rupture propagation velocity is close to the S wave speed, and it is difficult to identify the effect of further advancing the rupture front on the seismic wave radiation. Nevertheless, such an effect would be evident if the rupture velocity were slower, and in fact, we observed a slight difference between case 1 of Figure 5 and case (a) of Figure 9. However, the influence on the ground motion is limited (Figure 10). Ulrich and Aochi (2014) attempted to identify different sizes of patch by inversion, but any secondary small patches were difficult to identify whereas a large patch was able to be easily characterized. This is because the radiation from the small patches tends to be hidden by the waves from the greater area of a large patch. This effect should be important when the main rupture triggers remote small patches. In other cases, the rupture of small patches should be significantly delayed. However, we cannot rule out the possibility of delayed rupture of small patches for more realistic friction laws. For example, in the case when fault strength τp is time-dependent, the rupture may be delayed according to the accumulation of strain and some relaxation process. This is in effect the same mechanism as aftershock occurrence (e.g., Dieterich 1994; Helmstetter and Shaw 2009).

However, the differences in Δτ demonstrate the influence of small patches on ground motions. The mechanical features dynamically investigated in this study should be consistent with kinematic interpretations. In the original concept of SMGAs, Irikura and Miyake (2011) proposed a high stress drop on the SMGAs to distinguish it from background rupture propagation. This appears to make sense, given that Aochi and Dupros (2011) reconstructed a fault constitutive relationship from the SMGA source model (Irikura 2008) for the 2007 Mw 6.6 Niigata-Chuetsu-Oki earthquake in Japan and obtained a difference in Δτ (and Δτexcess) that was about double between SMGA area and the rest of the fault plane. The studies of Das and Aki (1977) and Mikumo and Miyatake (1978) have shown that the rupture process is very sensitive to stress parameters such as Δτ and Δτexcess (e.g., Madariaga and Olsen 2000; Gabriel et al. 2012). Even if our study were to take into account differences in fracture energy, the sensitivity to these parameters remains valid and coherent.

In the framework of our fractal patch models, it is more realistic that each patch is attributed different stress states on a probabilistic basis, mainly because each patch has a different history (e.g., Aochi and Ide 2009) and partly because of other uncertainties. In the case of the 2011 Tohoku-Oki earthquake, many small patches must have been involved during the rupture process, and only some of these should have had a marked Δτ that would allow them to be recognized as SMGAs. Figure 11 shows one simulation of rupture propagation in the case where we introduced a probabilistic distribution in Δτ for each patch. The distribution of Δτ is given by a Weibull distribution with a shape parameter of 2 and a scale parameter of 4. All the sets do not create the expected scenario (i.e., the foreshock remains the same and the main shock increases to the largest patch). In the illustrated example, some patches around the hypocenter were assumed to be homogeneous in order to assure rupture initiation. Compared with the homogeneous stress condition (Figure 1), we observed that the advancing rupture fronts on the small patches are more localized and can be identified. However, the stress condition is difficult to constrain prior to the earthquake. Therefore, to predict a probable rupture scenario, a probabilistic approach should be used in the dynamic rupture simulations.

Figure 11
figure 11

A simulation of dynamic rupture propagation introducing a randomly heterogeneous Δ τ (or τ r ) distribution of patches. Some patches around the hypocenter are kept uniform to assure that the initial rupture progresses. (a) Initial shear stress τ0, with the small circles corresponding to the initial rupture area of the foreshock and main shock. (b) Peak strength τp. (c) Residual stress τr. (d) Dc. (e) Snapshots of the dynamic rupture propagation. The hypocenter of the main shock is represented by a cross in the above panels and by a star at the bottom.

Conclusions

We have simulated the wave radiation process from multi-scale heterogeneous models of fracture energy for a mega-earthquake such as the 2011 Tohoku-Oki earthquake. Based on our previous research, we started by assuming heterogeneity only in the slip-weakening distance Dc under uniform stress conditions and frictional levels τp and τr (i.e., uniform stress excess Δτexcess and stress drop Δτ). However, in this case, once a rupture is initiated on a large patch, the influence of nearby small patches on the ground motions is insignificant and hidden by the dominant wave radiation from the large patch, which has a rupture propagation velocity similar to the S wave speed. We then introduced additional heterogeneity in Δτexcess and/or Δτ. Finally, we showed that heterogeneity in Δτ has a greater influence on ground motions than the same level of heterogeneity just in Δτexcess. In fact, small patches can be easily ruptured by the use of the multi-scale concept in fracture energy, but these are difficult to rupture behind a large patch even if high Δτexcess is assumed. Small patches exert an influence on SMGAs with low τr (i.e., those with a higher static stress drop Δτ). This effect can be identified from differences in the spatial distribution of peak ground velocities at different frequency ranges. Our models indicate that the distribution of large patches may be inferred from past seismicity when earthquakes recur frequently. However, uncertainties on the stress heterogeneity need to be treated stochastically. Our multi-scale heterogeneity earthquake model is consistent with inferred earthquake models, such as SMGA models, and may be a useful approach for probabilistically introducing stress state heterogeneity when predicting ground motion for quantitative seismic hazard studies.

References

  1. Andrews DJ: Rupture dynamics with energy loss outside the slip zone. J Geophys Res 2005., 110: B01307, doi:10.1029/2004JB003191 B01307

    Google Scholar 

  2. Aochi H, Douglas J: Testing the validity of simulated strong ground motion from the dynamic rupture of a fault system, by using empirical equations. Bull Earthq Eng 2006, 4: 211–229. 10.1007/s10518-006-0001-3

    Article  Google Scholar 

  3. Aochi H, Dupros F: MPI-OpenMP hybrid simulations using boundary integral equation and finite difference methods for earthquake dynamics and wave propagation: application to the 2007 Niigata Chuetsu-Oki earthquake (Mw6.6). Procedia Comput Sci 2011, 4: 1496–1505. doi:10.1016/j.procs.2011.04.162, 2011

    Article  Google Scholar 

  4. Aochi H, Ide S: Numerical study on multi-scaling earthquake rupture. Geophys Res Lett 2004. doi:10.1029/2003GL018708

    Google Scholar 

  5. Aochi H, Ide S: Complexity in earthquake sequences controlled by multi-scale heterogeneity in fault fracture energy. J Geophys Res 2009. doi:10.1029/2008JB006034

    Google Scholar 

  6. Aochi H, Ide S: Conceptual multi-scale dynamic rupture model for the 2011 Tohoku earthquake. Earth Planets Space 2011, 63: 761–765. doi:10.5047/eps2011.05.008

    Article  Google Scholar 

  7. Aochi H, Madariaga R: The 1999 Izmit, Turkey, earthquake: Non-planar fault structure, dynamic rupture process and strong ground motion. Bull Seism Soc Am 2003, 93: 1249–1266. doi:10.1785/0120020167

    Article  Google Scholar 

  8. Asano K, Iwata T: Source model for strong ground motion generation in the frequency range 0.1–10 Hz during the 2011 Tohoku earthquake. Earth Planets Space 2012, 64: 1111–1123.

    Article  Google Scholar 

  9. Bernard P, Herrero A, Berge C: Modeling directivity of heterogeneous earthquake ruptures. Bull Seism Soc Am 1996, 86: 1149–1160.

    Google Scholar 

  10. Campillo M, Favreau P, Ionescu IR, Voisin C: On the effective friction law of a heterogeneous fault. J Geophys Res 2001, 106: 16307–16322. 10.1029/2000JB900467

    Article  Google Scholar 

  11. Das S, Aki K: A numerical study of two-dimensional spontaneous rupture propagation. Geophys J R Astron Soc 1977, 50: 643–668. 10.1111/j.1365-246X.1977.tb01339.x

    Article  Google Scholar 

  12. Dieterich JH: A constitutive law for the rate of earthquake production and its application to earthquake clustering. J Geophys Res 1994, 99: 2601–2618. 10.1029/93JB02581

    Article  Google Scholar 

  13. Dupros F, Aochi H, Ducellier A, Komatitsch D, Roman J: Exploiting intensive multithreading for efficient simulation of seismic wave propagation. 11th International conference on computational science and engineering, Sao Paulo, 16–18 July 2008 2008. doi:10.1109/CSE.2008.51

    Google Scholar 

  14. Frankel A: High-frequency spectral falloff of earthquakes, fractal dimension of complex rupture, b value, and the scaling of strength on faults. J Geophys Res 1991, 96: 6291–6302. 10.1029/91JB00237

    Article  Google Scholar 

  15. Fukuyama E, Madariaga R: Integral equation method for a plane crack with arbitrary shape in 3D elastic media. Bull Seismol Soc Am 1995, 85: 614–628.

    Google Scholar 

  16. Fukuyama E, Madariaga R: Rupture dynamics of a planar fault in a 3D elastic medium: rate- and slip-weakening friction. Bull Seisim Soc Am 1998, 88: 1–17.

    Google Scholar 

  17. Gabriel AA, Ampuero J-P, Dalguer LA, Mai PM: The transition of dynamic rupture styles in elastic media under velocity-weakening friction. J Geophys Res 2012., 117: B09311, doi:10.1029/2012JB009468 B09311

    Google Scholar 

  18. Goto H, Yamamoto Y, Kita S: Dynamic rupture simulation of the 2011 off the Pacific coast of Tohoku earthquake: multi-event generation within dozens of seconds. Earth Planets Space 2012, 64: 1167–1175.

    Article  Google Scholar 

  19. Hartzell SH, Heaton TH: Inversion of strong ground motion and teleseismic waveform data for the fault rupture history of the 1979 Imperial Valley, California, earthquake. Bull Seism Soc Am 1983, 73: 1153–1184.

    Google Scholar 

  20. Hartzell S, Liu P, Mendoza C: The 1994 Northridge, California, earthquake: investigation of rupture velocity, risetime, and high-frequency radiation. J Geophys Res 1996, 101: 20091–20108. 10.1029/96JB01883

    Article  Google Scholar 

  21. Helmstetter A, Shaw BE: Afterslip and aftershocks in the rate-and-state friction law. J Geophys Res 2009., 114: B01308, doi:1029/2007JB005077 B01308

    Google Scholar 

  22. Ide S, Aochi H: Earthquakes as multiscale dynamic rupture with heterogeneous fracture surface energy. J Geophys Res 2005., 110: B11303, doi:10.1029/2004JB003591 B11303

    Google Scholar 

  23. Ide S, Aochi H: Historical seismicity and dynamic rupture process of the 2011 Tohoku-Oki earthquake. Tectonophysics 2013, 600: 1–13.

    Article  Google Scholar 

  24. Ide S, Baltay A, Beroza GC: Shallow dynamic overshoot and energetic deep rupture in the 2011 Mw 9.0 Tohoku-Oki earthquake. Science 2011, 332: 1426–1429. doi:10.1126/science.1207020

    Article  Google Scholar 

  25. Irikura K: Lesson from the 2007 Niigata-ken Chuetsu-oki earthquake on seismic safety for nuclear power plant. Bull Jap Asso Earthq Eng 2008, 7: 25–29. (in Japanese) (in Japanese)

    Google Scholar 

  26. Irikura K, Miyake H: Recipe for predicting strong ground motion from crustal earthquake scenarios. Pure Appl Geophys 2011, 168: 85–104. doi:10.1007/s00024–010–0150–9

    Article  Google Scholar 

  27. Kakehi Y, Irikura K, Hoshiba M: Estimation of high-frequency wave radiation areas on the fault plane of the 1995 Hyogo-ken Nanbu earthquake by the envelope inversion of acceleration seismograms. J Phys Earth 1996, 44: 505–517. 10.4294/jpe1952.44.505

    Article  Google Scholar 

  28. Kamae K, Irikura K, Pitarka A: A technique for simulating strong ground motion using hybrid Green’s function. Bull Seism Soc Am 1998, 88: 357–367.

    Google Scholar 

  29. Koketsu K, Yokota Y, Nishimura N, Yagi Y, Miyazaki S, Satake K, Fujii Y, Miyake H, Sakai S, Yamanaka Y, Okada T: A unified source model for the 2011 Tohoku earthquake. Earth Planet Sci Lett 2011, 310: 480–487. doi:10.1016/j.epsl.2011.09.009

    Article  Google Scholar 

  30. Kurahashi S, Irikura K: Source model for generating strong ground motions during the 2011 off the Pacific coast of Tohoku Earthquake. Earth Planets Space 2011, 63: 571–576. 10.5047/eps.2011.06.044

    Article  Google Scholar 

  31. Lay T, Kanamori H, Ammon CJ, Koper KD, Hutko AR, Ye L, Yue H, Rushing TM: Depth-varying rupture properties of subduction zone megathrust faults. J Geophys Res 2012., 117: B04311, doi:10.1029/2011JB009133 B04311

    Google Scholar 

  32. Madariaga R, Olsen KB: Criticality of rupture dynamics in 3-D. Pageoph 2000, 157: 1981–2001. 10.1007/PL00001071

    Article  Google Scholar 

  33. Matsu’ura M, Kataoka H, Shibazaki B: Slip-dependent friction law and nucleation process in earthquake rupture. Tectonophys 1992, 211: 135–148. 10.1016/0040-1951(92)90056-C

    Article  Google Scholar 

  34. Mikumo T, Miyatake T: Dynamical rupture process on a 3-dimensional fault with non-uniform frictions and near-field seismic-waves. Geophys J Royal Astro Soc 1978, 54: 417–438. 10.1111/j.1365-246X.1978.tb04267.x

    Article  Google Scholar 

  35. Nakahara H: Chapter 15 Seismogram envelope inversion for high-frequency seismic energy radiation from moderate-to-large earthquakes. Adv Geophys 2008, 50: 401–426. doi:10.1016/S0065–2687(08)00015–0

    Article  Google Scholar 

  36. Nakahara H, Sato H, Ohtake M, Nishimura T: Spatial distribution of high-frequency energy radiation on the fault of the 1995 Hyogo-Ken Nanbu, Japan, earthquake (Mw 6.9) on the basis of the seismogram envelope inversion. Bull Seism Soc Am 1999, 89: 22–35.

    Google Scholar 

  37. Ohnaka M: A constitutive scaling law and a unified comprehension for frictional slip failure, shear fracture of intact rock, and earthquake rupture. J Geophys Res 2003, 108(B2):2080. doi:10.1029/2000JB000123

    Article  Google Scholar 

  38. Ripperger J, Mai PM, Ampuero J-P: Variability of near-field ground motion from dynamic earthquake rupture simulations. Bull Seism Soc Am 2008, 98(3):1207–1228. doi:10.1785/0120070076

    Article  Google Scholar 

  39. Satoh T: Source modeling of the 2011 off the Pacific coast of Tohoku earthquake using empirical Green’s function method. J Struct Const Eng 2012, 77: 695–704. (in Japanese) (in Japanese) 10.3130/aijs.77.695

    Article  Google Scholar 

  40. Simons M, Minson SE, Sladen A, Ortega F, Jiang J, Owen SE, Meng L, Ampuero JP, Wei S, Chu R, Helmberger DV, Kanamori H, Hetland E, Moore AW, Webb FH: The 2011 magnitude 9.0 Tohoku-Oki earthquake: mosaicking the megathrust from seconds to centuries. Science 2011, 332: 1421. doi:10.1126/science.1206731

    Article  Google Scholar 

  41. Suzuki W, Aoi S, Sekiguchi H, Kunugi T: Rupture process of the 2011 Tohoku-Oki mega-thrust earthquake (M9.0) inverted from strong-motion data. Geophys Res Lett 2011, 38: L00G16. doi: 10.1029/2011GL049136

    Google Scholar 

  42. Ulrich T, Aochi H: Rapidness and robustness of finite-source inversion of the 2011 Mw9.0 Tohoku earthquake by an elliptical-patches method using continuous GPS and acceleration data. Pageoph 2014. doi:10.1007/s00024–014–0857–0

    Google Scholar 

  43. Wibberley CAJ, Shimamoto T: Earthquake slip weakening and asperities explained by thermal pressurization. Nature 2005, 436: 689–692. 10.1038/nature03901

    Article  Google Scholar 

  44. Zeng Y, Anderson JG, Guang Y: A composite source model for computing realistic synthetic strong ground motions. Geophys Res Lett 1994, 21: 725–728. 10.1029/94GL00367

    Article  Google Scholar 

Download references

Acknowledgements

We thank Hiroe Miyake and Nobuki Kame for inviting us to the working group on 'Dynamic Source Models for the Next Generation Ground Motion Prediction’ in 2013-2014, which was funded by the Earthquake Research Institute of the University of Tokyo. We also thank John Douglas for offering comments that improved this paper. The comments from Martin Mai and an anonymous reviewer were very helpful. This research began under the framework of the French-Japanese ANR-JST joint program DYNTOHOKU (2011-2013) and was then funded by the French national project S4 (Subduction: Slow and Standard Seismology, 2012-2014, ANR-2011-BS56-017) of the Agence National de la Recherche. Most of the calculations were made at the French national supercomputing center GENCI-CINES (grant 2013/2014-c46700). This work was also partially supported by JSPS KAKENHI (23244090) and MEXT KAKENHI (21107007) grants.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Hideo Aochi.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

HA carried out principally the numerical simulations. SI set up the parameter studies and brought the analysis scheme. Both authors read and approved the final manuscript.

Authors’ original submitted files for images

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 2.0 International License (https://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Aochi, H., Ide, S. Ground motions characterized by a multi-scale heterogeneous earthquake model. Earth Planet Sp 66, 42 (2014). https://doi.org/10.1186/1880-5981-66-42

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/1880-5981-66-42

Keywords