# Seismological Research Letters

- © Seismological Society of America

## ABSTRACT

The 24 August 2014 *M*_{w} 6.0 South Napa mainshock produced fewer aftershocks than expected for a California earthquake of its magnitude. In the first 4.5 days, only 59 *M*≥1.8 aftershocks occurred, the largest of which was an *M* 3.9 that happened a little over two days after the mainshock. We investigate the aftershock productivity of the South Napa sequence and compare it with other *M*≥5.5 California strike‐slip mainshock–aftershock sequences. While the productivity of the South Napa sequence is among the lowest, northern California mainshocks generally have fewer aftershocks than mainshocks further south, although the productivities vary widely in both regions. An epidemic‐type aftershock sequence (ETAS) model (Ogata, 1988) fit to Napa seismicity from 1980 to 23 August 2014 fits the sequence well and suggests that low‐productivity sequences are typical of this area. Utilizing regional variations in productivity could improve operational earthquake forecasting (OEF) by improving the model used immediately after the mainshock. We show this by comparing the daily rate of *M*≥2 aftershocks to forecasts made with the generic California model (Reasenberg and Jones, 1989; hereafter, RJ89), RJ89 models with productivity updated daily, a generic California ETAS model, an ETAS model based on premainshock seismicity, and ETAS models updated daily following the mainshock. RJ89 models for which only the productivity is updated provide better forecasts than the generic RJ89 California model, and the Napa‐specific ETAS models forecast the aftershock rates more accurately than either generic model. Therefore, forecasts that use localized initial parameters and that rapidly update the productivity may be better for OEF than using a generic model and/or updating all parameters.

## INTRODUCTION

On 24 August 2014, an *M*_{w} 6.0 strike‐slip earthquake occurred near the city of Napa, California. The South Napa earthquake occurred in the West Napa fault zone and was the largest earthquake in the San Francisco Bay area since the 1989 *M*_{w} 6.9 Loma Prieta earthquake (Brocher *et al.*, 2015). Within 18 min of the mainshock, the first aftershock probability report was released, based on generic California aftershock parameters in the Reasenberg and Jones (1989; hereafter, RJ89) model, currently used by the U.S. Geological Survey (USGS) to produce aftershock forecasts following *M*≥5 earthquakes in California, as well as globally on an *ad hoc* basis. However, in the days following the mainshock, many fewer aftershocks occurred than was expected for a California earthquake of its magnitude. During the first 4.5 days of the sequence, only 59 *M*≥1.8 aftershocks occurred, the largest of which was an *M* 3.9 that happened a little over two days after the mainshock. In contrast, during the same time period following the 2004 *M* 6.0 Parkfield earthquake, there were over 200 *M*≥1.8 aftershocks, 6 of which had magnitudes greater than 4.

Although the aftershock productivity of the South Napa earthquake seems unusually low, the coseismic and rapid postseismic surface deformation is unusually high for an earthquake of its magnitude. Rupture nucleated around a depth of 9–11 km, then propagated primarily up‐dip and to the northwest (e.g., Dreger *et al.*, 2015; Ji *et al.*, 2015; Wei *et al.*, 2015; Hardebeck and Shelly, 2016). The coseismic surface slip peaked in the northern segment of the fault rupture roughly 10 km away from the hypocenter, whereas the southern segment is dominated by rapid shallow postseismic slip (Floyd *et al.*, 2016; Lienkaemper *et al.*, 2016). Most of the aftershocks occur deeper than the coseismic and postseismic slip, and they occur to the west of the fault rupture, likely in a zone of secondary faulting (Hardebeck and Shelly, 2016). Hardebeck and Shelly (2016) suggest that the lack of on‐fault aftershocks in the southern segment of the fault where postseismic slip dominates implies a lack of stick‐slip patches, which could contribute to the low productivity of the aftershock sequence, although they note that this is unlike many creeping faults, such as the San Andreas fault at Parkfield.

