Aftershock rate and hazard modelling for the Mw7.8 14/11/2016 Kaikoura Earthquake

13 November 2017: this document will be updated as appropriate. For more information, please contact M. Gerstenberger.

To date, we have modelled aftershocks for the Mw 7.8 Kaikoura Earthquake using several different models. Here we describe the basic procedure and parameters we have used. We describe the procedures based on the forecasts that were issued at different time-periods.

Following the M7.8 Kaikōura earthquake on 14 November 2016, GNS Science revised the 2010 National Seismic Hazard Model (NSHM) published by (Stirling et al., 2012).  A new model for distributed seismicity sources was developed in a similar way to the hybrid forecasting model that was constructed following the 2010 Darfield and 2011 Christchurch earthquakes (Gerstenberger et al., 2014, Rhoades et al., 2016, Gerstenberger et al., 2016).  The new model has three components, short-term, medium term and long-term.  This new model is now used to calculate regular forecasts for the Kaikōura earthquake sequence. Here we provide an overview of the model, describe briefly the contributing model, explain the maps of expected shaking that are published on GeoNet’s website as part of the forecast and also outline how the forecasting model has changed during the Kaikōura sequence.

Overview of the Kaikōura hybrid model

This model is a hybrid of time-varying and time-invariant components – a short-term component and a medium-term component (time varying) and a long-term component (time invariant).  Each component takes the form of a regional earthquake likelihood model (RELM; Schorlemmer and Gerstenberger (2007)) and consists of a forecast of the number of shallow (hypocentral depth < 40 km) earthquakes expected to occur in spatial cells on a 0.05-degree square grid, and within magnitude bins of width 0.1 units. The lowest magnitude bin covers the magnitude range from 4.95 to 5.05, and the highest magnitude bin covers the range from 7.95 to 8.05.

The short-term component is constructed as a 50:50 weighted average of two standard aftershock models, the STEP and ETAS models (further described below). These model aim to only forecast earthquakes triggered by previous earthquakes, i.e., aftershocks. Aftershocks occur after almost all large earthquakes, and also after many smaller earthquakes.  The expected number of aftershocks is highest immediately after an earthquake and thereafter decays like a power-law in time.  The decay of aftershock rates is referred to as the Omori-Utsu law (Utsu et al., 1995) and forms the basis for modelling aftershock occurrence. 

The medium-term component uses the Every Earthquake a Precursor According to Scale (EEPAS) model (Rhoades and Evison, 2004). The EEPAS model is based on an increase in the rate of occurrence and magnitude of minor earthquakes observed prior to the occurrence of major earthquakes. Associated with this precursory scale increase are predictive scaling relations (Evison and Rhoades, 2002, 2004).  

The long-term component is a weighted average of three different long-term models, which are described in more detail below.

The hybrid model is constructed by taking the highest rate contribution from each model component in each spatial and magnitude cell.

The individual forecast models

The STEP model
The Omori law is the oldest empirical relationship in seismology and originally described how the number of felt aftershocks per day decayed with time t as 1/t for the 1891 Nobi, Japan earthquake (Omori, 1894). The relationship was found to be still ongoing 100 years after the Nobi earthquake (Utsu et al., 1995).  The STEP (Short-Term Earthquake Probability) model is an aftershock model based on the idea of superimposed Omori type sequences (Ogata, 1988) that expected behaviour of aftershock sequences over a given time period (Gerstenberger et al., 2004, Gerstenberger et al., 2005).

