Inverse analysis of seismic activity rate changes for severely incomplete sequences: comparison of aftershock activity patterns immediately following the 2023 M6.5 and 2024 M7.6 Noto Peninsula earthquakes

We estimate the true aftershock activity rate from the detected initial aftershocks of the M7.6 Noto Peninsula earth-quake on January 1, 2024, in addition to the preceding M6.5 Noto Peninsula earthquake on May 5, 2023. Estimating the initial aftershock activity rates has the serious problem of missing data, but in itself contains important clues to the physical structure associated with aftershocks and other external factors. We first estimate the detection rate models of the two aftershock sequences assuming the Gutenberg–Richter frequency laws


Graphical Abstract 1 Introduction
In the last century, there has been little seismic activity above M4 in and around the Noto Peninsula, except for the aftershocks of the 1993 M6.6 Noto Peninsula Offshore Earthquake (Ogata 2007(Ogata , 2011)).In 2018, seismic activity began to increase in the northern part of the peninsula, and by the end of 2020, microseisms (M ≥ 1) that had previously been observed in small numbers at shallow depths in the vicinity of Suzu City suddenly shifted to the depths of 14 km or more at the end of November 2020, causing intermittent swarm earthquakes.Then, with a delay of several months, they gradually spread clockwise to the surrounding area (Kumazawa and Ogata 2022a, b;2023a, b).Later, on May 5, 2023, an M6.5 earthquake occurred in the northern region, which cascaded into the M7.6 earthquake.
This paper compares the statistical characteristics of the aftershocks of the M6.5 Noto Peninsula earthquake of May 5, 2023, and the M7.6 Noto Peninsula earthquake of January 1, 2024.To discuss these in detail, we need to avoid the difficulties for analyzing significant missing data immediately after large earthquakes.In other words, applying the ETAS model to aftershock sequences immediately following large earthquakes is inherently a poor fit because of the severe lack of aftershock data, resulting in biased estimates specific to the aftershock sequences.This is because aftershocks even with relatively large magnitudes are not detected (e.g., Enescu et al. 2009;Sawazaki and Enescu 2014).To avoid estimation bias, setting the lower magnitude limit of the aftershock data used at a high value increases the estimation uncertainty because the sample size becomes smaller.To overcome this difficulty, Ogata (1983Ogata ( , 2001) ) has conventionally treated a main shock and its immediately following aftershocks as a history, and excluded them from the target interval for fitting the ETAS model.
In this paper, however, all aftershock data from the main shock are fitted to obtain the characteristics of the aftershocks without any delay.For this purpose, we propose different methods from Omi et al. (2014), Zhuang et al. (2017), Morikawa et al. (2021) and Hainzl et al. (2024) to estimate the detection probability of aftershocks and to correctly compensate the ETAS model accordingly.
Furthermore, to test the significance of time variation by external effects such as significant pore pressures or slow slips in addition to the static triggers, the result of the nonstationary ETAS model is compared against the (stationary) ETAS model.

Magnitude frequency distribution model with missing measurements
In this paper, we assume that earthquake magnitude series always follow the Gutenberg-Richter (GR) frequency law (exponential distribution), even if they have different scales (i.e., b-values).Large earthquakes are almost certainly detected in the earthquake catalog, while small earthquakes are almost always missed, and intermediate magnitudes are detected or missed by a certain percentage.Therefore, the detection rate of earthquakes is expressed as a cumulative function of the normal distribution of magnitudes, erf {・} (Fig. 1).That is to say, the parameter n represents the magnitude at which the detection rate is 50%, and s represents the range of the magnitudes where partial detection occurs.The GR distribution and the detection rate distribution are then multiplied together.Then, the shape of the density function equation, with the y-axis on a logarithmic scale, is very similar to the empirical histograms of earthquake catalogs, as seen for example in Ogata andKatsura (1993, 2006).Specifically, the detection distribution density function for each earthquake of magnitude M t , as a function of elapsed time t, is (see Fig. 1a) From this, the detection rate of earthquakes with M ≥ M c is (see Fig. 1b) The parameters b t , ν t , σ t that depend on the time t are represented by cubic spline functions of equally spaced time intervals.This smoothing is determined by a trade-off between minimizing the smoothness penalty, which is the integral of the weighted sum of squares of each of the first-and second-order derivatives of the spline function, and maximizing the log-likelihood function.The optimal method for determining the weights (hyperparameters) of the penalties is important and is obtained by minimizing the Akaike Bayesian Information Criterion (ABIC; Akaike 1980).
However, in this paper, in order to automatically relax the smoothing penalty in a short period immediately after a strong triggering effect of the ETAS model, a cubic spline function is defined on the equally spaced intervals in [0, N] where N is the total number of event data on [0, T], and the spline function of n = 1,2, …, N is transformed into the function of real time t n over [0, T].This provides a more effective representation than penalties applied in real time when earthquake events are densely clustered in real time, such as immediately after large earthquakes or bursting clusters of earthquake swarms.

