Space–Time Trade-Off of Precursory Seismicity in New Zealand and California Revealed by a Medium-Term Earthquake Forecasting Model

Abstract

The ‘Every Earthquake a Precursor According to Scale’ (EEPAS) medium-term earthquake forecasting model is based on the precursory scale increase (Ψ) phenomenon and associated scaling relations, in which the precursor magnitude $M_P$ is predictive of the mainshock magnitude $M_m$, precursor time $T_P$ and precursory area $A_P$. In early studies of Ψ, a relatively low correlation between $T_P$ and $A_P$ suggested the possibility of a trade-off between time and area as a second-order effect. Here, we investigate the trade-off by means of the EEPAS model. Existing versions of EEPAS in New Zealand and California forecast target earthquakes of magnitudes M > 4.95 from input catalogues with M > 2.95. We systematically vary one parameter each from the EEPAS distributions for time and location, thereby varying the temporal and spatial scales of these distributions by two orders of magnitude. As one of these parameters is varied, the other is refitted to a 20-year period of each catalogue. The resulting curves of the temporal scaling factor against the spatial scaling factor are consistent with an even trade-off between time and area, given the limited temporal and spatial extent of the input catalogue. Hybrid models are formed by mixing several EEPAS models, with parameter sets chosen from points on the trade-off line. These are tested against the original fitted EEPAS models on a subsequent period of the New Zealand catalogue. The resulting information gains suggest that the space–time trade-off can be exploited to improve forecasting.

1. Introduction

Medium-term earthquake forecasting with time-varying models is becoming increasingly important for operational earthquake forecasting and the development of seismic hazard models. For example, the New Zealand medium-term forecast model has end users interested in time-varying earthquake hazards and risk, including the land use planning and building sector, central and local government agencies and the insurance industry [1,2,3].

Empirical observations of precursory seismicity patterns have an important role in aiding the development of earthquake forecasting models [4,5,6,7,8,9,10]. One such pattern is the precursory scale increase (Ψ) phenomenon, which is an increase in the magnitude and rate of occurrence of small earthquakes [11,12]. Individual examples of Ψ were identified by examining the seismicity in arbitrary frames of space and time preceding the occurrence of a major earthquake, such as in Figure 1 for the 2014 Napa, California earthquake. A magnitude versus time plot (Figure 1b) and a cumulative magnitude anomaly (Cumag) plot (Figure 1c) were used to identify the onset of precursory seismicity [12]. The onset is marked by the minimum of the Cumag plot. The precursor time $T_P$ is then found as the time between the onset and origin time of the major earthquake. The space–time frame was chosen to informally maximize the increase in magnitude and seismicity rate at the time of the onset. Each example of Ψ provided a value of the mainshock magnitude $M_m$, precursory magnitude $M_P$, precursor time $T_P$ and precursory area $A_P$ (Figure 1c), within which the precursors, major earthquake and aftershocks all occurred.

From the combined identifications of Ψ from four well-catalogued regions, it was found that $M_m$, $M_P$, $T_P$ and $A_P$ were all positively correlated [12]. In particular, three scaling relations (Figure 2) allowed $M_m$, $T_P$ and $A_P$ to be predicted from $M_P$, defined as the average magnitude of the three largest precursory earthquakes. These three predictive relations became the basis for the ‘Every Earthquake a Precursor According to Scale’ (EEPAS) medium-term earthquake forecasting model [13].

Although $M_m$, $M_P$, $T_P$ and $A_P$ were all positively correlated, $A_P$ and $T_P$ were less correlated than the other pairs of variables, as shown by the low value of the coefficient of determination R² in Figure 3a compared with those in Figure 2a–c. In Figure 2, we highlighted the earthquakes for which $A_P$ was high and $T_P$ was low or vice versa relative to the fitted relations, a condition that is not uncommon. The same earthquakes are highlighted in Figure 3. Remarkably, the product of $T_P$ and $A_P$ was highly correlated with $M_m$, as seen in Figure 3b, with R² being higher than any of those values in Figure 2. These features pointed to a trade-off between $A_P$ and $T_P$. However, the origin of this trade-off was not clear. Could it have a physical origin related to, say, the tectonic setting or seismicity rate [14,15,16], or could it be a statistical side-effect? For example, in this case, if $\log T_P$ and $\log A_P$ were independently correlated with $M_m$, then their sum would be correlated even better, such as in Figure 3b.

