The “Every Earthquake a Precursor According to Scale” (EEPAS) model uses minor earthquakes in a catalogue to forecast major earthquakes, based on the precursory scale increase phenomenon and associated predictive scaling relations. It has previously been applied to regional earthquake catalogues for medium-term forecasting of earthquakes with target magnitudes exceeding thresholds in the range M 5.0–6.0 using previous earthquakes about two magnitude units below the target threshold, depending on the completeness of the catalogue. Here, we apply it to the ISC-GEM global earthquake catalogue and test its ability to forecast earthquakes with target magnitudes exceeding M 7.5, i.e., the largest and most important earthquakes for assessing earthquake hazard. We use background models based on both smoothed seismicity and spatial variation of strain rate estimates. The models are fitted to the period 1994–2013, with 85 target earthquakes, and tested on the period 2014–2019, with 19 target earthquakes. The EEPAS model shows larger information gains over a smoothed seismicity model than in most regional applications. Also, a hybrid composed of EEPAS, smoothed seismicity and the strain rate outperforms a hybrid consisting only of smoothed seismicity and strain rate components. These preliminary results suggest that EEPAS may be useful for forecasting the largest earthquakes globally but should be followed up by transparent prospective testing. There are several avenues for improvement of EEPAS modeling.

Earthquake forecasts estimate the likelihood of seismic activity within a specific region over a given timeframe, utilizing historical data and patterns from past earthquakes. In New Zealand, the GeoNet program within GNS Science is the main source of geological hazard information and has publicly provided earthquake forecasts since the Darfield earthquake in September 2010. The generation and provision of initial forecasts and subsequent updates have relied on extensive time commitments of experts. The growing use and the desire to make forecast delivery less dependent on personnel capacity have motivated the development of a robust software solution through a hybrid forecast tool (HFT). The HFT is composed of forecast models that cover several different timescales: short term (ranging from a few hours to several years, based on empirical relations for aftershock decay), medium term (spanning years to decades, utilizing the increased seismic activity preceding major earthquakes), and long term (covering decades to centuries, combining information from the spatial distribution of cataloged earthquake locations and slip rates of mapped faults and strain rates estimated from geodetic data). Originally, these models were developed over many years by individual researchers using various programming languages such as Fortran, Java, and R, operating on separate operating systems, with their features documented and published. The HFT unites these models under one umbrella, utilizing a Docker container to navigate disparate software library compatibility issues. Furthermore, the HFT offers user-friendly navigation through a graphical user interface and a command-line feature, facilitating the configuration of automatic and periodic forecast runs. The stability and integration provided by the HFT greatly improve the capability of GNS Science to provide forecasts that inform responses to significant regional seismic events and bring New Zealand closer to automated and operational earthquake forecasting. Although HFT is specifically designed for New Zealand’s earthquake forecasting, the framework, implementation, and containerization approach could also benefit forecasting efforts in other regions.

Abstract New Zealand entered a period of increased seismic activity in 2003, which, to date, has included around 20 large or damaging earthquakes. Building on decades of forecast model research and including research into formal model evaluation, GNS Science, the governmental institution tasked with providing natural hazards information, began issuing public earthquake forecasts following the 2009 Mw 7.2 Darfield earthquake. This article provides a review of the public earthquake forecasting methods and outcomes since that time. Initially the release of forecasts was motivated by the scientific team, but interest, use and understanding quickly rose in the public, government, and private sectors which led to regular refinement and additions to the type of forecast information provided. The basic tenet of forecasting has remained the same, with the forecasts being used to inform decisions requiring consideration of timeframes from days to decades. The operational models were all based on the previous 10–20 years of research in model development and testing; this included the development of methods for combining multiple models into a single hybrid model that typically provides increased statistical forecast skill. The hybrid models combine models from three classes of time horizons: (1) short-term aftershock clustering based on the epidemic-type aftershock sequence class of models, (2) medium-term clustering with time horizons of years to decades, and (3) long-term models with time horizons of decades to nominally time independent. Forecasts have been issued in response to 14 large or significant earthquakes with the forecasts and accompanying information regularly updated on the GeoNet webpages. An important component to the successful deployment and uptake of the models was the foundation laid in the previous decades in model development and testing. This provided understanding and confidence in the models’ ability to provide useful forecasts for response and recovery decisions.

