Visualization of Green’s function anomalies for megathrust source in Nankai Trough by reciprocity method

We used simulation by the reciprocity method to visualize the distribution of Green’s function amplitudes in the source of a megathrust earthquake in the Nankai Trough and considered the effects of various areas (asperities or strong motion generation areas) on the simulated long-period ground motions at Konohana in the Osaka basin. We employed a fault source model proposed for an anticipated M9-class event in the Nankai Trough and the 3D Japan Intergrated Velocity Structure Model developed for simulations of long-period ground motions in Japan. Green’s functions were calculated for about 1400 subsources by combining the finite-difference method and the reciprocity approach. Depths, strikes, and dips of subsources were adjusted to the shape of the upper boundary of the Philippine Sea plate. Ground motions with periods of 4–20 s were considered. The simulated distribution of peak amplitudes of Green’s functions identified two strongly anomalous areas: (1) a large along-strike elongated area just south of the Kii Peninsula and (2) a parallel area closer to the trench. The elongation of the anomalies corresponded well with depth isolines at the top of the Philippine Sea plate. Postulating that plate shape influences simulated ground motions, we investigated the effect on Green’s function amplitudes of phenomena related to plate shape: radiation pattern; variations of medium properties (e.g., velocity and density) at subsource depths; depth, strike, and dip; and the effect of soft sediments. We suggest that the cumulative effect on Green’s function amplitudes of subsource radiation patterns, medium properties at subsource depth, reflection from crustal interfaces, and passage through soft sedimentary layers plays a critical role in the formation of amplitude anomalies. Analysis of waveforms and the time delay of peak amplitude demonstrate that large-amplitude waves of Green’s functions in shallow parts of the plate boundary are composed mostly of surface waves.


Introduction
The 2011 Tohoku earthquake (M w 9) raised awareness of the dangers of future devastating earthquakes in nearby regions such as the Nankai Trough. The probability of a major earthquake occurring in the Nankai Trough has been estimated to be about 70 % within the next 30 years (Headquarters for Earthquake Research Promotion 2013). Figure 1 shows the distribution of strong motion generation areas (SMGAs) according to several recent source models. Tsurugi et al. (2006) constructed a source model for Nankai-Tonankai earthquakes by revising the source model of the Central Disaster Management Council of Japan (Central Disaster Management Council of Japan 2003), which was constructed by inversion of historical distributions of seismic intensity. Morikawa et al. (2013) proposed and simulated several source models with quasi-random distributions of SMGAs. Figure 1 also shows areas of large slip derived from waveform inversion for the 1944 Tonankai earthquake according to the Headquarters for Earthquake Research Promotion (Headquarters for Earthquake Research Promotion 2009) and for the 1946 Nankai earthquake .
Knowledge of the locations of SMGAs and asperities is important in earthquake source modeling. By waveform inversion, Yamanaka and Kikuchi (2004) identified asperities along the subduction zone off northeastern Japan and showed that there have been repeated ruptures on some asperities in the Tohoku region. They concluded that the asperities for future earthquakes are likely to be in areas where slip was large during previous earthquakes. From another side, it is important to know the largest likely ground motion in areas of critical infrastructure (e.g., nuclear power plants, high-rise buildings), which is dependent on the distribution of SMGAs.
Estimation of the null-space in the source inversion process is also important. This analysis usually occurs after inversion (frequently in the background): areas of large slip are individually tested to determine whether or not they produce waveforms of greater amplitude than background noise or inversion misfit (e.g., Sekiguchi et al. 2000;Kakehi 2004;Yoshida et al. 2011). A straightforward approach requires estimation of waveform amplitudes from every subsource to constrain slip for subsources that have small waveform amplitudes at all observation sites (e.g., Poiata et al. 2012).
Here we examined how the distribution of SMGAs can affect the amplitude of long-period ground motions at important sites (e.g., the megacity of Osaka). When the effect of seismic activity at SMGAs on ground motion at a particular location is tested by separate simulations for each SMGA, the effects of SMGA location and rupture propagation (e.g., the directivity effect) are mixed. In this study, we separated the effect of SMGA location from the effect of rupture propagation and calculated the distribution of peak amplitudes of Green's functions (GFs) at our target site, Konohana (OSKH02) in the center of the Osaka basin. We used the reciprocity theorem (e.g., Aki and Richards 1980) for our GF calculations in a complex 3D velocity structure model. This theorem states that when the locations and orientations of a seismic source and observation point are swapped, the same elastic response is observed. This procedure greatly reduced computational time, which allowed us to complete a detailed study in a reasonable time.

