Exceedance Probabilities for Large Earthquakes From DIY Local Earthquake Ensemble Nowcasting and Forecasting: Magnitude, Natural Time, and Calendar Time

1. Introduction

1.1. Earthquake Nowcasting

Previous papers have developed several techniques for earthquake nowcasting. The nowcasting approach is based on “natural time”, which is defined as counts of small earthquakes [1,2,3,4]. Natural time is the time scale that is relevant to the system dynamics, so earthquake nowcasting uses natural time to track the progression of the fault system through its cycle of large earthquake activity. The term “nowcasting “is used in the same sense as for weather and economic nowcasting [5,6], tracking the current state of a system in the recent past, current time, and the near future. Methods to produce earthquake nowcasts are the subject of many previous papers [5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39].

The nowcasting method is applied to a small circle of arbitrary radius surrounding a point of interest, typically a city. Computing the probability then involves a data mining phase using an expanding series or ensemble of larger rectangular regions surrounding the circle. These regions then comprise the ensemble for the method. Each member of the ensemble contains a number of earthquake “cycles” which begin and end on earthquakes larger than a target magnitude M_T. Each cycle contains a number of “small” earthquakes M_S < M_T. The similarity of the Gutenberg-Richter statistics for the local region and the surrounding ensemble regions is used to construct a time-scaling for the times in the ensemble cycles.

After time-scaling the calendar forecast times for the ensemble members, we then consider each cycle of activity in a given ensemble. Collections of time-scaled cycles of activity that have fewer small events than the number in the circle are considered to represent a set of might-have-been past histories for activity in the circle. Collections of time-scaled cycles of activity that have more small events than the number in the circle are considered to represent a set of possible future activity for the circle. Examples of the data for a circular region of radius 125 KM around Los Angeles, CA, and a selection of the ensembles is shown in Table 1.

The next step involves computation of the ROC curve for each member of the ensemble, a plot of the True Positive Rate (TPR) vs. the False Positive Rate (FPR) [40,41]. Once we are in possession of these curves, we compute the ensemble mean curve and its standard deviation. The area under the ROC curve (AUROC, or “skill”) can be interpreted as the ability of a model to discriminate between, or to correctly classify, different categories of events. In our case, an event is the occurrence of a major earthquake vs. non-occurrence [40,41]. The method then uses a simple data mining method on the longer earthquake cycles to compute the probability of future events.

This current paper is the third of a three-part series introducing methods for computing local earthquake forecasts and associated probabilities of future large earthquakes from local earthquake nowcasts. The first paper [15] detailed a method to compute local forecasts for a fixed future natural time interval, where natural time is the count of small earthquakes since the last large earthquake. The second paper [16] extended these methods to calendar time forecasts using an ensemble-of-regions approach. Here the local region of interest, a circle of radius r_C surrounding a central geographic point of interest, is enclosed by an expanding series (ensemble) of regions.

To summarize our results:

Local ensemble probabilities for large earthquakes can be computed using only counts of small earthquakes

Ensemble exceedance probabilities for magnitude and earthquake occurrence time can be computed for evaluating risk profiles

To test the basic assumption of uniformity of statistics, a nowcast transform can be defined to adjust the regional statistics to enforce uniformity within the ensembles of events

The method is then validated/benchmarked by comparison to the accepted current USGS forecasts.

2. Materials and Methods

2.1. Nowcast Function and the Ensemble Method

We begin by considering a circular region of radius 125 km (Figure 1a) surrounding a point of interest, in this case, the city of Los Angeles, CA [5,15,16]. Figure 1a also shows an ensemble of expanding rectangular regions. The last major earthquake within this circle was the magnitude M6.7 Northridge, CA earthquake on January 17, 1994. We wish to compute the probability of the next large earthquake of target magnitude M_T ≥ 6 occurring within the circular region, whose radius was chosen mainly because it is about twice the radial dimension of the aftershocks of the Northridge earthquake.

