Seven Thousand Felt Earthquakes in Oklahoma and Kansas Can Be Confidently Traced Back to Oil and Gas Activities

1. Introduction

Within Oklahoma and Kansas, the disposal of massive amounts of wastewater, co-produced during oil and gas production, started in the 1940s and began to scale significantly in the 1990s, and increased dramatically with the advent of hydraulic fracturing (HF) operations in the mid-2000s (Hough and Page, 2015). Between the year 2000 and 2021, 39.7 billion bbls (6.3 billion m³) of wastewater were disposed in the region. To put things into perspective, this is roughly equivalent to the total current volume of the Salton Sea lake. Wastewater disposal wells (SWD) inject into laterally extensive aquifers, which are at different shallow depths, with some (e.g. Arbuckle) directly overlying the Precambrian basement. Notably, little to no wellhead pressure is needed during injection (gravity-feed disposal) (Weingarten et al., 2015; Langenbruch et al., 2018), while the median disposal depth varies spatially but is around 1200m. The HF activities occur at shallow low-permeability formations (e.g. Woodford Shale) and at much higher wellhead pressures. The median HF depth is around 2200m, but some stimulations exceed 4000 m in depth. The most significant HF activities are taking in the SCOOP/STACK plays, the Arkoma Basin and the Anadarko Platform. Kansas has hosted very limited HF operations.

Historically, the seismicity rates in Oklahoma (OK) and Kansas (KN) were very low, at about two to four earthquakes above moment magnitude (M_w) 3 per year (Grigoratos et al., 2020a; Rubinstein et al., 2018). However, around 2009 these rates started to increase, then accelerated rapidly around the year 2014, and peaked at the end of 2015 at unprecedented levels, orders of magnitude above the historical levels (Figure 1a). Within OK, the largest event was a M_w 5.8 in Pawnee (Pollitz et al., 2017), with four events larger than M 5 in total (one in 2011 and three in 2016). The estimated statewide earthquake losses have exceeded US$ 85 million, mostly due to non-structural damages (Grigoratos et al., 2021). The frequency of felt seismicity caused unrest to the public as the regional infrastructure has been designed with little to no consideration of seismic demands. At the beginning, the epicenters were mostly concentrated in the center of OK, but then migrated further north and northwest. Since 2016, the seismicity rates have been steadily decreasing and are now much closer to the historical values. The majority of earthquake hypocenters is below the top of basement (Schoenball & Ellsworth, 2017; Choy et al., 2016), with a distinct overlap between the deepest disposal wells and the shallowest earthquakes (Cochran et al., 2020).

Scientific studies attributed this seismicity surge in OK to SWD (Keranen and Weingarten, 2018; McClure et al., 2017; Grigoratos et al., 2020b), including the four largest M5+ events (e.g. Goebel et al., 2017; Chen et al 2017). Most of these studies calibrated adhoc geomechanical models to past seismicity and injection data to argue that a triggering mechanism is plausible. Only McClure et al. (2017) and Grigoratos et al. (2020b) applied robust hypothesis testing protocols to support the alleged association. The SWD rates have been gradually increasing since at least the year 2000, with accelerated growth between 2011 and 2015 (Figure 1a). Since then, they have been steadily decreasing over time. Notably, even though disposal continuously increased between 2000 and 2015, the earthquake rates did not increase until about 2009 and then accelerated in about 2014. Thus, there is an apparent time lag of months to years between increases in the disposal rates and the onset of seismicity (Norbeck and Rubinstein, 2018; Keranen and Weingarten, 2018). That said, we should note that a similar time lag was not observed in the response of the seismicity after the disposal rates started their downward trend.

SWD is considered the primary but not the only source of elevated seismicity levels within OK. Holland (2013) and Skoumal et al. (2018) highlighted that some earthquake clusters within the state have been induced by HF, especially within the SCOOP/STACK zone. Ries et al. (2020) identified several geological and operational factors that increase the apparent probability of HF-induced seismicity within Oklahoma. Notably, in certain areas of the state, the SWD and HF operations overlap both in space and time (Figure 1b).

