Fluid‐Earthquake and Earthquake‐Earthquake Interactions in Southern Kansas, USA

An increase in injection activity associated with energy production in southern Kansas starting in 2013 has been linked to the occurrence of more than 130,000 earthquakes (M −1.5 to 4.9) between 2014 and 2017. Studies suggest that the dramatic increase in seismicity rate is related to wastewater injection into the highly permeable Arbuckle formation. Most of the seismicity is located in the underlying crystalline basement, for which hydrological properties and specific fault geometries are unknown. Additionally, some earthquake clusters occurred relatively far (tens of kilometers) from the main injection wells. Therefore, the effect of pore pressure diffusion may be insufficient to explain the relation between the volume of injected fluids and the spatiotemporal evolution of seismicity. Combining physical models (static stress and poroelasticity) and a statistical cluster analysis applied to a high‐resolution relocated catalog, we analyze the evolution of seismicity in southern Kansas. We find that pore pressure changes (Δp) and Coulomb stress changes (ΔCFS) due to fluid diffusion smaller than 0.1 MPa are enough to initiate seismic sequences, which then evolve depending on their distance from the major injection wells. However, we find that earthquake sequences have different seismogenic responses to Δp and ΔCFS in terms of triggering threshold. In regions located close to disposal wells (Harper area) our cluster analysis suggests that both earthquake interactions and fluid diffusion control the evolution of seismicity. On the other hand, at greater distances (Milan area), where clustering behavior suggests greater earthquake interactions, we find that coseismic ΔCFS are larger than Δp.

induced earthquakes (Zaliapin & Ben-Zion, 2016). In particular, the nearest-neighbor approach developed by Ben-Zion (2013a, 2013b) allows for differentiating between background events and clustered events. The former represents independent (Poissonian) seismicity due to a steady forcing such as tectonic stresses, while the latter represents spatially and temporally related events that are likely due to earthquake-earthquake interactions. For example, in an analysis of seismicity in geothermal fields, Schoenball et al. (2015) and Langenbruch et al. (2011) observed a higher rate of background events compared with tectonic seismicity, suggesting that earthquakes directly induced by fluid pressure changes behave as independent background events. Moreover, Schoenball et al. (2015) and Zaliapin and Ben-Zion (2016) also found evidence for the possible presence of repeating earthquakes in geothermal fields. Repeating earthquakes have been observed in Oklahoma and southern Kansas as well (Cochran et al., 2018;Schoenball & Ellsworth, 2017). However, in these regions, cluster analyses (Schoenball & Ellsworth, 2017;Trugman et al., 2017;Vasylkivska & Huerta, 2017) show the predominance of temporally and spatially clustered seismicity over background seismicity, suggesting that seismic sequences in Oklahoma and southern Kansas are initially triggered by pore pressure changes that subsequently evolve following earthquake interactions. Finally, using one month of data recorded by a seismic array in Oklahoma, Cochran et al. (2020) found a large proportion of independent events associated with high stressing rates from nearby disposal wells. This is in contrast with the results of Schoenball and Ellsworth (2017) who found a lower proportion of background events. This may due to the fact that Schoenball and Ellsworth (2017) applied the clustering method to a subset of events identified by a clustering algorithm, rather than the full catalog.
The works above have primarily focused on one particular aspect or methodology related to induced earthquakes. This study uses a novel approach that combines physical models (static stress and poroelasticity) and nearest-neighbor cluster analysis (Zaliapin et al., 2008) applied to a high-resolution relocated catalog (Cochran et al., 2018) to analyze the processes controlling the occurrence of seismicity in southern Kansas ( Figure 1).The benefit of the combined approach is that physical models will also help quantify the relative influence of injection-induced pore pressure and stress changes to coseismic static stress changes resulting VERDECCHIA ET AL. S Hmax from moderate earthquakes. With the statistical analysis, we calculate the proportion of background versus clustered seismicity, and among the latter, we discriminate between swarm-like sequences and mainshock-aftershock sequences.
In the next section, we describe the two largest events that occurred in our study region, and their related foreshocks and aftershocks. The data, methods, and results from the statistical cluster analysis, the poroelastic modeling, and the coseismic stress modeling are described in Sections 3-5, respectively. Finally, a discussion of the results from the three sets of analyses (Section 6) is followed by a summary of our conclusions (Section 7).