A study of the Ψ phenomenon in synthetic earthquake catalogues shed new light on the matter [17]. It was found that, in a synthetic catalogue generated by the earthquake simulator RSQSim [18,19], two or more equally plausible identifications of Ψ could be found for individual mainshocks. These identifications presented very different $T_P$ and $A_P$ values, consistent with a hypothetical space–time trade-off.

The evidence for the trade-off, whatever its origin, can also be strengthened through applications of the EEPAS model. One example was the EEPAS model fitted with different fixed lead times [20]. The lead time is defined as the time interval between the start of the catalogue and the origin time of a target earthquake. It was found that as the lead time increases, the mean of the EEPAS time distribution increases, and the variance of the location distribution decreases. The time and spatial scales involved varied by about a factor of two. Here, we aim to further understand the space–time trade-off by fitting the EEPAS time distribution with a fixed spatial distribution and the spatial distribution with a fixed time distribution.

In the next section, we review the defining equations of the EEPAS model and then describe the method and data for the present study. Our results show how the space–time trade-off is revealed through constrained fitting of the EEPAS model to the New Zealand and California catalogues. Finally, we indicate by way of a simple New Zealand example how the space–time trade-off might be exploited for improving the performance of medium-term earthquake forecasts.

2. EEPAS Forecasting Model

Although inspired by the Ψ predictive scaling relations (Figure 2), the EEPAS model does not involve the identification of precursory seismicity for individual major earthquakes. It treats every earthquake as a potential precursor of future larger earthquakes to follow in the medium term [13]. Depending on the magnitude, this period can range from months to decades. The model has a background component and a time-varying component. The background component is a smoothed seismicity model, with the spatial distribution depending on the proximity to the location of past earthquakes. It is, in principle, time-invariant, but it is updated at the origin time of each contributing earthquake. The time-varying component, based on the Ψ predictive relations, is obtained by summing the contributions from all past earthquakes after a starting time $t_0$ and exceeding an input magnitude threshold $m_0$. The expected earthquake occurrence rate density is a function of the time, magnitude and location denoted by λ. For times $t>t_0$, magnitudes m exceeding a target threshold $m_c$ and locations (x,y) within a region of surveillance R, the total rate density takes the following form:

$$\lambda \left(t,m,x,y\right)=\mu {\lambda }_0\left(t,m,x,y\right)+{\displaystyle \sum }_{t_i\ge t_0,m_i\ge m_0}\eta \left(m_i\right){\lambda }_i\left(t,m,x,y\right)$$

(1)

where μ is an adjustable mixing parameter representing the proportion of the forecast contributed by the background model component; ${\lambda }_0$ is the rate density of the background model; $\eta$ is a normalizing function and $t_i$ and $m_i$ are the origin time and magnitude of the ith contributing earthquake, respectively. The contributing earthquakes come from a larger search region, which needs to be big enough to include all earthquakes that might affect the rate density within R. The contribution from the ith earthquake to the rate density is given by

$${\lambda }_i\left(t,m,x,y\right)=w_{i}f\left(t|t_i,m_i\right)g\left(m|m_i\right)h\left(x,y|x_i,y_i,m_i\right),$$

(2)

in which w_i is a weighting factor and f, g and h are the densities of probability distributions which are based on the Ψ predictive scaling relations (Figure 2). These distributions depend on the magnitude m_i of the contributing earthquake. Following the notation in [20], the magnitude density g is a normal density of the following form:

$$g\left(m|m_i\right)=\frac{1}{{\sigma }_M\sqrt{2\pi }}\exp \left[-\frac{1}{2}{\left(\frac{m-a_M-b_M{m}_i}{{\sigma }_M}\right)}^2\right],$$

(3)

in which a_M, b_M and σ_M are adjustable parameters. The time density f is a lognormal density of the following form:

