Narrow Spatial Aftershock Zones for Induced Earthquake Sequences in Oklahoma

The current paradigm for estimating regional long‐term seismic hazard involves declustering whereby aftershocks are removed from an earthquake catalog to identify the underlying background seismicity rate. Inaccurate declustering can ultimately underestimate or overestimate the regional seismic hazard. In Oklahoma, estimating a background seismicity rate is complicated by the highly variable seismicity rate over the last decade. To improve conventional declustering methods used in hazard modeling in Oklahoma, we scrutinize the aftershock windows used for declustering and investigate how aftershocks decay in space and time to establish data‐driven parameters for aftershock windowing. We observe that the spatial decay of aftershocks is more rapid in Oklahoma than in Southern California, motivating the need for smaller spatial declustering windows in Oklahoma. Temporal aftershock decay is statistically indistinguishable between Oklahoma and Southern California, suggesting that temporal declustering windows derived for Southern California are likely sufficient for Oklahoma.


Introduction
A by-product of heightened unconventional oil and gas production in the Central United States over the past decade have been sharply rising and highly variable seismicity rates, mostly linked to the disposal of produced wastewater (Ellsworth, 2013). Most of the recent seismicity in the broader region has occurred in the state of Oklahoma, which has been characterized by high productivity swarms (Benz et al., 2015), elevated background rates (Walsh & Zoback, 2015), and four large mainshocks of M 5 and greater since 2011 Goebel et al., 2017;McGarr & Barbour, 2017;McMahon et al., 2017;Walter et al., 2017;Yeck et al., 2017). Although the earthquake rate has been declining since 2016 in Oklahoma due to both market-driven reductions in new production wells in central Oklahoma and mandated regional wastewater injection rate reductions (Baker, 2017), it remains well above pre-2009 levels. Given these conditions, various researchers have attempted to forecast future seismic hazard for the state Langenbruch et al., 2018;Langenbruch & Zoback, 2016;Norbeck & Rubinstein, 2018;Petersen et al., 2016).
Past locations and rates of earthquakes are a key component of earthquake hazard forecasts. Probabilistic seismic hazard models in the United States have historically been developed using long-term seismicity rates and patterns of tectonically driven background activity. Over long times, the occurrence of an event might be random within a given activity rate, such that the seismicity approaches a statistical Poissonian distribution (Cornell, 1968;Gardner & Knopoff, 1974;Petersen et al., 2016). The goal of such efforts are long-term or time-independent forecasts of future activity, requiring the removal of short-term rate bursts during aftershock sequences, where aftershocks are removed using different declustering techniques (van Stiphout et al., 2012). Recently published short-term U.S. Geological Survey (USGS) hazard forecasts for the Central United States use a space-time window declustering procedure from Gardner and Knopoff (1974), which was originally derived for Southern California (e.g., Petersen et al., 2016). In addition to mainshock-magnitude scaled space-time windows (Gardner & Knopoff, 1974;Uhrhammer, 1986), other commonly used declustering techniques include cluster linking using assumptions about postmainshock stress distributions (Reasenberg, 1985), stochastic declustering based on point processes (Zhuang et al., 2002), and nonparametric network-tree aftershock identification (Baiesi & Paczuski, 2004;Zaliapin et al., 2008).
In Oklahoma, both tectonic stresses and short-term variations in fluid injection activity are thought to influence earthquake rates. As a consequence, some hazard models for Oklahoma incorporate physical changes of fluid pressure (Langenbruch et al., 2018) and fault stressing conditions (Norbeck & Rubinstein, 2018) in an effort to link seismicity with injection rate changes. These studies produce different results, partly because the forecasted seismicity rates are compared with different types of catalogs that are either not declustered (Langenbruch et al., 2018) or declustered using parameters in the Reasenberg (1985) method, which were derived from a California catalog (Norbeck & Rubinstein, 2018). A properly declustered catalog could allow for meaningful comparisons between otherwise differently constructed models or forecasts. Additionally, a declustered catalog might allow for clarity in understanding how external effects such as pressure/stress changes drive background seismicity or lack thereof.
To date, declustering windows have not been defined specifically for Oklahoma. We derive mainshock magnitude-dependent aftershock identification windows for recent seismicity in Oklahoma, using techniques from statistical seismology to study the aftershock decay directly. For comparison, we also study Southern California aftershocks since this region has been used to derive multiple commonly applied declustering algorithms (e.g., Gardner & Knopoff, 1974;Reasenberg, 1985). We first examine the spatial decay of stacked aftershocks near mainshocks and define new spatial windows for different magnitude ranges. Then, we fit Omori-Utsu p values to the temporal decay of sets of individual sequences. We suggest declustering windows specific to Oklahoma seismicity and highlight the importance of constrained declustering parameters for understanding induced seismic hazard in Oklahoma and beyond.

