Electronic Supplement to
Prospective Earthquake Forecasts at the Himalayan Front after the 25 April 2015 M 7.8 Gorkha Mainshock

by Margarita Segou and Tom Parsons

Part I of this electronic supplement focuses on historic earthquake location, the implementation of epidemic-type aftershock sequence (ETAS), the conversion of stress changes in forecast earthquake rates, an example of San Francisco Bay area as related to the total stress method implementation (Fig. S1), and a short-term forecast evaluation of the predictability of the forecast models using such performance evaluation metrics as information gain (Fig. S2), together with a long-term (six-month) evaluation (Fig. S3). In part II of this electronic supplement, we present the original forecast as of 9 May 2015 (56 hours before the 12 May 2015 M 7.3 aftershock).

---

I. Historical Earthquakes and Development of Forecast Models

Historic Earthquake Location

We calculate historic earthquake locations from historical intensity observations (Bakun and Wentworth, 1997; Martin and Szelinga, 2010). The method performs a grid search for trial epicenters (5 km spacing in east and north directions) and, using an empirical intensity attenuation relation versus distance, locates the region where 95% of trial epicenters minimize root mean square misfits to observed intensity values.

Epidemic-Type Aftershock Sequence (ETAS) Implementation

Short-term statistical earthquake forecasts built on the idea that each earthquake triggers a number of offspring events, relative to its magnitude, and on empirical laws such as the Omori law decay and the Gutenberg–Richter magnitude frequency distribution. Here, we allow the ETAS (Ogata, 1988, 1998) model to use all available earthquakes to achieve optimum performance, although the Nepal catalog is complete to M ≥ 4.7. To estimate the necessary parameters, we use equations for estimating the apparent fraction n_a with respect to the real fraction of triggered events n (Sornette and Werner, 2005).

Analytically, following the ETAS model, the space–time seismicity rate λ (x,y,t) is given by $$\lambda (x,y,t)=\mu (x,y)+\sum _{i:t_{i<t}}^{}\frac{K{e}^{\alpha (M_i-M_{th})}}{{(t-t_i+c)}^p}{f}_i(x-x_i,y-y_i;M_i),$$with $$f(x,y;M)=\frac{q-1}{\pi D\left(M\right)}{(1+\frac{x^2+y^2}{D\left(M\right)})}^{-q},$$and D(M) = de^γ(M−Mmin). In our implementation, parameter values are as follows: M_max = 7.8, M_min = 4.7, m_d = 4.0, α = 0.7, b = 1.0, n = 0.67, n_a = 0.51, k = 0.2167, and c = 0.002 day. Parameter µ(x,y) is the background rate; c and p are Omori-law values governing the decay rate of aftershocks; α estimates the magnitude efficiency of an earthquake in generating its offspring; and d, γ, and q (d = 0.0071 deg², q = 1.96, γ = 0.7) are spatial fitting parameters in close agreement with parameters from similar seismotectonic environments (Parsons and Segou, 2014).

For the favorable case of lower magnitude scaling formulation with α < b = 1.0, the apparent branching ratio n_a is given by the equation $$n_a=\frac{kb\text{ln}(10)(M_{\text{max}}-m_d)}{1-{10}^{-b(M_{\text{max}}-M_{\text{min}})}},$$and the relation between n and n_a is $$n_a=n\left(\frac{{10}^{(b-a)(M_{\text{max}}-m_d)}-1}{{10}^{(b-a)(M_{\text{max}}-M_{\text{min}})}-1}\right),$$in which m_d is the detection threshold and M_min the minimum triggering threshold.

Conversion of Stress Changes to Forecast Earthquake Rates

Following the rate-and-state friction framework (Dieterich, 1994; Dieterich and Kilgore, 1996), the premainshock earthquake activity R is suppressed or enhanced by static stress changes; $R=\frac{r}{\gamma \dot{\tau }}$, in which ${\gamma }_n={\gamma }_{n-1}\text{exp}\left(\frac{-\text{ΔCFF}}{\mathrm{A\sigma }}\right)$.

The time-dependent seismicity rate R(t) is a function of state variable γ after a stress perturbation $${\gamma }_{n+1}=[{\gamma }_n-\frac{1}{\dot{\tau }}]\text{exp}\left[\frac{-\text{Δ}t\dot{\tau }}{\mathrm{A\sigma }}\right]+\frac{1}{\dot{\tau }},$$in which r is the steady-state seismicity rate, ∆CFF is the stress step, aftershock duration t_a is assumed to be 25 years, and daily forecast parameters are ασ = 0.5 (Toda and Enescu, 2011) and $\dot{\tau }=0.00239\mathrm{Mpa}/\mathrm{yr}$, based on regional thrust fault slip rates of ∼10 mm/yr (Ader et al., 2012). R(t) is related to earthquake probability over the interval Δt as $$P(t,\mathrm{\Delta }t)=1-\text{exp}[-\underset{t}{\overset{t+\mathrm{\Delta }t}{\int }}R(t)dt]=1-\text{exp}(-N(t)),$$in which N(t) is expressed as $$N(t)=r_p\{\mathrm{\Delta }t+t_{\alpha }\text{ln}[\frac{1+\left[\text{exp}\right(\frac{-\Delta \text{CFF}}{A\sigma })-1]\text{exp}\left[\frac{-\Delta t}{t_{\alpha }}\right]}{\text{exp}\left(\frac{-\Delta \text{CFF}}{A\sigma }\right)}\left]\right\},$$and $r_p=\left(\frac{-1}{\Delta t}\right)\text{ln}(1-P_c)$, in which P_c is a conditional probability.