$$f\left(t|t_i,m_i\right)=\frac{H\left(t-t_i\right)}{\left(t-t_i\right){\sigma }_T\ln \left(10\right)\sqrt{2\pi }}\exp \left[-\frac{1}{2}{\left(\frac{\log \left(t-t_i\right)-a_T-b_T{m}_i}{{\sigma }_T}\right)}^2\right],$$

(4)

in which H(s) = 1 if s > 0 and is 0 otherwise and $a_T$, $b_T$ and ${\sigma }_T$ are adjustable parameters. If all other parameters are fixed, the mean of the time distribution is proportional to ${10}^{a_T}$. Therefore,${10}^{a_T}$ can be regarded as a temporal scaling factor. The location density h is a bivariate normal density of the following form:

$$h\left(x,y|x_i,y_i,m_i\right)=\frac{1}{2\pi {\sigma }_A^2{10}^{b_A{m}_i}}\exp \left[-\frac{{\left(x-x_i\right)}^2+{\left(y-y_i\right)}^2}{2{\sigma }_A^2{10}^{b_A{m}_i}}\right],$$

(5)

in which σ_A and b_A are adjustable parameters. If all other parameters are fixed, the area occupied by the location distribution is proportional to ${\sigma }_A^2$. Therefore, ${\sigma }_A^2$ can be regarded as a spatial scaling factor.

The adjustable parameters are fitted to maximize the log likelihood of the target earthquakes in the region of surveillance over a fitting period ($t_s$, $t_f$) and a magnitude range ($m_c$, $m_{\max }$). If the target earthquakes have coordinates $\left[\left(t_{ij},m_{ij}, x_{ij},y_{ij}\right), j=1, \cdots ,N\right]$, the space–time point process log likelihood [21,22] is given by

$$\ln L={{\displaystyle \sum }}_{j=1}^N\ln \lambda (t_{ij},m_{ij}, x_{ij},y_{ij})+{{\displaystyle \iint }}_R{{\displaystyle \int }}_{m_c}^{m_{\max }}{{\displaystyle \int }}_{t_s}^{t_f}\lambda \left(t,m,x,y\right)dtdmdxdy.$$

(6)

Information gain statistics compare the performance of different models with the same data [23]. For models with the same number of fitted parameters or for testing pre-fitted models on an independent data set, the information gain per earthquake $I\left(X,Y\right)$ of one model $X$ over another model $Y$ is given by

$$I\left(X,Y\right)=(\ln L_X-\ln L_Y)/N.$$

(7)

where $\ln L_X$ is the log-likelihood of model $X$ and N is the number of target earthquakes [24].

3. Method

The EEPAS model is usually fitted with a time lag to prevent any influence on the parameters from short-term clustering. Here, a time lag of 50 days was applied. This means that no precursory earthquake contributed to the time-varying rate density until 50 days after its occurrence.

Two different weighting strategies are commonly adopted in applications of the EEPAS model: equal weighting and down-weighting of aftershocks. For down-weighting of aftershocks, the weight assigned to each earthquake depends on the ratio of the rate density of the background model to the rate density of an epidemic-type aftershock model at the time, magnitude and location of its occurrence. For details, see [13]. The down-weighted aftershocks strategy is preferable for investigating the space–time trade-off because it better respects the hierarchical nature of seismicity, as seen in aftershock occurrence as well as precursory seismicity [25,26,27].

We considered two versions of the down-weighted aftershocks EEPAS model, which we labeled EEPAS_1F. The models were called NZ EEPAS_1F and California EEPAS_1F. The model parameters are listed in

Table 1. EEPAS_1F model parameters for New Zealand (NZ) and California.

Parameter	Details	NZ EEPAS-1F	California EEPAS-1F
$m_0$	Minimum precursor magnitude	2.95 *	2.95 *
$m_c$	Minimum target magnitude	4.95 *	4.95 *
$m_u$	Maximum target magnitude	10.05 *,!	10.05 *!
$b_{GR}$	Gutenberg–Richter b-value	1.16 †	1.0 †
$a_M$	Equation (3)	1.10 †	1.74 †
$b_M$	Equation (3)	1.0 *	1.0 *
${\sigma }_M$	Equation (3)	0.35 †	0.60 †
$a_T$	Equation (4)	1.44 †	2.11 †
$b_T$	Equation (4)	0.43 †	0.40 †
${\sigma }_T$	Equation (4)	0.53 †	0.43 †
$b_A$	Equation (5)	0.37 †	0.35 †
${\sigma }_A$	Equation (5)	1.16 †	0.88 †
μ	Equation (8)	0.18 †	0.27 †