Data
We utilize earthquake catalogs from the Oklahoma Geological Survey for 1 January 2009 to 1 November 2018 and from Shearer et al. (2005) via the Southern California Earthquake Data Center for 1 January 1984 to 1 January 2003. This specific California catalog is used to validate our use of methodology from another study (Felzer & Brodsky, 2006), which used the same catalog. For both catalogs, most events below M 3.5 are reported as local magnitudes (M L ) and most events above M 4 are moment magnitudes (M W ). Due to the mix of reported magnitudes and magnitude calculations from each seismic network, we validate our overall approach by comparing Oklahoma catalogs from two different seismic networks and find that they both yield the same primary results of our study (Figures S1 and S2 in the supporting information).
In general, a well-defined catalog completeness, M c , is critical for statistical analysis of earthquake catalogs due to the heterogenous nature of data acquisition and processing within seismic networks (Gulia et al., 2012;Schorlemmer & Woessner, 2008). We estimate a single M c for each catalog by finding the point of maximum curvature of their frequency-magnitude distribution (Woessner & Wiemer, 2005) for 1,000 and 5,000 event-wide moving windows for Oklahoma and California, respectively. We then take the median value over all windows and find M c = 2.2 for Oklahoma and M c = 1.5 for Southern California, which is applied uniformly in all analyses for both regions ( Figure S3). Additionally, we test the sensitivity of M c on our subsequent analyses and find our results to be robust for M c = 1.5-2.5 in Oklahoma and M c = 0.9-2.2 in Southern California (Figures S4-S7).

Methods
To define clusters, we segregate catalogs for mainshocks of a certain magnitude range and identify neighboring earthquakes in time and space with fixed windows. We create composite catalogs of earthquakes associated with mainshocks by calculating epicentral distances from each windowed earthquake to its mainshock. Great-circle distances are calculated using the haversine formula as presented by Vasylkivska and Huerta (2017), originally derived by Mendoza y Ríos (1795). Cluster-specific subcatalogs with distances recorded between earthquakes and the mainshock are stacked with all others in a mainshock magnitude range and sorted by distance to the common mainshock. Linear density, the number of aftershocks per unit length, is calculated between neighboring stacked earthquakes using the nearest neighbor method (Silverman, 1986). The densities are calculated by taking the inverse of the differential distance between successive data points. We validate our overall stacking approach by reproducing the results of an earlier publication (Felzer & Brodsky, 2006; Figure S8).
In order to study the decay of primary aftershock sequences, we reduce contamination from other aftershock sequences by disqualifying potential mainshocks if larger earthquakes have occurred nearby in space and time. For Southern California, we require a distance separation of 100 km between mainshocks and for no larger earthquakes to occur within 3 days before or 0.5 days after potential mainshocks. For Oklahoma, we use the same time criteria but require a distance separation of only 25 km between mainshocks given the paucity of earthquakes greater than M 6. These parameters are unrelated to mainshock magnitude-scaled aftershock windows but rather are chosen to conservatively reduce contamination from larger mainshocks' aftershocks at regional distances for times around the occurrence of a candidate mainshock. The spatial decay observations for the two regions are insensitive to increases in all three parameters (Table S1).
We compile data sets of stacked earthquake linear density across space in Oklahoma and Southern California for multiple magnitude ranges and time windows. Earthquakes are selected within 250 km of each mainshock and for time windows short enough to minimize background seismicity. Theoretically, choosing short time windows reduces the effect of background seismicity and emphasizes earthquakes possibly linked with a given mainshock. However, in practice, data availability prohibits choosing time periods that are sufficiently short. Thus, we use time windows from 1 to 72 hr for comparing the two regions, which represents a good trade-off between statistical robustness within the smaller data set and the influence of background in Oklahoma (Table S2). The same time windows are used for both regions to minimize any biases due to window selection, since the number of stacked earthquakes, and thus, linear density values vary for different time windows in a given region.