Source model
We used the source model of Central Disaster Management Council of Japan (2012), which includes 1404 subsources on a 10-km mesh in area bounded by a dashed line in Fig. 1. To calculate GFs for this study, we placed the subsources at the upper boundary of the subducting Philippine Sea plate and used the 3D Japan Integrated Velocity Structure Model (JIVSM) of Koketsu et al. (2012).
In applications of the finite-difference method (FDM), it is difficult to place a large number of subsources exactly on a curved plate boundary. Moreover, because that boundary separates media with different properties, a point source represented by a cubic cell of 2 × 2× 2 grid sizes (Graves 1996) may have nodes unintentionally assigned to different media. This kind of source settings may need additional treatment and testing. To avoid possible ambiguity related to artificial rotational forces, for example, we ensured that subsources were within a homogeneous medium by placing them two grid spaces (about 1.5 km) above the plate boundary. Strike and dip angles were also adjusted according to the shape of the

Calculation of Green's functions
To calculate GFs, we used the 3D FDM. This is an accurate method, but time-consuming because of the large number of subsources used for source inversion and forward ground-motion simulation. For example, inversion of the 1946 Nankai earthquake by Kagawa et al. (2012) involved 144 subsources, in effect, 288 subsources considering the two perpendicular rake angles used for inversion. Thus, the Kagawa et al. (2012) and Central Disaster Management Council of Japan (2012) models would require 288 and 2808 runs, respectively, of the 3D FDM simulation. To reduce our simulation time, we employed the reciprocity method, described in detail by Eisner and Clayton (2001) and Graves and Wald (2001), which calculates responses to three single forces in the NS, EW, and depth directions at the target site. The calculated waveforms for each subsource are then manipulated using the moment tensor contributions at the desired source location for each of the reciprocal force calculation, to replicate the responses to the subsource moment tensor at the target site. Thus, application of the reciprocity method to the Central Disaster Management Council of Japan (2012) model means that only three FDM simulations (rather than 2808) per target site are necessary to calculate GFs from all subsources.
First, we verified the validity of the 3D reciprocity method for calculation of GFs in the Nankai Trough by simulating three point sources in different parts of the trough (locations shown in Fig. 3). The simulated depths, dips, and strikes at these three points were similar to those of the CDMC model (Fig. 2). We used a dummy value for seismic moment M 0 = 10 20 Nm, a short rise time of T rise = 2 s, and a triangular source-time function, as is widely used in multi-time-window source inversions. The rake angle used in the simulations was 90°. Note that the specification of dip, strike, and rake for reciprocity simulations is not required at the FDM simulation stage; instead they are specified at the stage of manipulation of the moment tensor, which increases processing speed and is an additional advantage of the reciprocity method. Considering the uncertainty of the 3D velocity structure model in comparison to the real structure, comparisons of the waveforms of forward and reciprocal simulations for EW and NS components (Fig. 3) show little difference.
For our 3D velocity structure model, we used the wellcalibrated and waveform tested ) JIVSM model of Koketsu et al. (2012). This model is designed to deal with both sedimentary basins and crustal structures. It incorporates the Osaka basin model of Kagawa et al. (2004) and the Philippine Sea subduction plate interface and accretionary prism models of Baba et al. (2006), which together make it well-suited for earthquake modeling.
We carried out our 3D simulations by using the conventional staggered grid, fourth-order FDM scheme (Graves 1996) with a nonuniform grid size (Pitarka 1999). This approach simplifies the simulation of a double-couple seismic source and reduces computation cost (memory and time) for models with low-velocity sedimentary layers. The shortest wave period considered in the simulations was 4.2 s. The lowest shear-wave velocity in the JIVSM model, that defines combination of grid size and wave period, was 350 m/s (in the uppermost soft sedimentary layer in the Osaka basin having 500 m depth). The horizontal grid size for FDM simulations was 300 m (five grids for the shortest wavelength). A nonuniform grid size was used in the depth direction. For sedimentary layers shallower than 4200 m (the deepest level of soft oceanic sediments), a vertical grid spacing of 220 m was used; from 4200 to 9200 m depth (hard oceanic sediments), grid spacing was 440 m; and at greater depths, it was 880 m. The maximum depth of the calculation volume was 70 km, and we used 10,350 time steps at intervals of 0.0145 s.