* Fixed. ^† Fitted. ^! Standard threshold used for CSEP models.

. The values were slightly different from the models previously tested in the New Zealand and California testing centers of the Collaboratory for the Study of Earthquake Predictability (CSEP) since 2008 and 2006, respectively [28,29,30]. The differences were due to looser constraints imposed in the fitting of µ.

The surveillance and search regions for New Zealand and California are shown in Figure 4 and Figure 5, respectively. Figure 4a and Figure 5a show the locations of earthquakes with magnitudes M > 2.95 contributing to their fitting between times $t_0$ and $t_f$. Time $t_0$ is the beginning of 1951 for New Zealand and 1932 for California, while $t_f$ is the end of 2006 for New Zealand and 2005 for California. Figure 4b and Figure 5b show the locations of the target earthquakes with M > 4.95 between times $t_s$ and $t_f$, where $t_s$ is the beginning of 1987 for New Zealand and 1986 for California.

To investigate the space–time trade-off, we varied the EEPAS model parameters in a controlled way. Starting with the parameter sets listed in Table 1, we separately changed the EEPAS_1F parameters ${\sigma }_A$ and $a_T$ while the other parameters, except the mixing parameter μ, remained fixed at their previously fitted values. We changed ${\sigma }_A$ in seven steps in either direction away from its optimal value (

Table 2. Controlled values of ${\sigma }_A$ in EEPAS_1F model for New Zealand (NZ) and California.

NZ EEPAS-1F	California EEPAS-1F
0.34	0.26
0.41	0.31
0.49	0.37
0.58	0.44
0.69	0.53
0.82	0.97
0.97	0.74
1.16 †	0.88 †
1.38	1.05
1.64	1.25
1.95	1.49
2.32	1.77
2.75	2.10
3.27	2.50
3.89	2.98

^† Fitted.

) and obtained the corresponding values of the temporal scaling factor ${\sigma }_A^2$. Subsequently, we changed the $a_T$ values in a similar manner (

Table 3. Controlled values of $a_T$ in EEPAS_1F models for NZ and California.

NZ EEPAS-1F	California EEPAS-1F
2.49	3.16
2.34	3.01
2.19	2.86
2.04	2.71
1.89	2.56
1.74	2.41
1.59	2.26
1.44 †	2.11 †
1.29	1.96
1.14	1.81
0.99	1.66
0.84	1.51
0.69	1.36
0.54	1.21
0.39	1.06

^† Fitted.

) and obtained the corresponding values of the temporal scaling factor ${10}^{aT}$. Over seven steps, each of the controlled scaling factors varied by an order of magnitude on either side of the optimal fit. For each controlled value of $a_T$ or ${\sigma }_A$, two free parameters, μ and either ${\sigma }_A$ or $a_T$, were refitted to maximize the likelihood of target earthquakes in the region of surveillance over time period ($t_s$, $t_f$).

4. Results

The likelihood of the refitted models declined with each step change in the controlled parameter away from its optimal value, as shown for New Zealand in Figure 6. The results for California were similar. The log-likelihood of the refitted model is plotted against the controlled spatial scaling factor in Figure 6a and against the temporal scaling factor in Figure 6b. An order of magnitude change in each scaling factor induced a modest reduction in the log-likelihood. The maximum reduction of about 34 units corresponded to an information loss per earthquake of about 0.2 relative to the overall optimal fit.

The refitted mixing parameter μ tended to increase as the controlled parameter shifted further away from its optimal value, as shown for New Zealand in Figure 7. Again, the results were similar for California. The variation of μ with the spatial scaling factor is shown in Figure 7a and against the temporal scaling factor in Figure 7b. The values of μ increased from about 0.15 at the optimal fit to greater than 0.5 when the temporal or spatial scaling factors were changed by an order of magnitude. The μ value represents the proportional contribution of the background model to the total EEPAS model rate density. Higher μ values thus indicate a greater contribution of the background component and a smaller contribution of the time-varying component. In other words, higher μ values indicate that there were fewer target earthquakes with precursors matching the changed spatial and temporal distributions.