Analysis and Results
For the stacked earthquake catalogs in Oklahoma and Southern California, the event density decay with distance from mainshocks allows for a qualitative separation of aftershocks and background. We fit an inverse power law using least squares to the aftershock decay portion of the data. Mainshocks of 3 < M < 5 show that Oklahoma has more rapid aftershock decay with distance when compared to Southern California for the same time windows, where those power law exponents clearly differ by~0.6-1.1 (Figures 1a and 1b). We assess the statistical significance of this observation by conducting a two-sample Kolmogorov-Smirnov test and confirm that the stacked aftershock data are indeed from different continuous distributions at the 1% significance level ( Figure S9). To further evaluate the robustness of this result, we vary the time windows used to create the stacked aftershock data sets. These tests show that spatial aftershock decay is consistently more rapid in Oklahoma than in Southern California, although the decay rates flatten with increasing time due to gradual inclusion of background seismicity ( Figure S10).
We also observe a relative difference in the distance range where aftershocks likely transition to background activity for the two regions. In general, mainshocks of 3 < M < 4 for Oklahoma clusters are contained within a 4-to 7-km radius of the mainshock and Southern California clusters are contained within 15-20 km. For 4 < M < 5, Oklahoma clusters are contained within 6-12 km and Southern California clusters are contained within 20-25 km. Since values of cluster window lengths are chosen only by visual inspection of event density changes with distance, we report a range of possible window lengths. This also allows us to factor uncertainty into the measurements, which in the following section are used to develop our spatial aftershock identification model.
Mainshocks of 5 < M < 6 show the same overall behavior in Oklahoma with rapid spatial aftershock decay and tightly confined aftershock clustering ( Figure S11). We have fewer mainshocks in this magnitude range to analyze within Oklahoma (Table S2). Since M c is important for evaluating statistics of earthquake catalogs, we compute M c for each of the four Oklahoma mainshocks individually. The rapid spatial aftershock decay of the larger events is qualitatively consistent with observations for the smaller mainshocks of 3 < M < 5. We also consider the effect of the uniformly applied M c on the aftershocks of all mainshocks of 3 Aftershock decay is fit to the median linear density values in log-spaced bins (red and blue dots). Green lines (G.&K.*) are spatial windows from Gardner and Knopoff (1974). Solid black lines (W.&C.*) are empirical subsurface rupture radii of mainshocks from Wells and Coppersmith (1994). < M < 5. We find that the difference in spatial decay rates between the two regions is unchanged for M c = 1.5-2.5 in Oklahoma and M c = 0.9-2.2 in Southern California (Figures S4 and S5).
At short distances close to the stacked common mainshock, we observe flattening of density values (Figure 1), because larger earthquakes rupture across a spatial dimension within the same order of magnitude of the distances over which the aftershocks occur. Thus, in Figures 1b and S11, we plot estimates for empirical subsurface rupture radii (Wells & Coppersmith, 1994) as a visual guide and observe power law decay beyond those rupture radii for both regions. Next, we test if higher location uncertainty in Oklahoma explains the observed differences in spatial decay using a high-resolution, waveform cross-correlated relocated Oklahoma catalog (Schoenball & Ellsworth, 2017b). We obtain consistent results as in Figures 1 and S10 using the relocated catalog, suggesting that the influence of relative location differences on spatial decay is insignificant in Oklahoma (Figures S12 and S13). Lastly, we test whether our spatial decay results are sensitive to the choice of cluster identification technique by utilizing the network-tree algorithm from Zaliapin et al. (2008) since it has often been applied in Oklahoma (Vasylkivska & Huerta, 2017;Zaliapin & Ben-Zion, 2016). We use the network-tree approach as applied in Oklahoma by Goebel et al. (2019) to identify clusters and find rapid spatial density decay of aftershocks that is, qualitatively, in close agreement with the results of Figures 1 and S11 (Figures S14-S16).