ETAS models with compensation
Consider a stationary ETAS model (Ogata 1988), which is applied to the aftershock data of magnitude M c and greater, where H t = {(t i , M i ); t i < t, M i ≥ M c } is the history of detected aftershocks.The missing data are accounted for by the following compensated conditional intensity function: using the detection probability function in (2).Then, the following log-likelihood function: (1) is maximized to obtain the maximum likelihood estimate (MLE) θ = ( μ, K , ĉ, α, p) of the original ETAS model in Eq. ( 3).The integration function in the second term of Eq. ( 5) is approximately calculated as a piecewise linear function.
In the nonstationary ETAS model (Kumazawa and Ogata 2013), Time-dependent functions µ(t) and K (t) are assumed to be smooth piecewise linear functions connecting the respective values µ i = µ(t i ) and K i = K (t i ) correspond- ing to earthquake i.That is to say, the log-likelihood in Eq. ( 5) is adjusted with a smoothing penalty to obtain the log posterior distribution function: Then, the ABIC is minimized by adjusting the smoothing weights w μ , w K , and the ETAS parameters (c, α, p) , to obtain the maximum a posteriori (MAP) estimation θMAP (see Appendix in Kumazawa and Ogata 2013, for   details).Note here that the constraint penalties in Eq. ( 7) are imposed on the order of the sequence of earthquakes to perform a detailed inverse analysis of the changes immediately after the main shock and major aftershocks.Namely, when the occurrence times are highly clustered or when large discontinuities are emphasized, the penalty applied to the sequence time advancing by one between two consecutive events, as in Eq. ( 7), may be better than ( 6) the penalty on the slope of the local linear function applied to the ordinary time (see Kumazawa and Ogata 2013, Appendix B).The superiority of this choice can be confirmed by comparing the ABIC values using a Bayesian method.
In the case of the maximum likelihood method, which maximizes the log-likelihood in (5), the goodness-of-fit comparison can be made by calculating the AIC (Akaike 1974) from the log-likelihood in (5) of the compensated stationary ETAS model in (4).However, when this is compared with a nonstationary ETAS model, the AIC values cannot be directly compared with the ABIC value because their respective common abbreviated terms are incompatible.
Instead, we examine whether the MAP solution of the nonstationary Bayesian model in ( 7) with sufficiently large weights w and w K , which make the functions m(t) and K(t) constant, is equivalent to the stationary ETAS model with the MLE.Then, we compare the ABIC value (= ABIC 0 ) of such a Bayesian model with the ABIC opt values of the optimal weights to test the significance of the fit.That is to say, their difference, ABIC = ABIC opt − ABIC 0 , the confidence probability of the nonstationary model is exp(−ABIC/2) times large (Akaike 1985; Ogata 2017, see Appendix A; Ogata et al. 2013).
When the superiority of the compensated nonstationary ETAS model is confirmed, we then obtain the inversion solution, the maximum a posteriori (MAP) estimate.In this way, the compensated nonstationary ETAS model

Aftershock detection rates immediately
following the M7.6 and M6.5 earthquakes We use all detected aftershocks in the catalog detected and compiled by the Japan Meteorological Agency (JMA).Specifically, we examined all detected aftershocks up to January 12, 2024, including the M7.6 earthquake on January 1 (Fig. 2b), and also detected aftershocks from May 1, 2023 to 12 days after the M6.5 earthquake (Fig. 2a).
The inversion results for the magnitude-frequency model in (1) are shown in Fig. 3. First, from Fig. 3a, the b-value of the aftershock activity for the M7.6 earthquake is constant.On the other hand, the change in σ values is significant, which may reflect the size of the aftershock area relative to the seismic station layout in the northern Noto Peninsula.Namely, there is a limited number of nearby stations in the interior, while the M7.6 aftershock zone extends to the distant coasts (Fig. 2b), where smaller events are better detected for the earthquakes occurring closer to the seismic station than for the more distant earthquakes.This is unlikely to be the case for the small aftershock area of M6.5 (Fig. 2a), so the σ value is constant throughout.
The difference in detection rate level between day and night is also adapted in the estimation of the 50% detection rate value (red curve) and also in the 97.5% detection rate (ν + 2σ value), represented by the green curve in Fig. 3b, d.
On the other hand, Fig. 3c shows that the σ value of the M6.5 aftershock activity is constant at 0.274 throughout the entire period, while the b-value varies from 0.9 to 1.3.It decreases immediately after the M6.5 main shock, increases rapidly after the largest aftershock, and then converges to an almost constant b-value around 0.9.
The blue curve in Fig. 3b, d shows the detection probabilities of M ≥ 2.0 for the M7.6 aftershock and M ≥ 1.0 for the M6.5 aftershock, respectively.These sensitively show decrease in the detection rates of secondary aftershocks immediately after a moderate or large aftershock in addition to some quasi-periodic fluctuations.