Results
The distribution of peak amplitudes of calculated GFs is shown in Fig. 4 for rake angles of 0°and 90°for transverse and radial components, respectively; the distributions of subsource location, depth, strike, and dip are shown in Fig. 2. Similarly to the validity test of the application of the reciprocity method in the Nankai Trough (discussed above), we used a triangular source-time function, a rise time of T rise = 2 s, and a seismic moment of M 0 = 10 20 Nm for all subsources. Ground-motion periods of 4.2-20 s were considered. For convenience, peak GF amplitudes were normalized to a maximum value of 20. The distribution of peak amplitudes of GFs was markedly irregular. Two areas showed anomalously large amplitudes (Fig. 4): a large area elongated along the plate strike just south of the Kii Peninsula for transverse component for rake = 0°(anomaly 1) and for radial component for rake = 90°(anomaly 2) and a large area off the Kii Peninsula, trenchward of anomalies 1 and 2, for radial component for rake = 90°(anomaly 3). Although GF amplitudes would normally be expected to decrease with increasing distance from the source, they anomalously increased with increasing distance perpendicular to the trench axis. Figure 5 shows example transverse and radial waveforms for a line perpendicular to the trench axis.
The elongations of both anomalies correspond well with iso-depth lines on the upper surface of the Philippine Sea , while target site Konohana is located on perpendicular to these isolines. For this reason, we speculate that plate shape may have a critical effect on simulated ground motions because of the cumulative effect of subsource radiation patterns and the distributions of strike and dip.

Effects of wave generation and propagation on GF amplitudes
There are six factors that may contribute to the observed GF amplitudes: source distance, source mechanism, generation of surface waves, seismic potency and seismic impedance at the source location, oceanic sediments in Fig. 3 Validation of the reciprocity approach. Comparison of reciprocal waveforms (thick lines) and forward simulated waveforms (thin lines) for EW (left) and NS (right) components for point sources Src1, Src2, and Src3, and for JIVSM model. Numbers at the left of each trace are normalized peak amplitudes the accretionary prism, and basin sediments (Osaka basin) at the target site. Below, we analyze each of these factors.

Source distance
GF amplitudes are expected to decrease with increasing distance from the source (approximately proportional to 1/R for body waves and 1/R 0.5 for surface waves, where R is the distance from the source). This is roughly the case for subsources southwest of the target site, along the strike of the subducting plate (Fig. 4). In contrast, for subsources progressively farther southeast of the target site, this regular decrease of amplitude due to geometric spreading is not apparent. Other effects, that are enumerated above and will be discussed later, overlap with distance attenuation effect and results in a complicated distance dependence of amplitudes.

