Earthquake forecasting is one of the defining problems of seismology. To provide useful solutions, forecasting models use a wide range of approaches: Coulomb rate-and-state (CRS) models (Cocco et al., 2010;Parsons et al., 2012;Mancini et al., 2019) calculate Coulomb stress changes and couple them with a lab-based constitutive friction law (Dieterich, 1994). On the other end of the spectrum are statistical models, with the epidemic-type aftershock sequence (ETAS) model being the best performing current statistical approach (Cattania et al., 2018;Taroni et al., 2018). First introduced by Ogata (1988), it models seismicity rate as the sum of background and aftershock events, where aftershocks are triggered according to regional empirical laws. In between the purely physics-based and purely statisticsbased approaches are models such as the short-term earthquake probability (STEP) model (Gerstenberger et al., 2005), the Inlabru model (Bayliss et al., 2020), and hybrid Coulomb and statistical models (Steacy et al., 2014). The STEP model combines clustering principles with fault information in a statistical model to produce time-dependent forecasts. The Inlabru model more generally allows the inclusion of diverse data sets as covariates to issue time-independent seismicity forecasts. A hybrid Coulomb/statistical model redistributes seismicity forecasted by STEP according to Coulomb stress changes.
While physics-based models aim to describe the processes and mechanisms underlying seismogenesis, statistical models are generally more empirical and data driven. Ultimately, "all models are wrong, but some are useful," to cite the famous statistician George Box (1979). Usefulness can be viewed from different perspectives. Different forecasting models can be useful for gaining new scientific insight, for producing the most accurate forecasts, or for producing forecasts that are most suited for operational earthquake forecasting (OEF), given the trade-off between accuracy and computational cost. Cattania et al. (2018) found in a pseudoprospective forecasting experiment for the 2010-2012 Canterbury, New Zealand, earthquake sequence that hybrid Coulomb/statistical models have a similar forecasting skill as CRS models, at a lower computational effort. Mancini et al. (2019Mancini et al. ( , 2020) ) conducted pseudoprospective experiments for the 2016 central Italy and the 2019 Ridgecrest, California, sequences, comparing CRS models of different complexity with ETAS forecasts. In both studies, the forecasting skill of CRS models increases with their complexity, with the most complex CRS model performing similarly to ETAS. Hardebeck (2021) investigated possible reasons for the general underperformance of the physics-based models relative to statistical models and suggested that understanding and incorporating heterogeneities in background conditions into physical forecasting models may be key in improving their skill.
Having been tested thoroughly and systematically (Woessner et al., 2011;Ogata et al., 2013;Strader et al., 2017;Taroni et al., 2018;Nandan et al., 2019b;Savran et al., 2020), ETAS models meanwhile remain the state-of-the art of earthquake forecasting and are being used or considered for OEF at various locations (Marzocchi et al., 2014;Rhoades et al., 2016;Field et al., 2017;Kamer et al., 2021;Nandan, Kamer, et al., 2021;van der Elst et al., 2022). Besides using the most basic formulation of ETAS, modelers also commonly refine the model. For instance, Bach and Hainzl (2012) enhanced ETAS with fault information, ShakeMaps, ground-motion models, or Coulomb stress changes. Seif et al. (2017) assessed the biasing effects of data incompleteness and model assumptions on the estimated ETAS parameters. Several techniques have been proposed to address the effects of short-term aftershock incompleteness (Mizrahi et al., 2021b;Grimm et al., 2022;Hainzl, 2022) or the assumption of isotropic aftershock triggering (Grimm et al., 2022;Page and van der Elst, 2022). Other studies focus on deriving spatial variations of ETAS parameters or background seismicity (Enescu et al., 2009;Nandan et al., 2017;Nandan, Ram, et al., 2021), also relating parameter variations with physical quantities such as heat flow. Others have refined the standard ETAS model with a relationship between magnitudes of triggered and triggering earthquakes and a magnitudedependent Omori kernel and found the resulting models to possess improved forecasting performance (Nandan et al., 2019;Nandan, Kamer, et al., 2021). A recent framework for modeling seismicity with an invariant Galton-Watson stochastic branching process provides a generalization of ETAS that is invariant with respect to various common deficiencies of earthquake catalogs (Kovchegov et al., 2022). However, this framework has not yet been used for forecasting seismicity.
A related forecasting topic, which has recently received attention, is ensemble modeling (Rhoades and Gerstenberger, 2009;Marzocchi et al., 2012;Taroni et al., 2014;Bird et al., 2015;Akinci et al., 2018;Llenos and Michael, 2019;Bayona et al., 2021). The idea, widely used for decades in the meteorological and climate forecasting community (Tracton and Kalnay, 1993;Leutbecher and Palmer, 2008;Eyring et al., 2016), is to combine different models in an overarching ensemble model to obtain more robust forecasts. Commonly, an ensemble is a linear or multiplicative combination of ingredient models (e.g., Bird et al., 2015), and the challenge is to optimize the weights given to each model. In a recent study, Bayona et al. (2021) found that the time-independent ensemble models WHEEL and GREAR1 (Bird et al., 2015) outperform the ingredient models of which they consist. Akinci et al. (2018) found that their time-independent ensemble model outperforms its ingredients and performs similarly to the best-performing time-independent model tested in the 2009 CSEP experiment (Schorlemmer, Zechar, et al., 2010;Zechar et al., 2010) for Italy. In the context of time-dependent models, Taroni et al. (2014) and Gerstenberger et al. (2014) used ensemble approaches, and Llenos and Michael (2019) found that ensembles of ETAS models perform best for the 2015 San Ramon, California, Swarm. Shebalin et al. (2014) proposed an iterative method to combine forecasting models and found the resulting models to have advantageous properties compared to the ingredient models or traditional linear combinations thereof. The emerging consensus across the mentioned studies is that ensemble modeling is a promising path to use for earthquake forecasting; this is also demonstrated by the fact that they are currently implemented in Italy's OEF system (Marzocchi et al., 2014). Yet, a breakthrough of ensemble models as established in the meteorological forecasting community is still pending.
For practical operational forecasting, especially in regions that are less studied due to a lack of data or resources, a balance must be achieved between model accuracy and simplicity. With this in mind, we relax some of the assumptions behind ETAS. We allow aftershock productivity and background seismicity to be described nonparametrically, providing event-specific productivity and background rates. This aims to better capture the real behavior of seismicity without making any choices on resolution, parametric form, and so on. Using pseudoprospective forecasting experiments in southern California and Italy, we evaluate whether these flexible ETAS (flETAS) models provide superior forecasts.
We also propose a novel approach for QDE modeling, fundamentally different from traditional ensemble modeling approaches. In the QDE approach, models are combined in the parameter space as opposed to the solution space. Several ETAS-like models are fit to the observed data, yielding an individual set of parameters for each model. A QDE model is then created by defining a new set of parameters based on a combination of the ingredient model parameters. The rules to combine parameters are defined by dividing the forecasting problem into several subproblems. Each subproblem addresses a question regarding the number of forecasted events or the spatiotemporal distribution of either background earthquakes or aftershocks. A QDE model can be viewed as a model that addresses different questions with different ingredient models. This approach allows the combination of ETAS variants but can be extended to combining more general types of seismicity models.
By including such QDE models in the forecasting experiments, we assess their forecasting capability in comparison with their ingredient models, standard ETAS, and flETAS. At the same time, the QDE approach helps to understand which ingredient models are best suited to solve different forecasting subproblems, thus, making it useful from the perspective of gaining new scientific insight.
The remainder of this article is structured as follows. We describe flETAS models and the QDE approach in the next section flETAS models. The setup for the forecasting experiments, the data analyzed, and the metrics used to evaluate forecasting performance are described in the Forecasting experiments section. We present and discuss our results in the Results and discussion section and finally provide our Conclusions.
The following subsections describe flETAS models and explain the QDE modeling. We begin by explaining the algorithm used to estimate the parameters of the ETAS model. Then, we describe how to relax some parametric assumptions of the ETAS model. Finally, we introduce a framework to create QDEs of flETAS models.
Consider an earthquake catalog: E Q -T A R G E T ; t e m p : i n t r a l i n k -; d f 1 ; 4 1 ; 3 5 3
consisting of events e i of magnitudes m i , which occur at times t i and locations x i ,y i .
The ETAS model describes earthquake rate as E Q -T A R G E T ; t e m p : i n t r a l i n k -; d f 2 ; 4 1 ; 2 7 6 λt,x,yjH t μ X i:t i <t gm i ,tt i ,xx i ,yy i : 2
That is, the sum of background rate μ and the rate of all aftershocks of previous events e i . The aftershock triggering rate gm,Δt,Δx,Δy describes the rate of aftershocks triggered by an event of magnitude m, at a time delay of Δt and a spatial distance (Δx,Δy) from the triggering event. We use here the definition:
E Q -T A R G E T ; t e m p : i n t r a l i n k -; d f 3 ; 4 1 ; 1 4 4 gm,Δt,Δx,Δy k 0 ×e am-m ref ×e -Δt=τ Δx 2 Δy 2 d×e γm-m ref 1ρ ×Δtc 1ω , 3 as in Nandan, Kamer, et al. (2021) and Mizrahi et al. (2021a). This formulation differs from other, more commonly used formulations of ETAS models in that it uses an Exponentially Tapered Omori Kernel (ETOK). In their article, Nandan, Kamer, et al. (2021) compare the ETAS model with ETOK to a more general version thereof, MDOK, which allows a magnitude dependency, finding that the more general version allows better forecasts. This indicates that including an exponential taper does lead to improved forecasts when compared to the commonly used Omori kernel. Besides allowing less heavy tails in the temporal distribution of aftershocks, this formulation of the Omori kernel makes it possible for the parameter ω to attain negative values, which is not possible in the traditional formulation. Also, our choice of this base model does not impact the main conclusions that can be drawn from comparing it to modified versions of itself.
To calibrate the ETAS model, the nine parameters to be optimized are the background rate μ and k 0 , a, c, ω, τ, d, γ, ρ, which parameterize the aftershock triggering rate g(m, t, x, y) given in equation ( 3). Implicitly, the model assumes that only earthquakes with magnitudes larger than or equal to m ref can trigger aftershocks. Most applications of the method define m ref as equal to the constant value of m c .
We build on the expectation maximization (EM) algorithm to estimate the ETAS parameters (Veen and Schoenberg, 2008). In this algorithm, the expected number of background events n and the expected number of directly triggered aftershocks li of each event e i are estimated in the expectation step (E step), along with the probabilities p ij that event e j was triggered by event e i , and the probability p ind j that event e j is independent. Following the E step, the nine parameters are optimized to maximize the complete data log likelihood in the maximization step (M step). E and M steps are repeated until convergence of the parameters. The usual formulation of the EM algorithm defines:
and st
with g kj gm k ,t jt k ,x jx k ,y jy k being the aftershock triggering rate of e k at location and time of event e j . For a given target event e j , equations ( 6) and ( 7) define p ij to be proportional to the aftershock occurrence rate g ij , and p ind j to be proportional to the background rate μ. As an event must be either independent or triggered by a previous event, the normalization factor Λ j : μ P k:t k <t j g kj in the denominator of equations ( 6) and ( 7) stipulates that p ind j P k:t k <t j p kj 1.
In the formulation of the ETAS model given in Equation ( 2), the rate of background earthquakes is described by the parameter μ, which does not vary with space nor time. During the maximization step of the EM algorithm, μ can be estimated independently from the other parameters as
in which A R and T denote the area of the study region and the length of the considered time window, respectively. In some approaches, the region of interest is divided into several subregions, which can have their own values for μ (Veen and Schoenberg, 2008). An iterative algorithm to estimate spatial variations of background rate based on maximum-likelihood estimation (Zhuang, 2012) uses a Gaussian kernel smoothing applied to the catalog event locations, weighted by their estimated independence probability, to obtain an estimate of μx,y. Here, we present a similar approach using EM, which has been shown to be more stable with respect to the initial conditions compared to maximum-likelihood approaches (Veen and Schoenberg, 2008). Our approach is similar yet not identical to the one described by Nandan, Ram, et al. (2021), which uses a regularized inverse power law for smoothing the locations. We define the background rate at a location (x, y) as
in which kΔx j ,Δy j is the Gaussian kernel with bandwidth σ applied to the distance Δx j ,Δy j of event e j to the location (x,y),
The bandwidth σ determines the smoothness of the background event density. In principle, σ could be calibrated itself, but we choose to fix it to 5 km for simplicity. Our next modification to the standard ETAS model is to allow flexibility of the aftershock probability. The number of directly triggered aftershocks lj is estimated during the expectation step of the EM algorithm as described in equation ( 5). We can, thus, replace the term k 0 × e am-m ref in equation ( 3) with κ j , in which κ j is stipulated to be proportional to lj . Instead of parameterizing aftershock productivity to be exponentially increasing with the magnitude of the triggering event, we allow each event to have its own productivity. This yields: E Q -T A R G E T ; t e m p : i n t r a l i n k -; d f 1 1 ; 3 2 0 ; 7 0 4 g j θ,κ j m,Δt,Δx,Δy
for given parameters θ c,ω,τ,d,γ,ρ and κ j . The EM algorithm is adapted as follows:
1. Define initial estimates of κ j as κ j e am j -m ref with a random guess for a; 2. define initial estimates of independence probability p ind j ≡ 0:1. The inversion result is not sensitive to this choice; 3. define random initial guesses for the parameters θ c,ω,τ,d,γ,ρ; 4. expectation step: calculate n, lj ,p ij ,p ind j using the current estimates of κ j ,θ, and p ind j . p ij ,p ind j are calculated using equations ( 6) and ( 7), but using the flexible definitions of g ij and μx,y of equations ( 9) and (11); 5. maximization step: optimize the parameters θ to minimize the complete data log likelihood (see Mizrahi et al., 2021a for details), given the current estimates of n, lj ,p ij ,p ind j ; 6. update κ new j to be κ old j × lj
is the expected total number of aftershocks of e j , given θ and κ old j . This ensures that lj G j θ,κ new j . We calculate G j θ,κ j as E Q -T A R G E T ; t e m p : i n t r a l i n k -; d f 1 2 ; 3 2 0 ; 3 7 7 G j θ,κ j ZZ R Z t end -t j 0 g j θ,κ j m j ,t,x,y dt dx dy, 12 in which t end is the end time of the considered time window, and we assume the spatial region R to extend infinitely in space, allowing a facilitated, asymptotically unbiased estimation of ETAS parameters (Schoenberg, 2013), and 7. repeat from step 4 until convergence of θ, that is, until
After the inversion, we calibrate an overall productivity law for the flETAS models with free productivity to avoid over fitting with event-wise productivity. From the individually estimated productivities κ j of magnitude m j events, we calibrate a law of the form:
by minimizing the sum of absolute residuals between the observed κm 1 nm P i:m i m κ i and the theoretical κm k 0 × e am-m ref , in which n(m) is the number of events with magnitude m.
Then, productivity is treated the same way as in the case of standard ETAS. In this way, the variability of productivity is only accounted for during the parameter inversion process and may lead to more accurate estimators of the productivity as well as the remaining ETAS parameters.
We propose a novel approach for QDE modeling, in which a forecast is created by combining model parameters of different ingredient models. The rules for how parameters can be combined are defined by questions that divide the forecasting problem into several subproblems: How many background events are expected? Where are they expected? When are they expected? How many aftershocks are expected? Where are they expected? When are they expected?
By answering each of these questions with different ingredient models, we create a suite of ensembles. The remainder of this section establishes rules to combine parameters based on the questions.
Consider a collection of ETAS or flETAS ingredient models, M i i0,…,n M . As they are sufficiently defined through their parameters, we can write:
In case M i is a flETAS model, μ i μ i x,y can vary with space. For simplicity, we denote with κ i the function that assigns to each event its appropriate value to replace the term κ j in equation ( 11). In our case, this means that we define κ i m k 0 i × e a i m-m ref , in which k 0 i and a i are either obtained during parameter inversion directly, or afterward in case M i is a flETAS model with free productivity. We chose the notation of κ i instead of (k 0 i ,a i ) to emphasize this possible distinction. We can then generally describe the aftershock triggering kernel g as E Q -T A R G E T ; t e m p : i n t r a l i n k -; d f 1 5 ; 4 1 ; 3 2 6 g i m,Δt,Δx,Δy
Let us now revisit the earlier questions.
1. How many background events are expected?
More precisely, what we want to ask here is how many background events do we expect in total in the region R and forecasting horizon T 0 ,T 1 we are issuing a forecast for.
The answer to this question, given out of the perspective of model M i , is E Q -T A R G E T ; t e m p : i n t r a l i n k -; d f 1 6 ; 4 1 ; 1 4 6
2. Where and when are they expected?
We address for now these two questions jointly. The spatiotemporal density of background events is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; d f 1 7 ; 3 0 8 ; 7 4 3 f B i x,y,t
which is effectively time independent due to our choice of a time independent μx,y.
Again, what we want to ask here is how many aftershocks do we expect in total in the region R and forecasting horizon T 0 ,T 1 we are issuing a forecast for. For an individual event e j , we expect it to have n A aftershocks, in which
T 0 g i m j ,tt j ,xx j ,yy j dt dx dy: 18
The total number of aftershocks N A i is then given as the sum of aftershocks of all events E Q -T A R G E T ; t e m p : i n t r a l i n k -; d f 1 9 ; 3 0 8 ; 5 3 5
4. Where and when are they expected?
We again answer these two questions jointly. If we define E Q -T A R G E T ; t e m p : i n t r a l i n k -; d f 2 0 ; 3 0 8 ; 4 5 7 G i x,y,t :
as the total rate of aftershocks at time t and location (x,y), consisting of the sum of aftershock rates of all events that occurred prior to the end T 1 of the forecasting horizon, the spatiotemporal density of aftershocks is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; d f 2 1 ; 3 0 8 ; 3 5 3 f A i x,y,t
We now construct a QDE model E klm as follows. The number questions (1) and ( 3) are answered with model M k , the background density question ( 2) is answered with model M l , and the aftershock density question ( 4) is answered with model M m . Questions (1) and (3) are addressed with the same model. This is a choice made to avoid unrealistic event numbers. If one model interprets the majority of events as background, and another model interprets the majority of events to be aftershocks, answering the two questions with two different models would lead to exceptionally high or low total event numbers, which is not intended by the two ingredient models.
In the earlier notation, which identifies a model with its parameters, this would give us:
Forecasting Experiments
To test whether flETAS models and QDE models, which consist of ETAS and flETAS models, provide better forecasts, we conduct pseudoprospective forecasting experiments for southern California and Italy.
In these experiments, we consider the following four competing ingredient models.
• M 0 : standard ETAS.
• M 1 : flETAS with free productivity and standard background.
• M 2 : flETAS with standard productivity and free background.
• M 3 : flETAS with free productivity and free background.
Out of these, 4 3 64 QDE models can be constructed.
M 2 is conceptually close to the models described by Zhuang (2012) and Nandan, Ram, et al. (2021).
We use interevent time horizons: whenever an event occurs, a forecast is issued, which is valid until the occurrence of the next event. A pseudoprospective model evaluation then aims to capture how well a forecast issued using data until event e j-1 can describe the occurrence of the next event e j .
An ETAS forecast always consists of the forecasted background seismicity rate plus the forecasted aftershock seismicity rate. With this flexible definition of forecasting horizon, our ETAS forecast can be calculated and evaluated analytically.
Consider λ i t,x,yjH t j-1 , the event rate under model M i as of time t j-1 of the (j-1)th earthquake. This formulation of λ i is valid for times t ∈ t j-1 ,t j between the occurrence of event e j-1 and event e j , and hence, this is the forecasting horizon we consider.
For the traditional experiment settings for which one is interested in the seismicity forecast of the next days, months, or years, such an analytical description of the forecasted seismicity is not possible. As soon as an event occurs during the forecasting period, its aftershocks are not part of the background seismicity, nor of the aftershock seismicity that was calculated at the start of the forecasting period. For this reason, ETAS forecasts for fixed forecasting horizons are usually produced through the simulation of a large number of possible continuations of the catalog.
In our case of flexible forecasting horizons, the log likelihood of observing e j under model M i is analytically defined (see Ogata et al., 2013;Daley and Vere-Jones, 2003) as E Q -T A R G E T ; t e m p : i n t r a l i n k -; d f 2 3 ; 5 3 ; 1 4 6 ln L i e j ln λ i t j ,x j ,t j jH t j-1 -ZZ R Z t j t j-1 λ i t j ,x j ,t j jH t j-1 dt dx dy:
23
We then define the information gain IG i 1 ,i 2 j of model i 1 over model i 2 during the jth forecasting period t j-1 ,t j as E Q -T A R G E T ; t e m p : i n t r a l i n k -; d f 2 4 ; 3 2 0 ; 7 4 3 IG i 1 ,i 2 j ln L i 1 e j L i 2 e j ln L i 1 e jln L i 2 e j : 24
The information gain per event (IGPE) over forecasting periods j 1 ,…,j K is defined as E Q -T A R G E T ; t e m p : i n t r a l i n k -; d f 2 5 ; 3 2 0 ; 6 7 9
the average of IGs over those testing periods.
Compared to evaluation techniques based on the simulation of large numbers of possible catalog continuations such as in Nandan et al. (2019a) and Mizrahi et al. (2021a), which are encouraged by CSEP (see Savran et al., 2022), this approach allows us to compare models much faster, accelerating the development and testing process. To apply these models operationally, in which forecasts are required for a fixed time horizon, simulations would still be required. This evaluation approach allows us to save time when developing and selecting the model to be used operationally and is especially useful for evaluating a large suite of QDE models.
For southern California, we consider the Advanced National Seismic System (ANSS) comprehensive earthquake catalog (ComCat), in the polygon given by the vertices in Table A1. We consider earthquakes of magnitude M ≥ 2.0 from 1 January 2010 until 1 January 2022. The first two years serve as an auxiliary period in the ETAS and flETAS parameter inversion, and thus, the start of the primary catalog is 1 January 2012. This means that the events between January 2010 and January 2012 can act as triggering events during the inversion but not as triggered events. Using the method described by Mizrahi et al. (2021b), we find that the overall catalog is complete at this threshold, although there are likely periods during which the catalog is incomplete due to shortterm aftershock incompleteness (STAI). Although Mizrahi et al. (2021a) have proposed a method to account for STAI in the ETAS model, we do not address this issue here.
For Italy, we consider the Italian Seismological Instrumental and Parametric Data-Base catalog (ISIDe, Group, 2007), in the area defined for the first CSEP experiment (Schorlemmer, Christophersen, et al., 2010, vertices given in Table A2). We consider earthquakes of magnitude M ≥ 2.5 from 16 April 2005 until 1 July 2021. This is the time horizon available to modelers in the upcoming prospective CSEP forecasting experiment in Italy, and the estimated magnitude of completeness provided in the experiment description. The start of the primary catalog is 1 January 2010.
For southern California, we consider 5 yr of testing, with the start of the first forecasting period at the occurrence of event e 0 , the first event on or after 1 January 2017. In Italy, we consider 3 yr of testing, starting at the occurrence of the first event on or after 1 July 2018. The idea of the pseudoprospective experiments is to only use data that would have been available at the time the forecast is issued to calibrate the models. One could, thus, recalibrate the model at the start of each forecasting period, whenever one more event becomes part of the catalog. To limit the number of computationally expensive parameter inversions for these experiments, we re-estimate the model parameters every 7 days in southern California, and every day in Italy, and use the latest available set of parameters at the start time of each forecasting interval. This does not mean that events between the calibration time and forecasting start are ignored. Their aftershocks are still considered in the calculated aftershock rate. We chose a shorter parameter updating interval for Italy to mimic the conditions of the CSEP experiment, and a longer one for southern California to limit computational cost.
We then calculate IG i 1 ,i 2 j for all j, and for all pairs of models M i 1 , M i 2 . If the IGPE over all forecasting periods of one model to another is positive, we consider the model to produce superior forecasts.
As one could argue that generating a large number of models and then selecting the best performing ones somewhat invalidates the pseudoprospective nature of our experiments, we consider the following additional model. At the start of the jth forecasting period, the total information gain of all QDE models during the last n forecasting periods, that is, periods j-(n + 1) to j-1, is compared. The model with the highest IG is selected to produce the forecast for the jth forecasting period. We call this model QDE-S n .
This type of model, if capable of producing a powerful forecast, would be well suited to be used in an OEF context.
The parameters that were obtained using the flETAS inversion algorithm are described in the Inverted parameters section. Here, we present the results of the Forecasting Experiments.
Figure 1 compares the IGPE over the standard ETAS null model (M 0 E 000 ) of all 64 QDE models in Italy and southern California. The IGPE varies between -0.64 and 0.45 in Italy, and between -0.13 and 0.12 in southern California. The best and worst performing QDE models are E 221 and E 112 , respectively, for both regions. The best performing model E 221 uses the free background model M 2 to answer the number and background density questions, and the free productivity model M 1 to answer the aftershock density question. Vice versa, the worst performing model E 112 uses M 1 to answer the number The symbol shape, fill color, and edge color describe the composition of the QDE. The shape, fill color, and edge color represent the ingredient model used to answer the background density (BG), number (N), and aftershock density (AS) questions, respectively. Box plots on top (for southern California) and to the right (for Italy) of the scatter plot: for N, BG, and AS questions, the four boxes represent the IGPE of four groups of QDE models. Each group contains the 16 QDE models that use a specific ingredient model (indicated by box color) to answer the indicated question.
and background density questions, and model M 2 to answer the aftershock density question. Generally, the models that perform well or poorly in Italy are also performing similarly in southern California.
The symbol shape, fill color, and edge color in the scatter plot of Figure 1 represent the ingredient model used to answer the background density (BG), number (N), and aftershock density (AS) questions, respectively. Models that perform well tend to answer the BG question with the free background ingredient model and the AS question with the free productivity model. Conversely, models that address the BG question with the free productivity model, and those that address the AS question with the free background model, tend to perform poorly. This is highlighted in the box plots of Figure 1. There, for each question, the distribution of IGPE of the 64 QDE models is given per possible answer. Although for the number questions, no clear trend can be inferred, it is evident that the free background model serves well at answering the BG question and the free productivity model serves well at answering the AS question. These trends are qualitatively very similar in southern California and Italy.
These results emphasize the added value generated by the flETAS approach, although most flETAS models individually do not outperform standard ETAS. Apparently, a model that gives full flexibility to the background rate during parameter inversion is more informative than others when addressing the background density question. And a model that is flexible at identifying aftershocks is more informative than others when answering the aftershock density question. These observations are made for both considered regions.
While conceptually it makes sense that a model that can more flexibly capture one particular aspect of seismicity is particularly successful at questions about this very aspect of seismicity, this is simultaneously a somewhat counterintuitive result. If flETAS with free background is more than other models at identifying background events, one would expect it, due to the self-consistent nature of parameter inversion, to also be more successful at identifying aftershocks, and thus at describing their occurrence times and locations.
A possible interpretation of the observation that E 221 , E 220 , and even E 223 can so clearly outperform E 222 , is the following. Compared to the null model M 0 , model M 2 E 222 allows the background seismicity to be free and, therefore, interprets a higher fraction of events in the training catalog to be background earthquakes, which manifests in a much higher background rate. M 2 can, thus, explain the spatial distribution of background events well, as well as the partitioning of seismicity into background events and aftershocks. Possibly, M 2 overestimates the background portion of the training catalog due to "too much freedom." The level of overestimation may be small enough so that M 2 still captures the fraction and locations of background earthquakes better than the other ingredient models do. Overestimation of the background seismicity comes with underestimation of the fraction of aftershocks in the training catalog. Although this underestimation may have a minor biasing effect on the number of background earthquakes and aftershocks, the spatiotemporal distribution of aftershocks can be affected in a more harmful way. Aftershocks that occur in the tails of the spatial or temporal distributions have higher chances to be falsely identified as background events compared to aftershocks that are close to their parent event. This leads to a distorted characterization of the aftershock triggering behavior of model M 2 , which can be fixed using the triggering parameters from models M 0 or M 1 , as indicated by the good performance of models E 221 and E 220 .
Another noteworthy observation is that model M 3 , which in principle has all the flexibility necessary to encompass the parameterization of model E 221 , is clearly outperformed by E 221 . We interpret this to be a consequence of the fact that the information that is optimized during model calibration and the information used for forecasting are not the same. This does not indicate a flaw in the method presented, but rather illustrates a complexity of the forecasting problem to which the QDE approach offers an apparently useful solution.
Figure 2a shows the cumulative information gain (CIG) over the standard ETAS model over time of the three flETAS ingredient models and the three best performing QDE models. The CIG of model i 1 over model i 2 at time t is given as the sum of IGs of all forecasting periods ending prior to time t: E Q -T A R G E T ; t e m p : i n t r a l i n k -; d f 2 6 ; 3 2 0 ; 3 6 6 X j:t j <t IG i 1 ,i 2 j : 26
In southern California, the flETAS ingredient models have a negative information gain following the Ridgecrest events in July 2019, meaning that during this time, the standard ETAS model (M 0 ) is better performing. The free background model M 2 outperforms M 0 immediately after the onset of the sequence and suffers from information loss later during the sequence. The other two ingredient models do not exhibit the initial information gain. Among the flETAS models, only M 2 can compensate for the information loss during the course of the 5 yr of testing and ends up with a positive overall information gain.
Among the QDE models presented, models E 221 and E 220 show an initial information gain after the onset of the Ridgecrest sequence, followed by a period of information loss. In contrast to the ingredient models, the information loss during the sequence is smaller than the gain at the beginning of the sequence, such that these models show positive information gain during the Ridgecrest sequence. The three QDE models in Figure 2a also show a rapidly accumulating information gain throughout the testing period, arriving at an overall IGPE of 0.12, 0.10, and 0.09.
From Figure 2b, it is clear that the IGPE is relatively close to zero in the Ridgecrest area, and the positive IG during the sequence must come from a few specific locations. In the rest of southern California, higher IGPE values are achieved, with a median grid-cell-wise IGPE of 0.66 for model E 221 shown in Figure 2b. Conversely, the median grid-cell-wise IGPE for the worst performing model E 112 shown in Figure 2c is -0.54. Generally, it performs poorly where E 221 performs well.
In Italy, all flETAS models have negative total information gain over M 0 . Nevertheless, two of the top three QDE models that perform best in southern California are also among the top three in Italy, with overall IGPE values of 0.
Figure 3 illustrates the composition and performance of QDE-S n models. The number n of past forecasting periods considered when selecting the forecasting model for the next period is in f1 2 0 , 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024 2 10 g for southern California, and n ∈ f1 2 0 ,…,512 2 9 g for Italy. We do not consider n = 1024 for Italy, as this would reduce the number of testing periods in which QDE-S n is defined by more than half compared to the QDE models. The top, middle, and bottom parts of Figure 3a,b show the ingredient model used by QDE-S n to answer the N, BG, and AS questions over time. Within each part, n increases from top to bottom. As expected, the composition of QDE-S n is more stable as n increases and is almost always defined via E 221 for large n, in both regions.
In southern California, a change in composition can be observed after the onset of the Ridgecrest sequence in July 2019. Specifically, the number questions are best answered by standard ETAS, free productivity flETAS, and free productivity and background flETAS, in this order, before moving back to answering with free background flETAS. The aftershock question is intermittently best answered by standard ETAS during the sequence. It is interesting to note here that the performance of E 221 and QDE-S 64 are almost identical throughout the 5 yr of testing, with the difference that QDE-S 64 does not show the information loss after the initial information gain after the onset of the sequence. This results in an overall IGPE of 0.13 and 0.12 for QDE-S 64 and E 221 , during the period in which both are defined, as is shown in Figure 3c. Thus, the QDE-S n model, which was originally designed to avoid a biased selection of the winning model after knowing the experiment outcome, is capable of outperforming the winning QDE model for good choices of n, and clearly outperforms all ingredient flETAS models for any tested choice of n.
In Italy, the best performing QDE-S n model is QDE-S 128 . It is almost always using E 221 to issue a forecast for the next period and, thus, unsurprisingly achieves the same IGPE. As in southern California, all tested choices of n yield a model that clearly outperforms all ingredient flETAS models. The simplest QDE-S n model, QDE-S 1 , which always selects the best QDE model of the previous forecasting period to issue the next forecast, already achieves a very high IGPE of 0.28.
We describe an adapted ETAS EM algorithm that allows a nonparametric inversion of aftershock productivity and/or background rate. Further, we introduce a novel approach of QDE modeling, which combines ingredient models by using them to answer different forecasting subproblems. In pseudoprospective forecasting experiments for southern California and Italy, we compare the forecasting skill of three flETAS models and a total of 60 nontrivial QDEs of flETAS and ETAS models, to that of the standard ETAS null model.
We find that the best models tend to use flETAS with free background to model the number of events and locations of background earthquakes and flETAS with free productivity to model the times and locations of aftershocks. The best model is the same in both regions and achieves an IGPE over standard ETAS of 0.12 in southern California and 0.45 in Italy. To address the possible concern of a biased selection of the winning model after knowing the experiment outcome, we also test the forecasting skill of a model that pseudoprospectively selects the currently best performing QDE model to issue the forecast for the next testing period. Depending on the criteria to identify the best QDE model, we find that the forecasting skill can be greater than that of the overall best QDE model. This approach thus provides a promising candidate for an operational earthquake forecast.
During the 2019 Ridgecrest sequence in southern California, different ingredient models are best suited to model the number of events during different stages of the sequence. The idea of operationally selecting different QDE models (i.e., selecting different ETAS model parameters) based on their recent performance is in this case related to the idea of Page et al. (2016). They considered sequence-specific parameters to be sampled from an underlying distribution and described a Bayesian approach to update this distribution as aftershock data become available.
Our results can also be viewed as a first step toward developing a potentially fruitful branch of earthquake forecasting research. Several key questions remain open and are to be addressed in future studies: Why do QDE models outperform ingredient models that were inverted in a self-consistent way? What drives the success of different QDE models during different phases of the Ridgecrest sequence? How does QDE performance increase when further ingredient models are considered? And what does all of this teach us about the dynamics of seismicity?
Downloaded from http://pubs.geoscienceworld.org/ssa/srl/article-pdf/94/2A/829/5793422/srl-2022230.1.pdf by University of Southern California user
Volume 94 • Number 2A • March 2023 • www.srl-online.org Seismological Research Letters