By 2013 the seismicity had migrated further north into southern KN, with the peak in rates later in 2015, coinciding with the peak in OK. Kansas had also experienced very little seismicity until that time, with only 15 events above M 3 since 1973 (Rubinstein et al., 2018). The largest recent earthquake in KN was the Mw 4.9 Milan event in 2014 close to the southern border (Choy et al., 2016). Wastewater disposal within the Arbuckle formation has been the suspected causal factor for southern KN (Rubinstein et al., 2018; Ansari & Bidgoli, 2020). However, this conclusion has been reached deterministically, demonstrating that such a causal link is physically possible. The uniqueness and statistical significance of this association remains unquantified. Figure 1b and Figure A1 demonstrate that the link in space between SWD and earthquake occurrence is not straightforward. In some regions, the seismicity is collocated with SWD wells, but in other regions there is no spatial correlation. There are vast areas within KN where massive SWD volumes have caused no detectable seismicity.

In this study, we apply an established probabilistic framework (McClure et al., 2017; Grigoratos et al., 2020ab; 2022) to examine the extent to which SWD and HF have contributed to the unprecedented seismicity levels observed in OK and southern KN over the last 15 years. The present investigation is the first one to examine the role both of these activities played throughout the two states, under the same framework. Our approach first hindcasts the seismicity rates after a given time on a spatial grid using either SWD or HF injection data as input, and then compares those rates against the null hypothesis of solely tectonic loading. In the end, each block is assigned a p-value, indicating the statistical confidence of its causal link with each oil and gas activity. Our analysis employed both declustered and non-declustered compilations of regional earthquake catalogs and established physics-based principles. This robust protocol has been previously successfully applied in Oklahoma (Grigoratos et al., 2020b) and West Texas (Grigoratos et al., 2022).

2. Data

First, we compiled a new unified earthquake catalog for the Central and Eastern United States with magnitudes above 1 over the time period of January 2000 to February 2023, merging multiple different sources (Table A1). Duplicate events among catalogs were removed following certain hierarchy rules, when more than one location or magnitude solution was available (Grigoratos et al., 2023). These rules, alongside the margins used during the duplicate-search are listed in Grigoratos and Wiemer (2023). Then, the unified catalog was split into various subcatalogs to account for spatio-temporal variations in the magnitude of completeness (M_c). The OK & KN area had its own subcatalog, containing 82 events above M 4 and 3298 above M 3. Most of our analyses focused on the bulk of the seismicity, below latitude 37.5. There, we estimate a magnitude of completeness (Mc) around 3 after 2006 in OK, around 2.5 after 2012 in OK, around 2.5 after mid-2014 in southern Kansas, and closer to 2 after 2016 for both. These estimates are based on the slope of the magnitude-frequency-distributions and the history of the networks (see Appendix A). The methods we used in this study are not sensitive to the M_c or to the b-value (Grigoratos et al., 2022). For this reason, a more quantitative approach for the estimation of M_c was not needed. A lack of need for accurate b-values is also why we did not homogenize the magnitudes into M_w, and we rather allowed all the various magnitude scales to co-exist. In very broad terms, most magnitudes below 3.5 are local magnitudes, while most magnitudes above 3.5 are moment magnitudes

Sources of wastewater disposal data were the B3 database (for OK), the Oklahoma Corporation Commission (OCC; for OK), the KGS (for KN), the Environmental Protection Agency (EPA; for a few counties in OK), Kyle Murray (Grigoratos et al., 2020a; for OK) and the studies by Weingarten et al. (2015; for all states), Barbour et al. (2017; for Pawnee county in OK) and Norbeck and Rubinsten (2018; for OK and KN). At the time of our analysis, the B3 wastewater disposal data are assumed complete till 2021. The OCC data were assumed complete for 2021, with only daily data for the Arbuckle wells being available for 2022 and early 2023. The KGS disposal were available only through 2021. The EPA data (relevant only for a few counties mostly within OK) were assumed largely incomplete after 2016 (given our sources).

Our SWD analyses required monthly injection rates, the coordinates of the well and is some basins the target-formation. The (true vertical) depth was also collected, although not used explicitly, because all our analyses were in 2D. We could not use the public pressure data because they are largely incomplete and not reliable (Murray, 2015). Only wastewater disposal wells were included in our calculations; Enhanced Oil Recovery (EOR) wells were not taken into account, given that EOR operations aim to stabilize, rather than increase, the pore pressure. Given the spatial extent of our investigation zone, we did not limit our analysis to the Arbuckle aquifer, but rather included all the disposal wells in the region, regardless of their target formation. Experience from the Delaware basin in West Texas has shown that SWD in shallow formations can also cause felt seismicity (Grigoratos et al., 2022). In the future, if precise focal depths become widely available, the analyses can be repeated with depth-specific filters applied to the seismicity and the injection data.