Source mechanism
To test the effect of source mechanism on GF amplitudes, we considered the amplitude distribution of S-waves from subsources in an infinite uniform medium according to the Aki and Richards (1980;their equation 4.33) as shown in Fig. 6. In contrast to the distribution of GF amplitudes shown in Fig. 4, the highest amplitudes are in areas of steeper dip closer to the target site, and low to medium amplitudes are in the areas of anomalies 1 and 3. We concluded from these results that the source mechanism alone cannot produce the simulated anomalies. Source depth distribution: generation of surface waves, seismic potency, and seismic impedance Shallow subsources generate large-amplitude surface waves that may result in GF amplitudes at the target site being larger because of the reduced effect of geometric spreading. Amplitudes of generated surface waves are proportional to amplitude of eigenfunctions, which in general increase both with decreasing of depth and with decreasing of velocity values at shallow depths. From another side, for a given seismic moment, the effects of seismic potency (Ben-Zion 2001) and seismic impedance at a source location are depth dependent. In order to evaluate the combined effect of these phenomena, we ran simulations with a 1D velocity structure (Figs. 7 and 8) that is equivalent to the onshore crustal velocity structure of the JIVSM model.
For the transverse component for rake = 0°, GF amplitudes were small for both deep and shallow subsources, but were clearly larger at intermediate depths (Fig. 7a), which is in good agreement with anomaly 1 (Fig. 4a). The larger amplitude waves of the transverse component should be surface Love waves or body S-waves. Amplitudes of Love waves are described by Aki and Richards (1980;their equations 7.148 and 7.149). For shallow subsources, Love wave amplitudes can be small because the terms related to small polar angles approach zero (true for subsources aligned on the perpendicular to the trench axis) due to the small horizontal moment tensor components for dip~0°and rake = 0°and the small depth derivatives of eigenfunctions of Love waves for lower modes at shallow depths.
For body waves, the effect of seismic potency at shallow depth, which increases amplitudes, may be smaller than the effect of the source radiation pattern, which decreases amplitudes for shallow subsources. Consequently, GF amplitudes for the transverse component were smaller at shallower depths and increased with increasing depth owing to the effect of the source radiation pattern (Fig. 7a). These interpretations are supported by example waveforms (Fig. 8); large-amplitude segment of transverse components correspond to large-amplitude S-waves (note waveforms for subsources numbered 6-12 in Fig. 8).
For the radial component for rake = 90° (Fig. 7b), at first glance, the distribution of peak GF amplitudes is consistent with our original assumption for surface-wave generation; that is, amplitudes are large for shallow subsources and decrease with increasing depth. On the contrary, if we refer to Aki and Richards (1980;their equations 7.150 and 7.151), we see that although the eigenfunctions of Rayleigh waves and their derivatives are not zero at zero depth, all terms of equations 7.150 and 7.151 vanish due to the vanishing of terms related to the small polar angle in the direction perpendicular to the trench axis and the small moment tensor components M xx , M xz , and M zz for dip~0°and rake = 90°. For this reason, we rejected the generation of Rayleigh waves as an explanation for anomaly 3 (Fig. 4).
On the other hand, the shape of the large-amplitude area in Fig. 7b corresponds very well with the shallow subsource area (normalized amplitudes around 4) in Fig. 6b, which reflects the effect of the radiation pattern of SV-waves. Although the effect of the radiation pattern for shallow subsources alone cannot explain the large shallow anomaly in Fig. 7b, it can explain it if combined with the effect of greater seismic potency at shallow depth. Example waveforms shown in Fig. 8 confirm this interpretation: large-amplitude waves for the radial component for rake = 90°from shallow subsources propagate at S-wave velocities.
However, despite the above investigation and analysis, an important question about the effect of seismic potency and seismic impedance at the source location remains unanswered. Namely, can these effects compensate the effect of the source radiation pattern or not? This question is considered in detail below, after our discussion of the effects on GF amplitudes of sedimentary layers in the accretionary prism and in the onshore basin.