Empirical Model for Spatial Aftershock Windows in Oklahoma
Using our observations of inferred separation between aftershocks and background for a range of mainshock magnitudes, we determine a model for aftershock identification in space for Oklahoma ( Figure 2). Given the range of distances expected to contain aftershocks for three magnitude bins in Figures 1 and S11 (red dashed lines), we plot spatial window radius versus mainshock magnitude. We populate the three boxes with a 0.1 km × 0.1-magnitude-unit mesh grid (n = 1,463) to obtain an unbiased distribution of possible aftershock windows. We fit, in a least squares sense, the distribution of data as an increasing exponential function: r ¼ 10 0:22M−0:02 ±2δ: (1) In equation (1), r is the circular window radius in kilometers around a mainshock, M is the mainshock magnitude, and 2δ = 2.56 km, which is the 95% prediction interval of one tail of the distribution as derived from the standard error.
The lower bound of this prediction interval allows for reasonable identification of aftershocks in Oklahoma with the lowest background contamination, which we apply in section 4 to study temporal decay for long time periods (months to years): r ¼ 10 0:22M−0:02 −2:56 km ð Þ: We also aim to constrain aftershock identification windows for Southern California to compare temporal aftershock decay with Oklahoma. However, given the order-of-magnitude more earthquakes in the California data set (Table S2), event rates are more sensitive to changes in the time interval used for aftershock selection than in Oklahoma. Thus, to better identify aftershocks and reduce background contamination for long time intervals in California, we let r = 10 km for 4 < M < 6 as seen in Figures 1 and S11 (vertical blue lines). Windowing models of Gardner and Knopoff (1974) and Uhrhammer (1986) are compared.

Methods
We analyze seismicity within 2 years after each mainshock, since temporal aftershock windows are on the order of months to years in models like Gardner and Knopoff (1974). To account for the longer time window and possible inclusion of background seismicity and secondary activity with time, we use conservatively small spatial windows to define aftershock zones for both regions to reduce the amount of unassociated seismicity in space. As described in section 3.3, we use equation (2) for Oklahoma and r = 10 km for Southern California. For each cluster, we systematically determine the Omori-Utsu parameters for aftershock decay with time (Utsu, 1969) following the modified Omori formula: In equation (3), dn/dt is the earthquake rate (in events per day), t is the time (in days) after the mainshock, K is the aftershock productivity, c is the completeness time of aftershock detection, and p is the decay rate of aftershocks with time. To estimate the Omori parameters for each sequence, we use the maximum likelihood method following Ogata (1999) with a constrained optimization algorithm for nonlinear, multivariate functions. This procedure allows us to find optimized parameters for the statistical model using bounded constraints on K (5-300), c (0.02-2), and p (0.2-2.7).

Analysis and Results
Using conservatively chosen spatial windows (section 3.3) and an identical procedure to estimate the Omori-Utsu parameters (section 4.1), we use p values as a relative metric of the temporal decay rate for sets of aftershock sequences. For mainshocks of 4.5 < M < 6, we find median p values of 0.78 and 0.83 for Oklahoma and California, respectively (Figure 3). These p values fall reasonably within a range of 0.6 to 2.5 as determined from a comprehensive, global study (Utsu et al., 1995). We determine 95% confidence intervals for the two distributions of p values using bootstrap resampling over 100 and 200 iterations for Oklahoma and California, respectively, and find that the distributions overlap ( Figure S17). We also consider the sensitivity of M c on the observed p value distributions and find the general trend to be unchanged for M c = 1.5-2.5 in Oklahoma and M c = 0.9-2.2 in Southern California (Figures S6 and S7). Overall, the temporal decay is statistically indistinguishable between the two regions. To assess the robustness of our p value results, we consider both data artifacts and properties of postmainshock seismicity, which may bias our results. Each subcatalog has 2 years of postmainshock seismicity and significantly more data per sequence than in the spatial decay analysis. We consider the effect of secondary activity larger in magnitude than the mainshock, as well as the amount of data to be used in the Omori parameter estimation, and find in both cases negligible effect on the p value results (Figures S18 and S19).
We also compare our p value results for different fixed c values and again find agreement with our initial result ( Figure S20). We show that p values for time windows from t = 3 days to t = 1 year are consistent with that of Figure 3 for the same conservatively small spatial windows ( Figures S21  and S22). However, upon changing the spatial window from our refined model to that derived by Gardner and Knopoff (1974), we find that the original median p values of 0.78 and 0.83 for Oklahoma and California decrease to 0.46 and 0.52, respectively ( Figure S23). This flattening is due to the inclusion of unassociated background seismicity outside of the inferred aftershock zone. Additionally, when we identify clusters with the network-tree approach (Goebel et al., 2019), we obtain a slightly higher median p value of 0.92 ( Figure S24). These findings highlight the sensitivity of temporal decay rates to the spatial windowing parameters and cluster identification approach used.

