Depth variation of seismic moment and recurrence interval in Japan

In repeating seismic event sequences within a specialized horizontal area, the moment magnitude is usually scaled with the recurrence interval. In addition to two horizontal dimensions, the vertical dimension plays a certain role in affecting the scaling law. However, whether and how the changing depth influences the scaling law remain enigmatic. Based on a large number of earthquake records with high-resolution epicenter locations in recent decades in Japan, we focus on a comparison between the 3-D seismic moment and seismic interval, which recognize the vertical dimension as the same dimension as the horizontal distances. The results show that (1) the seismic moment scaling law is applicable in the multiparameter 3-D models by visiting the 1.8 million events collected during a period of 15 years; (2) the vertical dimension of depth plays an important role in the Mo–SI relationship as well as in the variance in the 3-D seismic moment–interval magnitudes; and (3) the seismic moment rate, attributable to the plate convergence rate, varies with area and depth in influencing the regional earthquake recurrence frequency.


Introduction
The relationship between magnitude and cumulative earthquake numbers is understood, and the moment magnitude scale interprets the relationship between the seismic moment and the energy released by earthquakes (e.g., Kanamori 1983). Small repeating events are usually found after a large earthquake, such as the 1984 M6.2 Morgan Hill earthquake (e.g., Peng et al. 2005;Templeton et al. 2009), the 1989 M6.9 Loma Prieta event (e.g., Schaff and Beroza 2004), several subduction zone earthquakes in Japan (e.g., Uchida et al. 2004), and the 2004 M6.0 Parkfield earthquake (e.g., Taira et al. 2009;Lengline and Marsan 2009). Following the 2004 M6 Parkfield earthquake, many of the repeating events were associated with significantly reduced recurrence intervals that systematically increased with time (e.g., Chen et al. 2010). The observed shortened recurrence intervals of a repeating earthquake sequence follow an Omori's law distribution in time, implying accelerated fault traveling along some faults surrounding the rupture (e.g., Peng et al. 2005;Taira et al. 2009). The systematic temporal variations in the characteristics of repeating earthquake sequences in Parkfield through ~ 20 years of observation  indicate that they are mostly consistent with changes in fault strength (Taira et al. 2009). Half of a population of 17 sequences at Parkfield exhibited similar systematic temporal variations in seismic moment (Mo) and seismic interval (SI) after the Parkfield earthquake (Taira et al. 2009).
However, most previous studies focused on the seismic moment-seismic interval (Mo-SI) relationship (e.g., Eq. 1) are based on a horizontally prescribed area, in which the effect exerted by the third dimension (depth) in the Mo-SI law remains enigmatic. Considering that 1.8 million events of full-range magnitudes cataloged are available from the Japan Meteorological Agency (JMA) from October 1, 1997, to February 28, 2015, and given the strength of our research in 3-D numerical simulation (e.g., Ji et al 2016), we attempt to perform a seismicity statistical study of the dozens of 3-D spherical spatial bins of clustered earthquakes recorded in Japan to investigate whether and how the vertical dimension could play a role in the Mo-SI relationship and to unravel more earthquake recurrence properties varying with area and depth in Japan.

Data selection
In this study, with regard to the continuously increasing high-resolution seismic records of recent decades, we select our study region in Japan between longitudes 132°E and 147°E and latitudes 32°N and 46°N and employ a set of 1.8 million events of full-range magnitudes cataloged by the Japan Meteorological Agency (JMA) from October 1, 1997, to February 28, 2015, to enable a deeper exploration of the Mo-SI spatiotemporal relationship based on three-dimensional (3-D) seismicity distribution statistics. Considering that the new seismic travel time table JMA2001 (Ueno et al. 2002) and unification of data from universities, municipalities, and other institutes was put into use from October 1997, we select seismic records after October 1997 that have used the updated travel time table. The magnitude range of the events is full range spanning from −1.4 to 7.3 in this study (Table 1). Based on the method, we focus on the detailed 3-D distribution of 1,842,626 events cataloged by JMA.