The STEP model is comprised of two elements: 1) a background model, and; 2) a time-dependent clustering model.  The background model can consist of any model that is able to forecast a rate of events for the entire region of interest at all times. We have not applied a background model for this implementation.  The clustering model is based on the work of Reasenberg and Jones (1989), which defines aftershock forecasts based on the a- and the b- value from the Gutenberg-Richter relationship (Gutenberg and Richter, 1944) and the p-value from the modified Omori law (Utsu et al., 1995).  The STEP model defines these aftershock forecasts based on three model-components.  The first component is based on the average behaviour of aftershocks in New Zealand and uses the median Reasenberg and Jones parameter values for New Zealand aftershocks, with parameter estimates from (Pollock, 2007).  These are : a = -1.59; b = 1.03; p = 1.07; and c = 0.04. The second component uses the development of the ongoing aftershock sequence to attempt to improve the information in the forecast.  In this component, the Reasenberg and Jones parameters are automatically estimated for each individual aftershock sequence as it develops.  The third component allows for spatial heterogeneities within the sequence by calculating the Reasenberg and Jones parameters on the regional grid (i.e., 0.05˚×0.05˚).  The second and third components are only practical for large aftershock sequences with more than 100 aftershocks, therefore, in practice, the STEP forecasts are typically dominated by the first component.  In order to produce a single forecast based on the three component forecasts, each forecast is weighted by a score from the corrected Akaike Information Criterion (AIC) (Burnham and Anderson, 2002) and the three weighted forecasts are summed.  The AIC is based on a model's likelihood score, the amount of data and the number of free parameters that must be estimated, thus ensuring a smooth transition from the generic component to the spatially-heterogeneous component as more data becomes available.  Finally, for each grid node all aftershock forecasts are compared to the background forecast for the node and the highest forecast is retained.

The spatial distribution of the forecast aftershocks is partially decoupled from the rate calculation.  The rate forecast is isotropically distributed using a radius r based on the subsurface rupture length from Wells and Coppersmith (1994) and using a 1/r2 taper. If the radius is less than the grid spacing, a point source is assumed and the rates are smoothed in a circular fashion using the 1/r2 taper; otherwise a two segment fault is estimated based on the density of the aftershock distribution and the rates are smoothed using the 1/r2 taper from the fault. Importantly, because each aftershock can generate its own rate forecast as well as contributing to the overall rate forecast of the main shock, spatial heterogeneities are added by superimposing individual aftershock forecasts on the main shock forecast.  When more than one forecast is calculated for a single grid node, only the largest forecast is retained.

