Quantifying the expansion rates of aftershock zones for magnitude-7 class earthquakes around the Japanese archipelago

Earthquakes (mainshocks) trigger sequences of aftershocks, the frequency of which diminishes following a power-law decay, while the spatial domain of these aftershocks extends logarithmically over time. The delineation of the after-shock zone can be modulated by variables beyond the magnitude of the mainshock, encompassing the location of the fault (whether the fault is at a plate boundary), the depth at which the event occurs, and the prevailing local stress conditions. Here, we evaluate the expansion rate of aftershock zones by analyzing earthquakes of magnitude-7 class in the vicinity of the Japanese archipelago. Prior studies have offered approximate assessments of expansion rates; however, our approach involves the utilization of a straightforward algorithm for the automated estimation of this metric, facilitating the compilation of a catalog. Across the dataset, no pronounced correlations were discerned between the expansion rate and other examined parameters. Yet, an inverse relationship is identified between the expansion rate and the b value of aftershocks for mainshocks occurring at plate boundaries. This observation suggests that the expansion rate of aftershock zones predominantly mirrors the stress field following the main-shock. Such a pattern is not detected in mainshocks occurring within the plate’s interior. While the expansion rate of aftershock zones is likely influenced by various factors, aftershock zones may expand more rapidly with higher differential stress in areas surrounding hypocenters of major interplate earthquakes of magnitude 8 or 9.


Introduction
When an earthquake (mainshock) transpires, significant forces are exerted on the adjacent regions, leading to the occurrence of aftershocks.Initially, the force emanating from the mainshock is disseminated to the surrounding area through both dynamic and static stress changes (King et al. 1994;Felzer and Brodsky 2006;Toda et al. 2012).It is well established that the magnitude of the mainshock is directly proportional to the extent of the aftershock region (Tajima and Kanamori 1985;Yamanaka and Shimazaki 1990).Subsequently, processes that are delayed in time prompt alterations in the locations of aftershocks.Instances of such processes include postseismic deformation subsequent to the mainshock (Matsuzawa et al. 2004;Perfettini and Avouac 2004;Hsu et al. 2006) or the localized diminution of rock strength facilitated by pore fluids (Nur and Booker 1972;Yukutake and Iio 2017;Miller 2020), both of which can engender aftershocks.The spatial domain of aftershocks undergoes expansion over time (Chatelain et al. 1983;Henry and Das 2001), accompanied by a power-law decay in the frequency of aftershocks (Omori 1894;Utsu 1961).
The expansion rate of the aftershock zone is subject to influence by various factors, including the location of the fault (whether the fault is at a plate boundary), ambient temperature and pressure conditions, and local stress conditions.Nonetheless, a consensus has emerged in recent years around the general characteristic that expansion is proportional to the logarithm of time (Peng and Zhao 2009;Kato and Obara 2014;Hainzl et al. 2016;Ross et al. 2017).An underlying physical mechanism that shows promise is the concept that the aseismic afterslip following the mainshock serves as an effective trigger for aftershocks (Wesson 1987;Kato 2007;Ariyoshi et al. 2019).Beyond the afterslip, various transient (viscous) phenomena, such as power-law creep, are posited to play a critical role immediately following megathrust earthquakes (Agata et al. 2019;Morikami and Mitsui 2020).Independent of the presence of afterslip or other transient mechanisms, the aftershock zone is anticipated to undergo expansion to a certain degree as a result of a sequential occurrence of aftershocks (Helmstetter and Sornette 2002;Ozawa and Ando 2021).As mentioned above, a variety of cases have been reported and models have been proposed.However, there remains a notable scarcity of research that examines and compares actual seismicity in a comprehensive manner.
In this study, we examine the expansion rate of the aftershock zone and its determinants by analyzing magnitude-7 class earthquakes occurring in the vicinity of the Japanese archipelago, an area characterized by its advanced seismic observation networks.Initially, we perform an objective estimation of the aftershock zone expansion rate and compile a comprehensive catalog.Building on this groundwork, we proceed to investigate the controlling factors driving aftershock generation.Our exploration extends beyond the established characteristics of aftershocks, such as power-law decay and logarithmic spatial expansion, aiming to uncover new insights into the complex natural systems that govern seismic events.

