Improvements to the Third Uniform California Earthquake Rupture Forecast ETAS Model (UCERF3‐ETAS) Open Access

Modern earthquake hazard and risk analysis depends on two main scientific modeling components: (1) an earthquake rupture forecast, which gives a complete set of possible earthquake ruptures and their associated probability of occurrence, and (2) a ground‐motion model, which specifies the possible ground‐shaking levels produced by each rupture at a given site. This paper describes recent improvements to the Third Uniform California Earthquake Rupture Forecast ETAS Model (UCERF3‐ETAS) developed by the Working Group on California Earthquake Probabilities (Field, Jordan, et al., 2017; Field, Milner, et al., 2017). UCERF3‐ETAS is a hierarchical model in that it builds on a long‐term, time‐dependent model (UCERF3‐TD) that embodies elastic rebound effects (Reid, 1911), which in turn builds on a long‐term, time‐independent model (UCERF3‐TI) in which segmentation assumptions were relaxed and multifault ruptures were included (among other enhancements). Spatiotemporal clustering was included by adding an Epidemic Type Aftershock Sequence (ETAS; Ogata, 1988, 1998) component in UCERF3‐ETAS, which is simulation based in that it produces any number of synthetic catalogs for M ≥ 2.5 events for a specified timespan and conditioned on any or all prior M ≥ 2.5 events. Interestingly, obtaining realistic behavior requires the inclusion of elastic rebound on faults (otherwise the same event can immediately rerupture, which we do not see in nature) and also requires that we honor the non‐Gutenberg–Richter distributions implied by UCERF3‐TI near faults.

As such, UCERF3‐ETAS represents a candidate basis for operational earthquake forecasting (OEF, which provides real‐time information on evolving earthquake likelihoods; Jordan and Jones, 2010; Jordan et al., 2014), and it is the first such model to consider information on faults (proximity and associated recurrence interval and time since last event) over arbitrarily long sequences. It is also the only comprehensive forecast model that can represent complex events like the 2016 M 7.8 Kaikōura, New Zealand, earthquake, both in terms of it being an extreme multifault rupture (e.g., Hamling et al., 2017) and a possible continuation of the multiyear Canterbury earthquake sequence. A 100 yr time series of statewide event rates is shown in Figure 1, which reveals realistic looking aftershock sequences following larger events.

Similar to all models, however, UCERF3‐ETAS is a limited representation of the actual system in terms of embodying assumptions, approximations, and uncertainties. For example, ETAS is certainly not how earthquakes are physically triggered in the real world; what we hope is that it is an adequate statistical proxy in terms of making useful hazard and risk inferences. ETAS performs remarkably well at smaller magnitudes (e.g., Woessner et al., 2011; Ogata et al., 2013; Taroni et al., 2018; Nandan et al., 2019), but there is always some leap of faith at the larger, hazard‐dominating magnitudes due to limited observations. A number of potential issues and possible solutions have been discussed extensively in our previous UCERF3‐ETAS publications (e.g., Field, Milner, et al., 2017). However, because any improved model will also still be a limited approximation, we deliberately paused further improvements until the value of each could be more carefully ascertained.

A number of studies have examined potential UCERF3‐ETAS usefulness under the assumption that the model is correct. For example, Field, Porter, et al. (2017) presented a prototype operational loss modeling capability, which implies that expected statewide losses can vary by orders of magnitude due to earthquake triggering (and that this effect can be much more important than the epistemic uncertainties in the long‐term model). Similarly, Field and Milner (2018) showcased a variety of potential OEF products for the Haywired planning scenario (an M 7.1 mainshock on the Hayward fault; Detweiler and Wein, 2017), with a particular emphasis on the value of including faults. Both these studies emphasized two important considerations with respect to potential usefulness: (1) quantification of probability gains relative to “normal” or quiet times and (2) the large influence of forecast duration and messaging time delays (latency) due to the temporal decay of triggered‐event rates. More recently, Field et al. (2021) examined the influence of declustering and the Poisson assumption on various hazard metrics by comparing with fully time‐dependent results (based on UCERF3‐ETAS catalogs). Again, these valuation studies examined potential usefulness under the assumption that the model is correct, which is important because it is possible that a simpler model (e.g., one that ignores faults) could be more useful in terms of giving the same result with less computation time.