Seismically dense sphere method
The distances between every two arbitrary events are first calculated. For these events, all other events within a certain range of distances are also counted and understood. Then, every event is taken as the center of a sphere, around which we prescribed a radius of 10 km, to define a spatial range in which all of the adjacent events within 10 km (km) are counted. We call this sphere the seismically dense sphere because the sphere is a spatial range used to evaluate the seismicity activity density within it. A radius of 10 km is adopted because the range of error for hypocenter locations provided in the catalog is at a level of kilometers and the seismogenic fault rupture length for an M5.0 earthquake is approximately 4-5 km according to the seismic moment estimation. Every seismically dense sphere has an event at its center, and thus, we have 1.8 million seismically dense spheres. The seismically dense spheres with a significantly large event number are therefore counted and ranked. The so-called most seismically dense spheres (SDS) during the study period are shown in the map (Fig. 1), and the parameters are tabulated in Table 1. It is noted that the earthquakes' detection capability is heterogeneous in space and time, i.e., it is generally higher for the land region and gradually becomes lower seaward because most observation stations are located on land. Additionally, the detection capability becomes lower just after the occurrence of large earthquakes.
To avoid the presence of overlapping spheres, we set the lower limit of the distance between every two seismically dense spheres (from one center to another center) to be greater than 0.5 degrees (approximately 50 km). Provided two spheres are overlapping, the one with a larger event number will be kept, and the other will be ruled out. Using this method, we ultimately obtained 49 seismically dense spheres with the most clustered event distributions during the study period for this study region.

Maximum seismic interval method
The seismic interval (SI) used to discuss earthquake recurrence in this study is the time interval between two continuously recorded events, that is, the events with magnitudes of M ≥ M min (M min represents the minimum magnitude of events detectable in the catalog) in the seismically dense spheres. We focus on M < 4 (small-sized) earthquakes because there have been 71,789 (3.90%) M3 events, while the number of M4 events is only 12,553 (0.68%) and is easily contaminated by small events. We checked and found that 24.6% of the M3 events are uncontaminated within the interval, while 90% of M4 events are contaminated within the interval. To identify uncontaminated events, we observe that the maximum seismic interval (MSI) during a certain period (e.g., a day or a month) could potentially be used because the intervals of M3 events are at an average magnitude of days or 10 days. Due to the existence of pre-and postseismic repeating events (Uchida et al. 2004), the MSI within a month has proven to be a restraining factor regarding the noise of repeating events. Taking this into account, we use MSI as a proxy to measure the magnitude of the approximated real SI of M < 4 events. In addition, the M < 4 small-sized earthquakes compose 99.22% of all recorded events in this study involving the majority of seismicity, which potentially reflects the major properties of earthquake phenomena. Based on these assumptions, we calculate the seismic recurrence intervals and obtain the maximum interval during 1 day (day-based maximum interval, DMSI) or 1-month periods (month-based maximum interval, MMSI). The MMSI is filtered from the SI data using a time window threshold of 1 month (nearly 3.0 × 10 6 s) over the whole period, while the DMSI is defined within a time window of 1 day. Although both the DMSI and MMSI reflect the real SI, MMSI possesses higher resolution than DMSI. This enables us to further investigate the Mo-SI relation of regular repeating earthquakes without being affected by clustered seismic events and aftershocks.

Seismic moment and interval data source
The seismic moments of events are provided by JMA. Most of the recorded seismic magnitudes determined by JMA use the moment magnitude scale or the JMA magnitude scale, which is almost equivalent to a moment magnitude for a seismic event (Utsu 1983;Katsumata 1996). The intervals are calculated between two recorded events according to the JMA catalog. The large number of earthquake records observed within the high-density seismic observation network presents a good foundation for the study of the Mo-SI scaling law as well as further correlation analyses and investigations of the variance in Mo-SI magnitudes associated with changes in area and depth.