Probabilistic seismic hazard analysis requires a seismicity rate model, or in other words, a forecast of earthquake rates. In the New Zealand National Seismic Hazard Model 2022, the seismicity rate model is constructed through independent forecasts of earthquakes on mapped faults and earthquakes distributed over cells in a spatial grid. Here, we explore the seismicity rate model for upper plate (hypocenter >= 40 km) events, to investigate the shape of magnitude-frequency distributions (MFDs) considering events nucleating (or for which the hypocenters are located) within individual fault zone. We find that more than 80% of the fault zones have MFDs that are better described by a Gutenberg-Richter (GR) distribution, instead of a characteristic distribution (i.e., rates of larger magnitudes much higher than the GR trend). Furthermore, the MFD classifications are neither influenced by time -dependent (and time -independent) considerations nor directly affected by the size (or area) of the fault zones. Fault zones with faster slip rates ( > 20 mm/yr) exhibit characteristic MFDs, whereas those with slower slip rates may or may not. Although multifault ruptures are prevalent in the characteristic distributions, large maximum magnitude (M-w > 8.0) plays a pivotal role producing a characteristic MFD. On the other hand, physically unconnected multifault ruptures (i.e., involving rupture jumps >= 10 km) are mostly observed with GR distributions.

Recurrence intervals of ground-surface rupturing earthquakes are one of numerous datasets used to constrain the rates of fault ruptures in the 2022 revision of the New Zealand National Seismic Hazard Model (NZ NSHM 2022). Paleoearthquake timing and single-event displacement (SED) data in the New Zealand Paleoseismic Site Database version 1.0 alongside geologic and geodetic slip rates from the New Zealand Community Fault Model version 1.0 and NZ NSHM 2022 Geodetic Deformation Model were used to estimate recurrence intervals on faults across New Zealand for inclusion in the NZ NSHM 2022. Past earthquake timings were fit with lognormal, exponential, and Brownian Passage Time recurrence models to derive probability density functions (PDFs) of mean recurrence interval (MRI) in a Bayesian framework. At some sites, SED and slip-rate (SR) data were used to estimate PDFs of MRI; and at sites where timings, slip rate, and displacement data are available, the timings-based and slip-based PDFs were combined to develop tighter constraints on MRI. Using these approaches, we produce a database of maximum-likelihood MRIs and their uncertainties for 80 sites across New Zealand. The resulting recurrence interval dataset is publicly available and is the largest such dataset in New Zealand to date. It provides a valuable resource for future seismic hazard modeling and highlights areas that would benefit from future study.

ABSTRACT Using a new integrated earthquake catalog for Aotearoa New Zealand (described in a companion article), we estimate the magnitude–frequency distributions (MFDs) of earthquakes in the greater New Zealand region and along the Hikurangi–Kermadec and Puysegur subduction zones. These are key inputs into the seismicity rate model (SRM) component of the 2022 New Zealand National Seismic Hazard Model. The MFDs are parameterized by a b-value (describing the relative rates of small and large earthquakes) with its epistemic uncertainty expressed by three logic tree branches (low, central, and high), and by the annual rate of M ≥ 5 earthquakes, here called the N-value, which has a separate value conditioned on each b-value branch. The N-value has its own epistemic uncertainty besides the dependence on the b-value, and this is also estimated here and propagated through the SRM by scaling all event rates up and down by a “low” and a “high” scalar value on either side of 1.0, called “N scaling.” Adapting an approach used previously in California, we estimate these MFD parameters in the onshore and near-shore region incorporating data back to 1843, balanced with the better data in the more recent part of the instrumental catalog. We estimate the MFD parameters on the Hikurangi–Kermadec and Puysegur subduction zones using a slightly simplified version of this approach and more recent data. We then use a globally-based method to estimate the potential earthquake rate uncertainty on the Hikurangi–Kermadec subduction zone and an SRM-specific moment-rate-related argument to construct an appropriately wide rate uncertainty for the Puysegur subduction zone.

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.

The precursory scale increase (Psi) phenomenon describes the sudden increase in rate and magnitude in a precursory area AP, at precursor time TP, and with precursor magnitude MP prior to the upcoming large earthquake with magnitude Mm. Scaling relations between the Psi variables form the basis of the "Every Earthquake a Precursor According to Scale" (EEPAS) earthquake forecasting model. EEPAS is a well-established space-time point process model that forecasts large earthquakes in the medium term, that is, the coming months to decades, depending on Mm. In Aotearoa New Zealand, EEPAS contributes to hybrid models for public earthquake forecasting and to the source model of time-varying seismic hazard models, including the latest revision of the National Seismic Hazard Model. The Psi phenomenon was recently shown not to be unique for a given earthquake, with smaller precursory areas AP associated with larger precursor times TP and vice versa. This trade-off between AP and TP has also been found for the spatial and temporal distributions of the EEPAS models. Detailed analysis of the Psi phenomenon has so far been limited by the manual and labor-intensive procedure of identifying Psi in earthquake catalogs. Here, we introduce two algorithms to automatically detect Psi and apply them to real and simulated earthquake catalog data. By randomizing the catalog and removing aftershocks, we confirm that the Psi phenomenon is a feature of space-time earthquake clustering prior to major earthquakes. Multiple Psi identifications confirm the trade-off between AP and TP, and the scaling relations for both real and simulated catalogs are consistent with the original scaling relations on which EEPAS is based. We identify opportunities for future work to refine the algorithms and apply them to physics-based simulated catalogs to enhance the understanding of Psi. A better understanding of Psi has the potential to improve forecasting of large upcoming earthquakes.

