A positive answer on the existence of correlations between positive earthquake magnitude differences
Eugenio Lippiello Lucilla de Arcangelis Cataldo Godano



    Department of Mathematics and Physics, University of Campania “Luigi Vanvitelli”, Caserta, ItalyDepartment of Mathematics and Physics, University of Campania “Luigi Vanvitelli”, Caserta, ItalyDepartment of Mathematics and Physics, University of Campania “Luigi Vanvitelli”, Caserta, Italy
          The identification of patterns in space, time, and magnitude, which could potentially encode the subsequent earthquake magnitude, represents a significant challenge in earthquake forecasting. A pivotal aspect of this endeavor involves the search for correlations between earthquake magnitudes, a task greatly hindered by the incompleteness of instrumental catalogs.
        A novel strategy to address this challenge is provided by the groundbreaking observation by Van der Elst (Journal of Geophysical Research:
        Solid Earth, 2021), that positive magnitude differences, under certain conditions, remain unaffected by incompleteness. In this letter we adopt this strategy which provides a clear and unambiguous proof regarding the existence of correlations between subsequent positive magnitude differences. Our results are consistent with a time-dependent b-value in the Gutenberg-Richter law, significantly enhancing existing models for seismic forecasting.
          
        
        A positive answer on the existence of correlations between positive earthquake magnitude differences
            Cataldo Godano
            April 29, 2024
        ====================================================================================================
        
        
        
        Seismic catalogs suffer from significant incompleteness. Small earthquakes often go unreported due to their occurrence at distances beyond the detection range of monitoring stations <cit.>, or they may be overshadowed by the coda-wave of preceding larger earthquakes <cit.>. This latter mechanism, in particular, renders the catalog inherently incomplete <cit.> and has left fundamental questions about seismic activity unanswered.
        One such crucial inquiry involves seismic forecasting and the potential dependence of earthquake magnitudes on the organization of seismic events in time, space, and magnitude preceding its occurrence. Studies have convincingly demonstrated correlations between earthquake magnitudes, with these correlations depending on both temporal and spatial distances <cit.>. However, the challenge lies in discerning whether these observed correlations are genuine features or spurious artifacts of catalog incompleteness <cit.>. Consequently, the most commonly adopted statistical forecasting models, such as the ETAS model <cit.>, assume earthquake magnitudes to be independent and identically distributed (i.i.d.) variables. Within this framework, forecasting is severely limited, relying solely on the information that large earthquakes have a higher probability of occurring when and where the seismic rate is higher. Modifications of the ETAS model, incorporating magnitude correlations have been proposed in the literature <cit.> 
        
        Recently, Van der Elst <cit.> introduced a novel and remarkable tool designed to address issues associated with catalog incompleteness in a simple yet highly effective manner. This tool is grounded in the observation that, under specific conditions, positive magnitude differences remain largely unaffected by incompleteness <cit.>. Leveraging this approach, researchers have been able to derive the true magnitude distribution, the true occurrence rate <cit.>, and the correct parameters of the ETAS model <cit.>, even when working with incomplete catalogs. In this study, we employ this tool to unequivocally demonstrate the presence of genuine correlations between earthquake magnitudes.
        
        A seismic catalog is a multi-variate dataset C consisting of N elements, where the i-th element is the vector m_i, t_i, r⃗_i representing the earthquake magnitude, occurrence time, and epicentral coordinates, respectively, of the i-th earthquake. Specifically, we examine the relocated catalog for Southern California (SC) <cit.>, spanning from January 1981 to March 2022, encompassing N=800499 earthquakes with m ≥ -1. Additionally, we consider the relocated catalog for Northern California (NC) <cit.>, covering the period from January 1984 to December 2021, including N=879547 earthquakes with m ≥ -0.8.
        
        From the original catalog, we also construct a catalog of magnitude differences, denoted as D, where each element i is represented by the vector δ m_i, t_i, r⃗_i, with δ m_i = m_i+1 - m_i representing the magnitude difference between two successive earthquakes in the original catalog.
        
        
        
        To investigate correlations between earthquake magnitudes, we extend the method developed in <cit.>. We define the conditional probability as:
        
        
        P(δ m_i > M |δ y_i < Y) ≡N(M,Y)/N(Y).
        
        
        Here, δ y_i denotes either the epicentral distance δ r_i = |r⃗_i+1 - r⃗_i | or the temporal distance δ t_i = t_i+1 - t_i, depending on the context. In Eq.(<ref>), N(M,Y) represents the number of pairs of subsequent events where both δ m_i > M and δ y_i < Y, while N(Y) is the total number of pairs with δ y_i < Y.
        Starting from the initial catalog 𝒞, we generate a catalog with reshuffled magnitude differences 𝒟_ran, where the i-th element is the vector δ m^*_i, t_i, r⃗_i, with δ m^*_i = m_k - m_i, and k represents the index of an earthquake randomly chosen within the catalog. Specifically, the catalog 𝒟_ran retains the same occurrence time and epicentral coordinates as catalog 𝒟, but the magnitude differences are randomized. From a single instrumental catalog 𝒞, we can derive multiple reshuffled catalogs 𝒟_ran. In each reshuffled catalog, we evaluate the quantity P(δ m^*_i >  M |Δ y_i < Y), which exhibits a different value in each reshuffled catalog 𝒟_ran. This quantity follows a Gaussian distribution with a mean value P(M,Y) and a standard deviation Σ(M,Y).
        
        
        
        The key quantity of the study is defined as:
        Δ P(M,Y) = P(δ m_i> M |δ y_i<Y)-P(δ m^*_i> M |Δ y_i<Y).
        
        We then calculate the mean value E[Δ P(M,Y)], obtained by averaging over 1000 realizations of the reshuffled catalog 𝒟_ran. E[Δ P(M,Y)] for SC is plotted in Fig.<ref>a for the case where δ y_i=δ r_i and Y=10 km. The error bars represent twice the standard deviation Σ(δ M,Y). Values such as | E[ Δ P(M,Y)]|>2Σ(M,Y) indicate that it is highly improbable for a catalog with uncorrelated magnitudes, such as 𝒟_ran, to produce the same probability of observing a given magnitude difference P(δ m_i> M |Δ y_i<Y) of the real catalog. 
        Therefore, Fig.<ref>a demonstrates that magnitudes in the instrumental catalog 𝒞 are correlated.
        In Fig.<ref>a, we also examine the influence of a magnitude threshold M_th on E[Δ P(M,Y)]. By considering only earthquakes with magnitudes larger than M_th and increasing M_th, we observe a decrease in E[Δ P(M,Y)]. This decrease can be interpreted as a signature that magnitude correlations are primarily a spurious effect of incompleteness, as increasing M_th focuses on a catalog less affected by incompleteness. Nevertheless, the observation that even for M_th = 3, which exceeds the completeness magnitude M_c estimated from the magnitude distribution (Fig.<ref>b), | E[Δ P(M,Y)]|>2Σ(M,Y) suggests the persistence of genuine magnitude correlations beyond incompleteness. This scenario finds further support in the resilience of E[Δ P(M,Y)] to variations in the quality of the seismic network <cit.> and is reinforced by laboratory experiments on rock fracture <cit.>. These experimental findings demonstrate the existence of magnitude clustering irrespective of loading protocols, rock types, and observable magnitude ranges.
        
        Next, we will review these results in light of recent findings concerning positive magnitude difference statistics.
        We assume as null hypothesis that magnitudes are i.i.d variables which obey the Gutenberg-Richter (GR) law, p(m) ∝exp(-β (m-M_L))  for magnitudes larger than a lower limit M_L. Although a precise estimate of M_L is unavailable, it is reasonable to assume that M_L takes on a very negative value <cit.>. Additionally, we  assume that not all magnitudes are reported in the catalog due to detection issues. Therefore, the observed magnitude distribution is given by <cit.>
         p (m_i) =β e^-β (m_i-M_L)Φ(m_i-M_c),
         
         where M_c>M_L is the completeness magnitude and Φ(x) is a detection function which is a monotonic increasing function of x ranging from Φ(x)=0 for
         x ≲ -2σ to Φ(x)=1 for x ≳ 2σ. Essentially, earthquakes with magnitudes m_i<M_c-2σ are undetected, while those with m_i>M_c+2σ are always detected. A realistic functional form for Φ(x) is Φ(x)=(1+Erf(x/σ))/2, with typical estimates placing σ in the range σ∈ (0.2,0.3)<cit.>.
        Eq.(<ref>) is consistent with the magnitude distribution of instrumental catalogs (Fig.<ref>b), which conforms to the GR law only for magnitudes m_i>M_c, determined using the maximum curvature method  <cit.>, with M_c=2.6 for SC and M_c=2.7 for NC (Fig.<ref>b).
        
        
        Using Eq.(<ref>) in Eq.(<ref>), we obtain
            P(δ m_i>M |Δ y_i<Y)= β^2 ∫_m_L^∞dm_i ∫_m_i+M ^∞dm_i+1
            e^-β (m_i+1+m_i-2 m_L)Φ(m_i-M_c) Φ(m_i+1-M_c| m_i). 
        
        In the above equation, we explicitly use the notation Φ( m_i+1-M_c| m_i) to specify that the detection function must be evaluated under conditions where the previous earthquake m_i has been identified and reported in the catalog.
        It is worth noticing that it is exactly this conditioned detection function that introduces correlations between the magnitudes m_i and m_i+1. Indeed, the same expression Eq.(<ref>) holds for the probability P(δ m^*_i> M |Δ y_i<Y), with m_k instead of m_i+1.
        This makes a fundamental difference since, by construction, k is a random index, and therefore Φ(m_k-M_c| m_i)=Φ(m_k-M_c), independently of m_i.
        
        
        Accordingly, Δ P(M,Y) can be rewritten as
           Δ P(M | Y)= β^2 ∫_m_L^∞dm_i ∫_m_i+M ^∞dm_i+1 
         e^-β (m_i+1+m_i-2 m_L)
           Φ(m_i-M_c)
         [Φ(m_j-M_c| m_i) -Φ(m_j-M_c)].
        
        The origin of magnitude correlations in the null hypothesis scenario originates from the difference of the two detection functions in the term between square brackets in Eq.(<ref>). It is evident that, if M >0, m_j>m_i and Φ(m_j-M_c| m_i) is on average larger than Φ(m_j-M_c). Indeed, if a previous earthquake has been detected, the probability to detect a subsequent larger one is larger compared to the case where no information about previously detected earthquakes is available.
        At the same time, imposing m_j ≥ M_th and by considering increasing values of M_th, both detection functions Φ(m_j-M_c| m_i) and Φ(m_j-M_c) approach 1, and Δ P(M,Y) approaches zero. Accordingly, the behavior of E[Δ P(M,Y)] in Fig.<ref>a cannot exclude the null hypothesis scenario where magnitudes are i.i.d variables in the presence of detection problems.
        
        
        We now explore the prediction of the null hypothesis on the magnitude difference q_i = m_i+1 - m_i, restricting to  m_i+1 > m_i and  δ y_i=δ r_i <Y=d_0.
        The parameter d_0 is chosen sufficiently small to ensure that earthquake pairs are close enough in space, guaranteeing similarity in the distance to the stations necessary for their identification. Under this condition, since the earthquake m_i has been detected, it implies that m_i > M_c - 2σ, and therefore Φ(m_j + M - M_c | m_i) ≃ 1, for M > 2σ.
        Accordingly, from Eq.(<ref>), for M >2σ, P(q_i> M|Δ y_i<Y)=e^-β M K, where K is a constant given by K=∫_m_L^∞dm_i e^-2β m_iΦ(m_i-M_c). After derivation, we therefore obtain, for M>2σ
        P(q_i=M |Δ y_i<Y)=β K e^-β M,  
        
        showing that, even in the presence of detection problems, the distribution of magnitude differences P(q_i) follows an exponential GR decay for q_i>2 σ, independently of the previous seismic history.
        This behavior is observed in instrumental catalogs (Fig.<ref>b), where we find that P(q_i) follows a pure exponential decay when q_i is larger than a completeness threshold δ M_c=0.4, for both SC and NC, consistent with Eq.(<ref>), with the same value of the exponent (β=2.23±0.05 in SC, β=2.20 ±0.05 in NC)  extrapolated from the magnitude distribution P(m_i) (Eq.(<ref>)) when m_i>M_c (Fig.<ref>b).
        
        
        
        
        
        It is crucial to note that while detection issues introduce spurious correlations between m_i and m_i+1, positive magnitude differences remain unaffected by such detection problems. Consequently, when q_i > 2σ, the null hypothesis predicts no correlations between q_i and q_i+1. To delve deeper into this aspect, we select the subset of the catalog 𝒟, denoted as 𝒟', which contains only those earthquakes satisfying both constraints δ r_i < d_0 and δ m_i>0, and therefore the i-th element of the catalog 𝒟' will be the vector q_i,t_i,x_i,y_i.
        From the catalog 𝒟', we next obtain the catalog 𝒟𝒟, whose elements are the vectors δ q_i,t_i,r⃗_i with δ q_i=q_i+1-q_i.
        It is important to note that, in general, δ q_i ≠ m_i+2-m_i, and this equality only holds when m_i+2>m_i+1>m_i and dr_i and dr_i+1 are both smaller than d_0. We also generate several reshuffled catalogs 𝒟𝒟_ran, where the i-th element is the vector δ q^*_i,t_i,r⃗_i with δ q^*_i=q_k-q_i, and k is the index of a random earthquake in the catalog  𝒟'. We can define, P(δ q_i>Q |δ y_i<Y), P(Q,Y), and Σ (Q,Y) by simply replacing m_i with q_i and M with Q in all the definitions from Eq.(<ref>) to Eq.(<ref>).
        Accordingly, from Eq.(<ref>), under the null hypothesis Eq.(<ref>), we obtain
           Δ P_+(Q,Y)= P(δ q_i>Q |Δ y_i<Y)-P(δ q^*_i> Q |Δ y_i<Y) 
            =∫_0^∞dq_i ∫_q_i+Q^∞dq_j P(q_i)P(q_j) Φ(q_i)  [Φ_q(q_j| q_i) -Φ_q(q_j)],
        
        where Φ_q(q_j| q_i) and Φ_q(q_j) represent the conditioned and unconditioned probabilities, respectively, of detecting a magnitude difference q_j. Correlations between q_i and q_i+1 are induced by the term Φ_q(q_j| q_i). Nevertheless, in contrast to Φ(m_j-M_c| m_i) where the information that event m_i has been detected strongly constrains the detectability of m_i+1>m_i, in this case, the information that the previous magnitude difference q_i has been detected only weakly affects the probability Φ_q(q_j| q_i) to detect the next magnitude difference q_i+1. In fact, we expect that Φ_q(q_j| q_i) ≃Φ_q(q_j),
        and Δ P_+(Q,Y) is expected to be very small under the null-hypothesis scenario. In particular, by considering only magnitude differences q_i ≥δ M_th, as δ M_th increases, the detection probability Φ_q also increases and approaches 1 for δ M_th≥ 2σ. Consequently, we expect that Δ P_+(Q,Y) presents small values when δ M_th=0, and decreases as δ M_th increases, vanishing when magnitude differences are detected with probability 1 for δ M_th≥ 2σ.
        
        In the following we investigate the behavior of the mean value E[Δ P_+(Q,Y)] for the SC catalog. The main result is that E[Δ P_+(Q,Y)] exhibits a markedly different pattern (Fig.<ref>) compared to the one predicted by the null hypothesis (Eq.(<ref>)). While data for δ M_th=0 could be consitent with the null hypothesis Eq.(<ref>), as evidenced by E[Δ P_+(Q,Y)] being only slightly larger than 2Σ(Q,Y), the remarkable finding is that for larger values of δ M_th, E[Δ P_+(Q,Y)] not only fails to decrease but rather increases. Particularly notable is that for δ M_th≥ 0.5, significantly surpassing δ M_c ≃ 0.4, when the variables q_i remain unaffected by detection issues (Fig.<ref>b), the correlation between q_i and q_i+1 becomes even more pronounced. The significant growth of E[Δ P_+(Q,Y)] as a function of δ M_th  is notable until δ M_th=0.75, after which no substantial increase is observed up to δ M_th=1.
        Results plotted in Fig.<ref> unequivocally demonstrate the presence of non-trivial correlations among magnitude differences. Consequently, the result derived from the null hypothesis of i.i.d. magnitudes, as expressed in Eq.(<ref>), can be confidently dismissed.
        
        
        We note that the results presented in Fig.<ref> were derived by considering only earthquakes with epicentral distances smaller than d_0=20 km as elements of the catalog D'. Similar results were obtained for other choices of d_0 ≤ 50 km. Additionally, comparable outcomes were observed when analyzing the NC catalog (Fig.<ref>a). This suggests that correlations between magnitude differences appear to be a consistent characteristic of seismic catalogs, with patterns showing weak dependency on the geographic region.
        
        In order to gain further insights into the mechanisms responsible for correlations between q_i+1 and q_i, we examine the influence of Y on E[Δ P(Q,Y)] in Fig.<ref>. Specifically, in Fig.<ref>d, we consider δ y_i=δ t_i and investigate decreasing values of Y, ranging from Y=∞ to Y=0.1 hours. The results show that as Y decreases, imposing that earthquake pairs (q_i,q_i+1) are closer in time, their correlation strengthens. 
        In Fig.<ref>c, we conduct a similar analysis but with δ y_i=δ r_i, with Y decreasing from Y=∞ to 1 km. Here, we observe that magnitude correlations become more pronounced as we transition from Y=∞ to Y=10 km, but then stabilize, remaining roughly constant regardless of the specific value of Y. Consequently, magnitude correlations are stronger for earthquakes occurring closer in time, whereas they are only weakly affected by the spatial distance between earthquakes.
        
        
        Further insights are provided by examining the derivative p'(δ Q,Y)=-.∂ E[Δ P(Q,Y)]/∂ Q| _δ Q, which measures the difference in probability between the instrumental and reshuffled catalogs to observe q_i+1=q_i+δ Q. The results (Fig.<ref>b) reveal that it is more probable for a magnitude difference to be followed by a larger one. Given that the inverse of the average value of q_i coincides with β<cit.>, the origin of magnitude correlations could be linked to a decreasing β value over time during aftershock sequences. This observation aligns with the pattern identified by Gulia & Wiemer <cit.> from their analysis of 58 main-aftershock sequences.
        On the contrary, the scarcity of earthquakes preceding significant shocks poses a challenge in deriving conclusive insights from our analysis regarding the behavior of β for foreshocks, and in comparing them with the predictive framework proposed in references <cit.>.
        
        
        
        
        
        In summary, our findings unequivocally demonstrate the presence of non-trivial correlations between subsequent magnitude differences, suggesting potential modifications to be integrated into the ETAS model. One possible approach, as proposed by <cit.>, is to consider two distinct β values for aftershocks, β±δβ, depending on the magnitude of the triggering mainshock. Alternatively, a time-varying β value during the aftershock sequence could also be implemented <cit.>. This novel research direction holds promise for extracting valuable insights from the intricate patterns of seismicity, ultimately enhancing our ability to forecast seismic events and advancing our understanding of the mechanisms governing earthquake triggering.
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        apsrev4-2
         
          E.L. and L.de A. acknowledge support from the MIUR PRIN 2022 PNRR  P202247YKL,
          C.G. acknowledges support from the MIUR PRIN 2022 PNRR P20222B5P9
