Spatial variations in seismicity characteristics in and around the source region of the 2019 Yamagata-Oki Earthquake, Japan

The 2019 Mj\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${M}_{j}$$\end{document} 6.7 Yamagata-Oki earthquake occurred adjacent to the northeastern edge of the source region of the 1964 Mj\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${M}_{j}$$\end{document} 7.5 Niigata earthquake, offshore of Yamagata Prefecture, Japan. Few aftershocks occurred in the source region of the Yamagata-Oki earthquake immediately following the Niigata earthquake, and the recent seismicity rate in this region is low compared with the source region of the Niigata earthquake. This spatial variation in seismicity may allow us to elucidate plausible physical processes that shape the spatiotemporal evolution of these shallow-crustal environments. Here, we investigate the spatial variations in seismicity characteristics by applying the HIerarchical Space–Time Epidemic Type Aftershock Sequence (HIST-ETAS) model to an earthquake catalog compiled by the Japan Meteorological Agency for events in and around the Yamagata-Oki earthquake rupture region. We compare spatial variations in the background seismicity rate and aftershock productivity estimated from the HIST-ETAS model with the geophysical features in the study region. The background seismicity rate is high along the eastern margin of the Sea of Japan and correlates well with a previously identified zone that possesses a high geodetic E–W strain rate. The two major earthquakes occurred in and around a high E–W strain rate zone, suggesting that the background seismicity rate may serve as a key parameter for evaluating seismic hazard across the Japanese Archipelago. Furthermore, the source region of the Yamagata-Oki earthquake has a higher aftershock productivity and lower seismic-wave velocity than that of the Niigata earthquake. We interpret this low-velocity zone to be a well-developed damaged rock that resulted in an increase in aftershock productivity based on previous laboratory experiments and numerical results; this damage makes the rock more ductile at the macroscopic scale. The higher ductility in the source region of the Yamagata-Oki earthquake may have worked as a barrier against the propagation of dynamic rupture that occurred during the Niigata earthquake.


Introduction
The 2019 M j (Japan Meteorological Agency (JMA) magnitude) 6.7 Yamagata-Oki earthquake occurred off the coast of Yamagata Prefecture, Japan, at ~ 14 km depth, on 18 June 2019 (Fig. 1). The mainshock rupture caused strong shaking (JMA seismic intensity of 6 upper) in several towns in Yamagata and Niigata prefectures, and triggered a weak tsunami. The JMA-estimated centroid moment tensor of the mainshock rupture shows typical thrust faulting (Fig. 1) under an E-W horizontal compressive stress regime along the eastern margin of the Sea of Japan. The large historical earthquakes in this region possessed similar focal mechanisms (Abe 1975; Okamura et al. 2007). The Yamagata-Oki earthquake occurred adjacent to the northeastern edge of the source region of the 1964 M j 7.5 Niigata earthquake, which also showed thrust faulting (Abe 1975). There were few aftershocks Open Access *Correspondence: ueda@eri.u-tokyo.ac.jp in the source region of the Yamagata-Oki earthquake immediately following the 1964 Niigata earthquake (Kusano and Hamada 1991), and the recent seismicity rate has been low in this region compared with the source region of the 1964 Niigata earthquake (Fig. 1). However, the region with a low seismicity rate was filled by intense aftershocks after the Yamagata-Oki earthquake. These seismic observations have posed the following question: Which physical processes control the unique temporal and spatial evolution of seismicity in this region?
Previous studies have discussed the relationship between the observed seismicity and geophysical features in the seismogenic zone (e.g., Ogata 2004;Nandan et al. 2017;Zeng et al. 2018;Ben-Zion and Zaliapin 2019;Hainzl et al. 2019). For example, Zeng et al. (2018) demonstrated that spatial variations in the background seismicity rate in California correlate with the maximum shear strain rate, and Ogata (2004) suggested that aftershock productivity K (one of the parameters that define the modified Omori law) is high on the boundary of the source region hosting M 7-class earthquakes along the offshore Tohoku region. These previous studies imply that investigations of the spatial variations in given seismicity parameters (i.e., background seismicity rate, aftershock productivity, and aftershock decay rate) can allow us to infer the physical processes that shape the spatiotemporal evolution of these shallow-crustal environments. The unique setting of two adjacent large earthquakes, the 2019 Yamagata-Oki and 1964 Niigata earthquakes, provides us with an opportunity to understand the physical processes controlling the spatial variations in seismicity characteristics at crustal depths. Seismicity in the study region. a Epicenter distribution of earthquakes before the 2019 Yamagata-Oki earthquake. The yellow star is the epicenter of the 1964 Niigata earthquake. The black star is the epicenter of the M j 5.9 aftershock that occurred immediately after the 1964 Niigata earthquake (redetermined by Miyaoka et al. (2019)). Blue circles denote the aftershocks that occurred within 24 h of the 1964 Niigata earthquake. Gray circles are the earthquakes that occurred before the 2019 Yamagata-Oki earthquake . Gray lines are the surface traces of active Quaternary faults. b Epicenter distribution of the earthquakes that occurred before and within 24 h of the 2019 Yamagata-Oki earthquake. The black star, southern yellow star, blue and gray circles, and gray lines are the same as in a. The northern yellow star denotes the epicenter of the 2019 Yamagata-Oki earthquake. Red circles denote the earthquakes that occurred within 24 h of the 2019 Yamagata-Oki earthquake. Focal mechanisms were estimated by JMA (2019 Yamagata-Oki earthquake) and Abe (1975) (1964 Niigata earthquake). c Regional map of Japan, with the study region indicated by the blue square. d Magnitude-time plot of the earthquakes (1922-2019) that occurred in the study region. M j ≥ 1.8 earthquakes are shown in green. The two-way arrow indicates the time period used to estimate the HIST-ETAS parameters Here, we investigate spatial variations in a set of seismicity parameters by applying the HIerarchical Space-Time Epidemic Type Aftershock Sequence (HIST-ETAS) model (e.g., Ogata et al. 2003;Ogata 2004) to a JMA hypocenter catalog containing events in and around the Yamagata-Oki earthquake rupture area. The HIST-ETAS model enables us to model the observed spacetime changes in seismicity rate considering the location dependence of the seismicity parameters. We then compare the spatial variations in background seismicity rate and aftershock productivity estimated from the HIST-ETAS model with the geophysical features in the study region. We discuss spatial variations in the seismicity parameters, with a focus on the active tectonic zone and a relative abundance of ductility inferred from rock deformation mechanics.