The 2022 revision of Aotearoa New Zealand National Seismic Hazard Model (NZ NSHM 2022) has involved significant revision of all datasets and model components. In this article, we present a subset of many results from the model as well as an overview of the governance, scientific, and review processes followed by the NZ NSHM team. The calculated hazard from the NZ NSHM 2022 has increased for most of New Zealand when compared with the previous models. The NZ NSHM 2022 models and results are available online.

>Since the inaugural international collaboration under the framework of the Collaboratory for the Study of Earthquake Predictability(CSEP) in 2007, numerous forecast models have been developed and operated for earthquake forecasting experiments across CSEP testing centers(Schorlemmer et al., 2018). Over more than a decade, efforts to compare forecasts with observed earthquakes using numerous statistical test methods and insights into earthquake predictability, which have become a highlight of the CSEP platform.

We propose a method to estimate the uncertainty of the average rate of earthquakes exceeding a magnitude threshold in a future period of given length based on observed variability of the earthquake process in an existing catalog. We estimate the ratio R of the variability to that of a stationary Poisson process. R is estimated from subsets of the catalog over a wide range of timescales. The method combines the epistemic uncertainty in estimating the rate from the catalog and the aleatory variability of the rate in future time periods. If R is stable over many timescales, there is a solid basis for estimating the uncertainty of earthquake rate estimates. In the 2022 revision of the New Zealand National Seismic Hazard Model (NZ NSHM), estimation of the total shallow earthquake rate over the next 100 yr and its uncertainty is an important element. Using a 70 yr New Zealand catalog with hypocentral depths s 40 km and standardized magnitudes M z 4.95, we find stable estimates of R for timescales from 3 days to 2.4 yr. This gives a standard error of 0.95 on the estimated annual rate of M z 4.95, in the next 100 yr. R becomes unstable and has poor precision for longer subperiods. We investigate potential causes using synthetic catalogs with known inhomogeneities. Analysis of International Seismological Centre -Global Earthquake Model (ISC-GEM) catalog, to investigate the effect of higher magnitude thresholds, shows that R is lower for M z 6.95 than for M z 5.45. The ISC-GEM catalog restricted to New Zealand gives comparable stable estimates of R to the NZ NSHM 2022 catalog for M z 5.45 and lower estimates than the NZ NSHM 2022 catalog for M z 4.95. We also verify that magnitude standardization of the New Zealand GeoNet catalog has reduced the uncertainty of rate estimates by decreasing R throughout the entire range of timescales.

ABSTRACT A seismicity rate model (SRM) has been developed as part of the 2022 Aotearoa New Zealand National Seismic Hazard Model revision. The SRM consists of many component models, each of which falls into one of two classes: (1) inversion fault model (IFM); or (2) distributed seismicity model (DSM). Here we provide an overview of the SRM and a brief description of each of the component models. The upper plate IFM forecasts the occurrence rate for hundreds of thousands of potential ruptures derived from the New Zealand Community Fault Model version 1.0 and utilizing either geologic- or geodetic-based fault-slip rates. These ruptures are typically less than a couple of hundred kilometers long, but can exceed 1500 km and extend along most of the length of the country (albeit with very low probabilities of exceedance [PoE]). We have also applied the IFM method to the two subduction zones of New Zealand and forecast earthquake magnitudes of up to ∼Mw 9.4, again with very low PoE. The DSM combines a hybrid model developed using multiple datasets with a non-Poisson uniform rate zone model for lower seismicity regions of New Zealand. Forecasts for 100 yr are derived that account for overdispersion of the rate variability when compared with Poisson. Finally, the epistemic uncertainty has been modeled via the range of models and parameters implemented in an SRM logic tree. Results are presented, which indicate the sensitivity of hazard results to the logic tree branches and that were used to reduce the overall complexity of the logic tree.