Data
We utilize hypocenter catalog data from the Japan Meteorological Agency (JMA), selectively incorporating only those entries flagged with "high precision" in hypocenter determination.Our dataset comprises earthquakes of JMA magnitude 2.0 or greater, occurring within the timeframe of 2000-2021.During this period, major interplate earthquakes of magnitude 9 (the Tohoku earthquake in 2011) and magnitude 8 (the Tokachi-oki earthquake in 2003) occurred; however, due to the vastly different scale of aftershock activity associated with these events compared to smaller earthquakes, they are not categorized as "mainshocks" in this study.Additionally, aftershocks occurring within one day of these events are also excluded from being classified as "mainshocks." Among the remaining earthquakes, the largest recorded a magnitude of 7.4.
Earthquakes with JMA magnitudes ranging from 7.4 to 6.9, occurring around the Japanese archipelago between 2000 and 2021, are identified as potential "mainshocks" for the analysis of the aftershock zone expansion rate.From the perspective of balancing sufficient signal strength and a substantial population size, this study focuses on earthquakes down to a magnitude of 6.9.In this context, earthquakes in the vicinity of the Japanese archipelago can be categorized into subduction interplate earthquakes, which exhibit low-angle reverse faulting akin to the angle of the subducting oceanic plate, and intraplate earthquakes, including inland and intraslab events, as delineated by the centroid moment tensors (refer to Supplementary Table S1).In our categorization, we also refer to isobaths of the subducting plate and earthquake reports from the Headquarters for Earthquake Research Promotion.Due to the inferior quality (significant estimation errors) of the hypocenter catalog data, outer-rise earthquakes, which occur at a considerable distance from the seismic networks, are excluded from our study.
Consequently, we designate the 25 earthquakes depicted in Fig. 1 and Table 1 as "mainshocks." Among these, the magnitude 7.3 earthquake in 2011, cataloged as mainshock number 14 in Table 1, demands particular examination.Occurring approximately 2 days before the magnitude 9.0 Tohoku earthquake, this event is posited as a potential foreshock to the latter, more significant seismic occurrence (Hirose et al. 2011).In light of this, our analysis is confined to the examination of aftershocks that occurred prior to the Tohoku earthquake, with the aim of avoiding the influence of aftershock activity associated with the magnitude 9 event.In Table 1, mainshock numbers 3 and 4, with magnitudes of 7.1 and 7.4 respectively, occurred in close proximity to each other within 1 day, and their aftershock zones overlap (Park and Mori 2005).In our analysis for this study, both events are utilized; however, readers conducting their own analyses may consider excluding one of the two to avoid potential confounding effects due to the overlapping aftershock zones.In the case of the Kumamoto earthquake (magnitude 7.3) in 2016 (identified as mainshock number 22 in Table 1), prior indepth studies (e.g., Uchide et al. 2016;Yagi et al. 2016) have documented the active dynamic triggering of distant earthquakes at locales unassociated with the fault of the mainshock.This paper treats it similarly to other mainshocks, yet a cautious approach may be warranted when discussing it in detail.
We investigate the spatiotemporal dynamics of aftershocks, defined as earthquakes occurring in proximity to a "mainshock." Despite the spatial diversity observed in aftershock distributions, we utilize the three-dimensional distance from the "mainshock" hypocenter as an isotropic spatial coordinate to enable a cohesive analysis.This approach facilitates a standardized discussion across different seismic events.The spatiotemporal definition of aftershocks is not straightforward.In this study, we initially conduct our analysis using definition commonly employed in existing research.Subsequently, we also evaluate the implications of slightly altering the definition.

