Developing, Testing, and Communicating Earthquake Forecasts: Current Practices and Future Directions Leila Mizrahi1 , Irina Dallo1, Nicholas J. van der Elst2 , Annemarie Christophersen3 , Ilaria Spassiani4 , Maximilian J. Werner5 , Pablo Iturrieta6 , José A. Bayona5 , Iunio Iervolino7,8, Max Schneider9, Morgan T. Page2 , Jiancang Zhuang10 , Marcus Herrmann7 , Andrew J. Michael9 , Giuseppe Falcone4 , Warner Marzocchi7 , David Rhoades3, Matt Gerstenberger3 , Laura Gulia11, Danijel Schorlemmer6, Julia Becker12, Marta Han1 , Lorena Kuratle1, Michèle Marti1 , and Stefan Wiemer1 1Swiss Seismological Service at ETH Zurich, Zürich, Switzerland, 2US Geological Survey, Pasadena, CA, USA, 3GNS Science | Te Pū Ao, Lower Hutt, New Zealand, 4Istituto Nazionale di Geofisica e Vulcanologia (INGV), Rome, Italy, 5School of Earth Sciences, University of Bristol, Bristol, UK, 6GFZ German Research Centre for Geosciences, Potsdam, Germany, 7Università degli Studi di Napoli Federico II, Naples, Italy, 8IUSS – Scuola Universitaria Superiore di Pavia, Pavia, Italy, 9US Geological Survey, Moffett Field, CA, USA, 10The Institute of Statistical Mathematics, Tokyo, Japan, 11Dipartimento di Fisica e Astronomia, Università di Bologna, Bologna, Italy, 12Joint Centre for Disaster Research, Massey University, Wellington, New Zealand Abstract While deterministically predicting the time and location of earthquakes remains impossible, earthquake forecasting models can provide estimates of the probabilities of earthquakes occurring within some region over time. To enable informed decision‐making of civil protection, governmental agencies, or the public, Operational Earthquake Forecasting (OEF) systems aim to provide authoritative earthquake forecasts based on current earthquake activity in near‐real time. Establishing OEF systems involves several nontrivial choices. This review captures the current state of OEF worldwide and analyzes expert recommendations on the development, testing, and communication of earthquake forecasts. An introductory summary of OEF‐related research is followed by a description of OEF systems in Italy, New Zealand, and the United States. Combined, these two parts provide an informative and transparent snapshot of today's OEF landscape. In Section 4, we analyze the results of an expert elicitation that was conducted to seek guidance for the establishment of OEF systems. The elicitation identifies consensus and dissent on OEF issues among a non‐representative group of 20 international earthquake forecasting experts. While the experts agree that communication products should be developed in collaboration with the forecast user groups, they disagree on whether forecasting models and testing methods should be user‐dependent. No recommendations of strict model requirements could be elicited, but benchmark comparisons, prospective testing, reproducibility, and transparency are encouraged. Section 5 gives an outlook on the future of OEF. Besides covering recent research on earthquake forecasting model development and testing, upcoming OEF initiatives are described in the context of the expert elicitation findings. Plain Language Summary The exact location, time, and magnitude of future earthquakes cannot be predicted. However, based on past earthquake sequences, it is possible to assess probabilities for future earthquakes. This is called earthquake forecasting. Operational Earthquake Forecasting (OEF) systems are designed to provide near‐real‐time authoritative earthquake forecasts, based on current earthquake activity, to aid the decision‐making of various societal stakeholders. Setting up these systems is complex, involving decisions about which model to use, how to best test the model, and how to turn earthquake probability estimates into practical information. This review captures the current state of OEF worldwide and analyzes expert recommendations on the development, testing, and communication of earthquake forecasts. Section 2 provides an overview of OEF‐related research and the background knowledge required to understand the other parts. Section 3 describes existing OEF systems of Italy, New Zealand, and the United States in detail. Section 4 discusses an elicitation of expert views on modeling, testing, and communicating earthquake forecasts (Mizrahi, Dallo, & Kuratle, 2023, https://doi.org/10.3929/ethz‐b‐000637239). Data from the elicitation allow to identify consensus and dissent on OEF issues and provide guidance for future earthquake forecasting efforts. Finally, Section 5 gives an outlook on future OEF‐related research and planned OEF efforts at various institutions. REVIEW ARTICLE 10.1029/2023RG000823 Key Points: • We capture the state of earthquake forecasting systems in Italy, New Zealand, and the United States, and future plans in these and other countries • Experts encourage benchmark comparison, prospective testing, reproducibility and transparency, but avoid endorsing specific models or tests • Experts stress the need to co‐design forecast communication products with end‐users to ensure their societal rele- vance and usefulness Supporting Information: Supporting Information may be found in the online version of this article. Correspondence to: L. Mizrahi, leila.mizrahi@sed.ethz.ch Citation: Mizrahi, L., Dallo, I., van der Elst, N. J., Christophersen, A., Spassiani, I., Werner, M. J., et al. (2024). Developing, testing, and communicating earthquake forecasts: Current practices and future directions. Reviews of Geophysics, 62, e2023RG000823. https://doi.org/10.1029/ 2023RG000823 Received 3 NOV 2023 Accepted 18 JUL 2024 Author Contributions: Conceptualization: Leila Mizrahi, Irina Dallo, Stefan Wiemer Funding acquisition: Stefan Wiemer Methodology: Leila Mizrahi, Irina Dallo, Lorena Kuratle Project administration: Leila Mizrahi, Stefan Wiemer Supervision: Stefan Wiemer Writing – original draft: Leila Mizrahi, Irina Dallo, Nicholas J. van der Elst, Annemarie Christophersen, © 2024 The Author(s). This article has been contributed to by U.S. Government employees and their work is in the public domain in the USA. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. MIZRAHI ET AL. 1 of 70 https://orcid.org/0000-0002-5262-3168 https://orcid.org/0000-0002-3812-1153 https://orcid.org/0000-0003-1467-1414 https://orcid.org/0000-0003-2252-0202 https://orcid.org/0000-0002-2430-2631 https://orcid.org/0000-0002-4787-1343 https://orcid.org/0000-0002-3026-0327 https://orcid.org/0000-0001-9321-2990 https://orcid.org/0000-0002-9708-3871 https://orcid.org/0000-0002-2342-1970 https://orcid.org/0000-0002-2403-5019 https://orcid.org/0000-0002-2554-4421 https://orcid.org/0000-0002-9114-1516 https://orcid.org/0000-0002-0392-7114 https://orcid.org/0009-0006-8774-9170 https://orcid.org/0000-0002-4997-1153 https://orcid.org/0000-0002-4919-3283 https://doi.org/10.3929/ethz-b-000637239 mailto:leila.mizrahi@sed.ethz.ch https://doi.org/10.1029/2023RG000823 https://doi.org/10.1029/2023RG000823 http://creativecommons.org/licenses/by/4.0/ http://creativecommons.org/licenses/by/4.0/ 1. Introduction What is commonly referred to as earthquake prediction, namely the deterministic description of the location, time, and magnitude of a future earthquake “within narrow limits […], so that a planned evacuation can take place” (Main, 1999), is currently not possible. The term earthquake forecasting refers to a probabilistic assessment of earthquake occurrence, over a range of magnitudes, and within some specified time frame and region (e.g., Kagan & Jackson, 2000). On the spectrum between phenomena of complete randomness, such as a fair coin toss, and entirely deterministic processes, the occurrence of earthquakes falls somewhere in‐between, its most predictable components being the clustering behavior of groups of earthquakes in time and space (Kagan & Jackson, 1991; Omori, 1894; Zhuang et al., 2021), and the frequency distribution of earthquake magnitudes within a group (Gutenberg & Richter, 1944). The understanding of occurrence patterns of earthquakes and their implications for the short‐term increase in seismic hazard and risk can provide crucial information during an ongoing seismic crisis and enable informed decision‐making of civil protection, governmental agencies, the public, or other user groups. As first defined by Jordan et al. (2011), “Operational earthquake forecasting comprises procedures for gathering and disseminating authoritative information about the time dependence of seismic hazards to help communities prepare for potentially destructive earthquakes.” The label operational has been borrowed from meteorology, where it de- scribes the existence of automated and timely applications of forecasting models, which increases the action- ability of the information provided. To be considered authoritative, information must originate from an official, trusted source that is legally responsible for providing this information. Thus, the authoritativeness of an in- formation source is a quality that cannot be self‐assigned but must be granted by a higher authority, such as the government. For the remainder of this review, we will use a broader definition of operational earthquake forecasting (OEF) as authoritative near‐real‐time applications of earthquake forecasting. The concept of near‐ real‐time here refers to the delivery of information as soon as possible when an interest in updated information arises, acknowledging the possibility of some noticeable, unavoidable processing time. This is in contrast to the concept of true real‐time, which would refer to the delivery of information as soon as needed. Few governmental agencies worldwide have systems in place that fall under this OEF definition, and one could argue that only Italy has fully deployed an OEF system under the definition of Jordan et al. (2011), even though their OEF information has not yet been open to the public. The Italian OEF system was established in 2013 and has produced weekly forecasts since then (Marzocchi et al., 2014), including the expected seismic losses based on such forecasts (Operational Earthquake Loss Forecasting, OELF; Iervolino et al., 2015). The forecasts can be accessed only by authorized personnel and, only recently, the OEF information has been released to regional governments to trigger a discussion on how to use and communicate this information in different contexts. In New Zealand, earthquake forecasts that satisfy our broader OEF definition have been provided to the public since the 2010 M7.1 Darfield earthquake, and the content of these forecasts has evolved over time (Christophersen, Rhoades, & Colella, 2017; Christophersen, Rhoades, Gerstenberger, et al., 2017). The forecasts disseminated by the U.S. Geological Survey (USGS) since 1989 for California (Reasenberg & Jones, 1989) and since 2018 for the entire country (Michael et al., 2020), although fulfilling our OEF definition, are referred to as aftershock fore- casts, to emphasize that forecasts are only issued after the occurrence of large events. Here, we consider (operational) aftershock forecasting (OAF) to be a differently named special case of (operational) earthquake forecasting. There is currently no known difference in the physical mechanisms of what are often termed foreshocks, after- shocks, or mainshocks, that is, the earthquakes in a cluster, and their distinction can be considered a convention (Felzer et al., 2004). These labels can only be applied to earthquakes of a sequence in hindsight: the earthquake with the largest magnitude is generally referred to as the mainshock, while earthquakes that happened earlier than the mainshock are retrospectively defined to be foreshocks, and earthquakes that happened after the largest event are called aftershocks. However, because this distinction seems not to exist in nature, the classification cannot be done before the earthquake sequence has finished. During an ongoing earthquake sequence (the definition of which is itself ambiguous), it cannot be ruled out that a larger earthquake is yet to follow, rendering the mainshock‐aftershock terminology often misleading. Between 1980 and 2019, more than 10 percent of world- wide M ≥ 6.0 earthquakes were followed by a larger one within a space‐time window of 60 days and 100 km (Taroni, 2023). Aftershock forecasts published by the USGS clearly emphasize this possibility and it is a question of definition whether a larger earthquake following a first mainshock should be considered a second mainshock, a Ilaria Spassiani, Pablo Iturrieta, José A. Bayona, Iunio Iervolino, Jiancang Zhuang, Marcus Herrmann, David Rhoades, Laura Gulia Writing – review & editing: Leila Mizrahi, Irina Dallo, Nicholas J. van der Elst, Annemarie Christophersen, Ilaria Spassiani, Pablo Iturrieta, José A. Bayona, Iunio Iervolino, Max Schneider, Morgan T. Page, Jiancang Zhuang, Marcus Herrmann, Andrew J. Michael, Giuseppe Falcone, Warner Marzocchi, David Rhoades, Matt Gerstenberger, Danijel Schorlemmer, Julia Becker, Marta Han, Michèle Marti, Stefan Wiemer Reviews of Geophysics 10.1029/2023RG000823 MIZRAHI ET AL. 2 of 70 19449208, 2024, 3, D ow nloaded from https://agupubs.onlinelibrary.w iley.com /doi/10.1029/2023R G 000823 by M assey U niversity L ibrary, W iley O nline L ibrary on [23/10/2024]. See the T erm s and C onditions (https://onlinelibrary.w iley.com /term s-and-conditions) on W iley O nline L ibrary for rules of use; O A articles are governed by the applicable C reative C om m ons L icense https://agupubs.onlinelibrary.wiley.com/action/rightsLink?doi=10.1029%2F2023RG000823&mode= large aftershock, or whether the previous mainshock is now a foreshock to the only true mainshock. In Japan, aftershock forecasts ignoring the possibility of a strong earthquake being followed by a stronger one were pro- duced following 15 earthquakes between 1998 and 2016. After the 2016 Kumamoto sequence, in which a M6.5 event was followed by a M7.3 event 28 hr later, the Japan Meteorological Agency (JMA) stopped providing aftershock forecasts during the first week after a mainshock because of the “invalidity of the [previously applied] procedure” (Omi et al., 2019). The notion of foreshocks can likewise be misleading. Many mainshocks occur without detectably anomalous prior seismic activity, and vice‐versa, a temporal increase of seismic activity is not always followed by a large earthquake (Mogi, 1963). Series of earthquakes located close to each other in time and space that are not dominated by a clear mainshock are commonly referred to as earthquake swarms and are often associated with aseismic slip or slow‐slip events, fluid intrusions or magmatic processes (Hainzl, 2004; Hill et al., 1975; Hirose et al., 2014; Kisslinger, 1975; Passarelli et al., 2021; Ross & Cochran, 2021; K. D. Smith et al., 2004; Vidale & Shearer, 2006). Although the differentiation of swarms from other, mainshock‐dominated clusters is accompa- nied by proposed underlying physical mechanisms, the classification of earthquake clusters into mainshock‐ dominated clusters or swarms is, again, not free of ambiguity (Zaliapin & Ben‐Zion, 2016). Section 2 of this review provides an overview of the current state of research relevant for OEF. Besides an introduction to different model types currently applied in practice, it covers the history and necessary background information on forecast model testing as well as the state of the art of earthquake forecast communication. Thereafter, the current state of OEF systems in Italy, New Zealand, and the United States is described in Section 3 of this review. Currently, most countries, even in regions considered of high seismic risk, do not have OEF systems in place and only estimate the time‐invariant, long‐term earthquake probability as part of their seismic hazard and risk assessment, which serves different user groups involved in long‐term decision making (see Gerstenberger et al., 2020). In addition to long‐term seismic hazard programs, so‐called earthquake early warning (EEW) systems have been and are being established in various countries (e.g., Given et al., 2014; Nof & Kurzon, 2021; Orihuela et al., 2023; Suárez, 2022; see Cremen and Galasso (2020) for a review on EEW). These systems can disseminate warnings about earthquakes after they have already occurred, but, depending on the location and strength of the earthquake and the efficiency of the warning system, before the damaging seismic waves reach populated cities or critical infrastructure. The scarcity of OEF systems could be attributed to factors such as insufficient data, knowledge gaps, resource constraints, or a lack of guidance required for their establishment. Amidst a wealth of forecasting models and model tests proposed and published every year, it may be difficult to choose the most appropriate ones. On top of the more technical choices involved in the development and testing of forecasting models comes the additional challenge of how to communicate earthquake forecasts to different user groups. In New Zealand, in particular, much effort has gone into testing and evaluating forecast communication products with end‐users, providing initial guidance for other regions as well (e.g., Becker et al., 2020). In Europe and the United States, recent studies have similarly aimed to assess user needs, design preferences, and the impact of the communication format on people's risk perception (Dryhurst et al., 2022; Schneider, Wein, et al., 2023). It is crucial that earthquake forecasts serve a variety of stakeholders, are adapted over time as information needs evolve throughout the course of earthquake sequences, and are provided in different formats to consider stakeholders' knowledge, skills and requirements. In this process, two main challenges arise: (a) effectively communicating what can be considered low‐probability, high‐impact events and supporting stakeholders in interpreting probabilistic information and translating it into risk and risk mitigation actions, if needed; and (b) obtaining a thorough understanding of how earthquake forecasts influence people's emotions and, consequently, their behavior and well‐being. Section 4 of this review describes an elicitation of expert views on the three pillars Model Development, Model Testing, and Forecast Communication based on a Delphi study conducted among a group of experts in OEF (Mizrahi, Dallo, & Kuratle, 2023). The study involved two surveys and an online workshop, which facilitated meaningful discussions on the three pillars. The study aimed to establish areas of consensus among the experts and identify open questions within the research landscape to highlight future research directions. The results of the Delphi study provide a transparent and structured collection of expert opinions, offering valuable guidance for the future development of OEF systems. Reviews of Geophysics 10.1029/2023RG000823 MIZRAHI ET AL. 3 of 70 19449208, 2024, 3, D ow nloaded from https://agupubs.onlinelibrary.w iley.com /doi/10.1029/2023R G 000823 by M assey U niversity L ibrary, W iley O nline L ibrary on [23/10/2024]. See the T erm s and C onditions (https://onlinelibrary.w iley.com /term s-and-conditions) on W iley O nline L ibrary for rules of use; O A articles are governed by the applicable C reative C om m ons L icense https://agupubs.onlinelibrary.wiley.com/action/rightsLink?doi=10.1029%2F2023RG000823&mode= Section 5 of this review outlines emerging forecasting techniques potentially relevant for OEF in the future, an outlook on forecast testing efforts, as well as forthcoming developments in earthquake forecasting at various institutions represented by the authors. These developments are presented in the context of the guiding framework derived from the Delphi study, demonstrating the practical implications of the study's findings and emphasizing how they contribute to shaping the future of OEF systems. 2. The Current State of Research Relevant for Operational Earthquake Forecasting This Section 2 of the review introduces concepts, terms, and models that will be discussed in the subsequent sections: Examples of OEF Systems Worldwide, Elicitation of Expert Views, and Outlook. Readers who are familiar with the topic may wish to skip this section or parts of it. 2.1. Background on Earthquake Forecasting Models Used for OEF 2.1.1. Basics of Statistical Seismology The study of earthquake occurrence patterns goes back a long time. In 1756, Immanuel Kant noted that the frequency of earthquakes occurring in a region is more related to its proximity to mountain ranges and “fire‐ spitting mountains” than to the degree of Christianness of its inhabitants (Kant, 1756). Since then, descriptions of such empirical observations have become more precise and quantitative; the field of statistical seismology has emerged. One of its most central (empirical) laws, the Omori‐Utsu law (Omori, 1894; Utsu, 1961), describes a possible trend in the decrease of the expected number of aftershocks in time after a parent event as N(t) = k (t + c)p , (1) with parameters k, c, and p. The number of aftershocks (or earthquakes) per unit time, here N(t), is referred to as the aftershock rate (or earthquake rate). In short‐ or medium‐term forecasting, it is often provided in the unit of earthquakes per day. The parameter k in Equation 1 describes the overall number of aftershocks of a certain parent event. Note that we use the term parent event (rather than mainshock) to emphasize the possibility that the parent event can trigger aftershocks that are larger than itself. The exponent p of the denominator in Equation 1 quantifies how fast the aftershock rate decays with time—large values of p indicate a faster decay, and small values indicate a slower decay. Typically, p takes values of around 1 (Page et al., 2016), which corresponds to an aftershock rate decay that behaves like t− 1. The parameter c avoids a singularity at time t= 0, and interpretations of the parameter range from it being an artifact of more noisy data recordings right after a large earthquake (Hainzl, 2016a, 2016b; Kagan, 2004; Lolli & Gasperini, 2006) to it being a consequence of the physics of earthquakes (Lippiello et al., 2007; Narteau et al., 2009; Peng et al., 2006). Typically observed values of c are of the order of magnitude of several seconds to several minutes (Page et al., 2016). The operational provision of earthquake forecasts has practical value only if the earthquake probabilities provided are time‐variant. And the variability with time of most models currently applied for OEF in Italy, New Zealand, and the United States is centered around the Omori‐Utsu law. The EEPAS model (Every Earthquake a Precursor According to Scale; Rhoades & Evison, 2004; described later in this section) stands as an exception, although it, too, uses prior earthquakes as the sole precursory signals for impending earthquakes. Thus, one limitation common to all currently authoritatively applied forecasting models is their inability to predict (or forecast with probabilities considered high by decision‐makers) earthquakes. Similarly common to these models is their main strength of capturing the occurrence patterns of aftershocks following a moderate to large event. The possibility that an aftershock has a larger magnitude than the current mainshock can be captured in most short‐term forecasting models through another central empirical law of statistical seismology—the Gutenberg‐ Richter relationship (Gutenberg & Richter, 1944). It describes the relative frequency of earthquakes in relation to their magnitude as N(M) = 10a− b ·M, (2) Reviews of Geophysics 10.1029/2023RG000823 MIZRAHI ET AL. 4 of 70 19449208, 2024, 3, D ow nloaded from https://agupubs.onlinelibrary.w iley.com /doi/10.1029/2023R G 000823 by M assey U niversity L ibrary, W iley O nline L ibrary on [23/10/2024]. See the T erm s and C onditions (https://onlinelibrary.w iley.com /term s-and-conditions) on W iley O nline L ibrary for rules of use; O A articles are governed by the applicable C reative C om m ons L icense https://agupubs.onlinelibrary.wiley.com/action/rightsLink?doi=10.1029%2F2023RG000823&mode= where N(M) is the expected number of events with magnitudes greater than or equal to M, and a and b are pa- rameters generally referred to as a‐value and b‐value. In this formulation, the a‐value parametrizes the number of earthquakes with magnitude 0 or larger in the region and time horizon of interest. Note that using 0 as a threshold value is somewhat arbitrary and does not represent a lower limit on the magnitude. The a‐value naturally tends to vary in space; places close to boundaries of tectonic plates for example, tend to have more earthquakes and thus higher a‐values than places in the middle of tectonic plates. Spatial variations of the a‐value are merely a more abstract formulation of Kant's (1756) observations mentioned at the beginning of this section: the expected number of earthquakes may differ from place to place. The b‐value parametrizes the relative frequency of large magnitude events compared to small magnitude ones. In general, the b‐value for moment magnitudes is found to be around 1 (Kagan, 1999; Kamer & Hiemer, 2015; Utsu, 1972) for regions that are considered large compared to the seismic faults they contain, which implies that 10 times as many M ≥ 5 earthquakes are expected as M ≥ 6 earthquakes, 10 times as many M ≥ 4 earthquakes as M ≥ 5 earthquakes, and so on. Variations of the b‐value in time and space have been widely discussed and have been proposed to be related to stress, volcanic activity, rock heterogeneity, crack density, and other physical quantities (Main et al., 1992; Mori & Abercrombie, 1997; Murru et al., 2007; Scholz, 2015; Schorlemmer et al., 2004, 2005; Wiemer & Wyss, 1997, 2002). Combining the expected decay of the number of aftershocks given by the Omori‐Utsu law with the knowledge about the relative frequency of small and large magnitude events given by the Gutenberg‐Richter relationship enables an estimation of the expected probability that an aftershock larger than the current mainshock (or larger than any magnitude of interest) is yet to follow. In this sense, such probabilistic models provide earthquake probabilities, rather than aftershock probabilities. A crucial concept for the type of analysis of earthquake occurrence that has been so far described is that of an earthquake catalog, which is a list of earthquakes, where each earthquake is characterized through several at- tributes related to its source. Most catalogs contain each earthquake's epicenter (the point on the earth's surface above the point where the earthquake started), depth, origin time (the time when the earthquake started), and magnitude. Additional information such as details about the rupture (orientation of the fault plane, direction of slip, etc.), and other attributes or attribute uncertainties may be available too. An earthquake catalog is usually compiled for a specific region or area based on the recordings of seismic waves captured by a network of seismometers covering that region or area. Earthquake catalogs can be heterogeneous for various reasons. The seismic network only records signals at specific points, hence, naturally, not all locations are observed equally well. Some subregions may be more densely covered by the network than others. The quality of catalogs can vary with time, for example, because more instruments are deployed over time, or because newer earthquake detection methods become available (e.g., Hutton et al., 2010). In the shorter term, the ability to detect earthquakes can fluctuate after the occurrence of a relatively large earthquake, when the waves of smaller earthquakes are hidden within the large amplitude waves of the bigger one (Kagan, 2004). When discussing the issue of varying detection capability, the concept of the completeness or incompleteness of a catalog is often used. In general, large earthquakes are more easily detectable than small ones. The magnitude of completeness (also: completeness magnitude), which we will here denote with Mc, is the magnitude above which all earthquakes, in the region and time span of interest, are assumed to be detected and thus to be present in a catalog. Because it is impossible to know about the earthquakes that have not been detected, because they have not been detected, Mc has to be estimated and various methods have been proposed to do this (e.g., Schorlemmer & Woessner, 2008; Wiemer & Wyss, 2000). The magnitude itself is a quantity that deserves discussion. An earthquake catalog is a complex entity because the quantities it contains are inferred from seismograms through uncertain and non‐unique analyses. For the purposes of earthquake forecasting, the magnitude deserves particular attention because it does not have a clear physical meaning and forecasts are particularly sensitive to it. In contrast to the intensity of shaking, which varies between different locations based on their distance to the earthquake and other relevant factors, the magnitude is meant to be a unique quantity describing the size of the earthquake, independently of where it was measured. Several different definitions of magnitudes (i.e., magnitude scales) have been developed in the literature, and the way magnitudes are calculated varies strongly between the different scales, but also between different study regions. Two of the most widely used magnitude scales are local magnitudes (ML) and moment magnitudes (Mw). Local magnitudes are defined based on the maximum amplitudes of recorded transient displacement at a seismic station, and the effects of the location at which the recording was made on the recording itself are corrected for, to obtain a Reviews of Geophysics 10.1029/2023RG000823 MIZRAHI ET AL. 5 of 70 19449208, 2024, 3, D ow nloaded from https://agupubs.onlinelibrary.w iley.com /doi/10.1029/2023R G 000823 by M assey U niversity L ibrary, W iley O nline L ibrary on [23/10/2024]. See the T erm s and C onditions (https://onlinelibrary.w iley.com /term s-and-conditions) on W iley O nline L ibrary for rules of use; O A articles are governed by the applicable C reative C om m ons L icense https://agupubs.onlinelibrary.wiley.com/action/rightsLink?doi=10.1029%2F2023RG000823&mode= unique value which ideally does not depend on the location of the recording. This magnitude scale was originally introduced by Richter (1935) and is often called the Richter magnitude scale. Moment magnitudes (Hanks & Kanamori, 1979) are more tightly linked to earthquake source features: they are based on the energy (the seismic moment) released by the earthquake. Obtaining the moment released by an earthquake requires more complex processing of the recorded data. It is usually only done for a subset of earthquakes for which local magnitudes can be calculated. Both local and moment magnitudes, as well as other magnitude types not discussed here, scale logarithmically with the physical quantity they intend to measure. Although different magnitude types are usually defined such that the values of these different magnitude types of individual earthquakes coincide as much as possible, an exact match is not generally possible and the presence of different magnitude types within a catalog is another source of catalog heterogeneity which may need to be addressed by analysts. For example, using local or moment magnitudes for the same earthquake catalog can result in substantially different b‐value estimates (Deichmann, 2017) and some magnitude types adhere less closely to the Gutenberg‐Richter relationship than others (Herrmann & Marzocchi, 2021). In the following, we will refer to magnitudes mostly without specifying magnitude type, assuming what is discussed holds for different magnitude types. The macroseismic intensity scales, such as the Modified Mercalli Intensity scale (MMI, H. O. Wood & Neu- mann, 1931) on the other hand, are twelve‐level scales of Roman (ordinal) numerals which qualitatively indicate an increasing level of shaking based on observed effects, in vast areas, of an earthquake on people, animals, objects, infrastructure and landscape. Ultimately, these effects are what one is interested in when assessing the risk posed by earthquakes. While an earthquake ideally has a unique magnitude (ignoring the just discussed different possible magnitude scales), it can cause shaking of different MMI levels at different locations. However, because of the qualitative nature of macroseismic intensities, and the fact that they do not relate to specific sites of interest, quantitative risk assessment preferentially expresses earthquake hazard in terms of ground motion in- tensity measures, which are coupled with vulnerability and consequence models to compute (probabilistic) seismic risk metrics. 2.1.2. Description of Model Types Used for OEF 2.1.2.1. The Reasenberg and Jones Model Combining the Omori‐Utsu law with the Gutenberg‐Richter relationship yields the arguably simplest model currently used for OEF: the Reasenberg and Jones (R&J; 1989) model. It defines the rate λ(t,M) of earthquakes of magnitude M or above at time t after a parent event of magnitude Mm as λ(t, M) = 10a− b(M− Mm ) · (t + c)− p. (3) Here, the second factor is the denominator of the Omori‐Utsu law, and the first factor describes the expected number of aftershocks above a certain magnitude according to the Gutenberg‐Richter relationship. Note that the parameter a here has a different interpretation than in Equation 2; here it controls the number of aftershocks larger than the parent event. The R&J model is thus based on the assumption that the number of aftershocks larger than their parent event is independent of the size of the parent event. A b‐value of 1 would imply that, counting all aftershocks above a fixed magnitudeM, a magnitude 6 event has (on average) 10 times as many aftershocks as a magnitude 5 event. Using the same scaling for aftershock productivity and the Gutenberg‐Richter relationship produces a self‐similar process in which earthquakes in equal‐size magnitude bins produce the same total number of aftershocks. It also implies that the aftershock and foreshock processes are similar and share the same scaling (e.g., Felzer et al., 2004). Reading the previous paragraph, one may have noticed that it mentions parent events rather than mainshocks, in contrast to the originally used terminology of Reasenberg and Jones (1989). As noted earlier, this is to emphasize that a parent event can trigger aftershocks of any magnitude. In fact, if a large aftershock occurs, it is expected to trigger aftershocks of its own. It is a main limitation of the R&J model that it does not specifically model such higher‐order aftershocks (i.e., aftershocks of aftershocks, and so on), but considers all aftershocks as triggered by one event. This simplification allows the R&J model to be analytic and computationally efficient, advantages that were particularly important when it was introduced. Reviews of Geophysics 10.1029/2023RG000823 MIZRAHI ET AL. 6 of 70 19449208, 2024, 3, D ow nloaded from https://agupubs.onlinelibrary.w iley.com /doi/10.1029/2023R G 000823 by M assey U niversity L ibrary, W iley O nline L ibrary on [23/10/2024]. See the T erm s and C onditions (https://onlinelibrary.w iley.com /term s-and-conditions) on W iley O nline L ibrary for rules of use; O A articles are governed by the applicable C reative C om m ons L icense https://agupubs.onlinelibrary.wiley.com/action/rightsLink?doi=10.1029%2F2023RG000823&mode= 2.1.2.2. The Epidemic‐Type Aftershock Sequence Model Around the same time as when the R&J model was introduced, Ogata (1988) first described what is now widely known as the Epidemic‐Type Aftershock Sequence (ETAS) model, and it was first used as a forecasting tool by Console et al. (2003). As suggested by its name, it models seismicity to resemble an epidemic, where earthquakes recursively trigger more earthquakes, similarly to how infected beings can recursively infect others around them. Numerous variants of ETAS formulations have been and keep being proposed (Asayesh et al., 2023; Console et al., 2003; Guo et al., 2015; Harte, 2013; Helmstetter & Sornette, 2002; Nandan et al., 2017; Ogata, 1993, 1998, 2011; Ogata & Zhuang, 2006; Veen & Schoenberg, 2008; Z. Xiong & Zhuang, 2023). Generally, ETAS describes the rate λ(t,x,y) of earthquakes with magnitudes above a reference magnitude Mr at time t and location (x,y) as λ(t, x, y) = μ +∑ ti < t k(mi) · g(t − ti) · f (mi, x − xi, y − yi). (4) The key idea of the model is that earthquakes are partitioned into background earthquakes and triggered earth- quakes. All earthquakes not triggered by a previous one are considered background and are assumed to occur uniformly in time. The rate of background earthquakes with magnitudes aboveMr is captured by the parameter μ. In addition to the background earthquakes come the aftershocks with magnitudes above Mr of all prior earth- quakes. The summation in Equation 4 is over all earthquakes which occurred at times ti before time t and at locations (xi,yi) and had magnitudes mi. Each prior earthquake contributes to the rate of triggered events, and the rate of its aftershocks is defined through three components. The total number of expected aftershocks or after- shock productivity k(mi) depends on the magnitude mi of the parent earthquake; the temporal decrease of the aftershock rate is modeled through g(t − ti) and only depends on the elapsed time since the parent event. The spatial distribution of aftershocks, f(mi,x − xi,y − yi), can depend on the magnitude of the parent event and is often defined as isotropically decreasing with the distance to the parent event (Musmeci & Vere‐Jones, 1992; Ogata, 1998). Figure 1 schematically illustrates the concept of aftershock triggering within the ETAS model. The temporal aftershock decay is most often modeled according to the Omori‐Utsu law given in Equation 1. The aftershock productivity k(mi) is often described to depend on the magnitude of the parent event through an exponential relationship of the form k(m) = K0 · 10α · (m− Mr), (5) (Utsu, 1971). The parameter α parametrizes the increase of the aftershock number per magnitude unit of the parent event. If α equaled the b‐value of the Gutenberg‐Richter relationship, as is often suggested, observed, and imposed (Hainzl et al., 2008, 2013; van der Elst et al., 2022), the assumption made in the R&J model that the number of aftershocks larger than their parent event is independent of the parent event magnitude would be fulfilled. A variety of spatial aftershock decay kernels has been proposed and will not be further discussed here; instead, we refer to Zhuang et al. (2011) for a detailed technical overview of spatio‐temporal ETAS models. The description of the ETAS model given here is a simplistic one. Variants in which the background seismicity rate μ, but also other model parameters, vary in space and/or time are common. Generally, the estimation of ETAS parameters has been studied thoroughly, yielding improved parameter estimation techniques (e.g., Lippiello et al., 2014; Lombardi, 2015; Ross, 2021; Schneider & Guttorp, 2020; Veen & Schoenberg, 2008) and aiming to understand and address biases in parameter estimates due to data imperfections (e.g., Grimm et al., 2022; Harte, 2013; Seif et al., 2017). The important quality of the ETAS model that it accounts for higher‐order aftershocks is also the cause of a noteworthy drawback: even when the parameters describing background seismicity and aftershock triggering are known, the expected rate of future earthquakes cannot be calculated analytically. To obtain the expected rate of events at time t, the summation in Equation 4 is over all events that occurred prior to t. In practice, one may be interested in the expected rate of events hours, days, weeks or even years into the future. Thus, aftershocks of events that have not yet occurred when the forecast is calculated are not accounted for in Equation 4. Modelers address this by issuing many simulated stochastic event sets, containing simulated earthquakes and cascades of Reviews of Geophysics 10.1029/2023RG000823 MIZRAHI ET AL. 7 of 70 19449208, 2024, 3, D ow nloaded from https://agupubs.onlinelibrary.w iley.com /doi/10.1029/2023R G 000823 by M assey U niversity L ibrary, W iley O nline L ibrary on [23/10/2024]. See the T erm s and C onditions (https://onlinelibrary.w iley.com /term s-and-conditions) on W iley O nline L ibrary for rules of use; O A articles are governed by the applicable C reative C om m ons L icense https://agupubs.onlinelibrary.wiley.com/action/rightsLink?doi=10.1029%2F2023RG000823&mode= aftershocks, which yields more correct forecasts (Helmstetter & Sornette, 2003) but requires more computation time. 2.1.2.3. The Epidemic‐Type Earthquake Sequence Model Although some authors use the name Epidemic‐Type Earthquake Sequence (ETES) model as a synonym for the ETAS model (e.g., Helmstetter et al., 2006), we maintain its original use as a special case of the ETAS model based on the Omori‐Utsu law for the temporal aftershock decay, and a circular symmetry around the triggering event's epicenter for the spatial component of the rate (Console et al., 2003; Falcone et al., 2010). The main distinction between ETES and ETAS as defined in Equation 4 can be summarized in the following points. 1. The background component of the rate is still time‐invariant, but it depends on a failure rate, defined as the ratio between the expected number of background events and the total number of events. 2. In ETES, the expected spatial density of triggered events does not depend on the magnitude of the triggering event and hence is of the form f(x − xi,y − yi). 3. The parameter α in Equation 5 is set to 1.0. Physically, this means that the number of triggered events scales with the magnitude of the parent event in a similar way as its rupture area (Console et al., 2006; Hainzl et al., 2008; Murru et al., 2014). Note that assuming a b‐value of 1.0, setting α = 1.0 is equivalent to setting α equal to b, which was discussed earlier. The ETES model with α = 1 was applied to the aftershock sequence following the 2011 Mw 9.0 Tohoku‐Oki, Japan, earthquake and achieved a better performance at forecasting event numbers compared to other models submitted in the same test (Nanjo et al., 2012). 2.1.2.4. The Short‐Term Earthquake Probability Model The Short‐Term Earthquake Probability (STEP) model (Gerstenberger et al., 2004; Woessner et al., 2010) is based on the same core principles as the R&J and ETASmodels. It uses the R&Jmodel to describe aftershock rate, Figure 1. Schematic illustration of the ETASmodel. (a) Background earthquakes (turquoise) and aftershock cascades (black). Arrows indicate triggering. (b) Composition of the earthquake rate into a temporally invariant background rate and aftershock rate. The aftershock rate itself has three independent components: it decays with time and distance from the parent event, and increases with the magnitude of the parent event. All earthquakes (background events and triggered events) follow the same magnitude distribution. Reviews of Geophysics 10.1029/2023RG000823 MIZRAHI ET AL. 8 of 70 19449208, 2024, 3, D ow nloaded from https://agupubs.onlinelibrary.w iley.com /doi/10.1029/2023R G 000823 by M assey U niversity L ibrary, W iley O nline L ibrary on [23/10/2024]. See the T erm s and C onditions (https://onlinelibrary.w iley.com /term s-and-conditions) on W iley O nline L ibrary for rules of use; O A articles are governed by the applicable C reative C om m ons L icense https://agupubs.onlinelibrary.wiley.com/action/rightsLink?doi=10.1029%2F2023RG000823&mode= but differs from it in that it does, in addition, include a component that models the long‐term background seis- micity, does account for aftershocks of aftershocks, and does include a description of seismicity in space. It thus shares its main advantages over the R&J model with the ETAS model. The way in which these issues are addressed however differs between ETAS and STEP. While in ETAS, the background and the recursive space‐ time aftershock triggering are all self‐consistently described within a single model, the STEP model is a com- bination of several components that could themselves be viewed as individual models. Besides a time‐invariant component, a STEP model can have several time‐variant components: a generic one based on past seismicity, and sequence‐specific ones which can be derived in different ways. In these time‐variant components, recorded earthquakes are classified to be either primary events or secondary events. An event is considered secondary if it is an aftershock of a specific prior event, which can but does not need to be a primary event. Aftershocks are associated with their parent events through aftershock zones defined in space and time, and they can only be considered aftershocks of events of equal or larger magnitude. As in ETAS, all earthquakes are considered capable of triggering aftershocks. A distinction is made however between primary and secondary events, where a sequence‐specific set of parameters can describe the expected aftershock rate of primary events. In contrast, secondary events always use generic parameters. Within one component, the forecasted rate of events at a given location is given as the maximum of forecasted aftershocks at that location of all previous primary or secondary events. The different time‐dependent components are combined by weighting them according to how well they fit the observed data, and the combined rate forecast is added to the time‐independent rate forecast to obtain the final STEP forecast. An advantage of the STEP model is that, similar to the R&J model, it is analytic and computationally efficient, while partially including secondary triggering (unlike the R&J model). Its process of combining time‐dependent components with a long‐term hazard model allows it to easily integrate, and be consistent with, authoritative national seismic hazard models. STEP can accommodate different spatial kernels and can either fit a two‐segment fault model or take external finite fault input. A disadvantage of STEP, compared to ETAS, is that it does not fully model secondary triggering. 2.1.2.5. The Every Earthquake a Precursor According to Scale Model The Every Earthquake a Precursor According to Scale (EEPAS) model was introduced by Evison and Rhoades (2004) and Rhoades and Evison (2004). Although it shares the key principle of the ETAS model of describing earthquake rate as the sum of a background term and the sum of rate contributions of prior events (see Equation 4), the specific way in which prior events contribute to the forecasted rate differs substantially from the ETAS formulation. There is no suggestion that each earthquake may trigger its own offspring. Rather, each earthquake is regarded as evidence of a medium‐term seismogenic process taking place on the scale indicated by its magnitude (Rhoades, 2007). The probabilistic distribution of an earthquake's contribution is lognormal in time, isotropic bivariate normal in space and normal in magnitude (Rhoades & Evison, 2004). The EEPAS model is underpinned by the observation that prior to most large earthquakes the rate and magnitude of smaller earthquakes increases in the source region of the large earthquake. This precursory scale increase (Ψ‐ phenomenon) takes place within a precursor time TP and precursory area AP. Evison and Rhoades (2004) showed that TP, AP, and the mainshock magnitudeMm all scale with the precursor magnitudeMP (defined as the average of the three largest precursory earthquake magnitudes). Thus,MP is predictive of the time, magnitude and location of forthcoming mainshocks. The assumption behind EEPAS is that the three Ψ predictive relations are pervasive at all scales in the seismogenic process. The model regards every earthquake as a precursor, according to scale, of larger earthquakes expected to follow it in the coming months, years, or decades, depending on its magnitude and the regional earthquake rate. 2.1.2.6. Ensemble Models Rather than relying on a single model to provide a forecast, one can combine multiple models into one ensemble model. Ensemble modeling proved beneficial in various fields and applications such as weather/climate fore- casting (Hagedorn et al., 2005; Krishnamurti et al., 2000; Tebaldi & Knutti, 2007), flood forecasting (Cloke & Pappenberger, 2009), ecology (Buisson et al., 2010), computer security (Menahem et al., 2009), and long‐term Reviews of Geophysics 10.1029/2023RG000823 MIZRAHI ET AL. 9 of 70 19449208, 2024, 3, D ow nloaded from https://agupubs.onlinelibrary.w iley.com /doi/10.1029/2023R G 000823 by M assey U niversity L ibrary, W iley O nline L ibrary on [23/10/2024]. See the T erm s and C onditions (https://onlinelibrary.w iley.com /term s-and-conditions) on W iley O nline L ibrary for rules of use; O A articles are governed by the applicable C reative C om m ons L icense https://agupubs.onlinelibrary.wiley.com/action/rightsLink?doi=10.1029%2F2023RG000823&mode= seismic hazard assessment (Gerstenberger et al., 2020). Ensembles can emphasize the individual strengths of different types of models (e.g., physical, statistical, deterministic), aim to achieve better performance than the best single model, are more flexible than single models, and better capture (different kinds of) uncertainty. In earthquake forecasting research, several ensembles have been created using various methods. On short‐term (time‐varying) scales, weighted averaging of rate‐based statistical models is common (Helmstetter & Werner, 2014; Herrmann & Marzocchi, 2023; Llenos & Michael, 2019; Marzocchi, Zechar, & Jordan, 2012; Rhoades & Gerstenberger, 2009; Taroni et al., 2018), and some ensembles incorporate physical information (Király‐Proag et al., 2018; Steacy et al., 2013) or make an existing model more flexible (Gerstenberger et al., 2004; Mizrahi, Nandan, et al., 2023). These studies have shown promising results, even though ensembles have not always performed much better than the best individual model. On long‐term (time‐invariant) scales, a multiplicative procedure has been used to incorporate models and data other than earthquake rates, like strain and fault slip rates (Bird et al., 2015; Rhoades et al., 2014, 2017). Prospective testing of some of those has shown mixed results: although they can outperform their candidate models on a global scale (Strader et al., 2018), they may not do so on a regional scale (Bayona et al., 2022). Despite a wide range of possible ensemble methods, only weighted averaging allows for considering the dispersion of the single forecasts to inform about the epistemic uncertainty (Marzocchi & Jordan, 2014, 2017), that is, our lack of knowledge. A reasonable estimation of the epistemic uncertainty should be based on weights that account for the correlation of the forecasts (Marzocchi, Zechar, & Jordan, 2012) or maximize the forecasting skill of the ensemble (Herrmann & Marzocchi, 2023). Even if some ensembles have not greatly outperformed the best‐performing single model, they do not usually perform much worse. It remains to be seen whether the relatively limited skill improvements are due to sub- optimal ensemble methods, the limited skill of the single models, or because models did not complement each other well. A possible path for more skilled ensembles is to incorporate diverse models (e.g., physics‐based approaches, for short‐term earthquake forecasting like those proposed by Mancini et al. (2019), Cattania et al. (2018), Sharma et al. (2020); or Dahm and Hainzl (2022)), because diversity increases the predictive pool and enables exploiting well‐performing forecasting models. 2.2. Background on Earthquake Forecast Model Testing One decision‐making process in OEF involves the selection of one or more suitable forecasting models. While many factors might affect the selection process (including less or non‐scientific aspects such as required computational resources in real‐time, ability to explain the model to non‐experts, internal expertise to develop and maintain the model, maintaining some continuity with past model forecasts, etc., see also Marzocchi & Zechar, 2011), a key scientific question is howwell a model can forecast future seismicity, especially in relation to other available models. OEF itself is motivated by scientific studies concluding that some time‐dependent models forecast earthquake clustering better than time‐independent models (e.g., Cattania et al., 2018; Helmstetter et al., 2006; Werner et al., 2011), most obviously during aftershock sequences. Thus, a major community effort has focussed on the performance evaluation of available earthquake and aftershock forecasting models, including OEF candidate models. This subsection reviews some of the basic concepts and methods that the community has used to assess model performance. We will focus primarily on the major community effort organized by the global Collaboratory for the Study of Earthquake Predictability (CSEP; Jordan, 2006; Field, 2007; Michael & Werner, 2018; Schorlemmer et al., 2018; Savran et al., 2022). 2.2.1. Basic Concepts 2.2.1.1. Testing Model inference, evaluation, and selection are all well‐established and still evolving fields, with many approaches and applications from across many disciplines. The philosophical foundations of (probabilistic) model evaluation approaches continue to be debated (Harte & Vere‐Jones, 2005; Serafini et al., 2022), perhaps as a result of the different interpretations of probability (Marzocchi & Jordan, 2014). The seismological community has thus far largely focused on classical frequentist‐type testing of probabilistic forecasts, some of which are described below (see Basic Methods). A recognition of the shortcomings associated with frequentist testing and likelihood‐based model comparisons (e.g., Field, 2015; Nandan et al., 2019; Savran et al., 2022; Werner & Sornette, 2008) has recently led to more modern approaches, including simulation‐based inference approaches (e.g., Page & van der Elst, 2018; Savran et al., 2020, 2022), Bayesian (e.g., Marzocchi, Zechar, & Jordan, 2012) and unified approaches Reviews of Geophysics 10.1029/2023RG000823 MIZRAHI ET AL. 10 of 70 19449208, 2024, 3, D ow nloaded from https://agupubs.onlinelibrary.w iley.com /doi/10.1029/2023R G 000823 by M assey U niversity L ibrary, W iley O nline L ibrary on [23/10/2024]. See the T erm s and C onditions (https://onlinelibrary.w iley.com /term s-and-conditions) on W iley O nline L ibrary for rules of use; O A articles are governed by the applicable C reative C om m ons L icense https://agupubs.onlinelibrary.wiley.com/action/rightsLink?doi=10.1029%2F2023RG000823&mode= (Marzocchi & Jordan, 2014, 2017, 2018), the latter two of which are beyond our scope here but provide compelling frameworks. In the present context, we can define testing a forecasting model as formal procedures to identify discrepancies between forecasts and data, and interpreting these discrepancies in terms of their scientific and/or practical implications. Typically, models are also compared against each other in comparison tests. 2.2.1.2. Testability A study of the “tumultuous history of earthquake prediction” (Hough, 2010) shows that many past claims of successful earthquake predictions and forecasts based on models were not reproducible (see, e.g., National Academies of Sciences, Engineering, and Medicine, 2019), including those in peer‐reviewed publications. A range of reasons can explain this, but many can be summarized by a lack of testability of a model (e.g., Jack- son, 1996). To be testable, there must be (a) a clearly defined region and depth over which the forecast or pre- diction applies, a magnitude range, and a time window, (b) a clear definition of which characteristic of an earthquake is being forecast (e.g., epicenter or centroid, moment magnitude or surface wave magnitude) and (c) a pre‐specified authoritative, independent source of target/test data (e.g., Jackson & Kagan, 1999). 2.2.1.3. Authoritative, Independent Test Data Sources An ambiguity in testing involves the choice of the target or test data set, at times for the same earthquake. Earthquake catalogs are data products derived from seismic waveform observations, quantitative models of the Earth, fixed or adjusted parameters, and sometimes ad hoc human intervention. The same earthquake might be assigned different parameters by different networks, or even multiple Psource parameter estimates (e.g., magnitude) by the same network. To avoid this ambiguity and any potential bias due to a post‐hoc selection, the source and type of test data should be specified in advance, and originate from authoritative, independent net- works and agencies. 2.2.1.4. Prospective Testing The earthquake forecasting community places great emphasis on prospective testing, which comprises the evaluation of forecasts against future, yet‐to‐be‐collected seismic data. For example, a forecast created and archived in real‐time can be evaluated prospectively in the future against the earthquake data collected between when the forecast was issued and the end of the testing window (also called the forecast horizon). Prospective testing can also be conducted in a delayed manner, that is, with forecasts not created in real‐time, when the model is entirely isolated from any data beyond the forecast issue time, and the test data set is fully specified before the forecast issue time. This requires fully automated and transparent procedures to ensure a zero‐degrees‐of‐freedom environment. In its first phase, CSEP built and managed so‐called testing centers (Schorlemmer & Ger- stenberger, 2007; Zechar, Gerstenberger, & Rhoades, 2010; Zechar, Schorlemmer, et al., 2010) to implement delayed prospective testing by providing a controlled, secure and automated environment to test models with at least a 1‐month delay. In this way, delayed prospective testing can remain entirely out‐of‐sample (i.e., model fitting and testing data are strictly separate) and, importantly, testing is against future data that neither model nor modeler have seen. Prospective testing provides the most objective view of the forecast skill of a model due to excluding even unconscious biases of modelers. It is often assumed to provide the most reliable estimate of future performance, although it is not yet clear for how long models must be tested to obtain stable and robust per- formance estimates. This is however not a specific problem of prospective testing, but one of testing in general. 2.2.1.5. Retrospective Testing In contrast to prospective testing, retrospective testing is the testing of forecasts against past data. In retrospective testing, the test data were available to the modelers and thus consciously or unconsciously influenced modeling and/or testing choices (e.g., Schorlemmer et al., 2018). Test data may or may not have been explicitly used to fit or calibrate a model. To emphasize when test data were not known by the model, and when time‐dependent causality was preserved when separating the test data from the data used for model fitting, the term pseudo‐prospective testing, described below, is often used. Figure 2 visualizes the different testing modes in terms of the data split between test data, data known to the model, and data known to the modelers. An example of retrospective testing for which time‐dependent causality is not preserved is shown at the bottom of the figure, where the test data are not known by the model, but cover a period before the period used for model calibration. Reviews of Geophysics 10.1029/2023RG000823 MIZRAHI ET AL. 11 of 70 19449208, 2024, 3, D ow nloaded from https://agupubs.onlinelibrary.w iley.com /doi/10.1029/2023R G 000823 by M assey U niversity L ibrary, W iley O nline L ibrary on [23/10/2024]. See the T erm s and C onditions (https://onlinelibrary.w iley.com /term s-and-conditions) on W iley O nline L ibrary for rules of use; O A articles are governed by the applicable C reative C om m ons L icense https://agupubs.onlinelibrary.wiley.com/action/rightsLink?doi=10.1029%2F2023RG000823&mode= Because truly prospective testing is often associated with long waiting times to obtain the yet‐to‐be‐collected data against which forecasts are evaluated, most scientific publications report model evaluations conducted retro- spectively. The main drawback of retrospective testing is that modeling choices or testing choices (e.g., the choice of a particular data set for testing) may affect model performance, even when testing data were not explicitly used by the model. The extent of the resultant bias in performance is difficult to assess. However, a rule of thumb is that retrospective performance will provide an upper bound to truly prospective performance. Nevertheless, retro- spective tests provide added value by acting as sanity checks, ensuring that the models function as intended by the modelers. Figure 2. Different modes of testing, visualized in terms of the data split between test data, data known to the model, and data known to the modelers. Reviews of Geophysics 10.1029/2023RG000823 MIZRAHI ET AL. 12 of 70 19449208, 2024, 3, D ow nloaded from https://agupubs.onlinelibrary.w iley.com /doi/10.1029/2023R G 000823 by M assey U niversity L ibrary, W iley O nline L ibrary on [23/10/2024]. See the T erm s and C onditions (https://onlinelibrary.w iley.com /term s-and-conditions) on W iley O nline L ibrary for rules of use; O A articles are governed by the applicable C reative C om m ons L icense https://agupubs.onlinelibrary.wiley.com/action/rightsLink?doi=10.1029%2F2023RG000823&mode= 2.2.1.6. Pseudo‐Prospective Testing The subset of retrospective testing called pseudo‐prospective testing aims to mimic the circumstances of pro- spective testing and preserve time‐dependent causality. Models are being calibrated using data only until a certain time t0 in the past, and are then used to produce forecasts for the time after t0. In this way, models and modelers pretend not to know what occurs after time t0, but the modelers might be biased due to their knowledge of the testing data or certain features thereof. Thus, while pseudo‐prospective tests can formally be considered out‐of‐ sample, modeling and/or testing choices may still be consciously or unconsciously influenced by the testing data, rendering pseudo‐prospective tests less informative of a model's performance than truly prospective tests. Thus, while pseudo‐prospective tests can formally be considered out‐of‐sample, modeling and/or testing choices may still be consciously or unconsciously influenced by the testing data, rendering pseudo‐prospective tests less informative of a model's prospective performance than truly prospective tests. 2.2.1.7. Partially Prospective Testing In the literature, a testing approach that neither qualifies as prospective, nor can be considered fully retrospective and is thus not captured by the established terminology, can often be found. This approach, which we define as partially prospective testing, involves testing forecasts against partially future, partially past data. For instance, researchers might generate a forecast of aftershocks a few hours or days after a mainshock, for example, after quality‐controlled information about the mainshock source parameters have become available, such as magni- tude, faulting mechanism or the hosting fault, or researchers may have calibrated their aftershock model to a new region that experienced the mainshock (see e.g., Mancini et al., 2020; Milner et al., 2020, for forecasts issued during the 2019 Ridgecrest, California, sequence). During this time, aftershocks will already have occurred and been recorded, and indeed may influence the quality control of mainshock parameters, such as its rupture planes and spatial extent. The forecast horizon, however, might start in the past (at the time of the mainshock) and extend into the future, so as to be useful in some form. In addition, these forecasts might eventually be evaluated and published in a peer‐reviewed publication, but no time‐stamped archived record of the forecast might exist. As illustrated in Figure 2, we make an overall distinction between prospective and retrospective testing, where prospective testing can be done in real time or in a delayed manner, and retrospective testing can be done pseudo‐ prospectively (with a data split that does preserve time‐causality), or out‐of‐sample with a data split that does not preserve time causality, or in‐sample. In‐between these two categories lies partially prospective testing. 2.2.2. Basic Methods 2.2.2.1. Forecast Specification During the first CSEP experiment, the Regional Earthquake Likelihood Models (RELM) experiment of 5‐year forecast models in California (Field, 2007; Schorlemmer et al., 2007), modelers decided to define a forecast Λ as the expected number of events (mean rates) in a space‐magnitude discretization (cells/bins of a grid), such that Λ = {λi} for i = 1, ….,I bins. These long‐term forecasts were time‐independent, both in the sense of their mathematical formulations (statistical descriptions were not explicitly dependent on time) and the decision that models would not be updated with new data during the forecast time horizon. This implied that forecast models could be issued for testing as flat files of numerical data (Schorlemmer & Gerstenberger, 2007). Following this scheme, short‐term 1‐day RELM forecasts (relevant for OEF) were also specified as mean rates on a space‐ magnitude grid. However, it was required that the model source code be available for prospective testing, so the forecasts could be updated with new data after the 1‐day time window, while preventing the modelers from modifying the original model formulation. In practice, the first generation of short‐term forecasts handled by CSEP was tested in a similar fashion and with similar metrics as the long‐term forecasts, with the only difference of providing a temporal evolution of the testing results (e.g., Cattania et al., 2018; Taroni et al., 2018). 2.2.2.2. Likelihood‐Based Model Scoring, Comparison, and Ranking The probabilistic nature of forecasts allows them to be evaluated using likelihood functions. That is, the joint probability L(Ω,Λ) of the observations Ω (i.e., the testing catalog), interpreted as a function of the issued prob- abilistic forecast Λ, measures how well a forecast explains the observed data. This concept is illustrated in Figure 3. Forecast scoring usually employs the log‐likelihood score (or log‐score) LL(Ω,Λ) = lnL(Ω,Λ), with Reviews of Geophysics 10.1029/2023RG000823 MIZRAHI ET AL. 13 of 70 19449208, 2024, 3, D ow nloaded from https://agupubs.onlinelibrary.w iley.com /doi/10.1029/2023R G 000823 by M assey U niversity L ibrary, W iley O nline L ibrary on [23/10/2024]. See the T erm s and C onditions (https://onlinelibrary.w iley.com /term s-and-conditions) on W iley O nline L ibrary for rules of use; O A articles are governed by the applicable C reative C om m ons L icense https://agupubs.onlinelibrary.wiley.com/action/rightsLink?doi=10.1029%2F2023RG000823&mode= higher values implying a better performance. The original RELM testing formulation involved evaluating the log‐ likelihood function assuming a Poisson process, for multiple reasons: (a) it yields mathematical tractability for space‐magnitude discretizations due to the independence of the distributions of event counts in the resulting bins/ cells, (b) it can be derived from only one temporal parameter (the mean rate), and (c) its simple formulation enables confronting various types of models, even models that have no explicit likelihood functions. The multivariate Poisson log‐likelihood is given as LL(Ω, Λ) =∑ I i=1 ln ( λi) ωi ωi! e− λi , (6) where ωi is the number of events observed in the ith space‐magnitude bin. The log‐likelihood scores of different models can be compared to one another: a model with higher LL(Ω,Λ) than another provided a better forecast of the observationsΩ. A ranking of the models' performance can then be built by ordering the forecasts according to their scores. Differences in log‐likelihood scores of different models are often referred to as information gain (IG), or, when divided by the number of earthquakes contained in the testing catalog, as information gain per earthquake (IGPE). Hereby, a positive IG (or IGPE) suggests a better performance of the alternative model, while negative IG suggests poorer performance and an IG of 0 suggests equal performance in terms of the log‐likelihood score. 2.2.2.3. Frequentist (Consistency) Testing Under the frequentist approach, predictive statistical distributions of a forecast model can be submitted to standard hypothesis testing, which informs us how (in)consistent a model is with new data. In traditional CSEP tests (Kagan & Jackson, 1995; Schorlemmer et al., 2007; Zechar, Gerstenberger, & Rhoades, 2010; Zechar, Schorlemmer, et al., 2010), referred to as consistency tests, the likelihood function evaluated on observed data is compared to likelihood values generated by the forecast thought of as a random process. Specifically, consistency tests determine whether L(Ω,Λ) falls within or outside a given confidence interval of a likelihood distribution L̂ = {L(Λ j,Λ)} obtained from a set of synthetic catalogs {Λj} simulated using a Poisson process (e.g., Kagan & Jackson, 1995). A forecast is deemed inconsistent when a quantile score γ, defined as γ = |{L(Λ j, Λ) such that L(Λ j, Λ)< L(Ω,Λ)} ⃒ ⃒ ⃒ ⃒{Λ j} ⃒ ⃒ , (7) satisfies γ < p for a given critical value p, where |A| denotes the number of elements contained in a set A. The quantile score can be interpreted as the probability that a likelihood score sampled from L̂ is lower than that of the observationsΩ. This likelihood consistency test (L‐test) is complemented with the Number (N‐), magnitude (M‐) Figure 3. Schematic illustration of the concept of the log‐likelihood score. (a) Observed numbers of earthquakes in individual cells of a spatial grid. (b) Comparison of likelihoods of the observed numbers under two alternative models for a particular cell. To obtain the log‐likelihood score of a model, the logarithms of the likelihoods are summed over all grid cells (Figure adapted from Mizrahi et al. (2024)). Reviews of Geophysics 10.1029/2023RG000823 MIZRAHI ET AL. 14 of 70 19449208, 2024, 3, D ow nloaded from https://agupubs.onlinelibrary.w iley.com /doi/10.1029/2023R G 000823 by M assey U niversity L ibrary, W iley O nline L ibrary on [23/10/2024]. See the T erm s and C onditions (https://onlinelibrary.w iley.com /term s-and-conditions) on W iley O nline L ibrary for rules of use; O A articles are governed by the applicable C reative C om m ons L icense https://agupubs.onlinelibrary.wiley.com/action/rightsLink?doi=10.1029%2F2023RG000823&mode= and Spatial (S‐) tests, which follow the same principle, but the likelihoods are calculated from forecasts reduced to the corresponding N‐, M‐, or S‐space (e.g., only the total number of events in the region; only the magnitude distribution; only the spatial distribution), by summing the mean rate of Poisson distributions across the bins of the remaining dimensions. An extensive description, reformulation and implementation of these tests can be found in Schorlemmer et al. (2007), Zechar, Gerstenberger, and Rhoades (2010), Zechar, Schorlemmer, et al. (2010), Werner et al. (2011), Rhoades et al. (2011), and Savran et al. (2022). Due to the spatio‐temporal clustering and correlations in earthquake distributions, limitations of assuming a Poisson process when estimating L̂ quickly became evident in short‐term forecasting experiments (e.g., Lom- bardi, 2014; Lombardi & Marzocchi, 2010a, 2010b; Nandan et al., 2019; Werner & Sornette, 2008). But the Poisson and cell independence approximations can also break down in long‐term forecast experiments at higher magnitude thresholds when clustering persists. To reduce the undesired strong influence of clustering within cells on testing results, Bayona et al. (2022) used a binary likelihood of observing either zero or n ≥ 1 earthquakes in each cell. To avoid the unrealistically narrow Poisson distribution for the total number distribution, Werner et al. (2010) used an overdispersed negative‐binomial distribution to describe and assess the total number forecasts of models (see also e.g. Kagan, 2010). Nandan et al. (2019) smoothed ETASmodel simulations to derive ETAS‐specific cell‐wise likelihood functions, thus representing the marginals within each cell well, while maintaining the cell‐independence assumption. 2.2.2.4. Simulation‐Based Testing From a modeling perspective, dependencies between earthquakes are already captured in simulation‐based forecasting methods, such as ETAS, which provide synthetic catalogs directly through simulations, rather than from assumed statistical distributions of event counts in given cells. Savran et al. (2020) thus extended the traditional CSEP tests by relaxing the assumption that earthquakes follow independent Poisson distributions in discretized space‐time‐magnitude regions. Instead, CSEP catalog‐based tests are based on pseudo‐likelihoods (and other metrics) that are aggregated over the likelihood scores of earthquakes in the test set, rather than over the cells/bins of the forecast region (as in Equation 6), by using the log‐likelihood definition for a continuous point process, LL(Ω, Λ) =∑ N i=1 ln λ(ei|Ht) − ∫ R λ(e|Ht) dR, (8) where λ(Ht) is the conditional intensity function specifying the point process,Ht is the earthquake history until the time t of the forecast, and R is the (space‐time‐magnitude‐) region for which the forecast is issued, and the summation is over all events ei ∈ Ω. We direct the interested reader to Daley and Vere‐Jones (2003, 2008) for a thorough introduction to the theory of point processes. Because the analytical description of λ(Ht) may not be given, it is approximated from simulated catalogs as λ(Ht) = E[λ(e|Ht)|Rd] at a region discretization Rd. In this way, pseudo‐likelihood distributions can be built (and tested for consistency) for models that fully represent overdispersion and spatio‐temporal dependencies through simulations. 2.2.2.5. Other Methods Additional tests have been proposed in the literature, which do not depend on (pseudo‐) likelihood functions or have not been yet implemented in routine CSEP experiments. These include Turing‐style tests that determine if a set of simulations can imitate the behavior of observed catalogs (Page & van der Elst, 2018). For this, they devise a series of statistical descriptors that summarize the rate and magnitude distributions of a catalog, as well as foreshock/aftershock productivity and spatio‐temporal clustering. In the same fashion as consistency tests, Turing‐style tests observe if the observed catalog's descriptors lie within a confidence interval from those generated by a set of simulations. Furthermore, earthquake forecasts have been submitted to tests based on second‐order (spatial) statistics (e.g., Clements et al., 2011; Veen & Schoenberg, 2006), which are based on K‐ functions (Ripley, 2005) to determine the level of clustering or inhibition of a catalog. Clements et al. (2011) also describe a series of residual analysis for point‐process methods, which consist of transforming the points of a simulated/observed catalog (e.g., by rescaling, thinning, superpositioning), such that the resulting transformed process should be homogeneous‐Poisson if the original model were consistent with the observations. Reviews of Geophysics 10.1029/2023RG000823 MIZRAHI ET AL. 15 of 70 19449208, 2024, 3, D ow nloaded from https://agupubs.onlinelibrary.w iley.com /doi/10.1029/2023R G 000823 by M assey U niversity L ibrary, W iley O nline L ibrary on [23/10/2024]. See the T erm s and C onditions (https://onlinelibrary.w iley.com /term s-and-conditions) on W iley O nline L ibrary for rules of use; O A articles are governed by the applicable C reative C om m ons L icense https://agupubs.onlinelibrary.wiley.com/action/rightsLink?doi=10.1029%2F2023RG000823&mode= 2.3. State of Research in Earthquake Forecast Communication “What will happen next?” is one of the most frequently asked questions after an earthquake felt by the public. Short‐term earthquake forecasts allow us to reply to this question in a probabilistic way. Following the 2009 L’Aquila earthquake, a scientific discourse started on how (and if) earthquake forecasts should be made available to society (Jordan et al., 2011; Marzocchi, Iervolino, et al., 2015; Marzocchi, Jordan, & Woo, 2015). While the value of long‐term forecasts was clear (i.e., building codes, risk management plans (Marti et al., 2019)), the effective use of short‐term earthquake forecasts by societal stakeholders was still unexplored (Jordan et al., 2011). Through workshops with experts and stakeholders, key concepts and relevant issues to be considered when communicating forecasts to society were derived (Field et al., 2016; Marzocchi, Iervolino, et al., 2015; Mar- zocchi, Jordan, & Woo, 2015). Furthermore, New Zealand has evaluated earthquake forecast communication products (e.g., tables, maps) and provided advice for other regions aiming to communicate forecasts (Becker et al., 2018, 2019, 2020). Other research teams around the world have also been advancing the ways to best communicate earthquake forecasts (e.g., Dryhurst et al., 2022; Schneider, Wein, et al., 2023), which are sum- marized in the following paragraphs. Earthquake forecasts are useful for various stakeholders ranging from first responders and critical infrastructure owners to the general public (Field et al., 2016). Depending on their professions and experiences as well as the severity of the event, stakeholders have different information needs that evolve throughout the sequence (Becker et al., 2019). Becker et al. (2020) have shown that agencies and the public in New Zealand make use of earthquake forecasts for, among others, timing infrastructure repair and rebuilding, land‐use planning decisions, undertaking preparedness actions such as organizing essential supplies and emergency items, or deciding about safe access to buildings. This is in line with the purposes of earthquake forecasts in the United States (U.S.), where for example, used by the Federal Emergency Management Agency (FEMA) to set the temporal duration of disaster declara- tions (Michael et al., 2020; van der Elst et al., 2020). Furthermore, hospitals may use earthquake forecasts to postpone high‐risk treatment and elective surgery until the chance of aftershocks has significantly declined (Field et al., 2016). However, making sense of the probabilities of forecasts in ways that translate into actions and measures is not trivial and needs further investigation (Becker et al., 2018), especially in countries with moderate seismic risk where people only rarely experience earthquakes. Studies on stakeholders' information needs have shown that a wide range of information is requested, from basic facts about aftershocks through to more technical information. The general public, for example, desires infor- mation about the likely future earthquake locations, the likely future magnitudes/intensities and locations of aftershocks, the expected duration of the aftershock sequence, and guidance on what to do during the aftershock sequence (Becker et al., 2019). Communicating forecasts can also fill information vacuums that set the stage for amateur earthquake predictions and misinformation (Marzocchi, Iervolino, et al., 2015; Marzocchi, Jordan, & Woo, 2015). Furthermore, recently conducted workshops in the United States, Mexico, and El Salvador by a multinational team revealed that for numerous stakeholders (from emergency managers to engineers to science communicators), the spatial distribution of forecasted shaking was of top priority, though other types of infor- mation could support profession‐specific use cases (Schneider, Wein, et al., 2023). In general, it is crucial that the forecasts are put into a context relevant to the stakeholders, and the information is personalized (Becker et al., 2019). Moreover, one should make the purpose of the forecast clear, communicate the limitations, and explain the difference between a forecast and a warning (Freeman et al., 2023). Becker et al. (2020, p. 3343) thus recommend “developing a diversity of audience‐relevant OEF information for communication purposes, alongside advice on how such information could be utilized.” Forecasts should be provided in multiple formats (Becker et al., 2019) and with a layered approach, which means providing multiple pieces of information in a layered manner and, thus, allowing users to access the information they are interested in or need for their decision making (M. M. Wood et al., 2012, 2018). Studies have shown that maps were requested most frequently by users (e.g., Becker et al., 2018). If well‐designed, they can be an effective way to display the spatial distribution of natural hazards and risks (Stieb et al., 2019; Thompson Clive et al., 2021). So far, only a few studies have analyzed short‐term forecast maps (e.g., Becker et al., 2020; Schneider, McDowell, et al., 2022), but best practices from long‐term hazard maps can be used as a reference point (e.g., Dallo et al., 2023; Marti et al., 2019; Schneider, Cotton, & Schweizer, 2023). When using tables, it is recommended to put the forecasted number of shocks of a certain magnitude alongside the actual observed number for periods that have already passed because it can reassure information recipients of the accuracy of the Reviews of Geophysics 10.1029/2023RG000823 MIZRAHI ET AL. 16 of 70 19449208, 2024, 3, D ow nloaded from https://agupubs.onlinelibrary.w iley.com /doi/10.1029/2023R G 000823 by M assey U niversity L ibrary, W iley O nline L ibrary on [23/10/2024]. See the T erm s and C onditions (https://onlinelibrary.w iley.com /term s-and-conditions) on W iley O nline L ibrary for rules of use; O A articles are governed by the applicable C reative C om m ons L icense https://agupubs.onlinelibrary.wiley.com/action/rightsLink?doi=10.1029%2F2023RG000823&mode= forecast (Becker et al., 2019). Regarding textual information, Doyle et al. (2020) recommend to use the expression “within the next x” instead of “in the next x,” to avoid the propensity for people to think the event will happen at the end of this time window. Furthermore, Dryhurst et al. (2022) recommend presenting probabilities of events occurring as both percentages and expected frequencies. They conducted tests on diverse formulations of earthquake probabilities for online platforms, involving the general public across multiple countries. By considering best practices from various fields, they demonstrated that the following formulation enhances the public's comprehension of low percentages: “With current levels of seismic activity the chance of an earthquake of [insert magnitude or intensity information here] happening in [insert location information here] between [insert dates here] is: x%. Imagine 100,000 places with exactly the same chance of an earthquake as [insert location information here]. In the week of [insert dates here], with an x% chance, we would expect: An earthquake of [insert magnitude or intensity information here] to happen in y of these places. No earthquake of [insert magnitude or intensity information here] to happen in z of these places.” Given the implications that such framings can have on people's interpretation and use of forecasts, it is important to continue exploring the best ways of presenting such information. Studies have also identified various challenges. First, earthquake forecasts can cause anxiety, which can nega- tively affect people's psychosocial health and well‐being (Becker et al., 2019; Schneider, Wein, et al., 2023;Wein, Becker, et al., 2016; Wein, Potter, et al., 2016). Second, earthquake forecasts may negatively influence sectors such as tourism, if people perceive the earthquake risk as too high and delay their activities (Becker et al., 2018). Third, professionals and the general public struggle to interpret probabilistic information and translate it into actions (Becker et al., 2020), though effective methods for visualizing uncertainty have been found (Schneider, Freeman, et al., 2022), and there lies potential in providing advice alongside probabilities to guide people in making decisions (e.g., “check your preparedness”; M. M.Wood et al., 2012). Thereby, accounting for individual variations in social background, knowledge, and the associated response capabilities remains a challenge (Mileti & Sorensen, 1990). Fourth, people tend to compare earthquake forecasts with weather forecasts, thus expecting that the precise location and time of a future earthquake can be estimated (Dryhurst et al., 2022). Fifth, the absence of significant damage caused by a large event can lead to a normalization bias among non‐victims. This can limit their perception of the risk posed by subsequent damaging aftershocks and affect their responsiveness to forecasts (Mileti & O’Brien, 1992). Sixth, given that people have expressed willingness to know and understand the un- certainties associated with OEF (Becker et al., 2018), including stochastic (system variability) and epistemic (lack of knowledge) uncertainty, these warrant further study in terms of understanding how to communicate such concepts in ways that are meaningful to people and support action (Doyle et al., 2019; Hudson‐Doyle et al., 2018). Future research is needed to (a) test if people correctly interpret the forecasts and are able to apply them to practical actions (Doyle et al., 2020), (b) define thresholds, or impact‐based contexts, for when it is recommended to take certain actions (Becker et al., 2018), (c) develop interactive systems where each user can set their own thresholds since they are profession‐ and use‐case‐specific (Schneider, Wein, et al., 2023), (d) assess users' needs in countries aiming to communicate earthquake forecasts, following the approach of Schneider, Wein, et al. (2023), also enabling cross‐culture comparisons (Becker et al., 2018), and (e) evaluate how one could go from earthquake forecasts to earthquake loss forecasts (e.g., Iervolino et al., 2015). 3. Examples of OEF Systems Worldwide We here provide a comprehensive review of three examples of OEF systems in Italy, New Zealand, and the United States. This encompasses an overview of the model types, tests applied, and communication tools employed by these countries to provide earthquake and aftershock forecasts. The aim is to foster transparency in the interest of the public as well as modelers worldwide who intend to develop OEF systems. The countries covered here were selected because (a) they are among the few countries that do OEF according to the definition in this article, and (b) they are also represented by participants of the expert elicitation conducted by Mizrahi, Dallo, and Kuratle (2023) and presented in Section 4. For recent advancements in the development of OEF systems in China, we refer to Z. Liu et al. (2023), and to Omi et al. (2019) for the case of Japan. Furthermore, Kamer et al. (2021) and Nandan et al. (2021) describe the non‐authoritative “earthquake prediction platform” RichterX, which has provided global earthquake forecasts since 2019 and allows users to submit earthquake predictions. Reviews of Geophysics 10.1029/2023RG000823 MIZRAHI ET AL. 17 of 70 19449208, 2024, 3, D ow nloaded from https://agupubs.onlinelibrary.w iley.com /doi/10.1029/2023R G 000823 by M assey U niversity L ibrary, W iley O nline L ibrary on [23/10/2024]. See the T erm s and C onditions (https://onlinelibrary.w iley.com /term s-and-conditions) on W iley O nline L ibrary for rules of use; O A articles are governed by the applicable C reative C om m ons L icense https://agupubs.onlinelibrary.wiley.com/action/rightsLink?doi=10.1029%2F2023RG000823&mode= We divide the relevant issues into three pillars: the development, testing, and communication of earthquake forecasts. At the end of each pillar, a summary is provided in table format. After the description of the OEF systems of the three countries, the last part of this section of the review will describe the concept of Operational Earthquake Loss Forecasting (OELF): how to turn time‐variant earthquake probabilities into time‐variant esti- mates of seismic hazard and risk. 3.1. Model Development 3.1.1. Italy The probabilistic forecasts delivered by the OEF‐Italy system are obtained from an ensemble model, which is the weighted combination of three stochastic models typically used in statistical seismology: Epidemic‐Type Aftershock Sequence (ETAS) (Lombardi & Marzocchi, 2010a, 2010b), Epidemic‐Type Earthquake Sequence (ETES) (Falcone et al., 2010) and Short‐Term Earthquake Probability (STEP) (Woessner et al., 2010). Following the score model averaging (SMA) procedure, the weights are assigned depending on the forecasting skill that the single models have shown considering the catalog of the last week, and are updated daily (Marzocchi, Zechar, & Jordan, 2012). The system also stores the forecasts computed by the single components of the ensemble. Nevertheless, they are not considered in the output. The ensemble model is expected to be the most reliable, or never much worse than any (unknown) best‐performing model. More precisely, the ensemble model in OEF‐Italy derives from the merging model approach, which accounts for the uncertainties of the components (Vere‐Jones, 1995). Necessarily, the assignment of the weight is susceptible to subjectivity, but it still gives a gain in reliability when the selection satisfies the following three properties (Garthwaite & Mubwandarikwa, 2012): dilution (assign smaller weights to the models highly correlated to others), strong dilution (the weight of a model has to be split with a newly added model that is equal to the first one), monotonicity (if a new model is added, the others' weights should not increase). This is compiled by OEF‐ Italy, where ETAS, ETES and STEP are combined through the score model averaging, according to which the weights are updated every week, on non‐overlapping temporal windows, by considering the logarithm of the cumulative spatio‐temporal likelihood scored by every model until the starting forecasting time (Marzocchi, Douglas Zechar, & Jordan, 2012; Marzocchi et al., 2014). This procedure guarantees the creation of ensemble forecasts never much worse than those produced by each component. The models currently embedded in the OEF‐Italy system are probabilistic models based on differently param- eterized spatiotemporal kernels. Their forecasts are heavily dependent on past seismicity, as are all clustering models; they feed on progressively occurring events, and are particularly uncertain early in a sequence when limited observations poorly constrain the productivity. The magnitudes adopted to compute the probabilistic forecasts are all local magnitudes (ML), this being the type estimated for the great majority of events recorded by the Italian Seismic Monitoring Room of the Istituto Nazionale di Geofisica e Vulcanologia (National Institute of Geophysics and Volcanology, INGV). Only for the few strongest events a moment magnitude (MW) is also estimated, as theML may saturate. However, to keep the magnitudes homogeneous throughout the entire seismic catalog used to produce forecasts, only the ML values are considered for all events. As better specified in a later, the OEF‐Italy system operates in real‐time at a national scale by delivering weekly forecasts over a grid lattice of 10 km × 10 km cells covering the entire Italian territory, with no site‐specific parametrization. The models have been calibrated during a learning period from April 2005 to April 2009 for the whole national territory, without accounting for any spatio‐temporal variation. Productivity is not updated during specific seismic sequences. Incompleteness entailed by strong events is also not automatically corrected in the current version of the system; so far, corrections for incompleteness are applied by hand only immediately after a large earthquake. A first attempt to make this correction applied automatically has been proposed by Stallone and Falcone (2021), but these authors' method has not been introduced in the OEF‐Italy system yet. Although the current version of the OEF‐Italy model includes only clustering models, there are no limitations on the number of forecasting models to be added to the ensemble. Indeed, we argue that using models based on different ideas and thoughts is beneficial, and may allow the ensemble model to perform better. The only constraint imposed on the candidate forecasting models is that they must be submitted to one or more CSEP experiments to get an independent and transparent evaluation. Reviews of Geophysics 10.1029/2023RG000823 MIZRAHI ET AL. 18 of 70 19449208, 2024, 3, D ow nloaded from https://agupubs.onlinelibrary.w iley.com /doi/10.1029/2023R G 000823 by M assey U niversity L ibrary, W iley O nline L ibrary on [23/10/2024]. See the T erm s and C onditions (https://onlinelibrary.w iley.com /term s-and-conditions) on W iley O nline L ibrary for rules of use; O A articles are governed by the applicable C reative C om m ons L icense https://agupubs.onlinelibrary.wiley.com/action/rightsLink?doi=10.1029%2F2023RG000823&mode= 3.1.2. New Zealand In New Zealand, earthquake forecasting occurs on three different time scales: short‐term, medium‐term and long‐ term, which are each associated with different observed phenomena and are described by different models as outlined in Table 1. Researchers in New Zealand have long been engaged in the development of earthquake forecasting models, most notably, the Short‐Term Earthquake Probability (STEP) model (Gerstenberger et al., 2004, 2005) and the medium‐term Every Earthquake a Precursor According to Scale (EEPAS) model (e.g., Rhoades & Evison, 2004). The EEPAS model software package also includes an Epidemic‐Type Aftershock Sequence (ETAS) model, and there is a stand‐alone ETAS implementation for New Zealand, referred to here as ETAS‐Harte (Harte, 2013). New Zealand has a suite of long‐term models that contribute to earthquake fore- casting (e.g., Gerstenberger et al., 2023) in hybrid forecast models. The development, testing, and application of earthquake forecasting models is strongly dependent on the availability of earthquake catalog data. Prior to 2012, the earthquake location was not automated and there was often a delay in data processing, leading to delays and gaps in the earthquake catalog, particularly during aftershock sequences. Regional moment magnitudes were introduced in 2007 (Ristau, 2008) and highlighted that local magnitudes for moderate‐sized earthquakes systematically overestimated the size of earthquakes compared to moment magnitude (Ristau, 2009). In 2012, the earthquake location system changed to the automated system SeisComp3, which by default uses Californian attenuation relations for local magnitudes. Magnitudes ofM ≥ 5.0 were found to have a similar bias as previously (Rhoades et al., 2017), and are converted to the expected rate of moment magnitude for hazard es- timates (Rhoades & Christophersen, 2017). In contrast, there are many fewer earthquakes in the magnitude range of 3–5, particularly for the central North Island (Harte, 2019). To overcome the differences between the magnitude scales in the different time periods of the catalog, a new magnitude was derived, MLNZ20 (Rhoades et al., 2021). MLNZ20 uses attenuation relations and a depth dependence to be on average as consistent with the regional moment magnitude as possible. For the revision of the NSHM, a catalog was derived with magnitudes as consistent as possible with moment magnitude for the instrumentally located since 1931 (Christophersen et al., 2022). Work is ongoing to implement MLNZ20 into the routine earthquake processing, and to integrate the revised magnitudes into the standard earthquake catalog. At this stage, the models involved in public earthquake forecasting have their parameters derived prior to 2012. Therefore, there is hesitation to forecast earthquakes of magnitude M < 5.0 until the new magnitudes are fully implemented and available in real‐time, as well as all model parameters refitted. Below we describe each model in more detail, as well as the hybrid forecasts that combine short‐, medium and long‐term models for public earthquake forecasting. 3.1.2.1. The STEP Model The STEP model applied for public forecasting in New Zealand is an aftershock forecasting model purely based on seismicity. It applies the Reasenberg and Jones model (Reasenberg & Jones, 1989, 1990) with generic pa- rameters for New Zealand (Pollock, 2007) to each earthquake above a chosen cut‐off magnitude. Within a region, the earthquake with the highest forecast rate in a chosen time window contributes to the total forecast. Higher‐ order aftershocks are automatically included when updating the forecast during an ongoing sequence. Each rate contribution is initially circular in space. However, as the sequence evolves and rates from different regions dominate, the spatial pattern can become more complex. The model also has the option to fit a two‐segment fault line, or have a simplified fault line supplied, and then spatially decay the rates with one over distance from the fault. Table 1 Overview of the Time Scales, Observed Phenomena, and Models Used for Earthquake Forecasting in New Zealand Time scale Short‐term Medium‐term Long‐term Relevant time window Hours, days and weeks Months to years Decades Observed phenomenon Temporal clustering following large earthquakes Precursory seismicity increase prior to large earthquakes Long‐term spatial clustering Models used STEP, ETAS EEPAS PPE, NSHM background Reviews of Geophysics 10.1029/2023RG000823 MIZRAHI ET AL. 19 of 70 19449208, 2024, 3, D ow nloaded from https://agupubs.onlinelibrary.w iley.com /doi/10.1029/2023R G 000823 by M assey U niversity L ibrary, W iley O nline L ibrary on [23/10/2024]. See the T erm s and C onditions (https://onlinelibrary.w iley.com /term s-and-conditions) on W iley O nline L ibrary for rules of use; O A articles are governed by the applicable C reative C om m ons L icense https://agupubs.onlinelibrary.wiley.com/action/rightsLink?doi=10.1029%2F2023RG000823&mode= The start of the earthquake catalog is a model parameter, so the user can choose howmuch historical data to fit. To address incompleteness in an ongoing sequence, the model fits the magnitude of completeness in each run and then fits parameters to a particular earthquake sequence, once there ar