The ETAS model
The ETAS model is based on the modified Omori law for aftershock decay, just like the STEP model, but differs in the way model parameters are estimated and combined. ETAS forecasts have been calculated since the beginning of the Kaikoura sequence and are available separately here: (

The forecasts are based on the ETAS model, in particular, Model 6 in Harte (2013).  The forecasting method is further explained in Harte (2015).  To mix ETAS with STEP, we first forecast 2,000  simulated one-year earthquake catalogues using the ETAS model. From the simulated catalogues we calculate the expected rate within each grid cell over the forecast area. Finally, the total short-term forecast is calculated by combining 50% of the STEP forecast with 50% of the ETAS forecast in each space-magnitude bin.

The EEPAS model
The EEPAS model (Rhoades and Evison, 2004, Rhoades and Evison, 2005, Console et al., 2006, Rhoades and Evison, 2006, Rhoades, 2007) is concerned with the possibility that minor earthquakes may be precursory to major ones to follow them in the long run. It is based on the precursory scale increase phenomenon and associated predictive relations (Evison and Rhoades, 2002, Evison and Rhoades, 2004). The precursory scale increase is an increase in the magnitude and rate of occurrence of minor earthquakes, which has been observed to typically precede major earthquakes in the long term. The predictive scaling relations are simple linear regressions of mainshock magnitude, log precursor time and log precursor area on precursor magnitude. The model assumes that the precursory scale increase phenomenon occurs at all scales in the seismogenic process. It sets aside any attempt to distinguish precursory earthquakes from others, and simply regards every earthquake as a precursor of larger events to follow it. Each earthquake makes a transient contribution to the estimate of future earthquake occurrence in its surroundings. The magnitude distribution of that contribution is at about one unit higher than the magnitude of the earthquake. The mean of the time distribution increases with the magnitude of the precursor. For example, the contribution from a magnitude 4 earthquake peaks 2-3 years after its occurrence, but extends for more than a decade. The spatial distribution of the contribution occupies an area that increases with the magnitude of the earthquake.

The long-term models
The long-term component is a weighted average of three different long-term models. These are known as NSHMBG, PPE-SSR, PPE1950, and have associated weights of 0.2, 0.5 and 0.3 respectively:

NSHMBG: The NSHMBG model is the background seismicity model from the 2010 update of the existing New Zealand national seismic hazard model (2010 NSHM; Stirling et al. 2012). It is a smoothed seismicity model with a 50-km Gaussian smoothing kernel.

PPE-SSR: The PPE-SSR model is a multiplicative hybrid of three components: a baseline model which is spatially uniform, a smoothed seismicity model known as PPE, and a gridded map of shear strain rate in New Zealand computed from the Global Positioning System (GPS) observations over the period 1991-2011. Its development is described in detail by Rhoades et al. (2017), who fitted and retrospectively tested a large number of possible multiplicative models to the New Zealand region.

PPE1950: A PPE model relies on a complete catalogue above its minimum target magnitude threshold mc for its fitting. PPE1950 is a version of the PPE model that uses the early period of the earthquake catalogue up to 1950. The early catalogue from 1840-1950 contains much information on the locations of large earthquakes in New Zealand that is not included in our models fitted with mc = 4.95. 5.0. The PPE1950 model uses the New Zealand earthquake catalogue up to 1950 with mc = 5.95. For use at lower magnitudes, it is extrapolated down to magnitude 5.0 based on the Gutenberg-Richter frequency-magnitude relation. The PPE1950 model was included in the long-term component of the Canterbury 50-year hybrid forecasting model (Gerstenberger et al. 2014, 2016; Rhoades et al. 2016).

Prior to using the Kaikōura background model described above, the earthquake forecasts for the Kaikōura sequence were done with the PPE model only. The PPE (Proximity to Past Earthquakes) model is closely based on a model proposed by Jackson and Y. Kagan (1999). It is a smoothed seismicity model, in which earthquakes are forecast to occur close to the epicentres of past earthquakes above the magnitude threshold mc for the test region. The PPE model was described in detail by Rhoades and Evison (2004).  In this case we develop the model based on a non-declustered earthquake catalogue.

The PPE model has three parameters: a normalization parameter a, a spatial smoothing distance d, and a spatially uniform background rate s. The optimal values of the parameters ad, and s over the period 1965-2016, using earthquakes in the data-collection region from 1951 – 2016 the parameters are as follows: 

Parameter                  Value
a                                  0.346
d                                  3.02
g                                  8.095×10-12

Modified Mercalli Intensity Forecasts
Throughout the forecasting period we have calculated the probability of experiencing Modified Mercalli Intensity (MMI) 7 or greater. We have used the Dowrick and Rhoades (2005) intensity prediction equation (IPE) to calculate MMI. In this calculation 2000 earthquake catalogues were simulated based on the STEP 30-day earthquake rate forecast; for each simulated earthquake a MMI map was generated and that was sampled from the uncertainty distribution in the IPE. A preferred structural orientation, based on the New Zealand NSHM is applied in Dowrick and Rhoades (2005) and this orientation trends SW-NE in the region of the Kaikoura Earthquake.

For the MMI calculations, magnitudes from the aftershock forecast in local magnitude ML  have been converted to moment magnitude MW for compatibility with ground motion prediction equations using the formula of Rhoades and Christophersen (2017).

The development of the forecast model

14 - 20 November 2016
During this period we issued earthquake number forecasts for time lengths of 1-day, 7-days and 30-days. These forecasts included the mean forecast number, the Poisson probability of occurrence for one or more event, and the 95% Poisson confidence bounds on the mean forecast number.

Each of these forecasts was calculated based on the STEP model as described above. The one day forecasts were compared to the simulated (forecast) catalogues of an ETAS implementation (Harte, 2015). 

21 November - 19 December 2016
On November 21st we ceased calculating the daily forecasts, but added a one-year forecast to the disseminated information. The 7-day and 30-day forecasts continued to be calculated with the STEP model.

One-year Forecast Model
The one-year forecast is based on a combination of the STEP model, the EEPAS Model and the PPE model. Two forecast are calculated from each model: 1) a forecast from day 0 to day 30; and 2) a forecast from day 31 to day 365.  Within each of the two time periods, the maximum forecast is chosen in each space-magnitude bin. This two-step-maximum model method is used because STEP is anticipated to better forecast the short-term aftershock clustering and EEPAS to better forecast the longer-term clustering (Rhoades and Gerstenberger, 2009). PPE is used to provide a floor level, which the forecast should not drop below.