As the controlled parameter was changed, the refitted values of the other parameters changed in a way that was consistent with the notion of a space–time trade-off. The results are shown for New Zealand in Figure 8a and for California in Figure 8b.

In each plot, the pairs of scaling factors resulting from controlling ${\sigma }_A$ are shown as blue triangles, and those resulting from controlling $a_T$ are shown as black squares. The temporal scaling factor decreased as the controlled spatial scaling factor increased, and the spatial scaling factor decreased as the controlled temporal scaling factor increased. However, the curves had different slopes depending on whether ${\sigma }_A$ or $a_T$ was the controlled variable. An even trade-off line with a slope of −1 is drawn through the intersection of the two curves (straight blue line in Figure 8a,b). Its slope lies between the average slopes of the two controlled fitting curves.

5. Discussion

As seen in Figure 8, the controlled fits produced two curves which did not lie on the even trade-off line but instead had higher or lower slopes. This result can be explained by the limitations on the length of the catalogue and the size of the search region. The fitted parameters could only adjust to the precursors that were contained in the catalogue and not to those that were screened out by such limitations. We now consider in detail the trend of the fitted ${\sigma }_A$ value away from the even trade-off line for the controlled values of $a_T$. The trend of the other curve can be explained similarly.

As $a_T$ was stepped down to lower values (i.e., the time scale was shortened), fitting the trade-off required earthquakes at increasingly longer distances from the target earthquakes. However, at longer distances, more precursory events were screened out by the spatial limitation on the input catalogue. The precursors of the largest earthquakes in the target magnitude range would be most affected by the spatial limitation because they had larger precursory areas (Figure 2). The spatial limitation at small $a_T$ values forced the fitted values of ${\sigma }_A$ to increasingly fall below the even trade-off line. On the other hand, as $a_T$ was stepped up to higher values, the precursory time scale became longer and exceeded the available lead time. This temporal limitation most affected the largest earthquakes in the target magnitude range, which had the longest precursor times (Figure 2). Thus, more and more precursory earthquakes on the specified time scale were screened out by the limited time span of the catalogue. The remaining precursors for fitting ${\sigma }_A$ would be those at the lower end of the time distribution. Because of the space–time trade-off, these remaining precursors tended to be at longer distances than the screened-out events. This forced the fitted ${\sigma }_A$ to increasingly exceed the even trade-off line.

The space–time trade-off in the EEPAS model shows that as the mean of $f\left(t|m\right)$ in Equation (4) increased, the area of the fitted $h\left(x,y|m\right)$ in Equation (5) decreased and vice versa. This phenomenon can be interpreted in terms of the predictive scaling relations on which the EEPAS model is based (Figure 2). Figure 2 shows that $T_P$ and $A_P$ both increased with the precursory earthquake magnitude $M_P$. Similarly, in the EEPAS model, the mean of the time distribution $f\left(t|m\right)$ and the area of the location distribution $h\left(x,y|m\right)$ both increased with m. Now, the space–time trade-off observed in the EEPAS model can be interpreted in terms of the space–time distribution of precursors to an individual major earthquake; that is, the earliest precursors tend to occur very close to the source, and the later precursors to occupy a wide area around the source. This interpretation only applies to precursors occurring more than 50 days before the mainshock because of the time lag applied for EEPAS model fitting here.

The existence of this trade-off raises the question of how it can be exploited to improve the performance of the EEPAS model. The EEPAS model treats the time and location as independent variables, but the trade-off implies that they are correlated. We will illustrate how to improve forecasting by forming hybrid models. The hybrid models are mixtures of three EEPAS models with the values of $a_T$ and ${\sigma }_A$ chosen from points on the even trade-off line with a slope of −1. We constructed two models, Hybrid_1F and Hybrid_1R, starting from two different EEPAS models: EEPAS_1F and EEPAS_1R, respectively. EEPAS_1R was similar to EEPAS_1F in nearly all aspects, apart from having fewer optimized parameters. Its fixed and optimized parameters are given in