Results
The logarithms of the seismic moments and seismic intervals of successive events within the 49 seismically dense spheres are shown in Figs. 2, 3, 4, 5, 6, 7, 8. The distributions of the SI logarithm in most of the spheres (Figs. 2,3,4,5,6,7,8) have two typical but different stages (especially for the spheres located near the source region of the M9.0 Tohoku earthquake), which are possibly influenced by the most recent large earthquakes.
One is a sparsely distributed period that is mostly irregular (e.g., Fig. 2a before September 2001), and the other is a densely distributed period (e.g., Fig. 2a after September 2001 and lasting for almost 10 years until 2011). In the densely distributed period, the variation in the SI values in its upper part coincides with the variation in the recorded seismic moment (Figs. 2a,3,4,5,6,7,8).
The majority of the seismic frequencies that transition from a sparsely to densely distributed period sometimes appears to be synchronized with seismic events (M > 5, only for a comparison with the peripheral large earthquakes) within or outside the spheres within a distance of 500 km. However, it is sometimes difficult to determine whether they are correlated, as depicted by the pink open circles in the upper part of Figs. 2, 3, 4, 5, 6, 7, 8. The seismicity within a distance of 500 km peripherally can possibly influence the occurrence of earthquakes inside the spheres, and sometimes the minimum seismic interval drops to a level of minutes, but the maximum seismic interval is not influenced by seismic swarms or repeating earthquakes. The monthly maximum SI is plotted as the cyan curve in Fig. 2a. The results reveal a correlation between the seismic magnitude variation (pink curve) and the monthly maximum SI variation (cyan curve) in the sampled spheres (e.g., Fig. 2a). A similar correlation for moment and monthly maximum SI can also be found in the 49 seismically dense spheres spanning a wide range of regions from southwest Japan to Hokkaido (Fig. 1), with a comparatively acceptable confidence level (Table 1). This pattern indicates that a potential law likely governs this phenomenon because a great number of events have been included in the analysis. For example, the dense seismic sphere SDS1 with the largest event number 52965 during the period of October 1, 1997, to February 28, 2015, has its center located at 37.740°N, 139.971°E with a depth of 8.65 km (the radii of the spheres are all 10 km in Table 1), adjacent to the Bandaishan volcano in the inland Fukushima Prefecture (Fig. 1). The hypocenter is shallow in the upper crust, and all of the repeating events occur shallower than 18.65 km  1 The numbers of the seismic events from October 1, 1997 to February 28, 2015, in a full-magnitude range occurring within each SDS according to the catalog of the JMA. The full-range magnitude of events is adopted for each SDS 2 The completeness magnitude Mc is largely influenced by the station density, and the use of the earthquakes with magnitudes of 1.9 or larger above possibly covers a complete catalog (Nanjo et al. 2010). The b-value estimation in Table 1 includes M < Mc events and has an uncertainty due to the event detectability 3 The parameters of Eq. 1 in the main text. The variable k is a dimensionless value. The unit of c is (Nms −1 ) 4 The Pearson product-moment correlation coefficient (PPMCC) of the logarithm of the seismic moment and seismic interval for each SDS Fig. 1 Spatial distribution of the 49 seismically dense spheres (SDS). The background color is the surface elevation (km) from the Etopo digital elevation model (Smith and Sandwell 1997). The black and violet lines indicate the iso-depth (km) contours of the upper surface of the subducted Pacific (PAC) and Philippine Sea (PHS) plates (Nakajima and Hasegawa 2007). The overlapping colored circles marked by numbers 1 to 49 are the seismically dense spheres targeted in this study (8.65 + 10 km). The a-value is − 5.63, the b-value is 0.99, and the k-value is 1.52, which indicate that the seismic moment is proportional to the seismic interval to a power of 1.52. The Pearson product-moment correlation coefficient (PPMCC) of the logarithm of the seismic moment and seismic interval for SDS1 is 0.798, which indicates the confidence level. Most of the 49 seismically active zones feature b-values of approximately 0.82-1.08, except SDS3 (0.71, western Tottori), SDS32 (0.78, northern Hiroshima) and SDS36 (0.79, offshore Ibaraki). The b-values of these study regions are approximated to be close to 1.
In addition, nearly all of the centers of these seismically active zones are located in the upper crust (< 20 km) inland (intraplate earthquakes) or at a depth of 30-60 km offshore (interplate earthquakes). Exceptions are that SDS23 (65.47 km, Tokyo) and SDS34 (60.77 km, southern Ibaraki) are inland but influenced by overlapping subduction of the Pacific plate and Philippine Sea plate at a depth of 60-80 km, and SDS31 (37.28 km, eastern Aichi) and SDS44 (34.60 km, western Ehime) are also inland but in the lower crust, which is probably facilitated by upwelling fluids/ melts from the subducted Philippine Sea plate at a depth of approximately 60 km beneath. Figures 2, 3, 4, 5. 6, 7, 8 show that the seismic moment (red curves in the upper part of Fig. a) and seismic intervals (blue curves in the lower part of Fig. a) have an approximately scaling relationship in magnitude. The variation in the magnitude of the seismic moment is usually accompanied by a simultaneous similar variation in the magnitude of the seismic interval in all 49 study regions. This scaling phenomenon is found not only among the intraplate Fig. 2 a The relationship between the logarithm of the seismic moment (the green curve indicates the total catalog, the red curve indicates the DMSI catalog, and the pink curve indicates the MMSI catalog following the left vertical axis, log 10 [Nm]) and the logarithm of the seismic interval (the gray curve indicates the total catalog, the blue curve indicates the DMSI catalog, and the cyan curve indicates the MMSI catalog following the right vertical axis, log 10 [s]), with Mw < 5 during the time interval from October 1, 1997 to February 28, 2015. The pink curve correlates well with the cyan curve. The pink open circles in the upper part indicate contemporaneous peripheral earthquakes (Mw > 5) within a distance of 500 km and with a moment following the left vertical axis, log 10 [Nm]. The circle radii are proportional to the reciprocals of the hypocenter distances from the sphere center. b The relationship between the values of MMSI (red circles), DMSI (orange circles) and the seismic moments. The cyan line indicates the predicted Mo-SI distribution following Eq. 3, while the gray line indicates the predicted Mo-SI distribution following Eq. 4 earthquakes in the upper crust but also among the subduction earthquakes at depth. The b-value does not change much between the two types of earthquakes or at various depths. Since sphere-shaped regions with radii as small as 10 km are prescribed in this study, the results indicate that the scaling law for seismic moment and seismic interval stands for the third dimension of depth, which could be regarded as the same as the horizontal dimensions and distances. That is, the 3-D dimensions are applicable to the G-R law and Mo-SI scaling law. The seismically dense spheres within the same region or at the same depth added together still comply with the Mo-SI scaling law, which has been recognized in many previous studies, although the a-value and b-value vary with changes in area or depth. From this angle, we understand that the 3-D-based G-R law and Mo-SI scaling law are reasonable extensions of the traditional 2-D-based perspective, which reflects the essence of the G-R law and Mo-SI scaling law.