At the request of multiple engineering end-users in Wellington, we have also calculated aftershock hazard spectra for the Wellington region on site class B, C and D. For these hazard estimates we have modelled earthquakes as finite faults. We used a depth distribution of earthquakes based roughly on the observed aftershock depth-distribution. We have placed 35% of the events at 5km depth, 50% of the events at 12km depth, and 15% of the events at 18km depth.  Finite fault ruptures have strike orientations of 0°, 45°, 90° and 135° assuming a uniform distribution of probability (i.e. 20% weight). The dip of finite faults are set to 90°. Active fault sources from the National Seismic Hazard Model are also included in the hazard calculations. The annual occurrence rates of active fault sources are scaled to the forecast time period (e.g. 30 days).

For these calculations we have independently applied the McVerry et al. (2006) Ground Motion Prediction Equation (GMPE) and the Bradley (2010) GMPE for the largest horizontal component which is also used in the NZS1170 loading standards. These hazard calculations were done using the OpenQuake software engine. Similar to the MMI calculations, 100,000 simulated earthquake catalogues were created that sample the uncertainty in the earthquake rate distribution (Poisson) and in the GMPEs. 

19 January 2017 – 19 October 2017
From January 2017, we included the ETAS model in the forecast model. The medium-term were the same as in the Kaikōura hazard model and the long-term model was PPE only.


Bradley, B.A., 2010. NZ-specific pseudo-spectral acceleration ground motion prediction equations based on foreign models. in Research Report No.2010-03University of Canterbury.

Burnham, K.P. & Anderson, D.R., 2002. Model Selection and Multimodel inferences A Practical Information-Theoretic Approach, edn, Vol., pp. Pages, Springer, New York.

Console, R., Rhoades, D.A., Murru, M., Evison, F.F., Papadimitriou, E.E. & Karakostas, V.G., 2006. Comparative performance of time-variant, long-range and short-range forecasting models on the earthquake catalogue of Greece, J. Geophys. Res., 111(B9).

Dowrick, D.J. & Rhoades, D.A., 2005. Revised models for attenuation of Modified Mercalli intensity in New Zealand earthquakes, Bulletin of the New Zealand Society for Earthquake Engineering, 38, 185-214.

Evison, F. & Rhoades, D., 2002. Precursory scale increase and long-term seismogenesis in California and Northern Mexico, Annals of Geophysics, 45, 479-495.

Evison, F.F. & Rhoades, D.A., 2004. Demarcation and scaling of long-term seismogenesis, Pageoph, 161, 21-45.

Gerstenberger, M.C., McVerry, G.H., Rhoades, D.A. & Stirling, M., 2014. Seismic Hazard Modelling for the Recovery of Christchurch, New Zealand, Earthquake Spectra, 30, 17-29.

Gerstenberger, M.C., Rhoades, D.A. & McVerry, G.H., 2016. A Hybrid Time‐Dependent Probabilistic Seismic‐Hazard Model for Canterbury, New Zealand, Seism. Res. Lett., 87, 1311-1318.

Gerstenberger, M.C., Wiemer, S. & Jones, L., 2004. Real-time Forecasts of Tomorrow's Earthquakes in California: a New Mapping Tool, United States Geological Survey Open-File Report,, 2004-1390.

Gerstenberger, M.C., Wiemer, S., Jones, L.M. & Reasenberg, P.A., 2005. Real-time forecasts of tomorrow's earthquakes in California, Nature, 435, 328-331.