Estimation of the expansion rate of aftershock zone
The expansion rate of the aftershock zone has traditionally been deduced from hypocenter data through visual inspection.However, in the present study, this rate is quantitatively estimated through a regression analysis, incorporating two variables: the logarithm of the elapsed time since the mainshock and the three-dimensional (hypocentral) distance from the mainshock hypocenter.
The definition of the aftershock zone is conducted in accordance with the well-established methodology in previous studies (Gardner and Knopoff 1974;van Stiphout et al. 2012).Specifically, the temporal window of the range is defined by the elapsed time t [days] since the mainshock, and the spatial window is defined by the distance d [km] from the hypocenter of the mainshock, as delineated by the following equation: where M is the magnitude of the mainshock.In defining the aftershock zones for the 25 mainshocks examined in this study, the temporal parameter t is set between 911.4 and 945.6 days, and the spatial parameter d ranges from 68.7 to 79.3 km.Moreover, in light of the proximity between certain mainshocks (Table 1), the implications of utilizing shorter temporal interval ( t ) will be examined.
Figure 2 presents an illustration of the expansion of the aftershock zone and the calculation of the expansion rate for the mainshock numbered 10 in Table 1, namely the 2008 off Ibaraki earthquake (Kubo and Nishikawa 2020).Within Fig. 2a-d, it is evident that the aftershock zone expands over time.To quantify this expansion rate, a scatterplot depicting the logarithm of elapsed time versus the hypocentral distance (Fig. 2e) is utilized, followed by a simple linear regression analysis of the points along the outer boundary.The regression analysis focuses on aftershocks at the greatest hypocentral distances as time progresses.It is important to note that the points defining the outer boundary are determined based on the spatiotemporal distribution of the aftershock hypocenters, without the imposition of any artificial time window for their selection.We present more information for all the 25 mainshocks, in Supplementary Data S1.

Estimation of Gutenberg-Richter b value
We focus on the Gutenberg-Richter b value (Gutenberg and Richter 1944) as a potential indicator of local stress associated with the expansion of the aftershock zone.The Gutenberg-Richter law posits that, for any given magnitude M , the number of earthquakes exceeding this magnitude, denoted as N , adheres to the expression log 10 N ∝ −bM .This b value demonstrates a tendency to negatively correlate with differential stress (Scholz 1968;Wyss 1973), and has been extensively employed in the research about the stress conditions on faults (Schorlemmer et al. 2005;Nanjo et al. 2012;Spada et al. 2013;Chiba 2020).
(1) t = 10 0.032M+2.7389 (2) d = 10 0.1238M+0.983Fig. 1 Magnitude-7 class "mainshocks" around the Japanese archipelago.The beachballs display the double couple part of the centroid moment tensors estimated by JMA.Interplate and intraplate earthquakes, which can be determined from the centroid moment tensors (Supplementary Table S1), are shown in red and blue.The vague gray beachballs show the outer-rise earthquakes, which were excluded from our study.The light pink lines depict the trench and the isobaths of the subducting Pacific Plate, plotted at 20-km intervals up to a depth of 80 km (Iwasaki et al. 2015;Lindquist et al. 2004).The violet lines represent those of the subducting Philippine Sea Plate Addressing the magnitude completeness problem in earthquake catalogs (e.g., Ogata and Katsura 1993;Nanjo et al. 2010;Mignan 2012), the recent development of the b-positive method (van der Elst 2021) marks a significant advancement in accurately estimating b values.The original Gutenberg-Richter law can be expressed as a probability density function for earthquakes exceeding a threshold magnitude M c , which is influenced by the conditions of seismic observation, as follows: where x = M − M c and β = bln(10).M c is generally not uniform across space and time; for instance, early aftershocks tend to be obscured by overlapping seismic waves (Kato et al. 2012;Omi et al. 2013;Sawazaki and Enescu 2014).The b-positive method (van der Elst 2021) utilizes the difference in magnitude, denoted as M ′ , rather than the absolute magnitudes.Namely, it employs the magnitude differences ( M ′ ≥ 0.1) between successive (3) f (x) = βexp(−βx) for x ≥ 0 earthquakes where the magnitude of the second earthquake exceeds that of the first.This approach culminates in the following relationship: The b value for each aftershock sequence can be estimated based on the maximum likelihood method (Aki 1965), yielding the estimated b value as follows: where M ′ is the mean magnitude difference and 0.05 represents the resolution of the catalog (Guo and Ogata 1997).We present histograms of the magnitude differences for all the 25 mainshocks, in Supplementary Data S2.The 90% confidence interval is obtained through the analysis of 500 bootstrapping samples. (4)