The more influential studies with respect to the improvements described here have had a greater emphasis on model validation (testing how well the model agrees with reality). For example, Page and van der Elst (2018) defined a suite of “Turing Style” tests (inspired by Alan Turing’s question of whether communication with an artificial intelligence has become indistinguishable from that with a human); their goal was to test whether UCERF3‐ETAS catalogs are indistinguishable from real ones. This analysis led to two modifications implemented here. One is a relatively minor bug that caused simulated magnitude–frequency distributions to have a Gutenberg–Richter b‐value of 0.98 rather than 1.0 below M 3.3, which was caused by sampling events from a magnitude‐probability distribution rather than from a magnitude–frequency distribution. The second and much more influential modification was to add intrinsic variability in aftershock productivity from one sequence to the next (discussed in the Aleatory Variability of Aftershock Productivity section).

The 2019 M 7.1 Ridgecrest earthquake sequence provided a real‐time, real‐world application that produced two influential studies. Milner et al. (2020) describe having UCERF3‐ETAS simulations running on a high‐performance computing (HPC) cluster within 33 min of the M 6.4 foreshock (the 4 July Searles Valley earthquake). Although the model has never been fully operational (automated), the Ridgecrest deployment did benefit from previous work by Milner to enable others to launch simulations on HPC, which is exactly what had to happen following the M 7.1 mainshock on 5 July (Milner was unavailable, so coauthor Savran launched the latter simulations). This “on‐demand operationalization” enabled Milner et al. (2020) to evaluate several variables as the sequence progressed, including the influence of the evolving inferred Ridgecrest finite‐rupture surface, periodic updates to incorporate the aftershocks that had occurred up to that point (machine read from the U.S. Geological Survey [USGS] Comprehensive Catalog), and the influence of various model options. More importantly, these UCERF3‐ETAS simulations enabled a quantitative answer to perhaps the most frequently asked question: what is the likelihood of a large, triggered earthquake on the Garlock or San Andreas faults? However, the most important result with respect to the modifications described here is the fact that the observed aftershock productivity differed from the average value utilized in UCERF3‐ETAS. Recall that Page and van der Elst (2018) concluded that we needed intrinsic variability in aftershock productivity from one event to the next, so it is to be expected that the value for any specific event will differ from the mean. Following an actual event, the productivity uncertainty is epistemic in that it will become better known as more aftershocks occur, so the ability to set this for past events was added to UCERF3‐ETAS. Within a forecast simulation, however, the productivity is applied as intrinsic variability (aleatory uncertainty), as described in the Aleatory Variability of Aftershock Productivity section.

The second paper (Savran et al., 2020) used the Ridgecrest sequence to more formally evaluate UCERF3‐ETAS using a new set of tests developed in the Collaboratory for the Study of Earthquake Predictability (CSEP; e.g., Zechar et al., 2013; Michael and Werner, 2018). The main testing innovation is the direct use of synthetic catalogs (rather than using implied average event rates), which enables the relaxation of a limiting Poisson assumption utilized in previous CSEP tests (Savran et al., 2020). Their analysis identified a number of possible model improvements, with the most important being a reiteration of the need for sequence‐specific aftershock productivity. However, their analysis also highlighted an unsatisfactory visual manifestation of an approximation that causes the aftershock zones of large off‐fault events to be circular rather than clustering around a linear fault zone. Fixing this issue is the second UCERF3‐ETAS modification described here (in the Treating Off‐Fault Events with Finite Rupture Surface section).