In this article, we investigate the aftershock productivity of the South Napa earthquake and focus on its impact on operational earthquake forecasting (OEF). First, we determine whether the South Napa earthquake indeed had an anomalously low aftershock productivity, compared with other *M*∼6 earthquakes in California and compared with prior seismicity that occurred in the region. We show that aftershock productivity can vary widely on regional scales, and low productivity is, in fact, common for the Napa region. Then, we examine the impact that low aftershock productivity has on the seismicity‐rate models that are commonly used for OEF, in particular the RJ89 model, which is currently the basis for aftershock forecasts produced by the USGS, and the epidemic‐type aftershock sequence (ETAS) model (Ogata, 1988). We use the South Napa aftershock sequence as a test case to explore which strategies result in more‐accurate aftershock forecasts and demonstrate that accounting for regional variation in aftershock parameters can significantly improve forecasts.

## IS THE AFTERSHOCK PRODUCTIVITY OF THE SOUTH NAPA EARTHQUAKE ANOMALOUS?

To determine if the aftershock productivity of the South Napa earthquake is anomalously low, we compare it to the productivity of 16 other *M*≥5.5 strike‐slip earthquakes in California (Fig. 1a and Table 1). As a simple and robust gauge of productivity, we count the number of *M*≥3 aftershocks within 30 days and within a magnitude‐dependent spatial window around each mainshock as defined by Gardner and Knopoff (1974). The productivity of the South Napa earthquake is somewhat lower than other similar‐sized mainshocks (Fig. 2). Although there is some amount of variability in the aftershock counts, the expected scaling of aftershock productivity with magnitude is evident.

The mainshocks are classified into northern or southern California earthquakes, based on whether they occur in the authoritative regions covered by the Northern California Seismic Network (NCSN) or the Southern California Seismic Network (SCSN), respectively. Notably, northern California earthquakes in general seem to produce fewer aftershocks, on average, than southern California earthquakes. We perform an analysis of covariance (ANCOVA) to determine if the mean productivities of the two groups are significantly different, because ANCOVA allows us to separate the variation between the two groups from the variation in the common dependent variable shared by the two groups (magnitude). The linear regression and ANCOVA results are shown in Table 2. The best‐fitting lines to the two groups do not have significantly different slopes but do have significantly different intercepts (*F*‐test *p*‐value of 0.006), with the northern California aftershock productivities having an adjusted mean that is roughly a factor of 3 lower than the southern California productivities. It is possible that some of the difference may be due to differences in the magnitudes that are computed by the NCSN and SCSN. However, to account for the factor of 3 difference in mean productivity, the magnitudes would have to be off by roughly 0.5 units, and so it is unlikely that this can explain all of the signal that is observed. Therefore, the low productivity of the South Napa sequence may reflect regional variations in aftershock productivity.

Aftershock productivity is also parameterized in the earthquake‐rate models used in OEF, such as the RJ89 model. The RJ89 model combines the Gutenberg–Richter magnitude–frequency distribution (Gutenberg and Richter, 1944) with the modified Omori law describing the temporal decay of aftershock rate following a mainshock (Omori, 1894; Utsu, 1961, 1971). Therefore, the rate *λ* of earthquakes with a magnitude greater than or equal to *M* at time *t* following a mainshock of magnitude *M*_{m} is (1)in which *a*′ is a productivity parameter, *b* is the slope of the magnitude–frequency distribution, and *c* and *p* are the Omori–Utsu parameters. This equation can be rewritten as (2)(Michael, 2012), so that the product 10^{a′}10^{bMm} represents the aftershock productivity *κ* of the sequence, or productivity (3)To compute the productivity *κ* for each mainshock in our data set, we therefore need to first estimate *a*′ and *b*‐values for each earthquake. We obtain these parameters, shown in Table 1 and Figure 3a,b, using maximum‐likelihood estimation (Aki, 1965; Shi and Bolt, 1982; Reasenberg and Jones, 1989). The parameters *c* and *p* are held fixed at values of 0.05 and 1.08, respectively, which are the generic California model parameter values from Reasenberg and Jones (1989).