Results
Figure 3 and Table 2 present the estimated expansion rates of the aftershock zones for the 25 magnitude-7 class mainshocks, as depicted in Fig. 1 and Table 1.The spatial distribution and the triangle diagram of focal mechanisms (Frohlich 1992), as illustrated in Fig. 3, do not reveal any discernible patterns in the expansion rates of aftershock zones.Nonetheless, from a physical perspective, it is plausible to conjecture that the occurrence of aftershocks may be influenced by factors such as the location of the fault (whether the fault is at a plate boundary), the ambient temperature and pressure conditions, or the local stress conditions.Consequently, we undertake a statistical examination of these aspects.Initially, we examine the presence of significant disparities in the expansion rates of aftershock zones between interplate and intraplate earthquakes, as delineated in Table 1.The average expansion rate for interplate earthquakes is determined to be 12.46 [km/ decade], in contrast to 12.79 [km/decade] for intraplate earthquakes, with "decade" in this context referring to a single digit in the base-10 logarithm.No significant difference exists between the two types of earthquakes.
Subsequently, we shift our attention to the focal depth of the mainshock as a parameter approximately indicative of the ambient temperature and pressure conditions, and investigate its correlation with the expansion rate of the aftershock zones.The correlation coefficient is calculated to be − 0.28, implying that the depth at which the event occurs does not exert a controlling influence on the expansion of aftershock zones.
Finally, we establish a correlation between the seismic b value of the aftershocks (Table 2), as an indicator of local differential stress (Scholz 1968(Scholz , 2015)), and the expansion rate of the aftershock zones, depicted in  However, in Fig. 4, when distinguishing between mainshocks of interplate earthquakes and those of intraplate earthquakes, a certain characteristic becomes apparent.Namely, when isolating only the interplate earthquakes, a clear negative correlation is observed between the expansion rate of the aftershock zones and the b value.The correlation coefficient was found to be − 0.75, with a p value of 0.01.Mainshocks No. 13 and No. 14 hold significant importance.This negative correlation implies a positive correlation between the aftershock expansion rate and local stress after the mainshock, because the b value is negatively correlated with local differential stress (Schorlemmer et al. 2005;Nishikawa and Ide 2014;Scholz 2015).It would be mechanically plausible when considering the simpler aftershock propagation along the plane-like mature faults at plate boundaries.Additionally, aseismic afterslip or other transient mechanisms mentioned in the introduction section might contribute to this simpler propagation pattern.
In the relationship between the aftershock expansion rate and the b value for interplate earthquakes (Fig. 4), mainshocks No. 13 and No. 14 are characterized by their high expansion rates and low b values.Notably, the hypocenters of these mainshocks are located in close proximity to the larger magnitude 8 and 9 earthquakes, compared to other mainshocks, suggesting significant implications.For example, the 2008 event (No. 13) at longitude 144.152°, latitude 41.776°, depth 30.86 km, was near the magnitude 8.0 2003 Tokachi-oki earthquake (longitude 144.079°, latitude 41.779°, depth 45.07 km), and the 2011 event (No. 14), at longitude 143.280°, latitude 38.329°, depth 8.28 km, occurred near the magnitude 9.0 2011 Tohoku-oki earthquake (longitude 142.861°, latitude 38.104°, depth 23.74 km), potentially serving as a foreshock, as previously mentioned.This implies that the areas surrounding the hypocenters of major earthquakes of magnitude 8 or 9 tend to have higher differential stress compared to other areas.Such an effect is less apparent in the intraplate mainshocks, where their focal mechanisms and the plate boundary slip substantially differ.
The persistence of high stress in the area surrounding the hypocenter 5 years after the magnitude 8.0 2003 Tokachi-oki earthquake aligns with a previous study, suggesting that the stress release of this earthquake did not significantly alter the stress field (Tormann et al. 2015).Although numerous cases have been documented where b values increase in surrounding areas following large earthquakes (Gulia and Wiemer 2019), the 2003 Tokachioki earthquake appears not to conform to this pattern.To further explore this issue, additional investigation might be warranted using statistical models such as the epidemic-type aftershock sequence (ETAS) model, properly defining background seismicity distinct from aftershocks (Ueda and Kato 2023).