2.1.1. Scaled Similarity

These ensemble regions are assumed to posses Gutenberg-Richter magnitude-frequency (GR) statistics that are the same or similar to the statistics in the circular region. This assumption can be stated as an assumption of “scaled similarity”, that when the temporal statistics of the members of the ensemble are appropriately scaled in time, the members can be regarded as statistical “copies” of the earthquakes in the circular region.

The specific form of the GR law that is used here is:

$$N=e^{b\left(M_{T-}{M}_S\right)}$$

(1)

where N is the number of small earthquakes of magnitude M_S for every large target earthquake M_T, and b is a parameter that is usually near b ~1 but can typically vary by ± 0.05 or so.

2.1.2. Step 1: The Nowcast Function and Earthquake Cycles

The first step in the analysis is to introduce a nowcast function Φ(n):

$${\mathsf{\Phi }\left(n\right)=1-e}^{-(n/N)}$$

(2)

where n is the number of small earthquakes of magnitude M_S since the last large earthquake of target magnitude M_T [16]. Equation (2) was chosen as a consequence of the observation [42], since validated many times by others [6], that aside from aftershock clustering, earthquakes occur randomly in time with Poisson statistics. In addition, [] has shown that equation (2) generally characterizes the interval statistics (Figure 1b) between large earthquakes M_T. Note that in the figures, we plot the nowcast value as ranging between [0%,100%] rather than [0,1].

2.1.3. Step 2: Construct Ensemble of Cycles

Using equation (2), build a series of “large earthquake cycles” in the ensembles, shown schematically in Figure 2, surrounding the circular region and centered on Los Angeles (Figure 1a). Each of the cycles begins and ends at the time of a target earthquake. So an ensemble member with N_ENS large earthquakes will have J = N_ENS -1 cycles by this definition. These cycles of activity represent the data mining set for the simple machine learning application.

The ith ensemble member will have a set of J cycles n_i,j denoted by C_i:

$$C_i\left\{n_i\right\}=\{n_{i,1},n_{i,2},n_{i,3},...{,n}_{i,J}\}$$

(3)

Here n_i,j is the number of small earthquakes in the jth cycle of the ith ensemble.

2.1.4. Step 3: Scale the Forecast Time for Each Ensemble Member

Build conditional ROC curves for each member of the ensemble using the following method, then compute the mean ROC curve from the ensemble ROC curves. The curves are conditioned on the number of small events that have occurred prior to calculation of the ROC curve.

We choose a forecast time of interest T_F,C for the next target earthquake in the circular region. The number of small events in the circle will be designated as n_C, and the corresponding rate of small earthquake occurrence within the circle is designated as R_C.

Clearly each member of the expanding ensemble will have a larger number of small events than the number of small events in the circle. The ith ensemble member will have a total number n_i of small earthquakes, and a corresponding rate of occurrence R_i.

In Step 2, we scale the forecast time interval for each member of the ensemble, T_F,i , according to [15,16]:

$$T_{F,i}=\left({\displaystyle \frac{R_C}{R_i}}\right)T_{F,C}$$

(4)

2.1.5. Step 4: Build Conditional Receiver Operating Characteristic (ROC) Curves

To build the conditional ROC curves for each of the ensembles, we adopt a set of threshold values for the time series amplitudes [10,11,12,13,14]. For each ensemble member, we classify all points on the nowcast timeseries by sweeping the threshold values over all amplitudes of the time series on the interval [0,1].

We start by selecting an arbitrary threshold value τ. Given a small event occurring at time t, we ask if nowcast value of that event is above or below the threshold. We also ask if the next target earthquake M_T occurs after time t but within the time interval t + T_F,i , where again, T_F,i is the scaled forecast time for that ensemble member. Classification is then given by:

TP, if the nowcast value is above the threshold τ, and the next target earthquake occurs within t + T_F,i

FP, if the nowcast value is above the threshold τ, and no earthquake occurs within interval t + T_F,i