Comparison with Mo-SI relationship
Combining the moment magnitude scale law log 10 (M o ) ≈ 1.5M w + 9.09 (Hanks and Kanamori 1979) with the Gutenberg-Richter law log 10 N = a − bM (Gutenberg and Richter 1956), we have where M 0 is the seismic moment, c is the earthquake moment rate, c = 10 3a 2b +9.09 , T is the average seismic interval (unit: s), a and b are the coefficients in the Gutenberg-Richter law, and k = 3 2b (see Additional file 1: "Explanation for Eq. 1"). Equation 1 interprets the relationship we observed between the logarithm of the seismic moment and MSI in the observed seismically dense spheres (Figs. 2,3,4,5,6,7,8). In addition, assuming that k is 1 and c = M o T ≈ 10 4 −10 6 (Nms −1 ), which is equivalent to fewer than 100 earthquakes with Mw ≥ 1.0 recurring every day, we have (2) log 10 (M o ) = log 10 (T ) + log 10 (c) ≈ log 10 (T ) + 4 ∼ 6, which shows that if a giant earthquake occurs within 500 km, the value of log c may reach 7-8, resulting in a larger gap between the curves of SI and magnitude (e.g., seismically dense sphere 24).
If we consider the definition of the seismic moment M o = µDL 2 , where µ is the fault rigidity or fragility (Pa), D is the average slip amount (m), and L is the rupture length of a fault (m), and assume a seismic load rate v (m/s) during the period of a seismic interval SI (T), then D = vT , and the seismic moment can be described as: If the seismic load rate v and the fault rupture length L are constant, then the seismic moment M 0 will scale to T(SI) by c (c = µL 2 v) . The c-value is proportional to the seismic load rate v, the fault rigidity or fragility µ and the fault rupture length L, which are assumed to be correlated with the fault strength. This correlation can be used to interpret the observed Mo-SI relationship coincidentally from the slope of the seismic moment. Assuming that µ is 30 GPa, L is 30 m for Mw = 1 (i.e., logMo = 10.6, which is the rough average magnitude in seismically dense sphere 49, as shown in Fig. 2), and v is 0.05 m/year, we can infer c = µL 2 v = 4.28 × 10 4 (Nms −1 ), whose logarithm is 4.63, which complies with the empirical estimate in Eq. 2. The other seismically dense spheres also follow this estimate. For example, logMo ≈ 11 in SDS5 (Fig. 3), then L ≈ 41.1 m, c = µL 2 v = 8.03 × 10 4 (Nms −1 ), whose logarithm is 4.90; logMo ≈ 10 in SDS7 (Fig. 3), then L ≈ 19 m, c = µL 2 v = 1.73 × 10 4 (Nms −1 ), whose logarithm is 4.24.