Methods
We adopt the HIST-ETAS model (e.g., Ogata et al. 2003;Ogata 2004;Bansal and Ogata 2013) to fit the observed seismicity rate and investigate the spatial variations in seismicity characteristics. This model is a point-process model that includes the Omori-Utsu law (Utsu 1961;Utsu et al. 1995) which formulates the typical aftershock temporal decay, Utsu-Seki law (Utsu and Seki 1955) which represents the linear relationship between the logarithm of aftershock area and the mainshock's magnitude, and branching process, such that each earthquake, regardless of magnitude, has the ability to increase the probability of future seismic events (Iwata 2009).
The earthquake occurrence rate at time and location ( t, x, y) can be expressed by: The first term µ is the background seismicity rate, and the second term expresses the excitation rate of the aftershock occurrence by an earthquake at time and location ( t i , x i , y i ) with magnitude M i . K is the aftershock productivity; p is the aftershock temporal decay rate; α is the aftershock magnitude sensitivity; q is the aftershock spatial decay rate; The five seismicity parameters ( µ , K , α , p , and q ) are given as a function of space ( Fig. 2a-e), and are expressed as: where µ , K , α , p , and q correspond to the geometric mean value of each parameter averaged over the analysis region. We adapt each ϑ x, y value to the data by expressing each function using many coefficients, which locates each earthquake epicenter, and some additional points on the boundary of the analysis region. Each ϑ x, y value at an arbitrary location is linearly interpolated using the three values at the vertices of each Delaunay triangle (Fig. 2f). The two parameters c and d are constants for the sake of simplicity. S i is non-dimensional positive definite symmetric matrix for anisotropic clusters determined by identifying aftershock cluster following magnitude-based clustering algorithm Ogata 1998) and choosing best-fit ellipsoid which represents the cluster; H t is the history of occurrence times up to time t and their corresponding locations and magnitudes; M c is minimum magnitude used in the analysis. Usage of the unnormalized form of the Omori-Utsu law enables the p-value smaller than 1. The unknown parameters can be estimated via the maximum likelihood estimation method. The log likelihood is expressed as: where k is the number of events in the analysis, S is the analysis region, and [0, T] is the analysis interval. However, it is hard to estimate seismicity parameters stably because the number of unknown parameters is about five times larger than the number of events used in the analysis. We obtained stable solutions by penalizing the spatial gradient of the parameter functions under the assumption of using smoothed functions. We then estimated the seismicity parameters that maximized the penalized log likelihood function by objectively tuning the penalties of the optimal weights using a Bayesian procedure, which implies the optimal maximum posterior solution. The parameters and weights were alternately estimated. See Ogata et al. (2003) and Ogata (2004) for the details of the parameter estimation.
We used the JMA hypocenter catalog (the Preliminary Determination of Epicenters) as the earthquake catalog (which includes aftershocks of the Yamagata-Oki . We applied the HIST-ETAS model to the M j ≥ 1.8 earthquakes that occurred during the 1998-2019 period and which were located at ≤ 80 km depth and within the blue rectangular box in Fig. 1c. Here, we define the minimum magnitude to be 1.8 considering the goodness of fit proposed by Wiemer and Wyss (2000) and the temporary deficiency of aftershocks immediately following the 2019 Yamagata-Oki earthquake (Additional file 1: Fig.  S1). The earthquakes that occurred during the 1922-1997 period ( M j ≥ 1.8 , ≤ 80 km depth) were used as the precursory occurrence history of the HIST-ETAS model.

Results
The spatial distributions of the seismicity parameters that were estimated via the HIST-ETAS model in the present study are shown in Fig. 2. Both µ and K vary spatially ( Fig. 2a, b, respectively), whereas α , p , and q are almost constant ( Fig. 2c-e, respectively). It is clear that there is a local high in µ along the eastern margin of the Sea of Japan (Fig. 2a). The source region that hosted the Yamagata-Oki earthquake is located near the eastern boundary of a high-µ zone. Furthermore, the source region of the Yamagata-Oki earthquake possesses larger K-values than that of the Niigata earthquake (Fig. 2b). We estimated the seismicity parameter uncertainties at each location to evaluate the significance of the relative differences in each parameter between the source regions of the Yamagata-Oki and Niigata earthquakes. We first resampled the data 100 times by randomly extracting 90% of the earthquake data used in this analysis. We then applied the HIST-ETAS model to each resampled dataset and estimated the spatial patterns of the five seismicity parameters. We note that the use of randomly selected resampled data may introduce the possibility that the reference parameters ( µ , K , α , p , and q ) vary significantly among the 100 resampled datasets owing to the tradeoffs between each parameter, such that these variations may be larger than those from the spatial differences. Therefore, we normalized the model parameters by dividing the reference parameters estimated from each resampled dataset and then calculated the standard deviation of ϑ x, y /ln10 (common logarithm of the normalized parameter) at each location. The standard deviations of ϑ µ x, y /ln10 and ϑ K x, y /ln10 (Additional file 1: Fig. S2a and b, respectively) are quite small compared with the relative variations in log 10 µ and log 10 K between the source regions of the Yamagata-Oki and Niigata earthquakes (Fig. 2a, b, respectively). We obtain similar results for additional cases of the different random extraction rates (80% and 70%) (Additional file 1: Figs. S3a, b and S4a, b).
These results indicate that the relative variations in µ and K are significant. We do not discuss the other three parameters ( α , p and q ) in this paper because the relative variations in log 10 α , log 10 p , and log 10 q (Fig. 2c-e, respectively) are comparable or smaller than the standard deviations of ϑ α x, y /ln10 , ϑ p x, y /ln10 and ϑ q x, y /ln10 (Additional file 1: Fig. S2c-e, respectively), indicating that spatial variations in these three parameters are insignificant.

Discussion
We find significant spatial variations in background seismicity rate µ and aftershock productivity K (Fig. 2a, b, respectively). However, one may argue that the intensive seismicity associated with the Yamagata-Oki earthquake gives a bias for the spatial variations in aftershock productivity (especially large K-value in the source region of Yamagata-Oki earthquake). To verify the effect of Yamagata-Oki earthquake sequence, we estimate HIST-ETAS parameter using only earthquake catalog before the onset of the Yamagata-Oki earthquake (from January 1, 1998 to June 17, 2019). The spatial distributions of the seismicity parameters and their uncertainties are shown in Additional file 1: Figs. S5 and S6, respectively. The spatial variations in µ and K shown in Additional file 1: Fig. S5 are quite similar to those shown in Fig. 2, although absolute value of K is different owing to the trade-off between the parameters K and q . Therefore, we conclude spatial variations in µ and K are stable regardless of the occurrence of the 2019 Yamagata-Oki earthquake sequence.
In addition, we perform two kinds of synthetic tests to verify reliability of the spatial variations of the estimated parameters. A synthetic catalog is created based on the simulation algorithm proposed by Ogata (1998) by giving a spatial distribution of seismicity parameters. We use the same sequence of magnitudes and the same precursory occurrence history as the real earthquake catalog. Subsequently, the synthetic data are inverted applying the HIST-ETAS model. We then evaluate how much the estimated parameters are recovered, comparing with the initially given values.
First, we generate the synthetic catalog #1 using the seismicity parameters estimated from the real catalog (Fig. 2). Comparing the parameters inverted using the synthetic catalog #1 (Additional file 1: Fig. S7) with the initially given values (Fig. 2), the spatial variations of µ and K are well recovered, although absolute value of K is shifted owing to the trade-off between the parameters K and q.
Second, we created the synthetic catalog #2 assuming each parameter to be spatially uniform (reference parameters of our result). This aims to test whether the current method causes apparent spatial variation of each parameter during the inversion. From the distribution of the parameters inverted using the synthetic catalog #2, we recognize that the apparent variation in K is smaller than deduced from the real catalog (Fig. 2b). Based on the above two kinds of the synthetic tests, we conclude that the spatial variations in µ and K obtained from the present study is indeed reliable.
Previous studies have discussed the physical processes that control the background seismicity rate (e.g., Ide 2013; Zeng et al. 2018; Ben-Zion and Zaliapin 2019). Zeng et al. (2018) revealed that the maximum shear strain rate along the San Andreas Fault and Eastern California Shear Zone correlated with the distribution of the background seismicity rate in the shallow-crustal environment. Ide (2013) suggested that subduction zones around the world exhibit an approximately linear increase in the background seismicity rate with the relative convergence rate at each plate boundary. These previous studies suggest that the background seismicity rate can provide clues to understanding the evolution of deformation in the seismogenic zone. Therefore, we compared the spatial variations in background seismicity rate obtained by the current study with the spatial variations in geodetic E-W strain rate obtained by Meneses-Gutierrez and Sagiya (2016); the areas with high background seismicity rates correlate well with those possessing high E-W strain rates (contraction) (Fig. 3). The correlation coefficient between the E-W strain rate and the logarithm of the background seismicity rate is calculated to be − 0.49 ( Fig. 3c; the correlation coefficient between the E-W strain rate and the logarithm of the medium of µ (black dots in Fig. 3c) is calculated to be − 0.91).
The high E-W contraction area was a northern extension of the Niigata Kobe Tectonic Zone, stretching from Kobe to Niigata District, proposed by Sagiya et al. (2000). The source regions of both the Yamagata-Oki and Niigata earthquakes are located around and in the zone characterized by a high background seismicity rate and high E-W strain rate (Fig. 3). Nishimura (2017) suggested that the large shallow-crustal earthquakes in the Japanese Archipelago (M j ≥ 6) tend to occur in and around the areas with a high shear strain rate. Ogata (2017) suggested that many large earthquakes in California occurred in areas with high background seismic activities. Our results are consistent with these papers and imply that large shallow-crustal earthquakes are likely to occur in areas with high background seismicity rates. Therefore, capturing spatial variations in the background seismicity rate may assist in evaluating seismic hazard across the Japanese Archipelago. The background seismicity rate is especially effective in areas where a geodetic network has not been densely deployed, such as marine settings and some inland areas, to monitor the secular accumulation of elastic strain.
The relationship between the seismicity parameters that describe aftershock generation and the geophysical features in the seismogenic zone has been explored by previous studies (e.g., Mogi 1967;Ogata 2004;Marsan and Helmstetter 2017;Nandan et al. 2017;Zakharova et al. 2017;Hainzl et al. 2019). For example, Ogata (2004) suggested that aftershock productivity ( K ) is high on the boundary of each source region that hosts M 7-class earthquakes along the subduction zone, offshore NE Japan. The K-value in the source region of the Yamagata-Oki earthquake is larger than that in the source region of the Niigata earthquake, as shown in Fig. 2b. Here, we focus on spatial variations in the seismic-wave velocity and b-value to investigate the relationship between K and the observed geophysical features (Fig. 4). Figure 4a, b shows the compressional-wave (P-wave) velocity structures at depths of 10 km and 15 km (the depth layer close to the mainshock hypocenter depth, ~ 14 km) from a regional seismic tomography study (Matsubara et al. 2020). The P-wave velocity in the source region of the Yamagata-Oki earthquake is lower than that in the source region of the Niigata earthquake. It is noted that the shear-wave velocity structure at 15 km depth follows a similar trend to the P-wave velocity structure. But, the contrast of shear-wave velocity between the two source regions becomes weak at 10 km depth. This E-W seismic velocity contrast is also observed in the velocity model estimated by Zhao et al. (2015). We define two zones based on the seismic velocity structure: a low-velocity zone that includes the rupture area of the Yamagata-Oki earthquake (black rectangular region in Fig. 4a, b) and a high-velocity zone that includes the rupture area of the Niigata earthquake (blue rectangular region in Fig. 4a, b). We calculate the b-value in each region during the 1998-2019 periods (Fig. 4c) using the formula of binned magnitude (Tinti and Mulargia, 1987;Marzocchi and Sandri, 2003).We defined M1.9 as the cutoff magnitude considering the goodness of fit proposed by Wiemer and Wyss (2000), in addition to the temporary deficiency of immediate aftershocks following the 2019 Yamagata-Oki earthquake (Additional file 1: Fig. S1b). We evaluate the uncertainty of b-value using the bootstrapping method (1000 resampling). Figure 4c shows that the b-value in the low-velocity zone is lower than that in the high-velocity zone but given the uncertainty the difference in b-values is not so significant. It is quite difficult to compare b-value with P-wave velocity in higher spatial resolution, because the number of available earthquakes is not sufficient to obtain reliable b-value (high goodness of fit > 90%).
The source region of the Yamagata-Oki earthquake has a higher K-value and lower P-wave velocity than that of the Niigata earthquake (Fig. 5). Assuming that b-value in the source region of the Yamagata-Oki earthquake is lower than that of the Niigata earthquake, we could interpret the low-velocity zone of the Yamagata-Oki earthquake to be a well-developed damaged rock that contains many fractures and cracks over multiple scales. Laboratory experiments and numerical simulations of rock deformation have implied that the progressive damage associated with an increase in shear stress results in both a reduction in the b-value and an increase in aftershock productivity within a broad volume (e.g., Amitrano 2003). Furthermore, progressive damage makes the rock more ductile at the macroscopic scale because diffuse inelastic deformation is more dominant than localized deformation. Diffuse deformation provides a large spatial correlation dimension for damage, leading to a low bvalue and high productivity of brittle fractures within a broad volume (Amitrano et al. 1999;Amitrano 2003).
The ductility in the source region of the Yamagata-Oki earthquake is therefore considered to be higher than that of the Niigata earthquake. The relative increase in ductility may act as a barrier against the propagation of dynamic rupture during the Niigata earthquake. This is consistent with the case where the ductile flow of Fig. 4 Seismic-wave velocity structure and frequency-magnitude distribution. a Horizontal map of the P-wave velocity structure at 10 km depth, which was imaged by Matsubara et al. (2020). Black and yellow stars, blue circles, and gray lines are the same as in Fig. 1b. The black circles are the same as the red circles in Fig. 1b. b Horizontal map of the P-wave velocity structure at 15 km depth (Matsubara et al. 2020). See a for other details. c Cumulative frequency-magnitude distributions of the earthquakes that were recorded from 1998 and which occurred in the corresponding colored rectangular regions in a. The b-value is calculated using the formula of binned magnitude (Tinti and Mulargia 1987;Marzocchi and Sandri 2003). We define M1.9 as the cutoff magnitude considering the goodness of fit (Wiemer and Wyss 2000) and the temporary deficiency of aftershocks immediately after the 2019 Yamagata-Oki earthquake (Additional file 1: Fig. S1b). The uncertainty is the standard deviation estimated using the bootstrapping method thick sediments may have arrested the southwestward propagation of the dynamic rupture of the 2004 Niigata Chuetsu earthquake (Kato et al. 2009(Kato et al. , 2010. The crustal structure along the eastern margin of the Sea of Japan is quite complex owing to the evolution of a fold-and-thrust belt system under a W20°N-E20°S horizontal compressive stress regime (Okamura et al. 2007;Kato et al. 2009). It is therefore plausible that lateral heterogeneities in the crustal structure impacted the spatiotemporal pattern of the aftershock sequences.
The results presented suggest that different deformation styles play a key role in controlling seismicity characteristics. A systematic investigation of the spatial variations in seismicity characteristics and related geophysical features could provide new insights into the physics of earthquake generation.

Conclusions
We investigated the spatial variations in seismicity characteristics by applying the HIST-ETAS model to the JMA hypocenter catalog. The 1964 Niigata and 2019 Yamagata-Oki earthquakes occurred in and around a zone characterized by a high background seismicity rate and high E-W strain rate along the eastern margin of the Sea of Japan, suggesting that the background seismicity rate may be a valuable parameter for evaluating seismic hazard across the Japanese Archipelago. The source region of the Yamagata-Oki earthquake has a slower P-wave velocity and higher aftershock productivity than that of the Niigata earthquake. Differences in the macroscopic behavior of rock deformation may explain the spatial variations in seismicity characteristics.