FN, if the nowcast value is below the threshold τ, and the next target earthquake occurs within interval t + T_F,i

TN, if the nowcast value is below the threshold τ, and no target earthquake occurs within interval t + T_F,i

It should be emphasized that the conditional ROC diagram is conditioned on the current number of events in the circle. For example, if the current number of small earthquakes in the circle is, for example 100, then no cycles of length less than 100 are used in the classification. Examples of conditional ROC diagrams are shown in Figure 3. Cyan curves represent the ROC curves for the ensemble members, the red curve is the ensemble mean, and the dashed curves are the standard deviations.

Classifying all points on an ensemble timeseries will then produce a confusion matrix, or contingency table, composed of the quantities TP, TN, FP, FN . We should also explicitly note that these quantities are functions only of the threshold values, τ, so that we have TP(τ), TN(τ), FP(τ), FN(τ). Conditional ROC curves for each member of the ensemble are computed, and the mean and standard deviation of the curves are computed. From each confusion matrix, we compute the True Positive Rate (TPR) and the False Positive Rate (FPR):

TPR = TP/(TP + FN)

(5)

FPR = FP/(FP + TN)

As discussed, the basic assumption in this method is that the catalog statistics of the surrounding region are the same as the statistics of the circular region. Or more specifically, we make what amounts to an ergodic assumption, that time averages can be replaced by space (or ensemble) averages. In the present paper, we will introduce below a “nowcast transform” that enforces the equality of the GR statistics of the ensemble of regions with the statistics of the circular region.

2.1.6. Step 5: Compute Positive Predictive Value (PPV) from the ROC Curves

Once the conditional ROC curve is computed for a member of the ensemble, the Positive Predictive Value (PPV) is then computed:

PPV = TP/(TP + FP)

(6)

The PPV value is the probability of a future large earthquake of target magnitude M_T during the forecast interval T_F,C.

To compute the PPV values (probabilities) at the times of the small events following the last large event in the circle, we use as threshold values τ the nowcast values $\mathsf{\Phi }\left(n_C\right)$for the small events in the circle. So if there are 100 small earthquakes within the circle since the last large earthquake, we compute 100 values for $\mathsf{\Phi }\left(n_C\right)$ and use these as our threshold set [τ ($n_C$)] ≡ [Φ($n_C$)]. Note that each of these small earthquakes has a calendar time t.

Each of these small earthquakes in the circle then has a defined nowcast value, an index, and an occurrence time, all of which are associated with the small event. Again, note that the TP, FP, TN, FN are functions of the threshold only, which in this special case is the event nowcast value. Since the nowcast value is also a threshold, this association allows us to identify specific event sequence numbers and event times with a PPV value. Figure 4 shows examples of calendar forecasts of this type for both 1 year and 5 year time scales.

2.2. Conditional Exceedance Curves

Here we compute conditional exceedance curves (survivor distributions), conditioned on the number of small earthquakes that have occurred prior to computation of the probabilities.

2.2.1. Step 6: Compute Magnitude Exceedance

To build a magnitude exceedance probability, we note that all the cycles in an ensemble are terminated (and begin) with a magnitude M ≥ M_T. Thus we can build a set consisting of the next, or terminating, magnitudes M_i,j for each interval n_i,j:

$$C_i\left\{M_i\right\}=\{M_{i,1},M_{i,2},M_{i,3},...{,M}_{i,J}\}$$

(7)

Considering both $C_i\left\{M_i\right\}$ and $C_i\left\{n_i\right\}$, we then build a combined set of intervals and magnitudes:

$$C\{M\}={\bigcup }_{i,j}\{M_{i,1},M_{i,2},M_{i,3},...{,M}_{i,J}\}$$

(8)

$$C\{n\}={\bigcup }_{i,j}\{n_{i,1},n_{i,2},n_{i,3},...{,n}_{i,J}\}$$