The M4.6 Harper and the M4.9 Milan Earthquakes
The October 2, 2014 M4.6 Harper and the November 12, 2014 M4.9 Milan earthquakes are the two largest ever recorded in Kansas. As is true for the majority of the earthquakes in southern Kansas, the events in their respective sequences occurred predominantly in the basement at depths between 2.5 and 5 km for the Harper sequence (Figure 2), and between 3 and 7 km for the Milan sequence ( Figure 3).
The 2014 M4.6 Harper earthquake occurred on a NE-dipping normal fault as suggested by the relocated seismicity and the moment tensor solution (Figure 2). It was preceded by several events, the largest of which (M3.8) occurred 3 days prior, and very close to the mainshock location. Its largest aftershock (M4.1) occurred on June 5, 2015, possibly at the intersection of the main fault with an antithetic structure (Figures 2b and 2c). The moment tensor solutions and seismicity distribution are in agreement with a normal fault optimally oriented in the regional stress field (Alt & Zoback, 2017).
Besides being the largest recorded earthquake in Kansas history, the 2014 Milan M4.9 earthquake also caused structural damage to a few buildings in the town of Milan, located about 5 km from the epicenter (Choy et al., 2016). The mainshock occurred at about 7 km depth on a right-lateral strike-slip fault (Figure 3), possibly on the southern continuation of an already mapped structure which is part of the Nemaha fault system (Choy et al., 2016). The Milan event was preceded by very few events, with the largest being VERDECCHIA ET AL.   (Carr et al., 1986).
.9 located about 1.8 km south-west of the mainshock at a depth of 3.5 km ( Figure 3). Similar to the Harper sequence, the mainshock moment tensor solution and the seismicity distribution are in good agreement with right-lateral strike-slip faults optimally oriented in the regional stress field (Alt & Zoback, 2017).

Method
We apply a nearest-neighbor approach to the relocated catalog to identify seismicity clusters by combining spatial and temporal distance between events (Zaliapin & Ben-Zion, 2013a, 2013bZaliapin et al., 2008). The relocated earthquake catalog contains 24,704 events (Data Set S1) with magnitude of completeness of 0.75 and b value of 0.75 ( Figure S1) (Cochran et al., 2018  ) and the mainshock (M4.9) from the U.S. Geological Survey (USGS), and the Saint Louis University Earthquake Center (SLU). White triangles represent closest disposal wells. Orientation of maximum horizontal stress from Alt and Zoback (2017). (b, c) Cross-sections with earthquakes color-coded with time. Dashed horizontal lines represent the top of the Precambrian basement (Carr et al., 1986).  Figure S1).
This cluster analysis defines the distance η ij between an earthquake j and an earlier earthquake i as: where T ij and R ij are rescaled time and distance, respectively, t ij is the interevent time, r ij is the difference in epicentral or hypocentral distance, m i is the magnitude of the parent event, b is the Gutenberg-Richter frequency magnitude b value, q = 0.5 is the weight of the temporal and spatial distances, assumed to have equal contribution to the final results, and d f is the fractal dimension that we assume to be equal to 1.6, in agreement with the values used in other studies for the same region (Trugman et al., 2017;Vasylkivska & Huerta, 2017). Because our relocated catalog is characterized by similar horizontal and vertical errors, here we use hypocentral locations to estimate the spatial distance between events. The distribution of the nearest-neighbor distance η ij is usually bimodal, with the mode characterized by larger rescaled times and distances representing background, Poissonian-like seismicity, and the short-proximity mode representing the clustered seismicity (i.e., foreshocks and aftershocks). The division between the two can be drawn by defining a cluster threshold c, where the events with η ij > c represent the background seismicity and the events with η ij < c represent the clustered seismicity. The threshold is estimated using the 1D Gaussian mixture model (Martínez-Garzón et al., 2018;Zaliapin & Ben-Zion, 2016). Clusters characterized by multiple events are called families, while events which do not cluster with other events are called singles. The largest event in a family is called the mainshock, and all the events within the family which occurred before and after the mainshock are called foreshocks and aftershocks, respectively (Zaliapin & Ben-Zion, 2013a).
In order to consider the uncertainties connected with varying M c and b value, we examine three different cases ( Figure 4). The first case ( Figure 4a) uses a value of M c of 0.75 and b value of 0.75, the same values calculated for the relocated catalog. Because the original template catalog  and the nonrelocated enhanced catalog have a b value of 1, the second case ( Figure 4b) uses a b value of 1 and the relative M c of 1.5. This choice is justified by the fact that the relocation method has preferentially removed the smaller events causing the b value to vary between small and large events. The third case uses the minimum event magnitude (−0.42) rather than M c and b value of 1 ( Figure 4c). This is because including events below M c has almost no effect on the cluster structure of the events, but provides a larger number of families, and larger cluster sizes to study (Ruhl et al., 2016;Zaliapin & Ben-Zion, 2013a).
Finally, for each identified cluster, we calculate the average leaf depth d m , which is the average topological distance from the cluster leaves (earthquakes with no children) to the root (the first earthquake in the cluster) (Zaliapin & Ben-Zion, 2013b). Large values of d m are typical for swarm-like sequences, while small d m and large cluster sizes (number of events in a cluster) are associated with burst-like sequences.