Little to no wellhead pressure is needed for SWD, even for high-rate wells (Weingarten et al., 2015). This indicates that the large-scale, bulk permeability of the Arbuckle Group is likely towards the high end of the reported range (Langenbruch et al., 2018), which is around 0.1 to 2 m²/s (Grigoratos et al., 2020a). The permeability certainly varies across locations, and might increase over time due to secondary porosity from fractures (seismicity).

Sources of HF stimulation data were the FracFocus Chemical Disclosure Registry (for both states) and the OCC FracNotices (for OK). Reporting to FracFocus started in 2011 and became mandatory after January 2012. The OCC requires a Hydraulic Fracture Notice Form (FracNotice) to be filed 48 hours prior to any hydraulic fracturing operations within OK. FracNotices were first required to an OCC district office in 2012. Electronic filing became available starting in July of 2016 and was required after December of 2016. In the FracNotices the operators provide the well surface and bottom-hole location, anticipated number of stages, average fluid volume per stage, scheduled start and end dates of operation, and initiation date of well flowback. In summary, for OK, we used HF data through December 2022 (leveraging the FracNotices), while for KN through September 2022 (~5 month reporting lags are common within FracFocus).

Our post-processing of the HF and wastewater disposal data followed the work by Grigoratos et al. (2020a; 2022) and is described in Appendix A. The median duration of a HF stimulation was 4 days (Figure A2a), and the median daily volume of a HF well was 16890 bbls (Figure A2b). Between 2014 and 2022, the total injected HF volume was 4.5 billion bbls (715 million m³) and the total number of stimulation-days was 15700 (Figure 2). The gridded spatial distribution of cumulative HF and SWD volumes and the earthquake epicenters are mapped in Figure 3a and Figure A3 respectively.

3. Methods

3.1. Declustering

The earthquake recurrence models we use in this study can only model mainshocks. Thus, we removed foreshocks and aftershocks using the declustering algorithm of Aden-Antoniow et al. (2022), which utilizes the nearest-neighbor clustering algorithm of Zaliapin et al. (2008). This algorithm is robust for small samples and stable with Mc, yet sensitive to the d- and w-parameters that are defined by the user (we opted for 1.5 and 0.5 respectively, based on sensitivity tests). We post-processed the clusters to always keep the largest event in each one as the mainshock (instead of the earliest one).

3.2. Hypotheses Testing for Causal Factors of Induced Seismicity

Following Grigoratos et al. (2020b) and Grigoratos et al. (2022), we hindcasted the observed seismicity with a set of competing statistical models to derive likelihood ratios that are then converted to p-values via reshuffling tests (McClure et al. 2007). First, injection and seismicity data are aggregated into spatial blocks of 5 to 10km at fixed time intervals, on a monthly basis for SWD and on a daily basis for HF. Only blocks with nonzero volume for the corresponding operation and with at least 3 earthquakes above Mc are considered “active”. Two block-specific hypotheses are constructed: (i) a null hypothesis model that assumes no relationship between injection and seismicity (total likelihood L₀), and (ii) an alternative model that does assume a relationship between injection and seismicity (total likelihood L₁). The ratio of L₁ to L₀ is defined as R. Although R values greater than 1.0 (i.e., L₁ > L₀) indicate that the alternative hypothesis is more likely than the null hypothesis, this criterion is not sufficient to statistically reject the null hypothesis. To do the latter, we need a reference statistical distribution for the ratio R in which the null hypothesis is true (R_null). This reference distribution enables the comparison of the two hypotheses in a rigorous way, even if they differ in model-complexity. R_null was empirically generated by reshuffling the injection data (with subsequent recalculation of L₁ and L₀). The reshuffling procedure is described in detail in Grigoratos et al. (2020b) and Grigoratos et al. (2022). R_null essentially reflects how likely it is for the earthquake recurrence model to find purely coincidental correlation between the observed seismicity and random injection data. By comparing the generated R_null values in each block with the single R value obtained from the real injection and seismicity time-series (R_obs), we can statistically determine how confident we are that the improved correlation of the alternative hypothesis is not coincidental. The metric used to quantify this confidence level is the statistical p-value, defined as $p=(\eta +1)/\kappa$, where $\eta$ is the number of synthetic (reshuffled) datasets with R_null greater than the R_obs value from the real data, and $\kappa$ is the total number of synthetic datasets (in our case 200). A p-value of 0.05, for example, translates to a confidence level C of 95% that the null hypothesis can be rejected (given that C = 1 – p). For the conditions and assumptions under which these confidence levels can be interpreted as probabilities of causality, the readers are referred to McClure et al. (2017).