where the $\bigcup _{i,j}\left\{\right\}$ indicates the union of all sets of ensemble intervals and magnitudes into a single large set.

Once we have these two large sets, we can construct the conditional survivor distributions (exceedance probabilities) for terminating magnitudes. Survivor distributions for terminating magnitudes are shown in Figure 5a.

Conditional exceedance curves are shown for several examples of natural times (small event counts), for 0 small events (immediately after the last large earthquake), for 448 small events (the current count in the circular region), and for the current-plus-mean-projected number of small events in the circular region. The mean projected number is the average number expected until just before the next large earthquake. This mean projected number is computed by finding the mean of all intervals larger than the current number, 448, minus the current number. Note that a light smoothing has been applied to these curves.

2.1.2. Step 7: Compute Magnitude Exceedance for 25%, 50%, 75% Probability

Figure 6 analyzes this conditional exceedance data further by computing the 25%, 50% (median), and 75% exceedance probabilities as a function of the number of small events that have occurred. Note the blue vertical line, which indicates the number of small earthquakes that have occurred to date, which is 448. As more small earthquakes occur in the future, the magnitude exceedance probabilities, particularly for the median (50%) and 25% level will begin to increase sharply.

2.1.3. Calendar Time Exceedance

To compute the calendar exceedance, or waiting time until the next large earthquake M_T, we first create the set C(ΔT_i) of combined ensemble time intervals ΔT_i corresponding to the natural time intervals in equation (8):

$$C\{∆T\}={\bigcup }_{i,j}\{{∆T}_{i,1},{∆T}_{i,2},{∆T}_{i,3},...{,∆T}_{i,J}\}$$

(9)

Next we scale the time intervals ΔT_i -> ΔT’_i using the inverse of the scaling factor in equation (4):

$${\widetilde{∆T}}_{i,j}=\left({\displaystyle \frac{R_i}{R_C}}\right)∆T_{i,j}$$

(10)

leading to:

$$C\{\widetilde{∆T}\}={\bigcup }_{i,j}\{{\widetilde{∆T}}_{i,1},{\widetilde{∆T}}_{i,2},{\widetilde{∆T}}_{i,3},...{\widetilde{∆T}}_{i,J}\}$$

(11)

Using the sets (8) and (11), we can then construct the conditional exceedance curves for the waiting times, both in natural and calendar times. Results will be discussed below.

3. Results

Now we test the results of the ensemble method, and explore ways to build ensemble forecasts using two additional ideas. These involve the following:

• We introduce a “Nowcast Transform” to test the assumption of scaled similarity of the GR statistics in the ensemble.
• We filter cycle intervals at the 95% confidence level to test an assumption that outliers may determine the statistics.

3.1. “Nowcast Transform”

We define this simple method to test the assumption that the statistics of the ensemble members produce a stable forecast result, and to examine the assumption that the statistics of the ensembles being similar to the circle are reasonable.

Using equation (2) we compute the nowcast value for a typical interval in ensemble i:

$${\mathsf{\Phi }\left(n_{i,j}\right)=1-e}^{-(n_{i,j}/N_i)}$$

(12)

where the scale N_i for ensemble member i is computed from equation (1):

$$N_i=e^{b_i\left(M_{T-}{M}_S\right)}$$

(13)

and where b_i is the b-value for ensemble member i.

To compute the transformed intervals n’_i,j , we invert equation (2), by making the equivalence $\mathsf{\Phi }\left(n_{i,j}\right)=\mathsf{\Phi }\left({n\prime }_{i,j}\right)$, and using the statistics of the circular region:

$${n\prime }_{i,j}=-N_{C}log[1-\mathsf{\Phi }(n_{i,j}\left)\right]$$

(14)

and where b_C and N_C represent the statistics of the circle :

$$N_C={10}^{b_C(M_T-M_S)}$$

(15)

This process ensures that the transformed ensembles have the same b-value, b_C, as the circle. Thus the GR interval statistics of the ensembles are the same as those of the circular region.