In addition to the aforementioned studies, UCERF3‐ETAS forecasts were also considered in California Earthquake Prediction Evaluation Council deliberations following the 2019 Ridgecrest mainshocks. Although the council preferred the sequence‐specific aftershock probabilities given by the USGS’s online forecast for the sequence as a whole, the UCERF3‐ETAS results provided the only quantitative estimate of rupture of the Garlock fault, which was of particular concern to the council. UCERF3‐ETAS results were also used more directly in ensemble forecasting, along with conventional ETAS modeling, during the 2016 Bombay Beach swarm (McBride et al., 2019). This forecast was released to the public and emergency managers.

One final UCERF3 test is consistency with the observed century‐long hiatus in paleo‐observed events on the major faults in California. This lull was first noted by Jackson (2014), who found that the observed lack of ruptures at California paleoseismic sites precluded Poissonian or lognormal earthquake occurrence, provided that earthquake times between distant sites were not correlated. Biasi and Scharer (2019) did their own analysis and found only a 0.3% probability of the current 100 yr gap, assuming a subset of paleoseismic sites are independent. When examining the UCERF3 model for lulls in earthquake activity, there is no need to choose a subset of sites, because any correlations are explicitly represented by the finiteness of fault‐based ruptures. Page et al. (2019) therefore examined UCERF3 model predictions (synthetic catalogs) for a century‐long lull at all UCERF3 paleoseismic sites, with a number of UCERF3 characteristics turning out to be influential. As shown in Figure 2, although the chance of observing zero events over a generic 100 yr period is only 0.5% for the UCERF3‐TI, when elastic rebound effects are added (UCERF3‐TD), this increases to 1.7%. For the fully time‐dependent model with aftershock clustering, UCERF3‐ETAS, 3.5% of 100 yr periods have zero events (because adding clustering in some time periods requires others to be more quiet). Furthermore, the chance of a 100 yr hiatus sometime within a generic 1000 yr period is quite high (28%). Should we reject the UCERF3‐ETAS model based on the observed hiatus being 3.5% likely? Although this is not formally rejected at a two‐tailed 95% confidence, it may warrant more investigation in the future, especially with respect to the influence of the UCERF3‐ETAS modification described next.

In the ETAS model, every earthquake can trigger a set of primary aftershocks, the latter of which can trigger their own aftershocks (which are secondary with respect to the first one), and secondary events can trigger tertiary events, and so forth. The ETAS‐model parameter that controls primary aftershock productivity (the number of triggered events) is k, which we previously set to the statewide average inferred by Hardebeck (2013); this and other ETAS‐parameter values are listed in table 3 of the UCERF3‐ETAS main report (Field, Milner, et al., 2017). The expected number of M ≥ 2.5 aftershocks from a given event also depends on other ETAS parameters, the start and end time of the forecast, and the magnitude of the mainshock (see Field, Milner, et al., 2017, for equations and table 4 therein for example numbers for different mainshock magnitudes). For present purposes, we can quantify productivity as the expected number of primary aftershocks that have a magnitude greater than that of the mainshock (thereby removing the mainshock magnitude dependence, except at M > 6.5, in which values get reduced due to magnitude–frequency‐distribution rolloff); the value implied by the average productivity (k) is 0.053 within 10 yr, meaning M≤ 6.5 events have about a 5% chance triggering a primary aftershock that is bigger than themselves within 10 yr (and this percentage is approximately tripled when all generations of aftershocks are considered; table 4 of Field, Milner, et al., 2017).

We previously applied this average productivity to all events, with aleatory (intrinsic) variability represented by random sampling from a Poisson model using the expected number of M ≥ 2.5 primary aftershocks. Figure 3a, which represents one of the Turing tests applied by Page and van der Elst (2018), compares the cumulative distribution of the number of triggered events for mainshocks observed over the last 28 yr to those from 200 separate 28 yr catalogs from UCERF3‐ETAS simulations. This comparison implies that the observations exhibit both more quiet and more active sequences, meaning more intrinsic variability of productivity from event to event is needed in UCERF3‐ETAS.