Discussion
We examine the definition of the aftershock zone that has been established (Gardner and Knopoff 1974;van Stiphout et al. 2012).Given the high background seismicity surrounding the Japanese archipelago, the current definition of the time window (ranging from t =  (Frohlich 1992), to demonstrate the estimated expansion rates 911.4 to 945.6 days) from Eq. ( 1) does not preclude the possibility of incorporating earthquakes unrelated to aftershock activity within this timeframe.
Taking into account the results obtained thus far (Supplementary Data S1), where the time parameter ranged from t = 911.4 to 945.6 days, we modified this parameter to t = 10 days in order to conduct an analysis focused exclusively on the early aftershocks.The duration of 10 days corresponds to an order of magnitude greater than the one-day period previously examined as early aftershocks in previous studies (Omi et al. 2013(Omi et al. , 2014)).
Figure 5 presents a comparison between the estimated results of aftershock zone expansion rates and b values with the time parameter set to t = 10 days and the results from the original time parameters.The correlation coefficients are sufficiently high at 0.71 and 0.95, respectively, indicating that the estimated values do not significantly change even when the time parameter t is altered by two orders of magnitude.Sup- plementary Data S3 and S4 show results analogous to those in Supplementary Data S1 and S2, with the time parameter t adjusted to 10 days.
We examine the correlation with the various parameters upon adjusting the time parameter to t = 10 days.Despite the two orders of magnitude change in the value of t , there are no fundamental changes to the findings discussed in the results section; no significant correlation is observed between the aftershock zone expansion rate and either the depth of the mainshock or b value.When the analysis is limited to the interplate earthquakes, the correlation coefficient between the aftershock zone expansion rate and b value slightly decreases to -0.60 (from -0.75), which is still meaningful.
Furthermore, we address the results obtained when using another scaling relations (Uhrhammer 1986) that involve t [days] and d [km], which are different from those in Eqs. ( 1) and ( 2), in the definition of aftershocks.The equations are given by: For the 25 mainshocks in this study, the temporal parameter t is set between 284.7 and 527.9 days, and the spatial parameter d ranges from 92.2 km to 137.8 km.Supplementary Table S3 compiles the results for the      2), yet there are no fundamental changes to the findings.When focusing solely on interplate earthquakes, the correlation coefficient between the aftershock zone expansion rate and the b value stands at − 0.79.Finally, we check the impact of the method used to estimate the b value.We have estimated the b value using Eq. ( 5), which means that it has been estimated within the range of M′ ≥ 0.1 .The previous study (van der Elst 2021) recommended verifying a minimum magnitude difference greater than 0.1; therefore, we attempt to estimate the b value at M′ ≥ 0.2 .As a result, the estimated b values differed by a maximum of 0.04 and less than 0.01 in more than half of the mainshocks compared to the original values (Table 2).Therefore, this point does not affect the main conclusions of this paper.
While there is ongoing potential for refinement in our methodology for individual mainshocks, we propose that we have effectively discerned shared characteristics among aftershocks through the utilization of a straightforward approach.

Conclusions
We conducted a systematic investigation into the expansion rates, measured in logarithmic time, of aftershock zones associated with 25 mainshocks of magnitude-7 class (Mj 7.4-6.9)occurring in the vicinity of the Japanese archipelago.When limited to interplate earthquakes, a negative correlation was observed between the expansion rate of the aftershock zones and the b value of the aftershocks.This reflects a tendency for aftershock zones to expand more rapidly in areas with higher differential stress, such as areas surrounding the hypocenters of major interplate earthquakes of magnitude 8 or 9.