The lower the p-value, the better the correlation between the seismicity and the oil and gas activity in question. Very low p-values, that is, below 0.05, signal that the examined human activity is a key triggering factor of the seismicity in that block (Grigoratos et al., 2022). P-values between 0.05 and 0.10 might flag blocks where a subset of the clusters is induced. Blocks with low p-values in two different maps are likely affected by both types of human activity, with more sophisticated statistical analysis needed to properly quantify the exact ratio. Higher p-values above 0.20, and especially above 0.50, indicate that very little to no seismicity there was triggered by the analyzed type of human activity.

Here, we examined 2 causal factors related to energy production, namely HF stimulations and SWD. We ran the p-value analysis two times, assuming a single alternative causal factor at a time, each compared against the null hypothesis that assumes all earthquakes are of tectonic origin. Implementation details regarding spatial oversampling are provided in Grigoratos et al. (2022) and Grigoratos and Wiemer (2023).

Very low p-values, that is, below 0.05, signal that the examined human activity is a key triggering factor of the seismicity in that block. P-values between 0.05 and 0.10 might flag blocks where a subset of the clusters is induced. Blocks with low p-values from SWD and HF are likely affected by both activities, with more sophisticated statistical analysis needed to properly quantify the exact ratio for each. P-values between 0.10 and 0.20 are somewhat inconclusive, while p-values above 0.20, and especially above 0.50, indicate that very little to no seismicity there was triggered by the analyzed activity.

3.3. Generalized Seismogenic Index Model

To hindcast the spatio-temporal changes in the seismicity rates we need a physics-based earthquake recurrence model that takes into account the external driving forces, i.e. the oil and gas activities. Since most of the stress changes are caused either directly or indirectly by the injected volumes, the earthquake recurrence model should use as input the injection rates of HF or SWD. We employed the framework of Grigoratos et al. (2020a; 2022), who expanded the Seismogenic Index model (Shapiro et al., 2010) to large-scale HF operations and wastewater disposal activities. They took into account the background tectonic rate and the stressing-rate dependency of the time lag between injection and seismicity rate changes (Norbeck & Rubinstein, 2018). The original Seismogenic Index model is itself a modified version of the G-R relationship and predicts that the number of induced earthquakes is proportional to the pore pressure change, which can be approximated by the injected volume. It was first successfully used to hindcast induced seismicity related to hydraulic stimulations (Shapiro 2015; Dinske and Shapiro, 2016; Kwiatek et al., 2024). The governing equations behind our earthquake recurrence model, as well as details on how the time-lags were modeled are available in Appendix A.

Following Grigoratos et al. (2020a), we acknowledged that SWD in one block will affect the pore pressure and stress field in a neighboring block. HF volumes are injected into much tighter formations and thus this effect is extremely localized relative to the size of our blocks. The latter is 5 km for HF and 10 km for SWD (Grigoratos et al., 2022). Thus, we distributed the SWD volumes of each well in space and time following the Theis (1935) equation for transient, radial flow in nonleaky vertically confined aquifers of infinite areal extent (equation S4). The only parameter needed to compute the distribution factors is the large-scale diffusivity, $D$. We performed our computations for three different values of D, equal to 0.3, 1.0 and 2.0 m²/s (Figure 3b) and commented on the sensitivity of the results. More implementation details are provided in Grigoratos et al. (2020a) and in Appendix A.

3.4. Hydraulic Fracturing Radar

Following Grigoratos et al. (2022), we used the Hydraulic Fracturing Radar (HFR) as an an initial screening tool for HF within each region of interest. This metric is employed only as a sanity check to confirm or question results stemming from the p-value analyses. To create HFR, first we plotted a 5km base-grid with the percentage of (declustered) earthquakes above M_c that occurred during a HF stimulation and within 5km from it (${EQ}_{HF}^{\%})$. This distance threshold is meant to cover the epicentral uncertainty of the earthquakes and the horizonal extension of the HF wells (toe-to-heel extension). Given the uncertainty in the end-dates of the HF stimulations coming from the IHS database, we allowed for a 3-day buffer after the assigned end-date for the temporal window. Then, we adjusted ${EQ}_{HF}^{\%}$ to account for the maximum total duration of stimulations within a block over the investigated study period. The reason for that adjustment is the observed asymmetry in the intensity of HF activities across blocks.