Gutenberg, R. & Richter, C.F., 1944. Frequency of earthquakes in California, Bull. Seismol. Soc. Am., 34, 185-188.

Harte, D.S., 2013. Bias in fitting the ETAS model: a case study based on New Zealand seismicity, Geophys. J. Int., 192, 390-412.

Harte, D.S., 2015. Log-likelihood of earthquake models: evaluation of models and forecasts, Geophys. J. Int., 201, 711-723.

Jackson, D. & Y. Kagan, Y., 1999. Testable Earthquake Forecasts for 1999, Seism. Res. Lett.,, 70, 393-403.

McVerry, G.H., Zhao, J.X., Abrahamson, N.A. & Somerville, P.G., 2006. Response spectral attenuation relations for crustal and subduction zone earthquakes, Bulletin of the New Zealand National Society for Earthquake Engineering, 39, 1-58.

Ogata, Y., 1988. Statistical models for earthquake occurrences and residual analysis for point processes, Journal of American Statistical Association, 83, 9-27.

Omori, F., 1894. On aftershocks, Report of Imperial Earthquake Investigation Committee, 2, 103-109.

Pollock, D., 2007. Aspects of short-term and long-term seismic hazard assessment in New Zealand, Master Thesis Master Thesis, ETH Zurich.

Reasenberg, P.A. & Jones, L.M., 1989. Earthquake hazard after a mainshock in California, Science, 243, 1173-1176.

Rhoades, D.A., 2007. Application of the EEPAS Model to Forecasting Earthquakes of Moderate Magnitude in Southern California, Seism. Res. Lett., 78, 110-115.

Rhoades, D.A. & Christophersen, A., 2017. Magnitude conversion of earthquake rate forecasts, Bulletin of the Seimological Society of America.

Rhoades, D.A., Christophersen, A. & Gerstenberger, M.C., 2017. Multiplicative earthquake likelihood models incorporating strain rates, Geophys. J. Int., 208, 1764-1774.

Rhoades, D.A. & Evison, F.F., 2004. Long-range earthquake forecasting with every earthquake a precursor according to scale, Pageoph, 161, 47-72.

Rhoades, D.A. & Evison, F.F., 2005. Test of the EEPAS Forecasting Model on the Japan Earthquake Catalogue, Pageoph, 162, 1271-1290.

Rhoades, D.A. & Evison, F.F., 2006. The EEPAS forecasting model and the probability of moderate-to-large earthquakes in central Japan, Tectonophys., 417, 119-130.

Rhoades, D.A. & Gerstenberger, M.C., 2009. Mixture models for improved short-term earthquake forecasting, Bull. Seismol. Soc. Am., 99, 636-646.

Rhoades, D.A., Liukis, M., Christophersen, A. & Gerstenberger, M.C., 2016. Retrospective tests of hybrid operational earthquake forecasting models for Canterbury, Geophys. J. Int., 204, 440-456.

Schorlemmer, D. & Gerstenberger, M., 2007. RELM Testing Centre, Seism. Res. Lett.,, 78, 30-36.

Stirling, M., McVerry, G., Gerstenberger, M., Litchfield, N. & Van Dissen, R., Berryman, K., Barnes, P., Wallace, L., Villamor, P., Langridge, R., Lamarche, G., Nodder, S., Reyners, M., Bradley, B., Rhoades, D., Smith, W., Nicol, A., Pettinga, J., Clark, K., Jacobs, K., 2012. National seismic hazard model for New Zealand: 2010 update, Bull. Seismol. Soc. Am., 102, 1514-1542.

Utsu, T., Ogata, Y. & Matsu'ura, R.S., 1995. The centenary of the Omori formula for a decay law of aftershock activity, Journal of the Physics of the Earth, 43, 1-33.

Wells, D.L. & Coppersmith, K.J., 1994. New Empirical Relationships among Magnitude, Rupture Length, Rupture Width, Rupture Area, and Surface Displacement, Bull. Seismol. Soc. Am., 84, 974-1002.