Stress loading and the occurrence of normal-type earthquakes under Boso Peninsula, Japan

Boso Peninsula, Japan, was formed by the interaction of the Philippine Sea, Eurasian and Pacific plates around the trench–trench–trench Boso triple junction. Normal-type earthquakes are persistently observed in the subducting Philippine Sea slab under the peninsula at a depth of ~ 30 km, including a recent (2019) Mw 4.9 earthquake which caused shaking throughout the Kanto region (greater Tokyo). Such shallow intraplate earthquakes are potentially hazardous to this heavily populated region, yet their mechanism is poorly understood, especially in the context of a three-plate system. Here, we calculate stress rates in the Philippine Sea slab and the surrounding area, using a subduction model constructed in a previous study, to explain the generation of the regional stress field and its effect on earthquake occurrence. In general, the calculated stress rates under Boso Peninsula are horizontally extensional both above and below the Eurasian–Philippine Sea plate interface. We apply our calculated stress rates to the nodal planes of the observed earthquakes to calculate the Coulomb failure function (ΔCFF). These calculated ΔCFFs are generally positive on normal-type earthquakes under Boso. The ΔCFFs are also consistent with earthquakes in adjacent areas that are seismically active, for example, in the Philippine Sea plate to the south, in the collision zone around Izu Peninsula, and in the cluster in the Eurasian plate northeast of Boso Peninsula, which further supports our stress loading model. Calculation of the individual contributions of Philippine Sea plate and Pacific plate subduction shows that the development of the stress field around Boso is dependent upon contributions from both subducting plates. In contrast, the arc–arc collision at Izu Peninsula has little influence.