The Collaboratory for the Study of Earthquake Predictability (CSEP) is a global commu-nity dedicated to advancing earthquake predictability research by rigorously testingprobabilistic earthquake forecast models and prediction algorithms. At the heart of thismission is the recent introduction of pyCSEP, an open-source software tool designed toevaluate earthquake forecasts. pyCSEP integrates modules to access earthquake cata-logs, visualize forecast models, and perform statistical tests. Contributions from theCSEP community have reinforced the role of pyCSEP in offering a comprehensive suiteof tools to test earthquake forecast models. This article builds onSavran, Bayona,etal.(2022), in which pyCSEP was originally introduced, by describing new tests and recentupdates that have significantly enhanced the functionality and user experience ofpyCSEP. It showcases the integration of new features, including access to authoritativeearthquake catalogs from Italy (Bolletino Seismico Italiano), New Zealand (GeoNet), andthe world (Global Centroid Moment Tensor), the creation of multiresolution spatial fore-cast grids, the adoption of non-Poissonian testing methods, applying a global seismicitymodel to specific regions for benchmarking regional models and evaluating alarm-based models. We highlight the application of these recent advances in regional studies,specifically through the New Zealand case study, which showcases the ability of py CSEPto evaluate detailed, region-specific seismic forecasts using statistical functions. The enhancements in pyCSEP also facilitate the standardization of how the CSEP forecast experiments are conducted, improving the reliability, and comparability of the earth-quake forecasting models. As such, pyCSEP exemplifies collaborative research and innovation in earthquake predictability, supporting transparent scientific practices, and community-driven development approaches.

ABSTRACT We develop candidate hybrid models representing the spatial distribution of earthquake occurrence in New Zealand over the next 100 yr. These models are used within the onshore/near-shore, shallow component of the distributed seismicity model within the New Zealand National Seismic Hazard Model 2022. They combine a variety of spatially gridded covariates based on smoothed seismicity, strain rates, and proximity to mapped faults and plate boundaries in both multiplicative and additive hybrids. They were optimized against a standardized catalog of New Zealand earthquakes with magnitude M ≥ 4.95 and hypocentral depth ≤40 km from 1951 to 2020. We extract smoothed seismicity covariates using three different methods. The additive models are linear combinations of earthquake likelihood models derived from individual covariates. We choose three preferred hybrid models based on the information gain statistics, consideration of the ongoing Canterbury sequence and regions of low seismicity, and inclusion of the most informative covariates. Since the hazard model is designed for the next 100 yr, the preferred hybrid models are also combined with 20-year earthquake forecasts from the “Every Earthquake a Precursor According to Scale” model. Thus, in total, six hybrid spatial distribution candidates are advanced for sensitivity analyses and expert elicitation for inclusion in the final logic tree for the New Zealand National Seismic Hazard Model.

ABSTRACT The 2022 revision of the New Zealand National Seismic Hazard Model—Te Tauira Matapae Pūmate Rū i Aotearoa—requires an earthquake catalog that ideally measures earthquake size in moment magnitude. However, regional moment tensor solutions, which allow the calculation of moment magnitude MwNZ, were introduced in New Zealand only in 2007. The most reported magnitude in the national New Zealand earthquake catalog is a variation of local magnitude ML. In New Zealand, ML is systematically larger than MwNZ over a wide magnitude range. Furthermore, the introduction of the earthquake analysis system SeisComP in 2012 caused step changes in the catalog. We address the problems by converting magnitudes using regressions to define a standardized magnitude as a proxy for MwNZ. A new magnitude, MLNZ20, has an attenuation relation and station corrections consistent on average with MwNZ. We have calculated MLNZ20 for nearly 250,000 earthquakes between 2000 and 2020. MLNZ20 is a reasonable proxy for MwNZ for earthquakes with ML<5.5. For earthquakes with ML>4.6, MwNZ is reliably available. We have applied ordinary least squares (OLS) regression for MwNZ and MLNZ20 on ML before and after 2012. We argue that OLS is the most appropriate method to calculate a proxy for MwNZ from individual ML measurements. The slope of the OLS regression compares well to the slope from the method of moments, which accommodates equation error that is present when there is scatter beyond measurement error, as is the case for our magnitude data. We have defined as a proxy for MwNZ a standardized magnitude Mstd, which is Mw when available, MLNZ20 with some restrictions as a second choice, and otherwise the magnitude derived from regression. Standardization of the magnitudes reduces the total number of earthquakes with a magnitude of ≥4.95 by more than half and corrects step changes in the spatial distribution of earthquakes between 2011 and 2012.

