Inelastic strain rate in the seismogenic layer of Kyushu Island, Japan

Seismic activity is associated with crustal stress relaxation, creating inelastic strain in a medium due to faulting. Inelastic strain affects the stress field around a weak body and causes stress concentration around the body, because the body itself has already released stress. Therefore, the understanding of inelastic deformation is important as it generates earthquakes. We investigated average inelastic strain in a spatial bin of Kyushu Island, Japan, and obtained the inelastic strain rate distribution associated with crustal earthquakes, based on the analysis of fault plane solutions and seismic moments. Large inelastic strains (>10−7 year−1) were found in the Beppu–Shimabara area, located in the center of Kyushu Island. The strain rate tensor was similar to that of the stress tensor except the absolute value in the area, implying that the inelastic strain was controlled by the stress field. The 2016 Kumamoto earthquake sequence (maximum magnitude 7.3) occurred in the Beppu–Shimabara area, with the major earthquakes located around the high inelastic strain rate area. Inelastic strain in the volume released the stress. In addition, the inelastic strain created an increment of stress around the volume. This indicates that the spatial heterogeneity of inelastic strain might concentrate stress.Graphical abstract Inelastic strain rate in the seismogenic layer of Kyushu Island, Japan. Inelastic strain rate in the seismogenic layer of Kyushu Island, Japan.


Background
Stress and strain fields in the crust are important parameters for understanding earthquake activity and tectonic formation.Recent seismic ground and global navigation satellite system (GNSS) observations provide us with high-quality data to investigate stress and strain rates in the crust.
GNSS data provide the total strain rate at the ground surface, which consists of elastic and inelastic strain.Earthquake activity contributes to the inelastic deformation in a seismogenic zone.Kostrov (1974) proposed a method to determine the average strain rate in a region by summing the moment tensors released in the affected volume.Many researchers have applied this method to analyze seismic activity data in investigations of crustal deformation and tectonic evolution (e.g., Wesnousky et al. 1982;England and Molnar 2005).We also consider earthquakes to be inelastic phenomena in a medium, with a large number of earthquakes in a volume ΔV creating inelastic strain in that volume.In this study, we discuss inelastic strain produced by earthquakes and estimate the inelastic strain field of Kyushu Island, Japan, using focal mechanism and seismic moment tensor data.Kyushu Island is located in southwestern Japan and lies on the overriding plate above the subducting Philippine Sea slab (PHS) to the west of the Nankai Trough subduction zone (see inserted map in Fig. 1a).
A series of large earthquakes, with a maximum magnitude of 7.3, occurred in the center of Kyushu Island since April 14, 2016 (termed the 2016 Kumamoto earthquake by Shimizu et al. 2016).The earthquakes occurred in a seismically active region of Kyushu known as Beppu-Shimabara (Unzen), hereafter termed the B-U area, which is characterized by active volcanoes (Aso, Unzen, Kuju, and Tsurumi) and bounded by a shear zone (the western extension of the median tectonic line) to the south.The crustal structure in this region was found to be complex (Saiga et al. 2010), with a stress field different from other areas of Kyushu (Matsumoto et al. 2015).In this study, we discuss the relationship between the inelastic strain field and the Kumamoto earthquake sequence.