Figure 3c then shows the derived productivity *κ* for each earthquake. As in Figure 2, the northern and southern California earthquakes appear to be distinct populations. We once again perform an ANCOVA to determine if they are significantly different. For the northern California earthquakes, because the *M* 7.0 Loma Prieta earthquake is the only data point at the upper end of the magnitude scale, there is a great deal of uncertainty in the slope, such that we cannot rule out that the slopes of the northern and southern California earthquakes are the same. Therefore, consistent with our aftershock count results, the slopes of the regressions are again not significantly different, but the northern California earthquakes have a significantly lower intercept than the southern California earthquakes (Table 2). In fact, there appears to be even more of a separation between the northern California and southern California productivities. This can be seen visually in Figure 3c and quantitatively in the fact that the adjusted intercepts in Table 2 show that the adjusted mean productivity in the south is ∼10 times higher than in the north (compared to ∼3 times higher, from the aftershock counts). Using this parameterization, the South Napa aftershock productivity does not appear as anomalously low compared to other northern California earthquakes as it did using the aftershock counts (Fig. 2). It is possible that by holding the parameter *p* fixed at the same value for both northern and southern California earthquakes, actual differences in *p* could map into apparent differences in productivity. However, when we allow *p* to vary between sequences, there is no systematic difference in *p* between the northern and southern California earthquakes, but there is still a significant difference in productivity, so it is unlikely that changes in *p* are a major contributor to the difference in the productivity we observe.

To further examine changes in aftershock productivity in the Napa region, we can use the ETAS model (Ogata, 1988), a statistical model of seismicity rates that considers aftershock triggering using all the sequences in a region. Then, we can determine whether the productivity of the South Napa earthquake is anomalously low compared to past seismicity in the Napa region. We will fit an ETAS model to the seismicity that occurred in the region from 1980 through the day before the South Napa earthquake (23 August 2014) and extrapolate it forward in time to determine whether the South Napa sequence is well explained by the model.

In the ETAS model, the earthquake occurrence rate in a region is the summation of a background rate of independent events and an aftershock rate triggered by each earthquake in the catalog. Given the occurrence time *t*_{i} and magnitude *M*_{i} of each earthquake *i* in the catalog that occurs prior to time *t* and the completeness magnitude *M*_{c} of the catalog, then the earthquake occurrence rate *λ* is (4)in which *μ* is the background seismicity rate, *K* is the aftershock productivity, *α* establishes how efficiently an earthquake of a given magnitude triggers aftershocks, and *c* and *p* are the Omori–Utsu parameters that describe the aftershock decay after a mainshock (Omori, 1894; Utsu, 1961). The parameters *μ*, *K*, *c*, *α*, and *p* are optimized to fit an earthquake catalog from time *T*_{1} to *T*_{2} by maximizing a point‐process likelihood function: (5)(Daley and Vere‐Jones, 2002). We fit the ETAS model to a catalog of *M*≥1.8 earthquakes that occurred from 1 January 1980 to 23 August 2014 in a roughly 22 km×28 km box encompassing the South Napa earthquake (the area shown in Fig. 1b). The magnitude of completeness was estimated at 1.6 using the maximum curvature (MAXC) algorithm, although we use a more conservative value of 1.8, because the MAXC algorithm sometimes underestimates completeness magnitudes by 0.1–0.2 units (Wiemer *et al.*, 1998; Wiemer and Wyss, 2000; Woessner and Wiemer, 2005).

To test the fit of the ETAS model to the data, we convert the earthquake occurrence times *t*_{i} to transformed times *τ*_{i} using the theoretical cumulative function (6)(Ogata, 1992). The transformed time *τ*_{i} therefore represents the cumulative number of events that the optimized model predicts would occur in the time interval [0,*t*_{i}]. If the optimized model fits the observed data well, then the transformed times behave as a Poisson process with unit rate, and the transformed interevent times should be independent and drawn from an exponential distribution, which we can test with a suite of statistical tests (Lombardi *et al.*, 2010; Llenos and Michael, 2013). We use the Runs test to check for independence (Wald and Wolfowitz, 1940), the autocorrelation function (ACF) test to check for the presence of temporal correlation, and the Kolmogorov–Smirnov (K‐S) test to check the goodness of fit of the interevent times to an exponential distribution, using a correction based on Monte Carlo simulation to account for the mean of the distribution being estimated from the data (Lilliefors, 1969). All three tests (Runs, ACF, and K‐S) indicate that we cannot reject the null hypothesis that the transformed interevent times are independent and identically drawn from an exponential distribution, with significance‐testing *p*‐values of 0.4, 0.78, and 0.64, respectively. Therefore, the ETAS model fits this catalog well (Fig. 4a), and the best‐fitting parameters are *μ*=0.0038 events/day, *K*=0.0121 events/day, *c*=0.0197 days, *α*=1.56, and *p*=0.96. We will refer to these parameters as model ETAS‐0.