Following the 2016 M7.8 Kaikoura earthquake, a time-varying seismic hazard model (KSHM) was developed to inform decision-making for the reinstatement of road and rail networks in the northern South Island. The source model is the sum of a gridded 100-year earthquake clustering model and an updated fault source model. The gridded model comprises long-term, medium-term and short-term components. The 100-year gridded model is constructed as the sum of 100 annual forecasts. A discounting method trades off expected earthquake occurrences of the distant future against those of the near future. The fault source model includes updates to account for newly revealed faults that ruptured in the Kaikoura earthquake and other recently obtained new information, and new time-varying probabilities of rupture for four fault segments. Two different characterisations of the Hikurangi subduction interface are incorporated via a logic tree, with weights determined by an expert panel. A suite of ground motion prediction equations contribute to a logic tree in order to account for epistemic uncertainties in source modelling for each of four tectonic region types. Here, we compare the resulting hazard estimates with the 2010 National Seismic hazard Model and recorded motions in past New Zealand and global earthquakes.

SUMMARY The Every Earthquake a Precursor According to Scale (EEPAS) forecasting model is a space–time point-process model based on the precursory scale increase ($\psi $ ) phenomenon and associated predictive scaling relations. It has been previously applied to New Zealand, California and Japan earthquakes with target magnitude thresholds varying from about 5–7. In all previous application, computations were done using the computer code implemented in Fortran language by the model authors. In this work, we applied it to Italy using a suite of computing codes completely rewritten in Matlab. We first compared the two software codes to ensure the convergence and adequate coincidence between the estimated model parameters for a simple region capable of being analysed by both software codes. Then, using the rewritten codes, we optimized the parameters for a different and more complex polygon of analysis using the Homogenized Instrumental Seismic Catalogue data from 1990 to 2011. We then perform a pseudo-prospective forecasting experiment of Italian earthquakes from 2012 to 2021 with Mw ≥ 5.0 and compare the forecasting skill of EEPAS with those obtained by other time independent (Spatially Uniform Poisson, Spatially Variable Poisson and PPE: Proximity to Past Earthquakes) and time dependent [Epidemic Type Aftershock Sequence (ETAS)] forecasting models using the information gain per active cell. The preference goes to the ETAS model for short time intervals (3 months) and to the EEPAS model for longer time intervals (6 months to 10 yr).

Nearly 20 years ago, the observation that major earthquakes are generally preceded by an increase in the seismicity rate on a timescale from months to decades was embedded in the “Every Earthquake a Precursor According to Scale” (EEPAS) model. EEPAS has since been successfully applied to regional real-world and synthetic earthquake catalogues to forecast future earthquake occurrence rates with time horizons up to a few decades. When combined with aftershock models, its forecasting performance is improved for short time horizons. As a result, EEPAS has been included as the medium-term component in public earthquake forecasts in New Zealand. EEPAS has been modified to advance its forecasting performance despite data limitations. One modification is to compensate for missing precursory earthquakes. Precursory earthquakes can be missing because of the time-lag between the end of a catalogue and the time at which a forecast applies or the limited lead time from the start of the catalogue to a target earthquake. An observed space-time trade-off in precursory seismicity, which affects the EEPAS scaling parameters for area and time, also can be used to improve forecasting performance. Systematic analysis of EEPAS performance on synthetic catalogues suggests that regional variations in EEPAS parameters can be explained by regional variations in the long-term earthquake rate. Integration of all these developments is needed to meet the challenge of producing a global EEPAS model.

Strain rates have been included in multiplicative hybrid modelling of the long-term spatial distribution of earthquakes in New Zealand (NZ) since 2017. Previous modelling has shown a strain rate model to be the most informative input to explain earthquake locations over a fitting period from 1987 to 2006 and a testing period from 2012 to 2015. In the present study, three different shear strain rate models have been included separately as covariates in NZ multiplicative hybrid models, along with other covariates based on known fault locations, their associated slip rates, and proximity to the plate interface. Although the strain rate models differ in their details, there are similarities in their contributions to the performance of hybrid models in terms of information gain per earthquake (IGPE). The inclusion of each strain rate model improves the performance of hybrid models during the previously adopted fitting and testing periods. However, the hybrid models, including strain rates, perform poorly in a reverse testing period from 1951 to 1986. Molchan error diagrams show that the correlations of the strain rate models with earthquake locations are lower over the reverse testing period than from 1987 onwards. Smoothed scatter plots of the strain rate covariates associated with target earthquakes versus time confirm the relatively low correlations before 1987. Moreover, these analyses show that other covariates of the multiplicative models, such as proximity to the plate interface and proximity to mapped faults, were better correlated with earthquake locations prior to 1987. These results suggest that strain rate models based on only a few decades of available geodetic data from a limited network of GNSS stations may not be good indicators of where earthquakes occur over a long time frame.