Analyzing aftershocks using compensated ETAS models
We calculate the earthquake occurrence rates with the compensated ETAS model in ( 4), which takes into account the effects of the detection probabilities shown in Fig. 3. First, the aftershock activity immediately after the main shock of the M7.6 Noto Peninsula earthquake in January 2024 is analyzed.Equation ( 4), which multiplies the probability of detecting of an aftershock of M ≥ 2 shown in the green curve in Fig. 4a (the same as the blue curve in Fig. 3b) by the intensity function of the original stationary ETAS model in (3), defines the intensity function of the compensated ETAS model in (4) for the detected aftershocks, and its log-likelihood in ( 5) is maximized with respect to the original ETAS model parameters θ to obtain a maximum likelihood estimate (MLE) (Ogata 2007).The ETAS model for this MLE parameter value gives the true intensity function of the original ETAS model in (3), as indicated by the red curve of the intensity function in Fig. 4a.The occurrence of M ≥ 2 aftershocks by this ETAS model largely compensates for the missing data in the first few days and continues to compensate in detail for subsequent days.The μ value for this interval is constant (94.12 events per day).
We then applied the nonstationary ETAS model in ( 6) for the same aftershock sequence, with optimal smoothing by the ABIC method by adjusting penalty weights to the roughness of the time variation of μ and K values.Referring to Table 1a for each weight (hyper-parameter), the larger the weights of μ and K values, the smaller the ABIC values.As the weights of the smoothing constraints increase, the maximum a posteriori (MAP) solution becomes eventually equivalent to the stationary ETAS model since both the μ and K curves become almost constant, and its MAP solution is the same as the MLE solution in Eq. ( 3).
Similarly, for the aftershocks of the M 6.5 earthquake, the intensity function (4) of the detected aftershocks is defined by multiplying the intensity function ( 6 based on the log-likelihood function ( 7) with smoothing penalty is given in Table 1b.Based on this, the difference in ABIC between the model with the best weights ( ŵµ , ŵK ) = (10 4 , 10 5 ) and the model with the heaviest weights ( ŵµ , ŵK ) = (10 6 , 10 13 ) is −12.3; this difference is significant.In other words, the former model was exp{−ΔABIC/2} ≈ 468 times more probable than the latter, indicating that the compensated nonstationary model is far superior to the compensated stationary ETAS model.On the other hand, when the constraint of m value is further weakened (i.e., ŵµ ≤ 10 3 ) in Table 1, the performance becomes extremely poor.
In Fig. 4b, the optimal MAP solution of the original nonstationary ETAS is shown as a red curve, which appears to compensate for the apparent occurrence rate since the missing measurements of the first few days are substantially filled in and then the detection probability is close to 1.0 thereafter.The aftershock induction coefficient K (t) shows little variation.In particular, the back- ground intensity rate μ(t) increases sharply by a factor of 10 in the first few days after the main shock and then decays.Hardebeck et al. (2024) review physically hypothetical mechanisms in detailed aftershock generation scheme, including research problems of dynamic aftershock triggering related to pore fluid pressure changes, to improve physics-based aftershock forecasting.

Discussion
Regarding the swarm activity in Noto Peninsula, there are a number of papers (Nakajima 2022;Kato 2024;Yoshida et al. 2023a, b;Nishimura et al. 2023) suggesting that fluid was supplied from the deep source, and this fluid slowly slipped through the seismic-free zone and penetrated the surrounding vulnerable areas.
In addition, Kato (2024) indicates that the early aftershock front of the M6.5 mainshock initially migrated at a speed of ∼20 km/hr and suggests that this rapid migration of the immediate aftershocks could be driven by upwelling of crustal fluids along the intensely fractured and permeable fault zone by the dynamic rupture of the mainshock.
Figure 4b shows the time evolution of the background intensity, which might be presumably due to strong fluctuations in fluid pore pressure caused by the M6.5 mainshock throughout the aftershock region (Kumazawa and Ogata 2023b).Similar migrations were observed for the space-time aftershock pattern of the M6 class earthquakes referred by Ogata (2007Ogata ( , 2011) ) and Kumazawa and Ogata (2024).
On the other hand, the spatiotemporal aftershock distribution of M7.6 is characterized by the time-invariant, constant background intensity (Fig. 4a; Kumazawa and Ogata 2024).This indicates that there is no extension of the aftershock front as were seen in the first half examples of Ogata (2010) or Fig. 12 in Kumazawa and Ogata Table 1 Comparison of the ABIC values for the nonstationary ETAS model fit ΔABIC value for weights (w µ , w K ) = (10 x , 10 y ) ; for (a) the aftershocks of the 2024 M7.6 Noto Peninsula earthquake, and (b) the aftershocks of the 2023 M6.5 Noto Peninsula earthquake, respectively.2024).Such a pattern is typical of a sequence in which static stress changes due to the fault motion of the main shock dominate other factors, and is characterized by a very small constant background intensity, as seen in Fig. 4a, From Figs. 3 and 4, we can see that the b-value and μ value are constants for the M7.6 aftershocks.In addition, the change in b-value and μ value for the M6.5 aftershocks appears similar in shape.However, we cannot definitively explain the reason why b-value evolves significantly with time for M6.5 aftershocks while b remains constant for M7.6 aftershocks, but this may indicate the ratio of the pore pressure increase relative to the shear stress increase caused by the rupture.For the consistent examples, we can only refer to the two classic papers: first, Utsu (1962) showed that the magnitude distribution of aftershocks (b-values) did not change with time for three aftershock series, including the 1957 Mw9.1 Aleutian earthquake; and second, in Healy et al. (1968), the Denver earthquake (1962-68), a well-known waterinjection-induced earthquake, showed that well pressure (which dominates the earthquake rate) and b-value (0.6-0.9) are positively correlated.