Given these parameters, we then extrapolate the model through the South Napa earthquake and its aftershock sequence from 24 August through 30 September and find that the model fits the sequence well (Fig. 4b). The fit of the optimized ETAS model can be further assessed by plotting the observed cumulative number of events versus the transformed times, which should be linear with a slope of 1 if the optimized model fits the data well. Positive or negative deviations from this line indicate that the model is under‐ or overpredicting the number of earthquakes, respectively. If Λ(0,*T*) is the transformed time of the last event during the time period [0,*T*] over which the ETAS model is optimized, then the 2*σ* error bounds of the extrapolation can be found using *σ*=[*τ*−Λ(0,*T*)+{*τ*−Λ(0,*T*)}^{2}/Λ(0,*T*)]^{1/2}, based on the fact that the cumulative curves of the transformed times after Λ(0,*T*) should behave as a Brownian process (Ogata, 2005). Figure 4c shows that the transformed times from the South Napa sequence fall well within the 2*σ* error bounds of the prediction using model ETAS‐0. This suggests that the South Napa sequence is not anomalous compared to the past seismicity in the region. Therefore, a low aftershock productivity may be typical in this region.

## FORECASTING THE SOUTH NAPA AFTERSHOCK SEQUENCE

Aftershock productivity is a critical parameter for aftershock forecasting, and we now explore how it affects forecasts during the South Napa aftershock sequence. Both the RJ89 model and the ETAS model are used for aftershock forecasting, and so we can compare how well the two types of models forecast the daily rate of *M*≥2 earthquakes during the first week of the South Napa aftershock sequence. Forecasts with the RJ89 model are made using the equations developed in that paper; forecasts with the ETAS model are made by simulating 100,000 catalogs and using that distribution to estimate the number of earthquakes likely to occur and the uncertainty. We begin with the RJ89 model using the generic California model parameters from Reasenberg and Jones (1989) (*b*=0.91, *p*=1.08, *c*=0.05), with a corrected *a*′ value of −1.85 (Michael, 2012). We call this model RJ‐0.

Figure 5 compares the observed daily rates with the RJ‐0 forecast, starting from 0.01 days after the mainshock. The model does not do a good job of forecasting the daily rates, because it consistently overpredicts the number of *M*≥2 earthquakes throughout the sequence. This is because the generic California model parameters are unable to account for the low aftershock productivity of the region, which suggests that regional differences in aftershock parameters will be critical to keep in mind when producing aftershock forecasts. The ETAS‐0 model (that fits specifically to the long‐term seismicity of the region) results in better forecasts, because it does account for the low productivity of the region, as well as secondary triggering. To determine how much of the improvement comes from location‐specific parameters as opposed to secondary triggering, we also make forecasts using generic California ETAS parameters (Hardebeck, 2013), which we call ETAS‐CA. The ETAS‐CA model still overpredicts the seismicity, although not as much as the RJ‐0 model. To measure the success of a forecast model, we use the observed data in each temporal bin to compute the likelihood of the forecast model, with parameters determined solely from data prior to the start of that bin (Schorlemmer *et al.*, 2007; Zechar *et al.*, 2010). This metric was used in the Regional Earthquake Likelihood Models experiment in California and is currently used by the Collaboratory for the Study of Earthquake Predictability (CSEP). Taking the logarithm of the likelihood in each time bin and summing over all the time bins of the forecast results in the joint log likelihood. Forecast models with higher joint log likelihoods are thus more consistent with the observations. We, therefore, compute the joint log likelihood for each of the models, and the results bear out that the ETAS‐0 model (LL=−15.27) is indeed the better model when compared to the RJ‐0 model (LL=−206.91) and the ETAS‐CA model (LL=−75.47). A difference of 2 in log likelihood is considered significant. This demonstrates that incorporating location‐specific aftershock parameters can significantly improve forecasts.

