Optimizing Earthquake Nowcasting With Machine Learning: The Role of Strain Hardening in the Earthquake Cycle

Abstract

Nowcasting is a term originating from economics, finance, and meteorology. It refers to the process of determining the uncertain state of the economy, markets or the weather at the current time by indirect means. In this paper, we describe a simple two-parameter data analysis that reveals hidden order in otherwise seemingly chaotic earthquake seismicity. One of these parameters relates to a mechanism of seismic quiescence arising from the physics of strain-hardening of the crust prior to major events. We observe an earthquake cycle associated with major earthquakes in California, similar to what has long been postulated. An estimate of the earthquake hazard revealed by this state variable time series can be optimized by the use of machine learning in the form of the Receiver Operating Characteristic skill score. The ROC skill is used here as a loss function in a supervised learning mode. Our analysis is conducted in the region of 5° × 5° in latitude-longitude centered on Los Angeles, a region which we used in previous papers to build similar time series using more involved methods (Rundle & Donnellan, 2020, https://doi.org/10.1029/2020EA001097; Rundle, Donnellan et al., 2021, https://doi.org/10.1029/2021EA001757; Rundle, Stein et al., 2021, https://doi.org/10.1088/1361-6633/abf893). Here we show that not only does the state variable time series have forecast skill, the associated spatial probability densities have skill as well. In addition, use of the standard ROC and Precision (PPV) metrics allow probabilities of current earthquake hazard to be defined in a simple, straightforward, and rigorous way.

Keywords: earthquake cycle; earthquakes; machine learning; nowcasting; receiver operating characteristic; strain hardening.

Figure 1

(a) Map of the region, 5° × 5° area centered on Los Angeles. Smallest blackdots represent earthquakes M > 3.29. Blue circles represent earthquakes M ≥ 6, red circles represent earthquakes M ≥ 6.9. Yellow star is the location of GNSS and InSAR data that is shown in Figure 4a. Monthly seismicity M > 3.29 as a function of time for the region in panel (a and b). Blue dashed line is the Exponential Moving Average (EMA) of the monthly seismicity rate with N = 36.

Figure 2

Receiver Operating Characteristic (ROC, (a), state variable time series (b), and Precision (PPV, (c). The figure represents frame 380 of a movie showing the time development of the nowcast. The movie is available in the supplementary material. The red dot represents the date for the frame, which was 27 March 1999. The vertical dashed red lines in (a) are the dates of the large earthquakes M ≥ 6.9, the vertical black dotted lines are the earthquakes M ≥ 6.0. Panel (c) represents the chance that an earthquake M ≥ 6.75 will occur during the following 3 years. The chance of such an earthquake M ≥ 6.75 at that moment is 52.9% as shown, both by the number in the green box and the red triangle at bottom. Skill is the area under the ROC curve in (a), which we find as 0.708. The no skill line is the diagonal line from lower left to upper right. Since the no skill area is obviously 0.5, this implies that the nowcast time series has significant skill. The cyan lines in (a) are ROC curves for 200 random time series, whose mean is the diagonal no skill line. The dotted lines bracketing the no skill line from lower left to upper right represent the ±1σ value for the random time ROC series. More results are shown in Table 1.

Figure 3

(a) Spatial ROC diagrams (a) for two grid box sizes, 0.5° × 0.5° and 1° ×  1°. (b) Spatial probability densities (PDFs). For the PDFs, the future earthquakes during the 3 years after 27 March 1999 are shown as circles on the maps. Smallest circles are for 4.9 ≥ M ≥ 4.0, next larger circles are for 5.9 ≥ M ≥ 5.0, next larger circles (blue) are for 6.9 > M ≥ 6.0, and largest circles (red) are for M ≥ 6.9. In the ROC curves in for the two PDFs in (a), the cyan curves are 48 spatial ROC curves evaluated at 1‐year intervals beginning with 1970. The curves are evaluated using all earthquakes M > 3.29 for the following TW = 3 years. Thresholds are established for all values of the PDFs, and the values of TP, FP, FN, TN are computed, with the ROC curves being computed from those. Again, the dashed line from lower left to upper right is the no skill line. These and other data are listed in Table 1.

Figure 4

Surface deformation from Sentinel InSAR (a), and GNSS data (b) at a point having latitude 35.57, longitude −117.68, near the epicenter of the M7.1 Ridgecrest, California earthquake of 5 July 2019, over the years from 2015 to 2019 (yellow star in Figure 1, Cerro Coso Community College—CCCC). Using LiCSBAS, a spatial‐temporal filter was applied that is high‐pass in time and low‐pass in space with a Gaussian kernel that is identical to filters used in other time series InSAR software packages such as StaMPS. The spatial filter width used is 2 km and the temporal filter width used is 0.16 years. After applying the filter, the velocity computed is 4.9 mm/yr. The fact that the deformation over the years prior to the Ridgecrest event is linear suggests that small earthquakes did not stop or decrease with time. When compared to the decrease in small seismic earthquakes prior to the large earthquakes seen in Figure 2, it suggests that there was a transition from unstable stick‐slip events to stable sliding. This is because unstable events can only be seen with seismometers, whereas both unstable and stable slip can be observed with geodetic data.