Conclusions
The temporal variation of the b-values considering the missing aftershocks of the M7.6 Noto Peninsula earthquake is constant, but the aftershocks of the M6.5 earthquake show a large variation.Furthermore, the compensated ETAS models of the respective aftershocks lead to different feature inversion results between the M7.6 and M6.5 earthquakes.For the former sequence, the compensated stationary ETAS model fits better, while for the latter, the compensated nonstationary ETAS model fits the latter well.The small constant in the μ value for the background intensity rate after the M7.6 event may indicate that the effect of pore pressure was significantly smaller than the shear stress increments caused by the main rupture.In contrast, the increase in the background intensity rate after the M6.5 event reflects the possibility that the effect of pore pressure became comparable to the shear stress increments.
Since early aftershock, forecasting plays an important role in emergency risk mitigation and immediate response efforts, stable procedures for real-time forecasting immediately after a mainshock need to be implemented.

Fig. 1
Fig. 1 Schematic diagram of the detected magnitude distribution model.a Curved line with knots represents magnitude frequency model in Eq. (1), obtained by multiplying G-R law and cumulative normal distribution (error-function), where ν represents the magnitude of earthquakes with 50% detection, and σ represents a range of partially detected magnitudes.b Detected probability of earthquakes in Eq. (2), represent the ratio of theoretical cumulative frequency of magnitude M c (blue real segment with arrows) to cumulative G-R law (red dotted segment with arrows)

Fig. 3
Fig. 3 Temporal evolution of aftershock detection rate parameters.Panels (a, b) and (c, d) are b-values, σ-values, and detection rate changes of the aftershocks of the 2024 M7.6 and 2023 M6.5 Noto Peninsula earthquakes, respectively.The time unit is day.The gray dots in all panels represent the magnitudes of the detected earthquakes, as a function of time from the corresponding mainshocks.Panel (a, c): red curves are b-value changes (with onefold standard error) and blue curves are σ-values (with onefold standard error).Panels (b, d): red curves are the ν values for the 50% detection rate, green curves are the ν + 2σ values.Blue curves are the detection rates in (2) for the b M ≥ 2.0 and d M ≥ 1.0 aftershocks, respectively Fig. 4 The optimal ABIC MAP solution for the ETAS models.a and b aftershocks of the M7.6 and M6.5 events, respectively, where the time unit is day.Red curves are true activity rates for a M ≥ 2 and b M ≥ 1, respectively; black curves are observed (apparent) a M ≥ 2 and b M ≥ 1 aftershock rates; blue horizontal line and curve are true background activity rate μ(t) values (with onefold s.e.); purple curve in b is K(t) values (with onefold s.e.); green curves are detection rate changes in (2) for a M ≥ 2 and b M ≥ 1 aftershocks, respectively.Gray disks show magnitude versus time for all detected earthquakes The smallest ABIC values are 0.0