But in altering the values of the n_i,j intervals, we also need to consider the times at which these transformed events occur. There are two cases, one in which n_i,j > n’_i,j and another in which n’_i,j > n_i,j. We have therefore developed a simple algorithm to assign the small event times for the small events in the intervals:

For the case n_i,j > n’_i,j we choose event times in n_i,j at random, and remove a sufficient number of these. Continue until the necessary number of times have been removed.

For the case n’_i,j > n_i,j we pick two neighboring times randomly, find the mean of these two times, and insert that mean time between the two chosen times. Continue until the necessary number of times have been assigned.

A major feature of this simple algorithm is that the basic structure of the event times is preserved. Where times are densely clustered, as in aftershock intervals, proportionately more times are added or removed. Where times are sparse, as when quiescence prevails, fewer times are added or removed. Thus changes to the temporal density of small event times is proportional to the original density of times. Other algorithms are clearly also possible.

As a final point, since the nowcast transform alters the number of small events in an ensemble, the scaled forecast time must also be further scaled by modifying equation (4) as:

$$T_{F,i}=\left({\displaystyle \frac{R_C}{{R\prime }_i}}\right)T_{F,C}$$

(16)

4. Discussions

4.1. Choice of Circular Region

As before, we applied our methods to a circular region of radius 125 km surrounding Los Angeles, CA, using a target magnitude M_T = 6 and small earthquake M_S = 3.49. The most recent large earthquake in the circle was the Northridge earthquake, having magnitude M=6.7, occurring on January 17, 1994.

The choice of circle radius is of course arbitrary, but we use a rule of thumb that the radius should be at least twice the linear dimension of the aftershock zone of the previous large earthquake in the circle (see Figure 1). In addition, [38] has conducted a more quantitative analysis of appropriate circle radius for similar earthquake magnitudes in Greece, and also arrived at a radius of 125 km.

An additional criterion for the choice of radius is the “radius of significant ground shaking” for a given magnitude, for example the distance at which one might experience Mercalli Intensity VI (PGA ~ 0.1g) shaking. This distance is about 100-150 km for an M_T = 6 earthquake, which is our minimum size for the forecast [48,49].

4.2. Examples of Nowcast and Filtered Calculations

With these considerations, we can now calculate the ROC curves, the PPV calendar time forecasts, and exceedance curves for the original ensembles, the transformed ensembles and for the filtered ensembles.

Figure 7 shows plots of the PPV curves for the original intervals, the filtered intervals, and the nowcast transformed intervals. The plot for the original intervals (Figure 7 left) repeats a plot from Figure 4 for convenience of comparison. The results show that the final value of PPV for the original intervals are the smallest of the three, which are in order from left to right, 29%, 35%, and 43%. While generally consistent, the three plots should be viewed as alternative views of the forecast probabilities.

Figure 8 are the corresponding magnitude exceedance curves for the three types of intervals, and Figure 9 shows the same magnitude exceedances for the 25%, 50% (median), and 75% exceedance probabilities as a function of the number of small events that have occurred. Again, the three sets of plots are generally consistent.

We now compute calendar time exceedance probabilities, the waiting times or lifetimes ΔL(t), for large earthquake occurrence. To do this, we need to scale the conditional lifetimes since the last large earthquake M_T within each interval n_i,j. To that end, let the lifetime within an interval n_i,j in ensemble i be denoted by ${∆L\left(t\right)}_{i,j}$ . Then we scale the lifetimes by the inverse of the factor in equation (4):

$${∆L\prime \left(t\right)}_{i,j}=\left({\displaystyle \frac{R_i}{R_C}}\right){∆L\left(t\right)}_{i,j}$$

(17)

where now ${∆L\prime }_{i,j}$ is the scaled lifetime since the last large earthquake. We then use the scaled lifetimes to compute the exceedance probabilities.