Oceanic sediments (accretionary prism) and basin sediments at the target site (Osaka basin)
To test the effect of sedimentary layers on GF amplitudes, we ran simulations with a velocity structure model without sedimentary layers (V S ≤ 2.4 km/s) and compared the results with those obtained with the JIVSM model. Figure 9 shows the distributions of GF Dashed lines indicate arrival of waves propagating with velocities 5.8, 3.4, and 2.0 km/s, which correspond approximately to P-waves, S-waves, and surface waves, respectively amplitudes for simulations without sedimentary layers, and Fig. 10 shows some example waveforms. The cross section of the velocity structure shown in Fig. 10 can be compared with the cross section of the JIVSM model shown in Fig. 5. Because of the similarity of the 1D velocity structure model (cross section in Fig. 8) and the velocity structure without sedimentary layers (cross section in Fig. 10), the simulation results for the two  Close examination of the transverse component waveforms corresponding to anomaly 1 (subsources numbered 5-9 in Fig. 5) reveals noticeable dispersion, which is unexpected for waves propagating at S-wave velocities. Characteristic for surface waves dispersion in the Osaka basin may arise from the generation of basin waves by incident S-waves. Similar dispersion is not evident for the radial component. However, the predominant wave period for the radial component is 6.2-6.7 s, which is similar to that of the transverse component (6.0-6.7 s), and is also similar to the predominant wave period observed during the 2004 off Kii Peninsula earthquake (M w 7.4; e.g., Miyake and Koketsu, 2005). Miyakoshi et al. (2013) demonstrated that the generation of basin waves (Love and Rayleigh   Fig. 10 were simulated by FDM using the reciprocity method and the JIVSM velocity structure model, but with layers for which V S < 2.4 km/s were replaced by layers with V S = 2.4 km/s. S-wave velocity cross section is shown in Fig. 10. Other annotations as in Fig. 6 waves) is responsible for the predominantly long-period strong ground motions characteristic of earthquakes in the Osaka basin. This phenomenon may explain the amplification of waves from subsources associated with anomaly 2.
Moreover, Fig. 5 shows that for the radial component for rake = 90°there are noticeable time delays for the large-amplitude waves from the four shallowest subsources (numbered 13-16 in Fig. 5) for the JIVSM simulations, which were not observed for simulations with the velocity structure without sediments in the accretionary prism (Fig. 10). The time delays, as well as longer durations of the waves, indicate that generation of surface waves in the accretionary prism also increases the amplitudes of waves from subsources in the area of anomaly 3.
Detailed discussion of effects of seismic potency and seismic impedance of subsource One intriguing question remains unexplained: If the effect on GF amplitudes of body-wave propagation from the deep subsources closest to the target site is stronger than elsewhere, why are GF amplitudes there low, even though the effect of the radiation pattern shown in Fig. 6 Fig. 10 Example of waveforms simulated with velocity structure without sediments. Cross section of S-wave velocity structure and location of subsources (left), transverse component waveforms for rake = 0°(center), and radial component for rake = 90°(right). Other annotations as in Fig. 5 is strongest there? To investigate this question, we considered the distribution of subsources among the various layers of the velocity structure model; that is, we determined which of them were within each of the mantle wedge, lower crust, upper crust, and seismic basement layers of the model. Figure 11 shows a result of this analysis superimposed on the simulated distribution of GF amplitudes for the radial component for rake = 90°(as in Fig. 4b).
This analysis demonstrates that the lower boundary of anomaly 2 corresponds with the boundary between the lower and upper crust. GF amplitudes from subsources in the upper crust are about three times larger than those from subsources in the lower crust. There are two possible causes of the lower GF amplitudes from subsources in the lower crust and mantle: the effect of lower seismic potency at greater depths, as previously discussed, and the effect of waves reflected at the Moho and Conrad discontinuities. Reflection of upward incident waves would decrease GF amplitudes, whereas reflection of downward incident waves would increase GF amplitudes. Because it is difficult to isolate the effect of reflection of long-period waves in a complex 3D velocity structure, we instead evaluated the effects of properties of the source media. That is, we examined the effects of seismic potency and seismic impedance on GF amplitudes and assumed that wave reflections at the Moho and Conrad discontinuities accounted for all that we could not attribute to seismic potency or seismic impedance.
According to the generally accepted convention for earthquake source modeling, we used the same seismic moment for all subsources. Ben-Zion (2001; see also Aki and Richards 1980 second edition) considered that seismic wave generation in an earthquake source region may be better described by seismic potency, which is defined as the ratio of seismic moment and shear modulus (M 0 /μ). Another effect of source medium is related to changes of seismic impedance, which is the product of density and wave velocity (ρV S for S-waves, where ρ is density and V S is S-wave velocity). During propagation of waves from the source to target site, wave amplitude increases as the square root of the source-to-site impedance ratio. Because μ = ρ V S 2 , both effects can be described in terms of density and S-wave velocity near the source. Aki and Richards (1980;their equation 4.97) combined these and other effects of density and S-wave velocity into a single amplification coefficient (c): where V S and ρ are the S-wave velocity and density at the subsource, respectively, and V S0 and ρ 0 are the Swave velocity and density at the target site, respectively. The velocities and densities of the layers of the JIVSM model and values of c calculated for each layer according to Eq. (1) are shown in Table 1. GF amplitudes for subsources in the upper crust are larger by a factor of 1.33 than those in the lower crust (Table 1), which is smaller than the factor of 3.0 for corresponding GF amplitudes shown in Fig. 11. We assigned the difference (around 2.25) to the combined effect of amplitude losses due to reflection at the Conrad discontinuity of waves coming from subsources below the discontinuity, and amplitude increases due to constructive interference of direct waves and waves reflected from the discontinuity for subsources above it. Thus, in addition to seismic potency, wave reflections at the Conrad discontinuity may contribute to the formation of the lower boundaries of anomalies 1 and 2.