Physical and Geological Interpretations
In this study, we measure spatial and temporal aftershock decay in Oklahoma and Southern California using mainshock magnitudedependent space-time windows. We find that temporal aftershock decay is indistinguishable between Oklahoma and Southern California ( Figure 3). Our temporal decay result is consistent with other statistical studies of Oklahoma seismicity (Llenos & Michael, 2013;Schoenball & Ellsworth, 2017a;Walter et al., 2017), which found that aftershock decay rates for plausibly induced mainshocks are not distinguishable from tectonic mainshocks.
We also observe a significant difference in spatial decay between the two regions. Spatial aftershock decay is likely controlled by multiple factors, including the mainshock stress perturbation and fault network properties. Felzer and Brodsky (2006) found that the spatial decay of aftershocks in Northern California was more rapid than in Southern California, possibly indicating a smaller fractal dimension of the fault network in Northern California. This explanation may be consistent with our observations in Oklahoma. However, another study found that secondary aftershocks significantly decreased the aftershock decay rates in Southern California when a windowing approach is used (Marsan & Lengliné, 2010). We observe rapid spatial decay in Oklahoma despite using a window-based method, which may indicate that aftershock sequences in Oklahoma have less secondary triggering.

Implications for Declustering and Hazard Modeling
Given that declustering directly determines the rate used for rate-based seismicity forecasting, our study highlights improvements that could be made when declustering Oklahoma earthquake catalogs with window-based methods. We find that the spatial extent of aftershocks in Oklahoma is overestimated bỹ 20-30 km when using the windows of Gardner and Knopoff (1974), which causes seismicity independent of the mainshock to be treated as dependent and discarded from contributing to the Poissonian activity rate ( Figures 1 and S25). Windowing with this method, as done by the USGS (e.g., Petersen et al., 2016), produces catalogs that are "overly declustered" and likely reduces the amount of the background activity to be used in the hazard estimation. Our findings suggest that for the USGS 1-year hazard forecasts, the hazard estimate within the current framework might be underestimated in Oklahoma, regardless of the assumptions in the framework.
There is a growing body of evidence questioning whether seismicity follows a stationary Poisson process (Corral, 2004;Mulargia et al., 2017) as is assumed in probabilistic seismic hazard assessment (Cornell, 1968;Frankel et al., 1996). When we decluster the Oklahoma catalog using either Gardner and Knopoff (1974) or with our preferred parameters (equation (2) and Figure 2), we find relatively rapid changes in the background rate over the past decade, which cannot be described by a single, stationary Poisson process (Figure 4). We also compare our window-declustered catalogs with the network-tree-declustered catalog (Goebel et al., 2019) and find similarly varying background rates ( Figures S26 and S27). This observation is qualitatively similar to other studies (e.g., Ellsworth, 2013;Walsh & Zoback, 2015) that induced seismicity background rates are nonstationary and roughly follow changing regional injection volumes.

Conclusions
We compare spatial and temporal aftershock decay in Oklahoma and Southern California using mainshock magnitude-dependent, space-time windows. We find that spatial aftershock density decays more rapidly in Oklahoma than Southern California. This may be controlled by differences in the physical or geologic conditions of the fault networks in the two regions. Aftershock decay with respect to time is statistically indistinguishable between Oklahoma and Southern California. Understanding how rapidly aftershocks decay in Oklahoma allows us to constrain the parameters used to define aftershock identification windows. These results are important for any type of fixed window declustering that has been used in probabilistic seismic hazard assessment in the Central United States, either for induced or tectonic seismicity. Inclusion of our suggested parameters or overall approach to constraining fixed aftershock windows might improve future forecasting efforts. In Oklahoma, where seismicity rates remain well above historic levels, understanding changes in the background seismicity rate is vital for assessing the long-term induced seismic hazard.