Performance Evaluation Metrics

The modified N-test (Zechar et al., 2010) evaluates the consistency between the total number of predicted and observed events within the area of interest. This test is based on the equations, δ₁ = 1−F(N_obs−1|N_F) and δ₂ = F(N_obs|N_F), in which F(x|μ) is the right-continuous Poisson cumulative distribution function with expectation µ evaluated at x and N_F is the forecast number of events determined by the model. The quantiles δ₁ and δ₂ answer two questions under the assumption that the forecast is correct: (1) What is the probability of observing at least N_obs events? (2) What is the probability of observing at most N_obs earthquakes? These metrics share a complimentary role (${\delta }_1\simeq 1-{\delta }_2$), suggesting a forecast be rejected if either δ₁(t) < α_eff or δ₂(t) < α_eff, in which a_eff = 0.025 corresponds to the effective significance value. In Figure S2, we show the results of the N-test within 10-day time windows.

The S-test (Zechar et al., 2010) aims to compare the relative spatial performance of the forecast model using log-likelihood statistics estimated over 1000 simulations. The log likelihood L of observing ω events at a given expectation λ for a model j is defined by the logarithm of the probability p(ω|λ), expressed by $$L(\omega |{\lambda }^j)=\text{log}p(\omega \text{|}{\lambda }^j)=-{\lambda }^j+\omega \text{log}{\lambda }^j-\text{log}\omega !$$and $$L(\text{Ω|Λ})=\sum _{(i,jℰR)}\{-\lambda (i,j)+\omega (i,j)\text{log}[\lambda (i,j)]-\text{log}[\omega (i,j)!]\},$$in the case of the joint log likelihood, which represents the sum of log-likelihood values over all bins b_i. We present in Figure S3 maps of the log likelihood per spatial bin for the time intervals 0–10 and 10–20 days, respectively.

The T-test (Rhoades et al., 2011) evaluates the sample information gain per earthquake of a model A over model B defined by $$I_N(A,B)=\frac{1}{N}\sum _{i=1}^{\mathrm{N}}(X_i-{\text{Υ}}_i)-\frac{{\hat{N}}_A-{\hat{N}}_B}{N},$$in which I_N(A,B) is considered as the mean of a sample from a population with actual mean I(A,B), in which I_N(A,B) is the true information gain of model A over model B with X_i = logλ_A(i) and Y_i = logλ_B(i) as the log-likelihood value of a model A and B in the ith bin. Here, we use the statistical model as reference model due to the simplicity of this implementation. In Figure S4, we present the mean and the 95% confidence interval of the information gain per model for the time intervals used in our spatial mapping of log likelihood.

---

Figures

Figure S1. Effects of stress change on different rake and dip on the San Francisco Bay region stress changes following the 1906 earthquake. The dominant optimal plane strikes northwest, has a right-lateral rake, and dips vertically. Coulomb stress changes calculated on these planes are negative. However, because the San Andreas fault ruptured through a restraining bend, the same fault orientations with 45° dips and pure thrust rakes have positive coulomb failure stress.

Figure S2. Information gain. Mean and 95% confidence interval of the information gain of physics-based forecasts when the statistical model is taken as reference for time periods of (a) 0–10 days and (b) 10–20-days. Note the low standard deviation for the total-stress model (<0.004).

Figure S3. Long-term performance evaluation. (a) Prospective forecast update for the statistical benchmark model (cyan), parallel (green), optimal (red), and total stress (magenta) method, overlaid with observation (triangles) above M 4.7. In (b) and (c), respectively, the quantiles δ₁ and δ₂ for the N-test are shown as a function of time. Gray dotted line indicates the 0.05 significance level at which a forecast is rejected.

---

II. Prospective Forecast as of 9 May 2015 (56 Hours before the 12 May 2015 M 7.3 Aftershock)

The original rapid forecast was completed two weeks after the 25 April M 7.8 Nepal mainshock. It is followed by an e-mail confirming the submission time before the 12 May M 7.3 aftershock.

---

Other

Download: SRL-D-15-00195-R-ES-PART-B [Adobe portable document format; ∼3.2 MB]. Forecast and confirmation: impending earthquakes beneath Katmandu and the Himalayan front after the 25 April 2015 M 7.8 Nepal mainshock.

---

[ Back ]