Table 4. EEPAS_1R model parameters for New Zealand.

Parameter	Value
$m_0$	2.95 *
$m_c$	4.95 *
$m_u$	10.05 *
$b_{GR}$	1.16 †
$a_M$	1.00 †
$b_M$	1.0 *
${\sigma }_M$	0.32 *
$a_T$	1.40 †
$b_T$	0.40 *
${\sigma }_T$	0.23 *
$b_A$	0.35 *
${\sigma }_A$	1.74 †
μ	0.24 †

* Fixed. ^† Fitted.

. An important difference between the two models was that EEPAS_1F (Table 1) had a larger value of ${\sigma }_T$ than EEPAS_1R (Table 4). The parameter ${\sigma }_T$ was optimized in the fitting of EEPAS_1F but not in EEPAS_1R. In prospective testing over 10 years in the New Zealand CSEP testing center, EEPAS_1F significantly outperformed EEPAS_1R [29].

To construct the hybrid models, we replaced the time-varying component of each model’s rate density with the average rate density of the three models, with the values of $a_T$ and ${\sigma }_A$ chosen from the trade-off line. The three models were the original one and two others formed by an arbitrary increase and decrease in $a_T$ of $∆=0.5$. For an increase in $∆$ in $a_T$, the corresponding value of ${\sigma }_A$ on the trade-off line was found by multiplying the original ${\sigma }_A$ by ${10}^{-0.5\mathsf{\Delta }}$. The other parameters, including μ and ${\sigma }_T$, remained unchanged at their values in Table 1 and Table 4. Using information gain statistics, we compared the performance of the EEPAS_1F, EEPAS_0F, Hybrid_1F and Hybrid_1R models. For this, we used a test period from 2007 to 2017, during which there were 259 target earthquakes with magnitudes M > 4.95. Hybrid_1R outperformed all the other models, and EEPAS_1R was the weakest model (Figure 9). Figure 9a shows the information gain of EEPAS_1F, and Figure 9b shows that of Hybrid_1R over the other models. Both hybrid models and EEPAS_1F outperformed EEPAS_1R with 95% confidence according to the T-test [24].

This simple example of hybrid formation, even without fitting additional parameters, suggests that it might be possible to use the space–time trade-off to improve forecasting. However, much more work needs to be done to construct a formal method for optimal inclusion of the trade-off in the fitting of the EEPAS model. The temporal and spatial limitations of the catalogue are clearly among the issues to be considered. The spatial limitations can be resolved if a global catalogue is used, but then a higher threshold magnitude of completeness would apply. That in turn imposes further limitations. Additionally, there is evidence that the precursor time distribution is dependent on the strain rate in the vicinity of a target earthquake [17]. This dependence would have to be included in a global model. Temporal limitations can also be partly resolved by introducing a fixed lead time for all target earthquakes and then compensating for the lead time using the method described in [20].

6. Conclusions

A space–time trade-off of precursory seismicity has been investigated by repeated refitting of the EEPAS earthquake forecasting model to the catalogues of New Zealand and California. In a sequence of controlled fits, the temporal scaling parameter was constrained to vary in steps ranging over two orders of magnitude with the spatial scaling parameter before being refitted, and vice versa. The two resulting curves of the temporal scaling factor against the spatial scaling factor differed depending on which parameter was controlled and which was fitted. However, both curves were consistent with an even trade-off between space and time once the temporal and spatial limits of the contributing earthquake data were considered. As the controlled parameter deviated further from its optimal value, the likelihood of the refitted model decreased. In addition, the refitted model had an increasingly large background component and a diminishing time-varying component.

The trade-off implies that the earliest precursors to a major earthquake tend to occur very close to its source and that the later precursors occupy a wide area around the source. A simple example in which hybrid forecasts were created by mixing several EEPAS models with parameters chosen from the trade-off line suggests that it should be possible to exploit the trade-off for improved forecasting. However, more research is needed to develop a formal method for routinely incorporating the space–time trade-off into medium-term earthquake forecasts.