Seismogenic zone strain inferred from seismic moments
In this study, we consider a region with many faults displaced by earthquakes loaded with triaxial stress, with maximum (σ 1 ), moderate (σ 2 ), and minimum (σ 3 ) principal stresses.If there are many source faults in a region, and they are all smaller than the size of the region, deformation may be treated as inelastic from a macroscopic point of view.Recently, Matsumoto (2016) discussed the relationship between a deviatoric stress field and seismic moment tensors through inelastic deformation.We follow their formulation to demonstrate a method to estimate inelastic deformation in a volume ΔV.We consider a body with many earthquakes in ΔV, which is deformed by earthquake faulting that behaves as a volumetric source in the media surrounding ΔV.Moment tensors distributed in the volume express the inelastic deformation (Aki and Richards 1980;Ichihara et al. 2016;Kusakabe et al. 2016), and the relationship between inelastic strain ε p ij in the volume and the moment density tensor may be expressed using the compliance C ijkl of the volume (Noda and Matsu'ura 2010): where m kl is the moment density tensor.The term ε p ij corresponds to the increment of inelastic strain in the medium after the events occurring in ΔV (i.e., the increment of strain).This inelastic strain is identical to the plastic strain in a medium because earthquake faulting is macroscopically a plastic phenomenon.
Assuming a uniform distribution of earthquake point sources in the volume, we define the average moment density tensor in the volume as where the seismic moment tensor of the kth earthquake is M k ij and K is the number of earthquakes.Then, we obtain The basis of this formula was proposed by Kostrov (1974) to evaluate strain in a volume.Equation (3) reveals that the average strain increment may be estimated by summing the seismic moment tensors in a volume.The formula assumes that the distribution of inelastic strain in the volume is homogeneous after the occurrence of earthquakes for the target period.The homogeneous distribution of inelastic strain is equivalent to the homogeneous moment density in the volume.If only one large earthquake occurs in the volume, the high moment density tensor is distributed along the earthquake fault.In this case, the assumption of homogeneous strain in the volume might not be appropriate.Therefore, the spatial distribution of the moment tensor needs to be considered when evaluating the applicability of the method.

Inelastic strain in Kyushu, Japan
In this study, we analyzed a focal mechanism dataset from Matsumoto (2016) that comprised of selected data from Matsumoto et al. (2015) for shallow earthquakes with magnitudes greater than 1.0 that occurred at depths ranging from 0 to 30 km.The data cover the period from January 2000 to July 2013, during which the detectability would be expected to be constant.Saiga et al. (2010) and Hori et al. (2006) determined earthquake hypocenters in this area using three-dimensional velocity structures.The focal mechanisms were estimated from P-wave polarity data observed at eight or more stations using the HASH algorithm developed by Hardebeck and Shearer (2002).Only the focal mechanisms characterized by a low misfit in the algorithm were used (see Matsumoto et al. 2015). (1) (2) Figure 2 indicates the orientation of compression (P) and tension (T) axes for earthquakes at depths shallower than 30 km.
The method used in this study requires information regarding the magnitude of seismic moment tensors in addition to focal mechanisms.Matsumoto (2016) adopted seismic moment (M 0 ) data by Iio et al. (2006) and the moment tensor (MT) solution catalog from the National Institute of Earth Science and Disaster Resilience (NIED).For events not found in the moment tensor magnitude dataset, scalar seismic moments were inferred from empirical relationships between the scalar moment and magnitude estimated short-period velocity seismogram (M a ) as shown in Fig. 1b.A linear regression line fitted to the data gave the equation log(M 0 ) = 1.15M a + 10.548.
The strain rate tensors were estimated at grid points set at latitudinal and longitudinal intervals of 0.15°.The sum of the seismic moment tensor in the bin volume ΔV (i.e., 0.15° × 0.15° × seismogenic zone thickness) was calculated where more than ten focal mechanisms were obtained.The thickness of the seismogenic layer was estimated from the data used in Matsumoto et al. (2015), which incorporated 40,981 hypocenters.We adopted a depth range for the thickness of the seismic layer, based on the location where 90% of the events occurred; the layer is bounded by shallow and deep limits corresponding to 5% (D5) and 95% (D95) of the depth distribution of events.The D5 and D95 depth distributions for Kyushu Island are shown in Fig. 2. The D95 depth is shallow at Beppu and in southwestern Kyushu Island.Earthquakes occurred at crustal depths shallower than 2 km in the volcanic region, as shown in the D5 distribution; the D95 depth is also shallow around the volcanic region (Fig. 2).The general shallowing of the seismogenic zone around volcanoes might be due to weakening associated with the presence of magma and/or related volcanic fluids.The seismogenic zone is thin in the eastern part of the B-U zone.In contrast, seismic activity around Kumamoto extends to depths of 15-18 km and is scattered over a wide range, resembling an earthquake swarm.Thick seismogenic zones are also found in the aftershock area of the 2005 Fukuoka earthquake.In Fig. 3, the inelastic strain tensor is projected onto the lower hemisphere of the focal sphere beach ball diagrams, with the sphere size logarithmically proportional to the strain rate magnitude.The orientation of strain tensor contraction and extension axes are represented by focal mechanism P and T axes.The tension axis is generally oriented N-S in the center of the area, while it is NNW-SSE or NW-SE in the other regions.This tendency is close to the principal axis distribution of deviatoric stress obtained by Matsumoto et al. (2015).The stress field estimation used only fault plane solution data; therefore, the released scalar seismic moment could not be evaluated.The similarity between inelastic strain and stress fields reveals that large earthquake faulting occurs on the maximum shear plane under the regional stress field because large earthquakes strongly contribute to the inelastic strain estimation in this study.The estimated inelastic strain rates are greater than 10 −7 year −1 at several grid points (Fig. 3).The highest strain rate, found in northern Kyushu Island, was caused by coseismic deformation associated with the 2005 Fukuoka earthquake (M7.0).It is not appropriate to include this strain rate in the investigation, as a heterogeneous strain field with a scale length larger than the bin size would be expected around the main coseismic fault.The strain rate in the Kumamoto area, close to the hypocentral area of the 2016 Kumamoto earthquake, reveals a large strike slip deformation.The Beppu region is also categorized by a high strain rate, with inelastic deformation in a normal fault regime.These high values are comparable to the surface strain rates detected in GNSS studies (e.g., Nishimura and Hashimoto 2006;Wallace et al. 2009).
The inelastic strain pattern correlates with the behavior of active faults as shown in Fig. 3.The faults have been categorized as normal, left lateral, right lateral, and reverse by the National Institute of Advanced Industrial Science and Technology (AIST 2013).The Futagawa-Hinagu and Beppu-Haneyama fault zones in the Kumamoto and Beppu regions correspond to areas with high inelastic strain rates.Active faults also create inelastic deformation.Long-term inelastic deformation caused by active fault behavior is similar to the present inelastic deformation due to seismic activity.However, the inelastic deformation associated with small earthquakes in this study is not always located along fault planes but is distributed in a volume.This suggests that the deformation of the seismogenic zone occurs along active faults and also in other regions.
In the B-U area, the major inelastic strain tendency is for the tensor to have a non-double-couple component, which implies deformation by strike slip faulting as well as normal faulting.The amount of elastic strain in the normal fault component is almost the same as that in the strike slip component.This suggests that the formation in this area cannot be explained by the development of a simple graben structure.
According ratio ф, which is defined as (σ 2 -σ 3 )/(σ 1 -σ 3 ), in the B-U area reveals a uniaxial extensional stress.This requires a declining σ Hmax mechanism in the case of constant vertical principal stress.According to a geodetic study by Nakao et al. (2005), there is a large contraction in the E-W direction around the Aso and Kuju volcanoes that might be related to volcanic activity, indicating the possible existence of an anomalous mechanical structure as discussed by Matsumoto et al. (2015).This is a possible explanation of the declining σ Hmax , resulting in a uniaxial inelastic extension.

Discussion
We have demonstrated large inelastic strain distributed around the B-U area.In contrast, Nishimura and Hashimoto (2006) and Wallace et al. (2009) suggested that a shear zone, which might cause aseismic behavior, is required around the southern edge of the B-U area in Kyushu in order to explain the surface deformation detected by GNSS.However, Matsumoto et al. (2015) showed that earthquakes generated on the fault plane parallel to the shear zone, which would be indicative of a fault slip on the shear zone, are absent from the region.
One of the principal stress axes around the shear zone is oriented normal to it (Matsumoto et al. 2015), suggesting the shear zone behaves as a weak fault.These findings allow us to summarize the significant deformation characteristics of Kyushu Island: (1) a large inelastic deformation due to seismic activity in the B-U area and (2) a weak shear zone where aseismic slip might occur, located to the south of the B-U area, as shown in Fig. 4. The 2016 Kumamoto earthquake sequence involved several large earthquakes (Shimizu et al. 2016) in the B-U area, including along the Futagawa-Hinagu fault zone.In fact, the major slip of the largest event (M7.3, April 16, 2016) occurred on the Futagawa fault.However, other large earthquakes in the sequence were not restricted to this location.For example, an M6.5 earthquake occurred two days before the M7.3 event on a fault with a different fault strike.Large earthquakes also occurred just after the M7.3 event in the Beppu area about 100 km away.This type of complex seismic activity is a characteristic of the Kumamoto sequence.In order to evaluate the sequence with respect to the inelastic strain field, we calculated shear strain in the seismogenic zone.The apparent maximum shear strain rate is defined as (ε 1 − ε 3 )/2, where ε 1 and ε 3 are the maximum and minimum principal strain rates obtained from the estimated strain rate tensor, respectively.The apparent maximum shear strain rate distribution is plotted in Fig. 5.We also plotted the hypocenters of the events in the 2016 Kumamoto sequence of magnitude M5.5 or greater.High shear strain rates are seen in the B-U area, the northern part of Kyushu Island, and in both the western and eastern offshore areas.As described previously, it is not appropriate to use the strain rate calculated from the dataset incorporating large earthquakes to evaluate the results because it contravenes the assumption of the homogeneous moment tensor distribution.The major events of the Kumamoto sequence are located around the areas with high strain rates as shown in Fig. 5; specifically, their epicenters correspond to the edge of the higher shear strain rate areas.The strain rate is estimated from the moment tensor data acquired before the Kumamoto sequence commenced, so the strain rate distribution may relate to prior stress conditions.We can consider that the heterogeneous inelastic deformation creates stress concentration; a possible mechanism for this is shown schematically in Fig. 6.Stress in the area surrounding the volume, ΔV, is the sum of the regional stress and the perturbation caused by inelastic deformation in the volume.The stress tensor at position x, located outside of ΔV, is expressed as the sum of the regional stress tensor and the elastic stress tensor calculated by the moment tensor distribution in ΔV, through the representative theorem in elastodynamic theory (Aki and Richards 1980).The stress field becomes a function of the spatial position of the point, the inelastic strain in ΔV, and the regional stress loaded on the entire region.Stress concentration might appear at the edge of ΔV, as shown in Fig. 6.A large earthquake would more likely occur in the higher stress area than in other parts of the medium.Applying this to the case of the Kumamoto earthquakes provides an explanation for the pattern of large earthquake generation locations shown in Fig. 5.In addition, there is a correlation between the slip behaviors of the major events and inelastic deformation around the events.The M6.5 and M7.3 events were characterized by strike slip and normal faults, respectively.Inelastic strain rates around the M6.5 and M7.3 events also show deformations with left lateral and normal faulting.This implies that the behavior of major events in the Kumamoto sequence was affected by the orientation of the inelastic strain around the focal area prior to the event.

Conclusion
We estimated the inelastic strain rate in the seismogenic zone of Kyushu Island, Japan, using the focal mechanism dataset of shallow earthquakes.The result provided important characteristics for considering crustal deformation and development on Kyushu Island as follows: 1. High strain rates, greater than 10 −7 year −1 , were found in the Futagawa-Hinagu fault zone, Beppu area, and several other parts of Kyushu Island.2. Considering other factors of inelastic deformation like aseismic slip at the shear zone, we interpret crustal deformation in Kyushu as being composed of inelastic deformation due to seismic activity, aseismic slip, and elastic deformation.3. The major events of the 2016 Kumamoto earthquake sequence were located around an area with the highest spatial shear strain rate.This suggests that stress concentration occurred around the area prior to the commencement of the Kumamoto sequence.

Fig. 1 a
Fig. 1 a Distribution of P-and T-axis orientations for shallow earthquakes on Kyushu Island.Green and red indicate P-and T-axis dip angles of greater than 45°, respectively.Blue, green, and red axes denote strike, normal, and reverse faulting of the earthquakes, respectively.Purple lines indicate surface traces of the active fault.Yellow triangles indicate active volcanoes.The location map for Kyushu is shown in the lower left.b Relationship between the magnitudes of earthquake estimated from the amplitude of the velocity seismogram (M a ) and the scalar seismic moment (M 0 ).Blue symbols indicate values for earthquakes as determined by Iio et al. (2006) and the CMT catalog by NIED.Orange ovals represent the average moment value for each magnitude range at intervals of 0.2.The dotted line represents a linear regression line fitted to the data, and its formula is given in the text (after Matsumoto 2016)

Fig. 2
Fig. 2 Depth distributions of the Kyushu Island earthquakes, January 1993 to July 2013, showing a shallower (D5) and b deeper (D95) limits of the studied events.The depth scale is displayed at the bottom of each map.Solid circles indicate earthquake epicenters.The triangles denote active volcanoes

Fig. 3
Fig. 3 Inelastic deformation rate on Kyushu Island from January 2000 to July 2013.The tensor is plotted on the lower hemisphere of the focal sphere, and the symbol size is proportional to the logarithmic strain rate magnitude as shown in the legend.Active fault behavior, as determined by AIST, is denoted by colored lines: red, blue, green, and purple indicate normal, left lateral strike slip, right lateral strike slip, reverse, and unknown, respectively

Fig. 5
Fig. 5 Distribution of maximum shear strain rate on Kyushu Island.The rate is calculated from the data shown in Fig.3and is defined as (ε 1 − ε 3 )/2, where ε 1 and ε 3 are maximum and minimum principal strain rates, respectively.Red circles denote earthquakes of M5.5 or larger in the 2016 Kumamoto earthquake sequence.Blue circles indicate epicenters of events of M5.5 or larger that occurred in the strain rate estimation period.Triangles and purple lines show active volcanoes and active faults, respectively