Results
In all cases, the line separating the background and the cluster modes is given by log 10 η = −5, which allows us to differentiate the number of mainshocks, aftershocks, and foreshocks present in the catalog. All three cases suggest that clustered seismicity dominates over background seismicity, which represents between 25% and 35% of the relocated catalog. Running the same calculation exclusively on the Harper and Milan sequences (Figures 4d and 4e) demonstrates that the vast majority of the nearest-neighbor event pairs occur at small rescaled distances (spatially clustered seismicity), with the Harper sequence ( Figure 4d) showing a larger portion of events occurring at larger rescaled time and shorter rescaled distance when compared to the Milan sequence ( Figure 4e). Figure 5 shows the temporal evolution of mainshocks, aftershocks, and foreshocks calculated using M c of 1.5 and b-value of 1 in relation to injection rates and cumulative seismic moment. We choose to explore this VERDECCHIA ET AL.

10.1029/2020JB020384
5 of 17 particular case because we want to analyze the time-dependent relation between number of background events, foreshocks, and aftershocks and use a higher magnitude of completeness to minimize the potential influence of a changing station distribution. We note that the number of events may have been underestimated for the period of time between March 2014 and July 2014. During this time period, the template catalog  used for the catalog enhancement (Cochran et al., 2018) may be incomplete for M < 3.4 due to the limited number of available seismic stations, as the temporary array was being deployed.
In agreement with previous studies in the same region (Cochran et al., 2018;Rubinstein et al., 2018), we find that the number of earthquakes per day is highest while injection rates are highest, and decreases following the reduction in the injection rates. It is worth noting, however, that the highest peaks in seismicity also correspond with the largest steps in the cumulative seismic moment, at the time of, and immediately following, the two largest events in the catalog. Figure 5a shows how the number of aftershocks dominates over foreshocks and background seismicity during the periods of high seismicity rates (end of 2014 to end of 2015). To emphasize this behavior, in Figure 5b, we show the proportion of background seismicity, aftershocks, and foreshocks over the total calculated time with a moving window of 250 events. The percentage of aftershocks systematically increases with the increase in seismicity rate, and as a consequence, the percentage of background seismicity decreases. The foreshock percentage seems to be partially correlated with background seismicity.
Results from the average leaf depth calculations (d m ) show a general increase of d m with the cluster size in agreement with a mixed behavior of burst-like and swarm-like clustering (Figures 6a-6c). However, most of the largest events (M > 3.5) seem to instead have burst-like behavior. In particular, the family that includes the Milan mainshock (star in Figure 6) represents the classical burst-like behavior, with no foreshocks and a large number of aftershocks. We find that clusters last from a few seconds to several months, and the duration is not always proportional to magnitude (Figures 6d-6f). In fact, for all three choices of M c and b value, the largest event-the M4.9 Milan earthquake-is never the longest duration cluster.