The RJ89 model makes considerably better forecasts when its parameters are updated midsequence. While holding *b* and *c* fixed at the generic model values, we update the estimates for *a*′ and *p* daily during the sequence (Fig. 6a,b). These estimates are simply the best‐fit parameters to the updated data, rather than the Bayesian updates described in Reasenberg and Jones (1989) and Page *et al.* (2016). Aftershock activity significantly drops after the first ∼8 hrs of the sequence, and the low activity continues for the first two days, causing lower estimates for *a*′ and higher estimates for *p*. This new model, which we call RJ‐ap, does a much better job at matching the lower productivity of the sequence than RJ‐0, and the joint log likelihood is significantly improved (LL=−145.78). However, updating the parameters each day results in estimates that may not be stable and have large uncertainties, especially because of the relatively low number of aftershocks. Therefore, we can also choose to leave *p* fixed to its generic value and simply update the productivity *a*′ each day (model RJ‐a). This results in more‐stable parameter estimates, and a slightly improved joint log likelihood (LL=−143.32).

Similarly, we can update the ETAS model parameters midsequence. In model ETAS‐Kcp, the parameters *K*, *c*, and *p* are updated daily, whereas in model ETAS‐K, only the productivity *K* is updated (Fig. 6c–e). Unlike the RJ models with the updated parameters, ETAS‐K and ETAS‐Kcp do not result in significantly improved forecasts over ETAS‐0. This could be because the parameter changes are again not well constrained, due to the limited amount of new data at each time step or because the ETAS‐0 estimate for *K* is good enough that it does not need updating (e.g., Fig. 6c).

Notably, when comparing the ETAS models and the RJ models, all ETAS models have higher likelihoods than all RJ models. A large part of this, however, may be due to the fact that during the first time step (*t*=0.01 days), the RJ models are still relying on generic California model parameters, whereas the ETAS models, except for the generic ETAS‐CA model, are starting from ETAS‐0, which has taken the lower productivity of the region into account already. If we ignore the first time step and simply compute the joint log likelihoods based on the remaining time steps, we still see that the ETAS models are slightly better than the RJ models, but the improvement is more limited. For example, the log likelihood of the best RJ model is −14.05, only slightly lower than the log likelihood of the best ETAS model, −12.64.

## DISCUSSION AND CONCLUSION

The results of this study show that there is a significant difference in aftershock productivity between northern and southern California mainshocks, demonstrating that aftershock productivity can vary widely on regional scales. Previous studies have similarly found large spatial variations in productivity. Comparing aftershock parameters in different tectonic regions around the globe, Page *et al.* (2016) found that variations in mean aftershock productivity reached almost a factor of 10. Singh and Suarez (1988) examined the aftershock productivity of *M*_{w}≥7 megathrust earthquakes in circum‐Pacific subduction zones and found variations as large as a factor of 40, with western Pacific earthquakes generally triggering more aftershocks than eastern Pacific earthquakes of similar magnitude. They inferred that the variations correlated with the degree of seismic coupling on the plate interface, with the older more‐weakly coupled western Pacific subduction zones having more aftershocks than the younger more‐strongly‐coupled subduction zones. A more recent study by Wetzler *et al.* (2016) confirmed this observation and suggested that the variations in aftershock productivity may be due to differences in both regional tectonics and stress state, as well as differences in source parameters, such as stress drop. Examining whether there are changes in such parameters that correlate with the changes in productivity observed in our California data set is beyond the scope of this article but would be a logical next step.