Introduction
Off Boso Peninsula in central Japan, the Eurasian (EUR), Philippine Sea (PHS) and Pacific (PAC) plates meet to form a unique trench-trench-trench triple junction. The closest onshore expression of this triple junction is Boso Peninsula itself, which lies on the uppermost EUR plate. It is marked by a rapid Holocene uplift rate of up to 5 mm/year, estimated from the height of paleo-shorelines (Koike and Machida 2001).
Under the peninsula, normal-type earthquakes are persistently observed in the subducting PHS slab at a depth of ~ 30 km (Nakajima et al. 2011;Imanishi et al. 2019). A recent example is the 2019 May 25 M w 4.9 earthquake (Fig. 1), from which moderate shaking was felt throughout the Tokyo metropolitan area, in the Kanto Basin.
The occurrence of normal-type earthquakes in an overriding plate appears contradictory in a convergent plate setting. Nevertheless, normal-type earthquakes in subduction zones worldwide occur around major extensional structures, typically backarc basins such as the Mariana trough, the Lau Basin, and the Aegean Sea. In addition, earthquakes beneath the outer rise in subducting slabs are mostly normal-type, ascribed to slab bending (Christensen and Ruff, 1988). Below 70 km depth, unbending results in a double-seismic zone, with compressive earthquakes in the upper part of the slab and extensional below. (Engdahl and Scholz 1977;Hasegawa et al. 1978). Between depths of 20 and 70 km, shallow slabs are generally considered to be in a neutral stress state, although Seno and Yoshida (2004) found that earthquakes of M ≥ 7, mostly normal type, occur in this depth range in circum-Pacific slabs. They suggested that the character of intra-slab activity may be linked to a pattern of increasing compression from arc to trench, which indicates the importance of considering plate-toplate interactions in the context of the whole system on generating intraplate stress.
The Japanese islands are dominated by reverse or strike-slip type earthquakes except in the Okinawa trough or the Beppu-Shimabara graben, where extensional structures develop (e.g. Terakawa and Matsu'ura 2010). Upon closer look, however, there are minor clusters of normal-type earthquakes elsewhere. In the Fukushima coastal area, normal-type events, including the M w 6.6 Iwaki earthquake, were activated after the 2011 Tohoku earthquake, where previously the area had been calm with tiny extensional events (Imanishi et al. 2012). At Boso Peninsula, minor clusters of normal-type earthquakes occur, mixed with other type events (Imanishi et al. 2019).
Intraplate earthquakes are considered to release the stress which has been gradually accumulated over thousands of years through steady plate subduction. However, under Boso, the loading mechanism on these potentially hazardous, shallow, intra-slab earthquakes is less wellunderstood because of the additional complication of the three-plate setting.
To intuitively understand and quantitatively estimate a loading mechanism of this kind, modeling is a useful approach. Mechanical modeling of plate subduction has been developed to explain the dynamics of slab behaviors in a crust-mantle system, including slab rollback that leads to extension in the overriding plate (e.g., Funiciello et al. 2003;Enns et al. 2005;Holt et al. 2018). However, the resolution of this approach is too low for the purposes of this study. Sato and Matsu'ura (1988) and Matsu'ura and Sato (1989) developed a dislocation model for steady plate subduction where relative slip is kinematically assigned to the known plate interface. Internal deformation and stress are then solved as a boundary problem of an elastic-viscoelastic Earth. The dislocation model can be easily applied to a 3-D plate interface. Hashimoto et al. (2004) and Hashimoto and Matsu'ura (2006)  simulations for vertical displacements and stress rates around Japan, showing an overall match with observation. Takada and Matsu'ura (2004) and Hashimoto and Matsu'ura (2006) extended the model to treat continental collision by adding a slip-rate deficit to steady subduction. In contrast, Hashima et al. (2008a) showed that sliprate excess (forward slip) produces backarc extension.
Recently, Fukahata and Matsu'ura (2016) theoretically proved that simple 2-D subduction leads to extension in the overriding plate.
In a previous study, we developed a dislocation model for the subduction of both the PHS and PAC plates and the collision of the Izu-Bonin (Ogasawara) arc using onshore geological constraints for the last 1 Myr (Hashima et al. 2016). Here, we use this model to explain the development of the stress field around Boso Peninsula, and then discuss its relationship to seismicity.

Tectonic setting around Boso Peninsula
Boso Peninsula is part of the Kanto basin system (Fig. 1a), which developed as a forearc basin as a result of subduction of both the PAC and PHS plates (Matsuda 1991). It is now covered by thick Cenozoic sediments. An eastwest trending structural high cuts across southern Boso Peninsula at 35.2° N, bounding the forearc basin and the trench slope. The development of the Kanto Basin system has been also strongly affected from the west by the arcarc collision at Izu Peninsula, which was originally part of the Izu-Bonin arc (Sugimura 1972;Matsuda 1978).
Presently, the PHS plate subducts northwestward under the Sagami trough at a rate of ~ 3 cm/year, while the PAC plate subducts from the east at ~ 8 cm/year (Fig. 1). The Sagami trough is known to have hosted M8-class Kanto earthquakes west of Boso Peninsula in 1923 and 1703 (Fig. 1a). Based upon an earthquake cycle model with viscoelastic media, the recurrence interval of these events was estimated to be 350 years (Sato et al. 2016). To the east of Boso, seismic activity is of a different character, with ~ M6 slow-slip events (Fig. 1a) occurring approximately every 6 years (Ozawa et al. 2003;Sato et al. 2017). Further north, the Japan trench was ruptured by the 2011 M9 Tohoku earthquake (Fig. 1b), which significantly affected seismicity around the Kanto basin (Ishibe et al. 2011(Ishibe et al. , 2015. To characterize seismic activity around the Kanto Basin, we take shallow (≤ 40 km) moment tensors of magnitude M w > 3.5, with a variance reduction of ≥ 70%, for the period 1 January 2003 to 30 September 2019 from the F-net focal mechanism catalog of the National Research Institute for Earth Science and Disaster Resilience (NIED) (Okada et al. 2004). Using the geometry of the plate interface from the CAMP model of Hashimoto et al. (2004), we divide the earthquakes into PHS, EUR, and PAC earthquakes according to their location. Then, we exclude PAC earthquakes. The number of the PHS and EUR earthquakes is 196 and 202. It should be noted that we do not distinguish between inter-and intraplate earthquakes with this classification because it is difficult to identify a single EUR-PHS interface due to a scarcity of earthquakes with clear thrust-type mechanisms in the F-net catalog. Indeed, the position of the EUR-PHS interface under the Boso Peninsula varies widely among different studies (Nakajima et al. 2009 and references therein). It is possible that our classification which relies on a plate interface model may thus introduce a bias in interpreting the overall mechanism pattern, given the ambiguity of distinguishing inter-and intraplate earthquakes.
For the purposes of our plots, moment tensors τ ij are classified in terms of mean horizontal normal component τ m = (τ 11 + τ 22 )/2 , where directions 1, 2, and 3 denote east, north, and up, respectively. The moment tensors are normalized as the normalized mean horizontal moment tensor τ ′ m = τ ′ 11 + τ ′ 22 /2 varies from − 0.5 to 0.5. Positive values denote horizontal extension and negative values denote horizontal compression, which correspond to normal-and reverse-type events, with strike-slip type events falling in between. Figure 2 shows the PHS and EUR earthquakes in map and side view; the inset histogram is plotted with respect to τ ′ m . We find that the east coast of Boso Peninsula is seismically active in both the hanging wall and foot wall of the plate interface. Normal-type earthquakes occur on the PHS side ( Fig. 2a), showing northwest-southeast extension. Normal-type earthquakes are also found in the upper plate (EUR) (Fig. 2b), although the population is more mixed with reverse types. To show the detailed vertical distribution of fault types around Boso, we plot north-south cross-sections between 139.9° E and 140.5° E, which we separate into two periods, before and after the 2011 Tohoku earthquake (Fig. 2c, d). In terms of seismicity, the PHS plate is more seismically active than the EUR plate before the M9, but less active after, while the fault types remain basically the same (Ishibe et al. 2011(Ishibe et al. , 2015. Their frequency distribution with respect to τ ′ m shows the superiority of the compressional-type earthquakes especially after the 2011 Tohoku earthquake, but most of them occur close to the plate interface (Fig. 2c,  d). In both periods, we see the broad occurrence of normal-type earthquakes above and below the PHS upper surface. From this distribution, we infer that stress states in the both EUR and PHS plates favor normal-type earthquakes. In contrast, reverse-type earthquakes, which prevail around 35.4° N into the PHS plate, can be viewed as an anomalous localized activity. South of the Sagami trough (< 34.5° N), strike-slip type earthquakes dominate with an extension axis normal to the Sagami trough.
Subsequently, we re-classify the earthquakes around Boso (Fig. 2c, d) into three groups: (A) major intraplate earthquakes, (B) interplate earthquakes, and (C) anomalous reverse-type earthquakes around 35.4°N. Group B is defined by thrust events ( τ ′ m < −0.1 ) located within 5 km of the EUR-PHS plate interface of the CAMP model. Group C is defined by non-interplate events with τ ′ m < −0.1 located between 35.3° N and 35.5° N latitude. The rest fall into Group A. The number of events of Groups A, B, and C is 19, 13, and 4, respectively, before the 2011 Tohoku earthquake and 26, 15, and 5, respectively, after the Tohoku earthquake.
To the west, in the PHS plate, we see seismicity around Izu Peninsula related to the arc-arc collision (Fig. 2a). Seismic activity along the longitude of 139° E varies from reverse-type at around 35.5° N to strike-slip type at around 34.5° N but the compressive axis is commonly northwest-southeast. The southernmost cluster is mostly related to volcanic activities but may reflect the regional stress field. The cluster around ( M6.5 earthquake (Fig. 2a). On the EUR side of the Izu collision zone, seismicity is much sparser (Fig. 2b), but we can see that the compressive quadrants trend from E-W and NW-SE around 138° E to N-S around 139° E [the event at (139° E, 35.5° N) may be negligible because it is likely an interplate event]. This trend, presumably related to the Izu collision, is consistent with a radial stress pattern estimated from surface fault traces (Matsuda 1977). There is also a large cluster of seismicity northeast of Boso Peninsula around 141° E in the EUR plate (Fig. 2b), but most of this activity occurred after the 2011 Tohoku earthquake.

Model setting
Following Hashimoto et al. (2004)  Around Izu Peninsula, steady plate subduction is impeded by the buoyancy of the Izu-Bonin arc crust. According to Takada and Matsu'ura (2004) and Hashimoto and Matsu'ura (2006), this effect can be represented as a permanent slip deficit. In this study, we use the distribution of the slip-rate deficit in the final stage of the 1 Myr-long rotation of the slip-rate deficit zone (locked zone) presented in Fig. 8d of Hashima et al. (2016) (Fig. 1a), which is determined by geomorphological, geological, and thermochronological constraints. In this stage, the locked zone is mainly located to the northwest of Izu Peninsula with an east-west width of 110 km. The collision rate (slip-deficit rate/convergence rate) is assumed to be 1 (fully locked) in the locked zone.
Following Sato and Matsu'ura (1988), we calculate the stress rate σ ij (x) in the plate interior for steady slip rates w k (ξ) assigned over the plate interfaces S (ξ) as follows: where G ijk (x, t = ∞; ξ, 0) denotes the final state (t = ∞) of the stress response to a unit slip in the k-direction. The explicit expressions for G ijk (x, t = ∞; ξ, 0) are derived by Fukahata andMatsu'ura (2005, 2006), and Hashima et al. (2008bHashima et al. ( , 2014. The equivalence theorem of Fukahata and Matsu'ura (2006) is also applied to efficiently compute the final state of the viscoelastic response. The calculated stress rate is gradually accumulated within the plate to develop the regional stress field as steady plate subduction proceeds, while the effects of repeated locking and rupture at the plate interface are canceled in the long term (Sato and Matsu'ura 1988;Matsu'ura and Sato 1989).
There are several ways to examine the match between the calculated stress loading rates σ ij and the observed earthquakes (Fig. 2), which release stress. The easiest way might be to compare directly the calculated stress loading rates and the observed moment tensor, assuming long-term balance and statistical consistency. However, individual moment tensors are not strictly equal to the released stress itself and may also be affected by the factors like local weak planes or variable pore pressure values. In this study, we evaluate the Coulomb failure function (ΔCFF) for the nodal plane of each earthquake. ΔCFF denotes the preference of loaded stress for the rupture of a fault, which is simply defined by the difference between loaded shear stress and fault strength as where σ s and σ n are the shear and normal (positive in compression) stress on the fault plane, respectively, and where normal vector n i and slip direction vector ν i . n i and ν i can be obtained by the strike, dip, and rake of the nodal plane in the F-net catalog; μ′ is the apparent friction coefficient. We assume a standard value of 0.4 for μ′. Positive ΔCFF indicates that the fault shifts closer to rupture and negative ΔCFF indicates a shift away from rupture. Figure 3a, b shows the calculated stress rates at depths of 10 km and 30 km, which correspond to the stress above and below the PHS plate interface around Boso Peninsula. Figure 3a is the same result as Fig. 10a of Hashima et al. (2016) except that here we include the full stress components excluding the trace. Around Boso Peninsula, the calculated stress rates are largely extensional, ranging

Results
(2) �CFF = σ s − µ ′ σ n with σ s = σ ij n j ν i , σ n = −σ ij n j n i , from uniaxial extension near the Sagami trough to allaround horizontal extension (represented by a doughnut-like symbol) in the Kanto Basin, with normal-type mechanisms in between. This explains the occurrence of normal-type earthquakes under Boso Peninsula. Previous stress estimates inverted from focal mechanisms show various stress types in the area under Boso Peninsula, which might be due to the inclusion of interplate events (Terakawa and Matsu'ura 2010;Imanishi et al. 2019). To the west, the calculated stress rate field also appears consistent with earthquake mechanisms and previous stress estimates around Izu Peninsula, and further south, as discussed in Hashima et al. (2016).
In Fig. 3c, d, we decompose the stress rates at 30 km depth into the contribution of the subduction of the PAC and PHS plates, together, and the contribution of the Izu collision. We observe that the subduction effect creates normal-type to all-around horizontal extension (Fig. 3c) while the collision effect dominates the area above the locked zone with compression (Fig. 3d). Around the eastern part of Boso Peninsula, the collision effect is smaller and may be negligible. Figure 4 shows how the calculated stress rates match actual earthquake occurrences in terms of ΔCFF. Here, we take nodal plane 1 of the F-net catalog to compute ΔCFF. The results for nodal plane 2 show similar characteristics (Additional file 1: Fig. S1). The map-view representations (Fig. 4a, b) show generally positive ΔCFFs over our study area including those around Izu Peninsula and in Suruga Bay in the PHS plate, and in the cluster northeast of Boso Peninsula around 141° E on the EUR plate side. Focusing on Boso Peninsula, positive ΔCFFs appear a b c d less dominant because of the active thrust events, especially for nodal plane 2 (Fig. 4c, d, Additional file 1: Fig.  S1c, d). However, we can see the tendency of normal-type earthquakes both above and below the plate interface for positive ΔCFF and that of thrusts, which are more concentrated on the plate interface, for negative ΔCFF. To see this tendency clearly, we plotted the correlation between ΔCFF and fault type in terms of normalized mean horizontal moment tensor τ ′ m (Fig. 4e). Group designations A, B, and C are marked by different colors, while timing before or after the 2011 Tohoku earthquake is distinguished by circles and triangles. Figure 4e shows a preference (positive ΔCFF) for normal-type earthquakes (positive τ ′ m ). In particular, Group A (major intraplate) events have mostly positive ΔCFF, supporting our stress loading model for the occurrence of normal-type earthquakes within the plate. Events with a negative ΔCFF are mostly from Groups B and C, which represent interplate earthquakes and anomalous reverse-type activity, around 35.4° N, respectively. The loading mechanism of Group C is unclear, but it might be related to local stress perturbations or rheological heterogeneity. Events with a positive ΔCFF are also found in Groups B and C, but their ΔCFF values are smaller than those of Group A normal-type events. If we compare events from before and after the 2011 Tohoku earthquake, the correlation between positive ΔCFF and group designation becomes less clear after 2011, probably due to external perturbation of the stress field from the M9 Tohoku earthquake. The northern Kanto Basin around (140° E, 36° N), dominated by negative ΔCFFs (Fig. 4b), is an exception because the calculated horizontal extension is inconsistent with the observed reverse-type earthquakes. Our previous study also shows a discrepancy between calculated subsidence and observed uplift in this area (Hashima et al. 2016), which can be ascribed to the omission of the effect of compression tectonics in the northeastern Japan arc (Sato 1994).
We also examined the match between stress rates and moment tensors by taking the normalized tensor dot product (Additional file 1: Fig. S2). The result is similar to Fig. 4 and Additional file 1: Fig. S1, but the correlation between the normalized tensor dot product and τ ′ m is more clear (Additional file 1: Fig. S2e).

Discussion
In this study, we calculated stress rates accumulated around Boso Peninsula and the Kanto Basin, based on the dislocation model that we previously developed to explain the evolution of the geological deformation pattern in this area (Hashima et al. 2016). Though the model neglects the elastic slabs under the elastic-viscoelastic boundary at depth of 40 km, it successfully explains the observed occurrence of shallow intraplate earthquakes with normal-type mechanisms under Boso Peninsula, including the most recent M w 4.9 earthquake. This indicates the importance of steady plate subduction in stress accumulation.
Before the period covered by the F-net catalog, a moderate M 6.7 earthquake occurred in 1987 east of Boso Peninsula, resulting in injury, loss of life, and broad damage to the urban function. This earthquake is also considered an intraplate earthquake in the PHS plate. The focal mechanism is strike-slip (e.g. Okada and Kasahara 1990), but its extension axis directing northwest-southeast is similar to the normal-type earthquakes under Boso Peninsula. Nakajima et al. (2009) and Nakajima and Hasegawa (2010) pointed out that the hypocenter is located at the western boundary of the serpentinized front wedge of the PHS slab sandwiched by the EUR and PAC plates. We applied our model stress to this event, using data from the Harvard CMT catalog (Ekström et al. 2012), and obtain positive ΔCFFs for both nodal planes ( Fig. 4a and Additional file 1: Fig. S1a). This result seems promising but should be carefully examined. The depth of the Harvard CMT solution is estimated to be 20 km, much shallower than previous estimates of 40-60 km (e.g. Okada and Kasahara 1990) where our model would be inapplicable. If the stress field at this depth (40-60 km) is similar to our model results, then the 1987 earthquake could be explained with the stress loading. This should be tested in a model with elastic slabs.
The spatial range of the locked zone with the magnitude of the collision rate around Izu Peninsula is the primary control parameter of the model. In this study, we adopted Hashima et al. (2016)'s locked zone, which is determined by geomorphological, geological, and thermochronological uplift/subsidence data for the period of 1 Myr. Based on 3-D stress data (Terakawa and Matsu'ura 2010), Hashimoto and Terakawa (2018) developed a method of stress data inversion and applied it to obtain a detailed distribution of the collision rate around Izu Peninsula. They obtained a wider locked zone, closer to Boso Peninsula, with a low collision rate (≤ 0.4).
To test the effect of a different locked zone, we set a locked zone after Hashimoto and Terakawa (2018), with a uniform collision rate of 0.3, and calculated the stress rates (Additional file 1: Figs. S3-S6). Focusing on the east side of Boso Peninsula, the calculated stress pattern is allaround horizontal extension to normal type (Additional file 1: Fig. S3), showing only small differences with our model (Fig. 3). The effect of the collision is strong near the edge of the locked zone but rapidly decreases to the east. We also examine the match between the calculated stress rates and the observed moment tensors in terms of ΔCFF (Additional file 1: Figs. S4 and S5) and of normalized tensor dot product (Additional file 1: Fig. S6). The results are also similar to our results ( Fig. 4 and Additional file 1: Figs. S1, S2). Therefore, for both locked zone settings, earthquake occurrence under Boso Peninsula appears to be dominated by the subduction effect rather than the collision effect.
In Fig. 5, we attempt to show the separate contributions of PAC and PHS plate subduction to the regional stress. Figure 5a, b shows the map-view visualization of the PAC and PHS subduction effects at 30 km depth. Figure 5c, d shows cross-sectional views through Boso Peninsula normal to their respective trenches. We can see that normal-type earthquakes with east-west extension develop with PAC subduction (Fig. 5a, c). Similarly, north-south extension develops with PHS subduction except around Izu Peninsula, where the geometrical effect of the plate interface appears to dominate (Fig. 5b, d). Thus, the development of 'doughnut-type' all-around horizontal extension around Boso (Fig. 3c) is explained by the superposition of both PAC and PHS plate subduction.
Fukahata and Matsu'ura (2016) theoretically investigated the mechanism of stress development in subduction zones using the dislocation model. They found that extensional stress develops above the 2-D plate interface depending on its curvature, which is normally convex upward. Away from the plate interface, stress fields in their solution are antisymmetric about the mid-plane as elastic flexure problems (Watts 2001). In our results, however, PAC and PHS plate subduction is simply extensional at any depth away from the trench (Fig. 5c, d), which may be a result of the curved 3-D geometry of the plate interface. As Eq. (1) shows, the stress fields consist of the integrated effects of slip distributed over the plate interface. In purely 2-D plate geometry cases, the near-and far-field effects of slip over the plate interface cancel out, leaving a small vertical stress variation. In the case of 3-D plate geometry, however, the two effects do not cancel out, and the near-field effect of slips from nearby areas (within ~ 100 km), which are presumably extensional, may dominate. This lateral variation may be stronger than the subtle vertical variation that appears in 2-D cases, resulting in the stress fields depicted in Fig. 5c, d.