4. Results

Figure 4 shows the results of the declustering, for the time period between 2000 and 2022, for a magnitude cut off of M 2.5. In reality, we expect the M_c to be closer to 3 in early years, and reach 2.5 only after 2012 or so. The algorithm removed 34% of the events as foreshocks/aftershocks. The corresponding value for the Delaware basin (M ≥ 1.5) and the Midland basin (M ≥ 1.5) was much lower, at 14% and 19% respectively (Grigoratos and Wiemer, 2023). This could be due to the fact that a relatively larger percentage of seismicity in OK and southern KN occurs within the crystalline basement.

4.1. SWD

Regarding SWD, we chose to start the calibration in 2012 because the Mc decreased significantly in 2011 throughout OK (Appendix A). Figure 5and Figure A4 shows the statistical p-values for SWD (p_SWD) in the region for M ≥ 2.5 and with the large-scale diffusivity value D set to 2 m²/s. In particular, 46% of “active” blocks have p_SWD ≤ 0.05, and 55% of “active” blocks have p_SWD ≤ 0.10. Furthermore, 61% of earthquakes occurred within a block with p_SWD ≤ 0.05 and 71% within a block with p_SWD ≤ 0.10. Recall that a p-value of 0.05 translates to a 95% confidence interval for association. The two key calibrated parameters are shown in Figure A5. The results are insensitive to declustering (Figure A6) and only slightly sensitive to the selected D value (Figure A7, Figure A8). This demonstrates the stability of our statistical framework, which is based on physics-based principles.

Our observations are in agreement with previous studies that have identified SWD as the main driver of seismicity in central and north western OK and in southern KN. By not flagging SWD as a causal factor in the SCOOP/STACK play and in the Arkoma basin, we are also in agreement with Skoumal et al. (2018) who argued that HF is the main driver there. Notably, all Mw 4.9+ ruptures (Prague, Cushing, Fairview, Pawnee; Milan) occurred in blocks with p_SWD ≤ 0.05. Thus, we were able for the first time to confidently link all major ruptures to SWD, using a single model-formulation and parametrization, without having to employ complex poroelastic effects or earthquake-to-earthquake nucleations. The latter two mechanisms were required in previous sequence-specific deterministic studies (e.g. Goebel et al., 2017; Barbour et al., 2017; Chen et al., 2017). The generic nature of our framework makes it uniquely suitable for investigations that include multiple human activities at various distances and timeframes.

We should highlight that there are many zones in central and northern Kansas that hosted high-volume SWD operations (Figure 3) without significant levels of seismicity or low p_SWD values (Figure A4). In fact, we did not find any compelling evidence of induced seismicity in central and northern Kansas (Figure A4; Figure A6). Even though detected seismicity in the 1980s near the northern state-border occurred close to oil field operations, Evans and Steeples (1987) made a strong case for a natural tectonic origin. Furthermore, we believe that the evidence that another northern sequence between 1986 and 1992 is linked to SWD are inconclusive (Armbruster et al, 1989).

Hough and Page (2015; 2016) argued in favor of the El Reno 1952 Mw 5.7 event in central ok being linked to SWD. Our results do not support this claim, since the epicenter (Lat: 35.53; Lon: -97.85) is far from blocks with low p_SWD values (Figure 5, Figure A8). If fact all but one of the SWD-linked historical earthquakes outlined in Hough and Page (2015) occurred outside the perimeter of our SWD zone, as defined by our p_SWD maps.

We also investigated, on a block-by-block basis, the relationships between total diffused or injected disposal volume, p_SWD, maximum magnitude and seismicity rate (Figure A15). We limited out investigation to blocks with p_SWD below 0.05. The maximum magnitude was negatively correlated with p_SWD indicating that larger magnitudes tend to occur in blocks where SWD is the key triggering force. A very interesting finding was that the seismicity rate (above Mc) and the maximum magnitude were positively corelated only with the diffused disposal volume, and not with the injected (non-diffused) ones. This highlights how important it is to distribute the SWD volumes in both space and time following pore-pressure diffusion principles.

4.2. HF