Method
Next, we examine the pore pressure and poroelastic stresses caused by fluid injection in the study region. We use the finite element software Comsol Multiphysics ® following the linear poroelasticity theory (Biot, 1941) to simulate the coupled relation between pore pressure and solid matrix stress changes induced by wastewater injection in the Arbuckle formation from 102 Class II disposal wells located in southern Kansas and northern Oklahoma (Figure 1), which operated from January 2012 to December 2016. Injection data are from the Kansas fluid injection database (http://www.kgs.ku.edu/Magellan/Qualified/class2_db.html).
The governing equations of linear poroelasticity can be written as follows (Wang & Kümpel, 2003): where μ is the displacement vector, p is the excess pore pressure, ε = ∇ · μ is the volumetric strain, G is the shear modulus, ν is Poisson's ratio under drained conditions, α is the Biot-Willis coefficient, M is the Biot modulus, κ is the matrix permeability, η is dynamic viscosity of the fluid, f is the body force per unit volume acting on the solid matrix, and q is the fluid volume injection rate (fluid source density). Both f and q are functions of spatial position x and time t. The stress-strain relation of the solid matrix when pore fluid p is under pressure is given by: where  ij is the Kronecker delta.

10.1029/2020JB020384
8 of 17 (a) Our 3D model is a 140 × 90 × 15 km 3 volume ( Figure S2) that includes three horizontal layers with different elastic and hydrological properties. We set the boundary conditions such that the upper boundary of the model is a free surface, the lower boundary is fixed, and the four side surfaces can only move parallel to the surface, as dictated by the roller boundary condition. We model a no-flow condition for the boundary between the Arbuckle Group and the overlying sediments, and assume a constant thickness of 300 m for the Arbuckle formation (Carr et al., 1986).
We test five models with different permeability values for both the Arbuckle formation (10 −12 to 10 −13 m 2 ) (Carr et al., 1986;Morgan & Murray, 2015) and the crystalline basement (10 −14 m 2 to 10 −16 m 2 ) (Shapiro et al., 1997;Townend & Zoback, 2000) (Table 1). In the first model we assume an isotropic permeability for both the Arbuckle formation and the basement. In the four other models, we adopt an anisotropic permeability for the crystalline basement in order to explain the seismicity distribution at depths larger than 5 km, with higher permeability values parallel to the strike direction of the main mapped faults in the region (striking 020-049°) (Schwab et al., 2017). Table 1 contains the list of parameters for each of the five models.
We use the results from the poroelastic modeling to calculate the change in Coulomb failure stress (ΔCFS): where Δτ is the change in shear stress (positive in the direction of the receiver fault slip), μ is the friction coefficient, Δσ the change in normal stress (positive when the receiver fault is unclamped), and Δp is the pore pressure change. Here, we calculate ΔCFS on optimally oriented faults based on the geometry and kinematics of the faults responsible for the Harper and Milan earthquakes. We use both the relocated catalog and moment tensor solutions from the U.S. Geological Survey (USGS) and the Saint Louis University Earthquake Center (http://www.eas.slu.edu/eqc/eqc_mt/MECH.NA/) to reconstruct the geometry of the faults and to define the kinematics of the two mainshocks. Finally, we compare our models with the relocated seismicity distribution (Cochran et al., 2018)

Results
The largest values of Δp and ΔCFS are concentrated in the southwestern part of our study region, where several wells with the largest total amount of injected fluids are located (Figure 7). At the base of the Arbuckle Group (Figure 7c and 7d) Model 1 shows Δp ranging between 0.03 and 0.14 MPa, while ΔCFS varies between 0.02 MPa and 0.07 MPa. In the basement (Figure 7a and 7b), which is characterized by lower permeability (2 × 10 −14 m 2 ), Δp ranges between 0.02 and 0.09 MPa, and ΔCFS ranges between 0.01 and 0.06 MPa. In general, the contribution of poroelastic stresses to the overall ΔCFS is roughly an order of magnitude smaller than the contribution of Δp (Figure 7e and 7f). When looking at the seismicity distribution, the western part of the study region is characterized by a larger number of events, in agreement with the higher Δp and ΔCFS values calculated in the same area. Interestingly, the largest event in the catalog, the 2014 M4.9 Milan earthquake, is located in a region far from the high-volume injection wells. Modeling results calculated in the location of the Milan earthquake (Figures 8a, 8c, and (Figures 8b, 8d, and 8f).

Method
Here we calculate the ΔCFS due to the M4.6 Harper and the M4.9 Milan earthquakes to study the relation between the static stress changes resulting from two moderate events and the subsequent seismicity. In addition, we also explore the ΔCFS due to the M3.9 event preceding the Milan earthquake. The coseismic slip on the faults is modeled based on the aftershock distribution, and the empirical relationship among event magnitude, rupture length, width, and area (Wells & Coppersmith, 1994), and earthquake stress drop VERDECCHIA ET AL.

10.1029/2020JB020384
10 of 17  (Trugman et al., 2017). We infer the geographic extent of the rupture from the aftershock distribution assuming that the part of the fault that slipped is the one where the aftershock density is lowest. We then adjust the slip distribution within the patches by matching the magnitude and the stress drop.

Results
The coseismic ΔCFS distribution from both the M4.6 Harper earthquake (Figure 9a) and the M4.  (Cochran et al., 2018). (c, d) Pore pressure changes (Δp) evolution (colored lines) calculated at the hypocentral depth of the two events considering different models (Table 1). (e, f) Coulomb stress changes (ΔCFS) evolution calculated at the hypocentral depth of the two events considering different models (Table 1). Dashed lines represent the monthly injections rates of the closest disposal wells. In southern Kansas, the absolute number of background events is clearly higher than before injection started, and it increases concurrently with injection rates ( Figure 5). However, the relative proportion of clustered events is higher than the background events (Figurex 4 and 5a), showing that both an absolute increase in background events and relatively lower proportion of background events are not mutually exclusive. We find that only 25%-35% of earthquakes in the analyzed catalog are independent background events, lower than the 83% calculated for The Geysers, and the 56% for the Coso geothermal field (Zaliapin & Ben-Zion, 2016), regions with very low or no known tectonic seismicity. The Salton Sea geothermal area, however, shows percentages of background seismicity (Zaliapin & Ben-Zion, 2016) similar to our results (31%). The similarity may be due to the fact that the Salton Sea geothermal field is in a region characterized by a complex network of active faults, and is very often struck by events not related to geothermal operations (Martínez-Garzón et al., 2018). In southern Kansas, earthquake-earthquake interactions seem to play a big role in the temporal and spatial evolution of seismicity (Schoenball & Ellsworth, 2017). While southern Kansas shows aftershock percentages comparable to tectonically active regions (Zaliapin & Ben-Zion, 2016), the pore pressure changes may result in overall higher foreshock percentages (∼20%) (Figure 4). The relation between foreshock sequences and fluids has also been observed for tectonic sequences (Jansen et al., 2019;Ruhl et al., 2016). Ruhl et al. (2016) hypothesized the strong involvement of fluids in the 2008 Mogul earthquake sequence based on the behavior of foreshock clusters. This conclusion was also supported by stress inversion and numerical modeling from Jansen et al. (2019).
Repeaters are thought to occur on creeping faults because of the repeated rupture of asperities embedded in an aseismically creeping fault. Their association with induced seismicity is thought to be because Δp are imposed repeatedly at the same location. Using a matched filter technique to analyze the enhanced catalog also used in this work, Cochran et al. (2018) found that colocated groups of events (near-repeaters) have VERDECCHIA ET AL.

10.1029/2020JB020384
12 of 17 longer interevent times in regions closer to injection wells (Harper area), and very short interevent times (temporal clustering) in regions far from the major injection wells (Milan area). This is confirmed by our results (Figure 4d and 4e), which show a larger concentration of events at large rescaled times and short rescaled distances (repeaters) for the Harper sequence when compared to the Milan sequence. The same characteristic is also observed when running the cluster analysis with Mc = 1.5 ( Figure S3).
Our cluster analyses agree with Cochran et al. (2018), who noticed that groups of near-repeating events farther from high-volume injection wells show stronger clustering of interevent times and shorter sequence durations. Their supposition, which is supported by the modeling conducted in this work, was that fluid diffusion plays a bigger role in controlling seismicity clusters located closer to major injection wells, while earthquake interactions may be the dominant factor in the evolution of seismic sequences at greater distances. This same conclusion can also be made by examining the earthquake families we identify through clustering analysis ( Figure 6). The mixed burst-like and swarm-like behavior of most of the earthquake clusters (where d m increases with cluster size) (Figures 6a-6c) suggests the involvement of both fluid diffusion and earthquake interaction as possible triggering mechanisms. However, the cluster that contains the Milan mainshock (star) is a pure "burst-like" sequence, having no associated foreshocks in the same family (Figures 6a-6c) and has a shorter duration for both total sequence (Figures 6d-6f) and the aftershock sequence (Figures 6g-6i) when compared to the two other large events that are located closer to high-rate injection wells. This may further suggest the weaker involvement of fluids in the aftershock evolution of the Milan sequence when compared to other aftershock sequences of M > 4 events located closer to major injection wells. In addition to pore pressure diffusion, the structural characteristics of the two fault zones may explain the differences in behavior between the Harper and the Milan sequences. The temporal (Figure 8a) and the spatial (Figure 3) evolution of the Milan sequence show a more "episodic" behavior, with several earthquake clusters activating different small structures. The Harper sequence, on the other hand, shows a continuous behavior (Figure 8b), with most events occurring on a normal fault and its associated antithetic structure (Figure 2). When we look at the temporal evolution of seismicity in southern Kansas, we find a correlation between the monthly number of events and the monthly injection rates from the wells used in this work ( Figure 5). The peak in seismicity also corresponds with the occurrence of the two largest events in the catalog, the M4.6 Harper earthquake and the M4.9 Milan earthquake. Aftershocks from these two events seem to dominate in the period between late 2014 and early 2015, either when considering the absolute number (Figure 5a) or the relative number of aftershocks over the total (Figure 5b). In particular, Figure 5b shows a recurring negative correlation between background proportion (blue line) and aftershock proportion (black line); when the percentage of aftershocks increases, the percentage of background seismicity decreases. The percentage of foreshocks is fairly stable (∼20%) over the study period.

Significance of the Poroelastic Modeling Results
Because of the high-permeability values (10 −12 to 10 −13 m 2 ) assumed for the Arbuckle formation, our models show relatively small values of cumulative Δp and ΔCFS even in the proximity of the injection wells (∼0.1 MPa), and a large spatial extent of the pressure changes (Figures 7 and S2). The calculated low-poroelastic stress changes when compared to Δp from fluid diffusion are also another consequence of adopting high-permeability values for the Arbuckle formation. This larger contribution of Δp to the final ΔCFS has also been found in other regional-scale hydromechanical models (Zhai et al., 2019(Zhai et al., , 2020. Interestingly, Goebel et al. (2017), modeling the possible triggering mechanisms for the 2016 M5.1 Fairview earthquake, concluded that at distances larger than 20 km, poroelastic stress changes are greater than pore pressure changes. The Fairview region may represent a specific case where several high-rate injectors (up to an average of 100,000 m 3 /mo) are concentrated in a relatively small area, increasing as a consequence the effect of poroelastic stress changes in the far-field. In our model, on the contrary, disposal wells are more sparsely distributed, and only three of them have an average injection rate up to 78,000 m 3 /mo. However, regardless of the mechanism involved, both pore pressure changes and poroelastic stress changes should be taken into consideration when studying induced earthquakes located far (>20 km) from major injection wells (Goebel et al., 2017;Zhai et al., 2019).
Here, we look in detail at the two largest earthquake sequences, the 2014 Harper sequence and the 2014 Milan sequence, the first relatively close to the major injection wells and the second at the north-eastern VERDECCHIA ET AL.

10.1029/2020JB020384
corner of our study region (Figure 8). While adopting different poroelastic models does not significantly change the calculated Δp and ΔCFS values for the Milan sequence, the results for the Harper sequence depend on the permeability value assumed for the Arbuckle formation, with higher Δp and ΔCFS predicted using lower permeability values. As expected, regions surrounding the major injection wells are affected by higher Δp and ΔCFS when we adopt lower permeability values for the targeted formation. It is also worth noting that the major injections wells located in the southwestern portion of our study region strongly influence the Δp and ΔCFS evolution in the area where the Milan sequence occurred (about 50 km from major disposal wells). As a consequence, disposal wells in northern Oklahoma could additionally increase Δp and ΔCFS in southern Kansas as already suggested by Zhai et al. (2020). Therefore, because here we do not consider several high-volume disposal wells in northern Oklahoma, our models may underestimate the cumulative Δp and ΔCFS in southern Kansas.
The vast majority of earthquakes in southern Kansas occur in the crystalline basement (2-8 km of depth), for which hydrological properties and specific fault geometry are debated. We found that whether we assume an isotropic permeability for the basement (Model 1, Figure 7) or instead a regional anisotropy (Models 2-5, Figure S4) following the orientation of the main structures in the region (Schwab et al., 2017), our final results do not significantly change. It is possible instead that Δp in the basement are controlled by local anisotropy (e.g., faults). For instance, Hearn et al. (2018), using a groundwater flow model, suggested that significant Δp (0.06 MPa) at the Milan earthquake hypocenter occurs only if a narrow conductive channel is adopted (fault zone) rather than a conductive 3D volume. However, that study only included disposal wells relatively close to the Milan earthquake in their model , while as suggested in this work and by other authors (Langenbruch et al., 2018;Zhai et al., 2020), higher volume disposal wells at larger distances could also have affected the Milan area. Zhai et al. (2019Zhai et al. ( , 2020 also suggested that the Arbuckle formation and the basement are hydraulically connected by localized fault zones so that Δp and ΔCFS in the basement can be approximated by those at the Arbuckle-basement boundary (Figures 7b, 7d, and 7f).
Here, similar to Zhai et al. (2020), we compare our modeling results with the well-bottom pressure measured in a Class I well (KS-01-077-002) in the northwestern part of our study area ( Figure S5) (Ansari et al., 2019;Peterie et al., 2018). We show that all our models largely underestimate the pore pressure measured in the Class I well. However, when we adopt hydrological parameters similar to Zhai et al. (2020), who assumed a basement permeability two orders of magnitude smaller than our models, the modeled pore pressure values significantly increase. Still, a difference (0.2 MPa) between modeled and measured pore pressure exists, in contrast with the good correlation showed by Zhai et al. (2020). We interpret this difference as being due to the number of injection wells considered in the modeling. While we included 102 disposal wells mostly located in southern Kansas, Zhai et al. (2020) modeled the effect of 668 disposal wells located in both southern Kansas and western Oklahoma, concluding that the high-volume wells in Oklahoma strongly affect Δp and ΔCFS in southern Kansas. However, these differences do not negate the general conclusion of our modeling that relatively small Δp and ΔCFS (≤0.1 MPa) are enough to trigger seismicity in southern Kansas.
Looking in detail at the Harper and Milan areas (Figure 8), we find that these two sequences have different seismogenic responses to Δp and ΔCFS in terms of triggering threshold (between 0.05 and 0.09 MPa for the Harper sequence, 0.01 MPa for the Milan sequence). For tectonic earthquakes, a ΔCFS threshold of 0.02 MPa is often used as the lower limit of triggering stress (Reasenberg & Simpson, 1992;Stein, 1999), although smaller values (<0.01 MPa) have been used to explain aftershock triggering in induced seismicity contexts (Goebel et al., 2004), and earthquakes triggered by strong tides (Cochran et al., 2019). The values cited above suggest that both faults responsible for these two moderate events were close to failure, in particular the fault responsible for the Milan earthquake. Interestingly, when we look at the temporal evolution of Δp and ΔCFS rate of change ( Figure S6), our results show that the largest events (the two mainshocks, and a M4.1 aftershock of the Harper event) always occur slightly after large peaks of rate of change. This, as also concluded by Barbour et al. (2017), indicates that seismicity may be sensitive to the rate of change of ΔCFS rather than cumulative ΔCFS. Moreover, as suggested by the cluster analysis (Figure 6), while the Milan earthquake may have been triggered exclusively by fluid diffusion, the Harper mainshock may have occurred as a consequence of both fluid diffusion and earthquake-earthquake interactions (foreshocks). While we find no single triggering threshold for all faults, the same assumption of near criticality possibly applies to all the structures responsible for events in southern Kansas. However, Figure 7 shows that some areas are characterized by very low or no seismicity (e.g., within the doughnut-shaped pattern of seismicity at about latitude 37.1°N and longitude 97.9°W), even though the pressure increase is likely higher due to the presence of several nearby disposal wells. This heterogeneity in seismicity is consistent with the hypothesis that the distribution of faults and/or the seismogenic state of the faults in Oklahoma and southern Kansas are spatially variable (Langenbruch et al., 2018;Rubinstein et al., 2018)

Coseismic ΔCFS due to Light-to-Moderate Earthquakes
Several studies have highlighted the importance of coseismic ΔCFS in the spatiotemporal evolution of induced seismicity (Brown & Ge, 2018;Chen et al., 2017;Sumy et al., 2014). For both the M4.6 Harper earthquake and the M4.9 Milan earthquake we find that most of the aftershocks and subsequent nearby events are in regions of positive ΔCFS (Figure 9), confirming the "burst-like" behavior of these two families detected in the cluster analysis. However, ΔCFS very close to the fault plane is difficult to interpret, because it strictly depends on the detailed fault geometry and slip model assumed. Here, the slip models are constructed from the aftershock distributions and empirical relationships that relate magnitude to fault length (Wells & Coppersmith, 1994), and not on geodetic (GPS, InSAR) or seismological data (e.g., waveforms). Therefore, our results may be reliable only for the off-fault regions at more than one fault length from the source faults. If we consider the 0.02 MPa ΔCFS threshold mentioned in the previous paragraph (Reasenberg & Simpson, 1992;Stein, 1999), the Harper mainshock likely affected the timing of a group of events located a few kilometers to the south-west (Figure 9a), while the Milan mainshock may have contributed to the triggering of seismicity up to 7 km from the mainshock (Figure 9b). In both cases, the magnitude of the ΔCFS produced by the two events is larger than the Δp and ΔCFS due to wastewater injections within distances of 2-5 km from the mainshock (Figure 7), in agreement with the findings of Brown and Ge (2018). This points out the importance of coseismic ΔCFS not only in triggering aftershocks but also in affecting nearby seismicity. Note that the ΔCFS showed in Figure 7 are calculated at 4 km depth for the Harper earthquake and at 5.5 km depth for the Milan earthquake; these values represent the average depths of the events located in each of the two areas.
Finally, both events were preceded by nearby, minor earthquakes: a M3.8 event 2 days before the Harper event, and a M3.9 earthquake about two months before the Milan event. Due to the close temporal and spatial proximity between the Harper earthquake and its foreshock, we hypothesize that the M4.6 Harper mainshock may have been encouraged by the occurrence of the M3.8 event. This is also suggested by the cluster analysis where both events are included in the same family. The two events appear to have occurred on the same normal fault (Figure 2b and 2c). The nearest-neighbor analysis results ( Figure 6) also shows that the family that includes the Harper mainshock included several foreshocks. On the other hand, according to our ΔCFS modeling ( Figure S7), the M4.9 Milan mainshock does not appear to have been encouraged by the preceding M3.9 event (ΔCFS < 0.01 MPa). Our cluster analysis also shows no foreshocks in the family that contains the Milan mainshock ( Figure 6). Therefore, we speculate that the Milan earthquake may have been triggered only by the effect of fluid diffusion and related poroelastic stress changes, despite the fact that the model fluid-related stress changes are small. Finally, the large distance between the two events suggests that there is no relation between the occurrence of the Harper earthquake and the Milan earthquake.

Conclusions
Combining physical models (static stress and poroelasticity) and a statistical cluster analysis applied to a high-resolution relocated catalog, we analyze the processes controlling the occurrence of seismicity in southern Kansas. Our results show that: 1. About 35%, 45%, and 20% of the events in the catalog are identified as mainshocks, aftershocks, and foreshocks, respectively, with most of the earthquake families having a mixed "swarm-like" and "burstlike" behavior 2. The absence of foreshocks in the Milan cluster and the short duration of its aftershock sequence in comparison with that of the Harper cluster suggest that the Milan earthquake was possibly triggered by fluid diffusion but subsequent events were driven by earthquake-earthquake interactions 3. The Harper mainshock was possibly triggered by a combination of Δp due to fluid diffusion and static stress changes from a nearby M3.8 foreshock 4. The effect of fluid diffusion on the spatial and temporal evolution of seismic sequences is more prominent on earthquake clusters located closer to major injection wells 5. Δp due to fluid diffusion smaller than 0.1 MPa are enough to initiate seismic sequences when we consider permeability values for the Arbuckle formation of 10 −12 to 10 −13 m 2 , and a conductive 3D volume. However, it is not possible to identify a single triggering threshold value that is consistent within the entire catalog 6. Following each of the two largest events in the catalog (the M4.6 Harper earthquake and the M4.9 Milan earthquake), the coseismic ΔCFS strongly controlled the temporal and spatial distribution of aftershocks as well as the occurrence of nearby clusters

Data Availability Statement
The relocated and the nonrelocated catalogs used in this work are available as supplementary materials (Data Set S1, Data Set S2) or in the Zenodo data repository (https://doi.org/10.5281/zenodo.4318041). Moment tensor solution from the Saint Louis University Earthquake Center is available in http://www.eas.slu. edu/eqc/eqc_mt/MECH.NA/. Injection data are from the Kansas fluid injection database (http://www.kgs. ku.edu/Magellan/Qualified/class2_db.html).