Fig. 2
Fig. 2 Expansion of the aftershock zone and estimation of the expansion rate for the 2008 off Ibaraki interplate earthquake (magnitude 7.0).Panels a-d show the spatiotemporal distribution of the aftershocks by elapsed time from the mainshock.Panel e exhibits the scatterplot of the logarithm of the elapsed time versus the hypocentral distance.The pink circles denote the temporal progression of the locations of aftershocks at the greatest hypocentral distances (outer boundary).The slope of the dashed red line represents the estimated expansion rate of the aftershock zone, whereas the range of error is indicative of the 90% confidence interval

Fig. 4 .
Fig. 4. The correlation coefficient is -0.19, indicating no significant correlation.However, in Fig.4, when distinguishing between mainshocks of interplate earthquakes and those of intraplate earthquakes, a certain characteristic becomes apparent.Namely, when isolating only the interplate earthquakes, a clear negative correlation is observed between the expansion rate of the aftershock zones and the b value.The correlation coefficient was found to be − 0.75, with a p value of 0.01.Mainshocks No. 13 and No. 14 hold significant importance.This negative correlation implies a positive correlation between the aftershock expansion rate and local stress after the mainshock, because the b value is negatively correlated with local differential stress(Schorlemmer et al. 2005;Nishikawa and Ide 2014;Scholz 2015).It would be mechanically plausible when considering the simpler aftershock propagation along the plane-like mature faults at plate boundaries.Additionally, aseismic afterslip or other transient mechanisms mentioned in the introduction section might contribute to this simpler propagation pattern.In the relationship between the aftershock expansion rate and the b value for interplate earthquakes (Fig.4), mainshocks No. 13 and No. 14 are characterized by their high expansion rates and low b values.Notably, the hypocenters of these mainshocks are located in close proximity to the larger magnitude 8 and 9 earthquakes, compared to other mainshocks, suggesting significant implications.For example, the 2008 event (No. 13) at longitude 144.152°, latitude 41.776°, depth 30.86 km, was near the magnitude 8.0 2003 Tokachi-oki earthquake (longitude 144.079°, latitude 41.779°, depth

Fig. 3
Fig. 3 Estimated aftershock zone expansion rates for the magnitude-7 class mainshocks (Fig. 1).The left figure represents the spatial distribution, while the right figure shows the triangle diagram of focal mechanisms(Frohlich 1992), to demonstrate the estimated expansion rates d = exp(0.804M− 1.024) ) a In proximity to each other within one day with overlapped aftershock zones b With only aftershocks up to about two days before the magnitude 9 Tohoku earthquake in 2011 c With active dynamic triggering of distant earthquakes d With only aftershocks up to the end of 2021 (due to catalog limitation)The expansion rate is derived from the regression analysis on the logarithm of the elapsed time and the hypocentral distance.The b value (aftershocks) represents the b-positive estimation results for aftershocks.The error ranges represent the 90% confidence intervals

Fig. 4
Fig. 4 Scatterplot of the aftershock zone expansion rate versus the seismic b value of the aftershocks.The range of error is indicative of the 90% confidence interval

Fig. 5
Fig. 5 A comparison between the estimated results of aftershock zone expansion rates and b values with the time parameter set to t = 10 days and the results from the original time parameters.The left figure displays the results for the aftershock zone expansion rates (correlation coefficient = 0.71), while the right figure represents the outcomes for the b values (correlation coefficient = 0.95).The range of error is indicative of the 90% confidence interval

Table 1
Magnitude-7 class "mainshocks" around the Japanese archipelago in this studyThe longitude, latitude, and depth of hypocenters and the magnitude are based on the earthquake catalog provided by JMA.The type of earthquakes, interplate ("inter") or intraplate ("intra"), is determined from the centroid moment tensors estimated by JMA

Table 2
Estimated parameters for each mainshock (Table1