Comparison of the different dimensions and b-values of the fault systems
Different fault systems loaded along various fault planes account for various fault slip behaviors. The direction of the load rate v is dependent more on the strike and dip of the fault system than on the relative motion of the plate direction. More importantly, the slip behavior probably changes its direction along a curved fault plane during the rupture occurrence. For a 3-D curved fault plane, the rupture length and the area of the ruptured zone likely change during fault slip. In this case, L (fault rupture length) is a variable similar to D (slip), and we must consider whether it should be divided by the seismic interval (T). Therefore, the seismic moment rate could also be two-or three-dimensional, involving full ranges of dip, rake and strike. Compared with M 0 = cT k , if we consider the dimension to be k = 3 (corresponding to b = 0.5) and the fault rupture length L to be approximately proportional to the seismic slip amount D, then we would have a scaling law related to two-and three-dimensional slip behaviors: where τ = D/L = E/M 0 is the ratio of the seismic energy and seismic moment, which is proposed to be a constant of 5.0 × 10 −5 (Kanamori 1983), and α is a proportional coefficient between the rupture length and characteristic length (e.g., Leonard 2010) and is approximately a constant for M < 6 earthquakes (e.g., Weng and Yang 2017). αµτ −2 v 3 is the seismic moment rate (c-value), and v is the seismic load rate. Equation 4 has an advantage in that it does not involve the characteristic fault length L, which further suggests that the key factor influencing the relationship between the seismic moment and SI is not simply the moment rate but is also the seismic energymoment ratio over a long period. For the MMSI dataset (red circles), both equations are applicable to the DMSI catalog (orange circles) with a lower resolution. Earthquakes with M 0 > 10 15 Nm (M > 4) have a smaller MMSI than the values predicted because the M > 4 events are partially contaminated by the clustered small events.
Equations 3 and 4 correspond to the empirical maximum and minimum b-values of 1.5 and 0.5 (e.g., Wiemer and Benoit 1996), respectively, which indicates that the dimensions of slip and rupture are complex and that the dimensions probably play an important role in determining the MO-SI coefficient. The recorded events are likely accompanied by two-or three-dimensional slips dominated by Eq. 3 or 4. In addition, M < 4 (small size) earthquakes are likely more two-and three-dimensional than M > 4 events due to their smaller b-values. The range of b-values from 0.5-1.5 for earthquakes has also been estimated in other regions, such as Alaska (Freymueller et al. 2008), South Asia (Kayal 2008) and Yellowstone (Farrell et al. 2009).