In this section we will analyze the potential link between the observed seismicity and HF. The largest volumes are in the southern and western part of OK (Figure 2 and Figure A3). There we also find the highest HFR values (Figure 6), in agreement with Skoumal et al. (2018). Most, if not all, of the seismicity within the SCOOP/STACK plays, the Arkoma Basin and the Anadarko Platform appears highly correlated with HF operations (Figure 6). According to our window-based associations, 7% of declustered earthquakes are linked to only 3% of the HF stimulations. We should note that 83% of HF stimulations that eventually induced seismicity (according to HFR), started doing so during the stimulation period, while 71% of “active” blocks have a fitted time-lag (equation A3; Figure A9a) of 3 days or less. Thus, we caution against using larger time-lags in well-to-earthquake association filters. Lastly, the median Σ_HF value was -2.6 (Figure A9b).

Figure 7 shows the statistical p-values for HF (p_HF) in the region, for the period between 2014 and 2022. In particular, 22% of “active” blocks have p_HF ≤ 0.05, and 36% of “active” blocks have p_HF ≤ 0.10. Furthermore, 8% of declustered earthquakes (above Mc) occurred within a block with p_HF ≤ 0.05. When we combine the HFR windows with p_HF, 3% of declustered earthquakes above M 2.5 are both flagged by HFR and occurred within a block with p_HF ≤ 0.05 (Figure 8). These earthquakes are almost certainly induced by HF. Rare exceptions to this rule may be events that also have a low p_SWD value, hosted in blocks with relatively stable seismicity rates through time (Figure A14).

The p_HF values are stable over time and were able to detect zones where HF is a causal factor early on (Figure A10). They are also almost identical when we use the constant time-lag ${\mathrm{t}}_{lag}^c$ function (Figure A11) instead of the time-dependent one (equation A3).

Crucially, the vast majority of blocks with low HFR values also have low p_SWD values, indirectly validating the consensus that the seismicity in central OK and sKN is mostly triggered by SWD. That said, there are several blocks in OK with HFR between 0.10 and 0.25 and a few with HFR above 0.25 within zones linked to SWD (Figure 5 and Figure 6). Those blocks also have low p_HF values in Figure 7, since HFR and p_HF values are highly correlated (Figure A12), as expected. Therefore, there are areas within central and northern OK where the seismicity is likely occurring at different depths, with shallow earthquakes triggered by HF and deeper ones (within the basement) triggered by SWD. This is a novel finding of this study. The poor hindcasting performance of the earthquake recurrence model for HF (Figure A13) within blocks with p_HF ≤ 0.05 indicates that in blocks where both causal factors are in play, the majority of the seismicity is linked to SWD.

As far as the HF wells are concerned, only 1% of stimulations are both linked to seismicity and occurred within a block with p_HF ≤ 0.05. On the other hand, 92% of stimulations are not linked to seismicity and occurred within a block with p_HF ≥ 0.05.

We also investigated, on a block-by-block basis, whether either the maximum magnitude was correlated with the p-value or with the total volume (provided that p_HF ≤ 0.05) (Figure A15). The maximum magnitude was weakly correlated with p_HF (Figure A15a), indicating that larger magnitudes tend to be linked to SWD, and not HF. Surprisingly, the maximum magnitude was not positively correlated with the total HF volume per block, implying that geologic conditions are more important than the injected volume as far as the size of the ruptures is concerned.

5. Conclusions

The seismicity levels in Oklahoma and southern Kansas have increased dramatically over the 15 years. Several past studies have identified the disposal of wastewater co-produced during oil and gas extraction as the driving force behind this surge, with a small percentage being attributed directly to hydraulic fracturing (HF) stimulations. While the present investigation validates that SWD is the main causal factor, it identifies several previously undocumented clusters of seismicity within Oklahoma that appear triggered by HF stimulations. We did not find evidence of induced seismicity in central and northern Kansas, nor close to epicenters of old historical earthquakes within Oklahoma.

When we expand our findings from both examined human activities (Figure 5 and Figure 7) to the entire catalog between 2009 and 2022, oil and gas operations are probabilistically linked at the 95% confidence level with 65% of the felt (M≥2.5) earthquakes (7731 events), including 46 events above magnitude 4 and all 4 events above magnitude 5. Our results are essentially grid-independent and are not particularly sensitive to the declustering of the earthquake catalog, to the magnitude of completeness, to the investigation period, and to the large-scale hydraulic diffusivity of the aquifers.

Notably, the seismicity is often unevenly correlated with the two oil and gas activities, with vast areas that have undergone intense operations having hosted no detectable seismicity. Thus, geomechanical factors dominate over operational ones, both for HF and SWD. This was the first time a study examined the role both of these activities played throughout the two states, under the same physical and statistical framework.