Regardless of their physical cause, variations in aftershock productivity can have a significant impact on aftershock forecasting models. We showed that the ETAS‐0 model, tuned to seismicity specifically around Napa, produced more‐accurate forecasts during the South Napa sequence than the generic California model parameters used in the RJ‐0 model and the generic California ETAS parameters used in the ETAS‐CA model (Fig. 5). Clearly, as Page *et al.* (2016) suggest, OEF can be greatly improved by having regional aftershock parameters already available for use, especially early on in aftershock sequences when not enough events have occurred yet to constrain sequence‐specific parameters or, in the case where the aftershock productivity is particularly low, making it difficult to constrain parameters. The current OEF practice in other parts of the world, such as Italy and New Zealand, generally relies on initial aftershock parameters that are spatially constant (e.g., Gerstenberger *et al.*, 2014; Marzocchi *et al.*, 2014). One model that does incorporate spatially varying aftershock parameters is the hierarchical space–time ETAS model (Ogata *et al.*, 2003; Ogata, 2011), which is a currently undergoing prospective CSEP testing in Japan and showing some success (Nanjo *et al.*, 2011; Tsuruoka *et al.*, 2012). Our results similarly suggest that ETAS‐based aftershock probabilities with localized parameters would offer additional improvements in forecast accuracy over what is currently done in the United States, where forecasts are still based on the RJ89 model, but regional parameters are starting to be implemented.

An important issue that could affect estimates of aftershock productivity and other triggering parameters is missing early aftershocks following large‐magnitude and/or high‐productivity mainshocks (e.g., Kagan, 2004; Enescu *et al.*, 2009; Omi *et al.*, 2014). This was not an issue for our South Napa study, because of our conservative choice for completeness magnitude and the relatively small mainshock magnitude. We did a sensitivity analysis fitting parameters to aftershock data starting 15 min to 2.4 hrs after the mainshock and found that our results are robust even if some early aftershocks are missing. However, in general, missing early aftershock data can have a large effect on estimates of earthquake rate following a mainshock, and therefore rate‐ and time‐dependent catalog incompleteness models are becoming an important component of earthquake forecasts (e.g., Omi *et al.*, 2013; Hainzl, 2016). Steps are now being taken to account for time‐dependent catalog incompleteness in USGS OEF efforts globally (Page *et al.*, 2016) and in the United States (Hardebeck *et al.*, 2017). Our results support the importance of these efforts by showing that having a set of regional parameters available that capture the lower productivity of a region would provide a better initial model, which could then be improved upon as more aftershock data are collected.

Our estimates of productivity derived from the RJ89 model illustrate another complication in comparing aftershock sequences; in the RJ89 model, the sequence‐specific productivity *κ* is dependent not only on an overall productivity *a*′ but also on the *b*‐value, because productivity scales with 10^{bMm} (equation 3). In a more general formulation, the productivity would scale with 10^{αMm}, as in the ETAS model, with *α* not necessarily equal to *b*. However, with *α* constrained to *b* in the RJ89 model, one can match a lower productivity simply by lowering the *b*‐value, and vice versa. This can be seen, for example, in our California mainshock data set. The lower productivity for the northern California mainshocks relative to the southern California mainshocks is not obvious from the *a*′‐values (Fig. 3b), but the northern mainshocks do have lower *b*‐values, on average (Fig. 3a), which contributes to the distinction between the two groups in Figure 3c. Therefore, some caution is necessary when interpreting differences in either *b*‐values or productivities between earthquakes, because variations in one parameter may be influenced by variations in the other.

Finally, in this article we utilized daily forecasts during the first week of the South Napa aftershock sequence primarily because after one week there were too few events occurring daily to analyze. However, we found that updating model parameters on a daily basis can lead to issues with stability, especially in a low‐productivity sequence like the South Napa. A critical issue for OEF moving forward is to determine both a time window over which forecasts are meaningful and have some skill at estimating future earthquake rates, as well as an optimal look‐back time window over which to estimate parameter changes.

## DATA AND RESOURCES

Earthquake catalogs used in this analysis were obtained from the Northern California Earthquake Data Center at www.ncedc.org (last accessed October 2014) and the Southern California Earthquake Data Center at scedc.caltech.edu (last accessed April 2015).

## ACKNOWLEDGMENTS

The authors are grateful to Morgan Page and Jeanne Hardebeck for their thoughtful comments and reviews, as well as an anonymous reviewer, the Associate Editor, and Editor‐in‐Chief Zhigang Peng for their comments and suggestions.