To address this, and based on the results of Page et al. (2016), we assume that the event‐to‐event variability of k follows a lognormal distribution with the same mean as used previously and with a coefficient of variation (standard deviation divided by the mean, written as kCOV here) being an adjustable parameter. Setting kCOV = 0 produces the same results as obtained previously, whereas increasing values imply more event‐to‐event variability. A productivity value for each simulated event is sampled randomly from the lognormal distribution. After trying a range of values, we settled on kCOV = 1.5, the results of which are shown in Figure 3b with respect to the Turing test (for comparison with the kCOV = 0 case in Fig. 3a). The random productivity gets applied only to primary aftershocks of the given event (secondary, tertiary, etc., events get their own random value). In addition, kCOV = 1.5 is not a formal best‐fit value but is rather obtained via more limited trial and error given the computational demands of simulations (each trial consisted of 1000 500 yr simulations and used approximately 1400 central processing unit hours).

Figure 4 shows the broad range of productivity values implied by kCOV = 1.5. The most likely productivity value (the mode of the distribution) is a factor of 6.7 below the mean, and the median value is about 50% of the mean, consistent with the need for more quiet aftershock sequences. There is a 28% chance a productivity value greater than the mean will be chosen and a 5% chance a productivity 3.7 times greater than the mean will be selected.

UCERF3‐ETAS makes a distinction between earthquakes on explicitly modeled faults versus those that occur elsewhere, in which the latter are represented with a gridded seismicity model based on ∼2 km sized cells in three dimensions. We previously treated such “off‐fault” events as point sources, even at larger magnitudes, but in which a finite‐source correction was made in any downstream hazard or risk calculations (consistent with USGS National Seismic Hazard Model standard practice; Moschetti et al., 2015). The problem is that this produces a circular distribution of aftershocks around large, gridded‐seismicity events, which does not look realistic, as noted by Savran et al. (2020). For example, Figure 4a shows one of many UCERF3‐ETAS simulations following the M 6.4 Ridgecrest foreshock (the 4 July Searles Valley earthquake); this particular simulation was chosen because it triggered an M ∼7.1 gridded‐seismicity event (similar to the 5 July Ridgecrest mainshock), for which the aftershock distribution is circular rather than revealing a fault surface.

To address this issue, we added the option to assign a randomly sampled finite‐rupture surface to all events that exceed a user‐specified magnitude threshold (M 6.0 was used here). The magnitude, rake, and average dip of gridded seismicity ruptures are predefined (in the underlying UCEFT3‐TI model). We first compute the rupture area from the event magnitude using the scaling relationship (also specified in the underlying UCERF3‐TI model). If the average dip is 90°, the hypocentral depth is less than 13 km (the average seismogenic thickness in UCERF3), and the square root of the rupture area is also less than 13 km, we assume a square rupture shape (length is equal to down‐dip width) and center the surface at 6.5 km; if the hypocenter is above the resultant top of rupture then the surface is translated up to just include the hypocenter (or likewise translated downward if the hypocenter is below the implied rupture bottom). If the square root of rupture area is greater than 13 km, then the rupture top and bottom are set to 0 and 13 km depth, respectively, and the length is equal to the rupture area divided by 13 km. The rupture strike is chosen at random (a uniform distribution between 0° and 360°) and the finite surface is centered on the hypocenter (with the trace assumed to be a straight line). For ruptures with a hypocentral depth greater than 13 km, the rupture bottom and maximum down‐dip width are assumed to be equal the hypocenter depth (rather than 13 km). The algorithm is similar for nonvertical dipping ruptures, but in which maximum down‐dip widths are greater (by 1 over the sine of the dip), and the surface is also translated laterally to ensure that the hypocenter intersects the rupture surface (see Data and Resources for a link to the code for full details).