Analysis of the time delay of peak amplitudes
We next considered the spatial distribution of the time delay of peak amplitude (peak delay time). Figure 12 shows the distribution of reduced peak delay time for the radial component for rake = 90°, calculated as where T peak is the peak delay time and R is the hypocentral distance; we assumed V S = 3.4 km/s, the JIVSM velocity for the upper crust. Figure 12 confirms that the propagation times of peak amplitudes for anomaly 2 are compatible with body S-waves, whereas those for anomaly 3 are compatible with the surface waves.

Discussion
All effects considered above, as well as their influence on GF amplitudes, are compiled in Table 2 for further comparison. Our analysis supposes that effects of source mechanism and accretion prism may be strong for the formation of anomaly 3, effect of basin sediments may be strong for the formation of anomaly 2, while effect of reflection from Conrad discontinuity may be strong for the formation of anomalies 1 and 2.
Comparison of the locations of previously defined SMGAs ( Fig. 1) with the locations of anomalies 1, 2, and 3 gives an indication of which of the SMGAs may be primarily responsible for ground motion at our Konohana target site. Simulation parameters for these SMGAs need to be fine-tuned to provide accurate predictions of strong ground motion there.
The velocity model we used for our simulations did not include an oceanic water layer. Inclusion of a water layer and seafloor topography reduces the amplitude of surface waves because of the scattering effect of topography and may also reduce the amplitude of surface Raleigh waves because of the suppression effect of the water mass (Petukhin et al. 2010). Both of these effects may reduce the GF amplitudes of anomaly 3.
In a near trench area, a simplified fault plane parallel to the subducting plate is assumed in this study. A real fault system consists of deep mega-splay fault (detachment fault) and seaward decollement; both are plateparallel, and the mega-splay fault that is branching to the surface and having a dip angle is considerably larger than assumed here (e.g., Tsuji et al. 2014). Ruptures at subsources on the mega-splay fault would result in large GF amplitudes due to intensive generation of surface waves.
Although the effect of the generation of surface waves in the accretionary prism is considered to be important to explain anomaly 3, this is a badly constrained part of the velocity structure model. Recently developed network of the ocean-bottom seismometers DONET (e.g., Nakano et al. 2013) should be helpful to improve accretion prism structure model and to validate effect of generation of surface (basin?) waves by direct observation.
It is important to note that the GF anomalies discussed here are specific to the Konohana target site, which lies at the center of the large Osaka basin hosting a megacity. For other megacities, such as Tokyo, Nagoya, and Fukuoka in Kyushu, the distribution of anomalies differs from those presented here. If our conclusions here are valid, GF anomalies relevant to Nagoya would lie east of the anomalies discussed here, those relevant to Hiroshima would be off Shikoku, and those relevant to Fukuoka would be in the Hyuga-Nada region. Tokyo, however, lies in a nodal direction from the megathrust source, so small GF anomalies may be expected there, although the influence of the directivity effect due to rupture propagation is strongest among all megacity cases.
Wave interference within a basin can cause large variations of GF amplitudes that are related to the location of the target site inside the basin, the target period and component orientation of ground motion. Sekiguchi et al. (2008) and Kawabe and Kamae (2008) provide examples of these for the Osaka basin. Direction of wave incidence also affects responses within a basin. For all these reasons, GF amplitudes may differ even for sites that are close together inside the basin. As concluded above, reflection of waves from deep or shallow subsources at the Moho and Conrad discontinuities has a considerable effect on attenuation or amplification of those waves; this can be attributed in part to the relative simplicity of the JIVSM model, which consists of layers of constant internal velocity and with large impedance contrasts at layer boundaries. For a model such as this, Graves and Pitarka (2010) considered that reflections from the Moho discontinuity make an important contribution to GF amplification because they arrive simultaneously with the direct waves. This was supported during validation of the JIVSM model by Petukhin et al. (2012) and T. Masuda (personal communication) using records of crustal earthquakes; they noted that the amplitudes of waves reflected at the Conrad discontinuity were larger than the observed waves. It may be necessary to revise the structure of the velocity model by introducing velocity gradients within layers and reducing impedance contrasts between them to improve simulations and strong ground-motion predictions.