**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 governmen<sup>t</sup> agencies and the insurance industry [1–3].

Empirical observations of precursory seismicity patterns have an important role in aiding the development of earthquake forecasting models [4–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 *TP* 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 *Mm*,

**Citation:** Rastin, S.J.; Rhoades, D.A.; Christophersen, A. Space-Time Trade-Off of Precursory Seismicity in New Zealand and California Revealed by a Medium-Term Earthquake Forecasting Model. *Appl. Sci.* **2021**, *11*, 10215. https://doi.org/ 10.3390/app112110215

Academic Editor: Robert Shcherbakov

Received: 3 October 2021 Accepted: 20 October 2021 Published: 31 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

precursory magnitude *MP*, precursor time *TP* and precursory area *AP* (Figure 1c), within which the precursors, major earthquake and aftershocks all occurred.

**Figure 1.** Identification of Ψ phenomenon for the August 2014 *M*6.0 South Napa, California earthquake. (**a**) The precursory area *AP* (dashed rectangle) with the epicenters of the precursory seismicity, mainshocks and aftershocks. (**b**) Magnitude versus time of prior and precursory earthquakes. Dashed lines show the precursory increase in magnitude level. *Mm* is the main shock magnitude, and *MP* is the precursor magnitude. (**c**) Changes in the cumulative magnitude anomaly (Cumag) over time; see [12] for the definition. Dashed lines show the precursory increase in the seismicity rate in 1998. The protractor translates the Cumag slope into the seismicity rate in magnitude units per year (M.U. yr<sup>−</sup>1). *TP* is the precursor time.

> From the combined identifications of Ψ from four well-catalogued regions, it was found that *Mm*, *MP*, *TP* and *AP* were all positively correlated [12]. In particular, three scaling relations (Figure 2) allowed *Mm*, *TP* and *AP* to be predicted from *MP*, 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) mediumterm earthquake forecasting model [13].

> Although *Mm*, *MP*, *TP* and *AP* were all positively correlated, *AP* and *TP* were less correlated than the other pairs of variables, as shown by the low value of the coefficient of determination *R*<sup>2</sup> in Figure 3a compared with those in Figure 2a–c. In Figure 2, we highlighted the earthquakes for which *AP* was high and *TP* 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 *TP* and *AP* was highly correlated with *Mm*, as seen in Figure 3b, with *R*<sup>2</sup> being higher than any of those values in Figure 2. These features pointed to a trade-off between *AP* and *TP*. 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–16], or could it be a statistical side-effect? For example, in this case, if log *TP* and log *AP* were independently correlated with *Mm*, then their sum would be correlated even better, such as in Figure 3b.

**Figure 2.** Predictive scaling relations and 95% tolerance limits derived from 47 examples of ψ from four regional earthquake catalogues, taken after [12]. (**a**) Mainshock magnitude *Mm* versus precursor magnitude *MP* (coefficient of determination *R*<sup>2</sup> = 71%). (**b**) Precursor time *TP* versus *MP* (*R*<sup>2</sup> = 65%). (**c**) Precursory area *AP* versus *MP* (*R*<sup>2</sup> = 48%). Enlarged and colored points are for 1990 Weber (blue square), 1968 Puysegur Bank (red square), 1969 E. Hokkaido (blue circle), 2000 W. Tottori (red circle), 1948 Karpathos (blue triangle), 1983 Kefallonia (red triangle), 1966 Colorado D. (blue cross) and 1980 S. Cascadia (red cross).

**Figure 3.** Scaling relations and 95% tolerance limits derived from 47 examples of Ψ from four regional earthquake catalogues, taken after [12]. (**a**) Precursor time *TP* versus precursory area *AP* (*R*<sup>2</sup> = 34%). (**b**) Product of *AP* and *TP* versus mainshock magnitude *Mm* (*R*<sup>2</sup> = 75%). Symbols are enlarged and colored as in Figure 2.

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 *TP* and *AP* 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 *mc* and locations (*<sup>x</sup>*,*y*) within a region of surveillance *R*, the total rate density takes the following form:

$$
\lambda(t, m, \mathbf{x}, y) = \mu \lambda\_0(t, m, \mathbf{x}, y) + \sum\_{t\_i \ge t\_0, m\_i \ge m\_0} \eta(m\_i) \lambda\_i(t, m\_i \ge y) \tag{1}
$$

where *μ* is an adjustable mixing parameter representing the proportion of the forecast contributed by the background model component; *λ*0 is the rate density of the background model; *η* is a normalizing function and *ti* and *mi* are the origin time and magnitude of the *i*th 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 *i*th earthquake to the rate density is given by

$$
\lambda\_i(t, m, x, y) = w\_i f(t | t\_i, m\_i) g(m | m\_i) h(x, y | \mathbf{x}\_i, y\_i, m\_i), \tag{2}
$$

in which *wi* 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 *mi* of the contributing earthquake. Following the notation in [20], the magnitude density *g* is a normal density of the following form:

$$\lg(m|m\_i) = \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],\tag{3}$$

in which *aM*, *bM* and *σM* are adjustable parameters. The time density *f* is a lognormal density of the following form:

$$f(t|t\_i, m\_i) = \frac{H(t - t\_i)}{(t - t\_i)\sigma\_T \ln(10)\sqrt{2\pi}} \exp\left[-\frac{1}{2} \left(\frac{\log(t - t\_i) - a\_T - b\_T m\_i}{\sigma\_T}\right)^2\right],\tag{4}$$

in which *H*(*s*) = 1 if *s* > 0 and is 0 otherwise and *aT*, *bT* and *σT* are adjustable parameters. If all other parameters are fixed, the mean of the time distribution is proportional to 10*aT* .

Therefore, 10*aT* can be regarded as a temporal scaling factor. The location density *h* is a bivariate normal density of the following form:

$$h(\mathbf{x}, y | \mathbf{x}\_i, y\_i, m\_i) = \frac{1}{2\pi\sigma\_A^2 10^{b\_A m\_i}} \exp\left[ -\frac{(\mathbf{x} - \mathbf{x}\_i)^2 + (y - y\_i)^2}{2\sigma\_A^2 10^{b\_A m\_i}} \right],\tag{5}$$

in which *σA* and *bA* are adjustable parameters. If all other parameters are fixed, the area occupied by the location distribution is proportional to *<sup>σ</sup>*2*A*. Therefore, *σ*2*A* 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 (*ts*, *tf*) and a magnitude range (*mc*, *m*max). If the target earthquakes have coordinates *tij*, *mij*, *xij*, *yij*, *j* = 1, ··· , *<sup>N</sup>* , the space–time point process log likelihood [21,22] is given by

$$\ln L = \sum\_{j=1}^{N} \ln \lambda \left( t\_{\text{ij}}, m\_{\text{ij}}, \, x\_{\text{ij}}, y\_{\text{ij}} \right) + \iint\limits\_{R} \int\_{m\_{\text{c}}}^{m\_{\text{max}}} \int\_{t\_{\text{s}}}^{t\_{f}} \lambda \left( t, m, x, y \right) dt dx dy. \tag{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 *<sup>I</sup>*(*<sup>X</sup>*,*<sup>Y</sup>*) of one model *X* over another model *Y* is given by

$$I(X,Y) = (\ln L\_X - \ln L\_Y) / N. \tag{7}$$

where ln *LX* is the log-likelihood of model *X* and *N* is the number of target earthquakes [24].