Figure 5b shows the consequent result for another UCERF3‐ETAS simulation following the M 6.4 Ridgecrest foreshock, which again was chosen because an M ∼ 7.1 gridded‐seismicity event was triggered (like the actual Ridgecrest sequence and the case in Fig. 5a). Here, the randomly chosen rupture surface of the M ∼ 7.1 “mainshock” is clearly revealed by its aftershocks, representing a much more realistic aftershock spatial distribution than that depicted in Figure 5a. Although this improvement may be strictly cosmetic with respect to hazard and risk inferences (depending on the specific application), it should nevertheless improve model performance in CSEP tests.

Adding aleatory variability to primary aftershock productivity and random finite surfaces for large gridded‐seismicity events are both worthy UCERF3‐ETAS improvements. That said, the match obtained in going from a kCOV of 0–1.5 is still far from perfect (Fig. 3b), and it is not presently clear how to improve it further. It could be that catalog completeness issues are at play, other ETAS parameters need aleatory variability, or that we are reaching the limits of what the ETAS proxy model can provide given the real world is certainly more physics based. Similarly, nothing is currently stopping one of the randomly generated, finite ruptures from crossing an explicitly modeled fault (as actually occurred with respect to the Garlock fault in Fig. 5b), and the strike direction could be constrained by nearby observed focal mechanisms. In other words, UCERF3‐ETAS is still an approximation of the actual system. Although these and other limitations could undoubtedly be improved upon further, it is not clear whether doing so would make the model any more useful, so once again we pause until more validation and valuation analyses can be conducted.

Ongoing and near‐term validation efforts include retrospectively testing the model against observed seismicity using the new CSEP tests (Savran et al., 2020) mentioned in the UCERF3‐ETAS Evaluations section. With respect to valuation, we also continue to explore the use of UCERF3‐ETAS for operational earthquake forecasting, including in the definition and pricing of catastrophic bonds (e.g., Franco et al., 2020). UCERF3‐ETAS could also be used to quantify uncertainties on catalog‐based magnitude–frequency distributions inferred at different regional scales (to quantify sampling errors associated with our single observed catalog). Similarly, we can test the efficacy of our declustering and smoothed‐seismicity algorithms in terms of inferring the long‐term spatial distribution of earthquake rates; that is, doing this with synthetic catalogs should reproduce the long‐term rates assumed in the model. In fact, analogous tests could be conducted for any metrics that are inferred from limited datasets, such as fault‐slip rates from a finite number of offset dates. We previously demonstrated that statewide average annual losses fluctuate by up to an order of magnitude due to spatiotemporal clustering, but to quantify ultimate usefulness we also need to explore the confidence bounds implied by the aleatory variability (i.e., are the mean fluctuations actionable given the wide range of possible outcomes?). Longer term, we need to develop and utilize multicycle physics‐based simulators (e.g., Tullis, 2012, and papers in that special issue) to evaluate UCERF3‐ETAS assumptions, because our data collection rate will be inadequate to do so empirically, especially at larger magnitudes.

All simulation data presented in this article are available at http://www.WGCEP.org/UCERF3-ETAS, and all calculations were made using OpenSHA (http://www. OpenSHA.org), which in turn utilizes JFree‐Chart (http://www.jfree.org/jfreechart) for making plots. Scripts and documentation for configuring and running the model described in this article are available at http://github.com/opensha/ucerf3-etas-launcher. All websites were last accessed in May 2021.

The authors acknowledge that there are no conflicts of interest recorded.

Third uniform California earthquake rupture forecast model (UCERF3) development was supported by the California Earthquake Authority, U.S. Geological Survey (USGS), and the Southern California Earthquake Center (SCEC) (contribution Number 11677). SCEC is funded by National Science Foundation (NSF) Cooperative Agreement EAR‐1600087 and USGS Cooperative Agreement G17AC00047. Calculations were performed at the University of Southern California (USC) Center for Advanced Research Computing. The authors thank the USGS Powell Center for Analysis and Synthesis for supporting workshops on UCERF3‐epidemic‐type aftershock sequence model (ETAS) development. The authors also thank two anonymous TSR reviewers and Thomas Jordan for constructive comments.

Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.