Figure 10 shows the comparison for the calendar time exceedance probabilities. In all three cases, we show the Poisson probability curve (dashed line) that is based on the average interval between large earthquakes. Similar to Figure 5, we show the curves for several calendar times in the earthquake recurrence cycles. These exceedance curves are for immediately after the last large earthquake, at the current time (Today); at the current time plus 15 years; and at the current time plus 30 years. Note that a light smoothing has been applied to these curves.

A feature of these calendar time curves is that once the initial large earthquake clustering period has passed, the lifetime until the next large earthquake becomes longer and longer as time passes, until the next large earthquake does eventually occur, as noted in [43,44,45,46,47]. This feature, can be inferred by comparing the zero time curve with the Poisson curve. Initially the zero time curve lies below the Poisson curve until a cutoff is reached (shorter waiting times more probable). Subsequently, the zero time curve lies above the Poisson curve (longer waiting times more probable).

We can then repeat the exceedance curve calculation for natural time (small event counts) using equation (8) directly (Figure 11). Conditional exceedance curves are shown for several examples of natural times (small event counts): for 0 small events (immediately after the last large earthquake); for 448 small events (the current count in the circular region); and for the current-plus-mean-projected number of small events in the circular region. The mean projected number is the average number of small events expected until just before the next large earthquake. This mean projected number is computed by finding the mean of all intervals larger than the current number, 448, minus the current number. Note that a light smoothing has been applied to these curves.

In this case however, one finds the more intuitively expected result that the longer it has been in natural time since the last large earthquake M_T, the shorter the expected natural time until the next large earthquake. Or stated another way, as time since the last large earthquake increases, the probability increases for a shorter lifetime, or time to the next large earthquake.

The reason for the different implications of the calendar and natural time curves is the increasing, nonlinear small earthquake quiescence as time passes. This observation has been discussed in [12]. There it was concluded that the reason for this anomalous slowing down in calendar time (relative to the Poisson curve) is due to an anomalous stiffening/strengthening of the crustal rigidity as the next earthquake approaches, due to the closure of small cracks that, when open, weaken the crust. And as the rigidity of the crust increase, a transition occurs from sliding via unstable stick slip to sliding via stable slip [6,50], as is well known from many laboratory experiments.

4.3. Validating/Benchmarking the Forecasts

An important consideration is to validate the forecasts by comparison to a recognized, accepted forecast or benchmark. The obvious choices are the UCERF2 time-independent forecast, and the UCERF3 time-dependent forecast [51,52]. Figure 12 shows the comparison. The UCERF2 forecast for magnitude 6.7, produced in 2007, computed a 30 year forecast of 67% for the Los Angeles box, a rectangular region that we show in the new Figure 5, and a 63% for the San Francisco Bay area for that rectangular box. In 2016, the UCERF3 forecast published the 30 year forecast of 60% for the same Los Angeles box. The UCERF2 30-yr forecast is especially relevant because we are past the midpoint of the 30-year time window. The UCERF3 time-dependent forecast is a valid comparison because it is the official USGS forecast.

The UCERF3 forecast was produced in 2014, whereas the ensemble forecast was produced now. But the numbers are similar, and within the quoted error bounds for the ensemble method. This similarity provides a useful benchmark for evaluating the quality of the ensemble forecast.

5. Conclusions

We have developed a new method for earthquake forecasting for arbitrary forecast times and for small geographic regions encompassing multiple earthquake faults. It is generally not as appropriate for individual faults, where quiescence may be the normal mode of behavior. Since the forecasts are conducted in regions, rather than on individual faults, stability of the Gutenberg-Richter statistics and the associated b-value are observed using least squares fits to the magnitude-frequency curves [53]. The assumption of similar Gutenberg-Richter statistics also amounts to the assumption that we confine analysis to a single tectonic province [53].

The method is appropriate for use when a forecast for such small regions is desired on a frequent basis, and can be carried out rapidly and easily on a laptop computer in a few minutes of computation. For that reason, the ensemble method constitutes a convenient addition to present forecast methods.