Comparison with the seismic sequences
The deviation of the seismic slip rate from the average seismic slip rate may account for the long-term variation in the value of c, as seen from the gap size between the pink and cyan curves of M 0 and MSI, respectively, in Figs. 2,3,4,5,6,7,8. This deviation possibly suggests coseismic slip or a slip deficit along a seismogenic (4) fault. A diverging trend of the two curves appears after the 2011 Tohoku-oki great earthquakes at the seismically dense spheres close to Tohoku, such as seismically dense spheres 20, 24 and 40. This trend might suggest a decreasing seismic slip rate, from the coseismic peak to a postseismic normal rate, possibly representing fault recoupling. This kind of aftershock sequence has also been found to follow the modified Omori law n(t) = K(t + c) −p (Utsu et al. 1995) or epidemic-type aftershock sequence model (Ogata et al. 2003). Most seismic sequences experience an immediate increase in Mo and an immediate decrease in SI and a subsequent decay as Mo and SI approach prequake durations, which show that the degree of seismic variation in Mo and SI is a function of the radius (r) and nucleation zone size (h*) of the velocity-weakening patch (e.g., Chen and Lapusta 2009;Chen et al. 2010;Uchida et al. 2015). If we use T = 1/n(t) to replace T in Eq. 3, then we easily have n(t)M 0 (t) = µL 2 v , and log n(t)+logM 0 (t) = logµL 2 v , which reasonably interprets the simultaneous decay of Mo and SI (T) for most of the seismic sequences. The sparsely distributed stage of seismicity usually accompanies a narrow gap between the pink and cyan curves and is interpreted as a slip deficit (e.g., Fig. 2, before September 2001), indicating a-value of logµL 2 v.

Comparison with the seismic rate and global seismic moment
Based on the calculated background seismicity rate (Mw > 4.5) for 117 subduction zones globally, a proportionality relationship, in which the relative plate motion velocity correlates with the seismicity rate, has been demonstrated (Ide 2013). If the total seismic energy release rate is hypothesized to be less variable, then Eq. 4 is equivalent to where R = 1/T is the seismic rate (s −1 ), V is the plate seismic load rate (m/s), and M ot is the total seismic moment (Nm) during a long period in the subduction zone or the global system. The advantage of Eq. 5 is a potential proportionality of the seismic rate with the relative plate motion velocity.

Calculation of the scaling coefficient in Japan
We further enlarged the modeled region to include all of Japan using the same catalog (Fig. 9). If a radius of 20 km and a threshold of ten thousand events for each sphere are applied, the resulting distribution of the c-values of the earthquake moment rates is shown in Fig. 9. Along offshore Tohoku and the Japan Trench, high c-values are identified despite other parts of inland Japan showing (5) R ∝ µ −1 τ 2 M ot − 1 3 V , a predominantly low c-value. This difference might be caused by the different plate seismic load rates and the different rigidities or fragilities of the continental crust, according to Eq. 3. In contrast, along with the Izu-Bonin island arc chain, the high dominant c-value is possibly caused by the thinner Philippine Sea (PHS) oceanic crust over the subducted Pacific (PAC) plate. Figure 10 depicts the calculated extrapolated distribution minimum c-values for every sphere from 2009 to 2010, implying a larger c-value (shorter seismic interval) for the regions abundant in M5 + earthquakes. This provides a potentially important way to study the details of the featured regional seismicity. Although the seismic load rates and c-values vary greatly in seismogenic zones, the regional seismicity is probably controlled by fault interaction properties such as fault strength and patch size (e.g., Taira et al. 2009;Chen et al. 2010) that are influenced by a tectonic background, which has the potential to provide insight into earthquake prediction provided that a sufficient number of recorded earthquakes enable such a study. Fig. 9 Distribution of the logarithm of the seismic moment rate (log c ) obtained by the MMSI method and compared with the hypocenters of recorded M5 + earthquakes from 1997 to 2015 at depths shallower than 60 km using a threshold of 10,000 events for each seismically dense sphere. The color indicates the seismic moment rate (log c ) value in and around Japan in this study. The squares and triangles indicate that the seismically dense spheres are unable to calculate the MMSI due to insufficient foreshocks before (squares) or after (triangles) December 31, 1998, while the circles indicate the seismically dense spheres with sufficient foreshocks to calculate the MMSI. The number inside the squares, triangles, and circles indicates the magnitude of M5+ earthquakes. The size of the squares and triangles is fixed. The radii of the circles indicate the ratio value of the smallest seismic moment rate of the foreshocks versus the seismic moment rate of the M5+ mainshock. Small radii of circles indicate that the M5+ earthquake (usually landward) occurred immediately after a large seismic interval or the foreshocks are associated with a small seismic interval