New Estimates of Magnitude‐Frequency Distribution and b‐Value Using Relative Magnitudes for the 2011 Prague, Oklahoma Earthquake Sequence

The magnitude‐frequency distribution (MFD) describes the relative proportion of earthquake magnitudes and provides vital information for seismic hazard assessment. The b‐value, derived from the MFD, is commonly used to estimate the probability that a future earthquake will exceed a specified magnitude threshold. Improved MFD and b‐value estimates are of great importance in the central and eastern United States where high volumes of fluid injection have contributed to a significant rise in seismicity over the last decade. In this study, we recalculate the magnitudes of 8,775 events for the 2011 Prague, Oklahoma sequence using a relative magnitude approach that depends only on waveform data to calculate magnitudes. We also compare the distribution of successive magnitude differences to the MFD and show that a combination of the magnitude difference distribution (MDFD) and relative magnitudes yields a reliable estimate of b‐value. Using the MDFD and relative magnitudes, we examine the temporal and spatial variations in the b‐value and show that b‐value ranges between ∼0.6 and 0.85 during the aftershock sequence for at least 5 months after the M 5.7 mainshock, though areas surrounding the northeast part of the sequence experience higher b‐values (0.7–0.85) than the southwestern part of the Meeker‐Prague fault where b‐value is the lowest (0.6–0.7). We also identify a cluster of off‐fault events with the highest b‐values in the catalog (0.85). These new estimates of MFD and b‐value will contribute to understanding of the relations between induced and tectonic earthquake sequences and promote discussion regarding the use of b‐value in induced seismic hazard estimation.

may be recognized and mitigated during the injection process.Probabilistic seismic hazard assessments (PSHA) are widely used to determine hazard for tectonically active regions.PSHA is a multistep process that combines seismic source characterization, estimations of magnitude recurrence based on magnitude-frequency statistics such as the b-value (detailed in the next section), and estimates of ground motion given a particular magnitude to determine probabilities of ground motion exceedance.Some studies have attempted to apply these methods to induced earthquakes (Petersen et al., 2016(Petersen et al., , 2017) ) but have faced challenges particularly in the case of magnitude recurrence intervals.For induced earthquakes, historical models of seismicity may no longer be valid for use in PSHA due to variations in injection practices and active well locations.Thus, the seismic hazard maps for induced earthquakes can only consider short-term seismicity to be accurate.Therefore, the central United States and similar regions will require methods of hazard assessment that incorporate variability of injection practices.
One solution that attempts to address this need is the implementation of "traffic-light systems" which prescribe thresholds for "acceptable" earthquake magnitudes based on community decision of tolerance levels.When those thresholds are passed, certain action is taken in an attempt to mitigate seismicity.In response to heightened seismicity, the Oklahoma Corporation Commission instituted several hazard mitigation procedures (Murray et al., 2023) including a "traffic-light style" well permitting system (Baker, 2015;Zhang et al., 2020).However, despite multiple mitigation efforts, three large magnitude earthquake sequences occurred in 2016 (M 5.1 Fairview, M 5.0 Cushing, and M 5.8 Pawnee).The occurrence of large magnitude sequences in previously quiescent regions highlights the need for new approaches to understanding induced seismic hazard that may utilize shortterm variations in earthquake behavior and statistical products such as the magnitude-frequency distribution.

Magnitude-Frequency Distribution and the b-Value
The MFD describes the number of earthquakes that occur for each magnitude unit and is commonly modeled by the Gutenberg-Richter (GR) Law (Equation 1) where N m is the number of earthquakes of magnitude M (or all earthquakes with magnitude ≥M for the cumulative distribution) above the magnitude of completeness M c , and a and b are constants that express the total earthquake productivity and the proportion of small to large earthquakes in the given catalog, respectively (Gutenberg & Richter, 1954).The b-value is routinely determined using a maximum-likelihood fit (Aki, 1965) to the MFD.A b-value of one indicates that for every increase in magnitude unit, there is a 10-fold decrease in the number of earthquakes.Deviation from b = 1 implies that the catalog in question experiences a skewed proportion of relatively large or small earthquakes.If the b-value is a steady state estimate of the distribution of magnitude occurrence, then a lower b-value suggests that a region is more likely to experience relatively large earthquakes.
While the global MFD fits the GR Law well with b-value close to 1 (Frohlich & Davis, 1993;Kagan, 1999), local catalogs and individual sequences have shown variations of b-value in a range of approximately 0.5-2.5 (El-Isa & Eaton, 2013;Mousavi et al., 2017, and references therein).
Besides the role of b-value in PSHA, there are various observations that link statistically significant deviations in the b-value to subsurface pore-pressure perturbations (Bachmann et al., 2012;Scholz, 1968;Shelly et al., 2016), shear stress regime or faulting style (Amitrano, 2003;Schorlemmer et al., 2005;Gulia & Weimer, 2010), material properties (Mogi, 1962(Mogi, , 1967)), or varying subsurface pressure and temperature conditions (Wiens, 2001).These observations and their physical interpretations suggest that b-value is inversely correlated with stress.Further, implications of b-value as a proxy for the state of stress on a fault has led others to investigate whether the b-value may exhibit characteristic spatial or temporal variations before and/or after the occurrence of a mainshock (Gulia et al., 2018;Gulia & Weimer, 2019;Main et al., 1989;Smith, 1981;Zuniga & Wyss, 2001).However, most of this body of knowledge is based on observations of naturally occurring tectonic sequences or laboratory experiments.
The variations in b-value due to subsurface fluid injection and their relation to induced seismic hazard are still not well understood.

Magnitude Discrepancies
The variations in MFDs and b-value of induced seismic sequences may provide vital information about the seismic hazard.However, estimation of the MFDs for small earthquakes suffers from uncertainties in earthquake 10.1029/2023JB026455 3 of 17 magnitudes, which can be nontrivial for small events.As a result, any uncertainty in magnitudes used to calculate the MFD will propagate to the b-value and may produce heavily biased results (Herrmann & Marzocchi, 2021).Different methods can be used to determine earthquake magnitude depending on the size of the earthquake, frequency content, and available geologic data.Moment magnitude (M w ) is usually the preferred method to calculate earthquake magnitude as it is an estimate of the seismic moment and tied directly to the source parameters of the earthquake.However, it is challenging to determine M w for small events often due to low signal to noise ratios, thus M w is not routinely calculated for events below M ∼ 3 (Atkinson et al., 2014).
Local magnitude (M L ) is widely used, and is based on the maximum amplitude of the observed waveform after correcting for local attenuation and distance.Ideally M w = M L , but in practice there are often discrepancies between the two methods due to reasons including (but not limited to) rupture complexity, radiation pattern or rupture directivity effects, and assumptions made regarding attenuation and distance correction terms (Deichmann, 2006;Shelly et al., 2016).In the central United States, magnitude differences of more than 0.5 magnitude units have been found during comparison of various magnitude estimates of M 2-4 earthquakes (Shelly et al., 2021).Conversion relationships between different magnitude scales have been developed in an attempt to provide M w estimates for small events and to correct discrepancies between M L and M w , but they can also be biased due to the lack of M w estimates for small events (Deichmann, 2017;Shelly et al., 2021).For this reason, we use relative magnitudes because they do not rely on distance or attenuation correction terms and attempt to reduce radiation pattern and directivity effects by comparing two events measured at the same station with high waveform cross correlation.
In this study, we aim to improve the quality of MFD and b-value estimates for the Prague, OK earthquake sequence.We first reevaluate magnitudes using the relative magnitude method based solely on the amplitude ratios of seismic waveforms.Then, we determine the temporal and spatial variations in the b-value throughout the sequence.We find that the Prague aftershock sequence is characterized by very low b-values between ∼0.6 and 0.8 with temporal and spatial variation.Finally, we discuss the implications of these new b-value estimates and how they may contribute to our understanding of the Prague, OK sequence and the greater picture of induced seismicity and seismic hazard.

Data Sources
The 2011 M 5.7 Prague, Oklahoma earthquake (Figure 1) was, at the time, both the largest earthquake recorded in the state of Oklahoma and the largest earthquake to be attributed to wastewater disposal until the M 5.8 Pawnee, OK sequence in 2016.Fluid injection into the Hunton Limestone and Arbuckle Group from several active wells nearest to the sequence began in 1993.However, nearby seismicity was not instrumentally detected until February of 2010 (Walter et al., 2020).Keranen et al. (2013) asserts that ∼18 years of injection from nearby wells increased the subsurface pore pressure, lowering the effective stress enough to trigger the first large event in the Prague sequence, a M 4.8 foreshock (Figure 1).The stress changes from this event led to cascading ruptures along the previously unmapped Meeker-Prague fault (Cochran et al., 2020;Sumy et al., 2014), including the M 5.7 mainshock which ruptured less than 24 hr after the M 4.8 foreshock.This was followed closely by a M 4.8 aftershock 2 days later (Figure 1).
The Prague sequence was recorded by a temporary deployment of 31 seismometers installed shortly after the M 4.8 foreshock, 7 existing US Transportable Array broadband seismometers, and 9 USGS NetQuakes accelerometers (Sumy et al., 2014(Sumy et al., , 2017)).Cochran et al. (2020) utilized these stations to enhance the spatial resolution and detect small events in the Prague sequence between 6 November 2011 and January 2012.In addition, Skoumal et al. (2020) used multistation template matching to develop a catalog of events throughout the state of Oklahoma from 2009 to 2016.Utilizing the detections from both catalogs allows us to take advantage of the long time period of the Skoumal catalog and high resolution of the Cochran catalog.The combined catalog contains 12,718 unique events occurring from mid-2009 to 2016, though any further analysis in this paper focuses on the 10,812 events that occur between March 2010 and March 2012.Of these events, 1,100 occur in both catalogs, 3,907 are unique to the Skoumal catalog, and 7,711 are unique to the Cochran catalog.For events recorded in both catalogs we preferentially select the Cochran catalog parameters (location, magnitude, etc.) because they are derived from stations closer to the sources.

Magnitude Recalculation
One issue that arises from the use of both catalogs is the difference between M L reported by the Skoumal catalog and M w reported by the Cochran catalog.Skoumal et al. (2020) estimates magnitude using a single-station Richter-scale approach (M L = log 10 [A/A 0 ]) for the template events and a relative amplitude measurement for the detected events.Cochran et al. (2020) also applies template matching to a prior catalog of 900 events from Sumy et al. (2017), where M w was measured using a method for small earthquakes based on pseudospectral acceleration (Atkinson et al., 2014;Novakovic & Atkinson, 2015).We identify 350 events that are recorded in both catalogs with |M L − M w | > 0.5 (Figure 2), the majority of which are events with M < ∼3.
Rather than searching for biases in the existing magnitudes, we choose to recalculate magnitude for all events using a relative magnitude method.A version of this method is commonly used in template matching where the magnitude of a template waveform (M 1 ) and the amplitude ratio (α) between the detected and template waveforms gives the detected event's magnitude (M 2 ) (Equation 2) (Chen et al., 2018;Huang & Beroza, 2015;Lee et al., 2020;Li et al., 2021;Ross et al., 2019).
where c is a scaling constant.This method allows one to calculate magnitudes for detected events, but robust magnitude estimates are still required for the template events.However, if we consider Equation 2 for every pair of events with similar waveforms and not just template-detection pairs, magnitudes may be recalculated by solving a system of linear equations without any prior magnitude estimates and using only the waveform data (Cleveland & Ammon, 2015).
We utilize 10 stations in the nearby region (five rapid deployment stations, four US Transportable Array stations, and one US Geological Survey station) based on location, recording time period, and signal-to-noise ratio (Figure 1; Table S1 in Supporting Information S1).We note that there are more stations available in this region, but we find that while data from more stations would contribute to more robust measurements, the increased computational load outweighs the benefit of adding more stations to our analysis.(Sumy et al., 2017).Focal mechanisms for the three largest events are retrieved from the global CMT catalog (Dziewonski et al., 1981;Ekström et al., 2012) and fault lines are retrieved from the Oklahoma Fault Database (Holland, 2015). 10.1029/2023JB026455 5 of 17 For a particular station, we download waveforms for each event from the Incorporated Research Institutions for Seismology (IRIS) Data Services and apply a bandpass filter of 0.5-15 Hz to the waveforms based on the expected frequency content for events in this magnitude range.We then assess the data quality using two threshold tests that examine the signal to noise ratios (SNRs) of individual waveforms and cross-correlation (CC) coefficients between waveform pairs (Figure 3).The first threshold test discards any waveforms with SNR < 3. Waveforms are then paired with every other waveform recorded on the same channel type.These waveform pairs must then have a CC coefficient above a predetermined threshold.A high CC threshold (>0.8) ensures that only very similar waveform pairs will be used for magnitude calculation, but much of the data will be eliminated.A lower CC threshold will allow for many events to be recalculated, but lowering the CC threshold too much may include event pairs with strongly varying path effects or source-receiver distances.We investigate the effects of varying cross-correlation thresholds on the resulting distribution and find that the MFD is stable for CC thresholds between 0.55 and 0.75 (Figure S1 in Supporting Information S1).Therefore, we choose to use CC = 0.6 as it allows for recalculation of a large number of events while still maintaining the same distribution shape as it would have with slightly higher CC values.
We then calculate the amplitude ratio for every passing waveform pair.Some methods estimate amplitude ratio from peak waveform amplitudes (Huang & Beroza, 2015) or the ratio of L2 norms (Schaff & Richards, 2014), which are both easily implemented but may not perform well for waveforms with low SNR.Singular value decomposition (Rubinstein & Ellsworth, 2010) is another avenue for amplitude ratio calculation but requires very coherent waveforms.
We choose to use a principal component fit as this method can accommodate waveforms with low cross correlation (Shelly et al., 2016) and has been used successfully in similar relative magnitude methods (Chen et al., 2018).The principal component method uses the covariance matrix between two waveforms to determine the largest eigenvalue and its associated eigenvector.The elements of the largest eigenvector (v 1 and v 2 ) then give the amplitude ratio (α) (Equation 3).
For a given pair of events, we calculate the amplitude ratio for each channel of the three-component seismogram that have passed the SNR and CC tests and then take the mean.Once we obtain amplitude ratios from a single station, we repeat this process at other stations and take the mean value for event pairs that have results from more than one station.We then use the amplitude ratios to estimate relative magnitude by expanding Equation 2 into a set of linear equations representing each pair of highly correlated events (Equation 4).We assume that every event's magnitude is unknown and determine the "best fit" magnitudes for the linear system using a least squares approach where α i,j is the amplitude ratio between events i and j, and M n is the best fit relative magnitude for the nth event.
Figure 3. Relative magnitude method workflow schematic.W i represents waveform i that must have SNR > 3 to pass through the first data quality threshold.Then pair P j must have a cross-correlation coefficient above a set threshold to pass the second data quality test.Amplitude ratio is then calculated for every passing pair and aggregated with data from other stations.Calibration events may also be considered for absolute magnitude calculation during the least squares inversion. 10.1029/2023JB026455 6 of 17 An optimal value must also be determined for the magnitude scaling constant c.Some prior studies choose c = 1 (Cleveland & Ammon, 2015;Huang & Beroza, 2015) assuming that the logarithm of waveform amplitude scales directly with local magnitude, that is, a 10-fold increase in peak amplitude corresponds to a magnitude unit increase of 1.However, for small events it has been suggested that the amplitude ratio should follow the scaling of seismic moment to magnitude wherein c = ⅔ (Hanks & Boore, 1984).Previous studies have also suggested the use of an intermediate value between ⅔ and 1 when calculating relative magnitudes of different types (Bakun, 1984;Hainzl & Fischer, 2002;Shelly et al., 2016).Since our method relies purely on comparing waveform amplitudes, we chose c = 1 for this analysis, though we note that this is an avenue for future investigation.
Equation 4 yields a set of values where the differences between them reflect the relative differences in magnitude and the average of the values ≈ 0. In some applications it may be necessary to calibrate the relative magnitudes to an absolute scale.Previous methods use the total moment of the originally cataloged magnitudes to calibrate the relative magnitudes (Chen et al., 2018;Cleveland & Ammon, 2015), but since we assume here that all the original magnitudes are unknown, we must also assume that the total moment derived from those magnitudes is unknown as well.In this study, we choose not to calibrate our magnitudes to an absolute, standard magnitude scale because the absolute value does not affect the shape of the MFD or b-value, though we modify Equation 4to include individual calibration events and discuss the process of scaling relative magnitudes to a standard scale in the supplementary material (Text S3 and Figure S2 in Supporting Information S1).It is worth noting that when using the relative magnitude method, it is possible for clusters of events to form with high CC within the cluster but no high CC links between clusters, which would cause each to be scaled so that the individual cluster's average magnitude is zero.While this is less likely to occur in the case of an individual sequence where events are located close together, we also set the cross-correlation threshold low enough to capture as many linked pairs as possible and ensure that there are no suspicious peaks of events centered near zero.

MFD and the b-Value
Once the relative magnitudes have been determined, we can estimate the b-value for the Prague sequence.We compare two different methods for estimating b-value for our relative magnitudes as well as the originally cataloged magnitudes.The traditional b-value method uses a maximum-likelihood fit to the noncumulative MFD (Equation 5) (Aki, 1965;Bender, 1983): where   is the sample mean with M ≥ M c , ∆M is the binning width, and M c is the magnitude of completeness.We determine M c from the maximum curvature method (Weimer & Wyss, 2000) and add 0.2 as recommended by Woessner and Weimer (2005).
The second method follows van der Elst (2021) who shows that the distribution of magnitude differences between successive earthquakes follows a Laplace distribution that can also be used to estimate the b-value (Equation 6): where   ′ is the sample mean of the positive magnitude differences where   ′ ≥  ′  and   ′  is a minimum magnitude difference threshold that we set to be 0.2 in accordance with van der Elst (2021).However, we do note that an absence of small magnitude events may affect the stability of the b + estimate.This is discussed further in supplementary Section 4 (Text S4 and Figures S3-S5 in Supporting Information S1).While the b GR method is more widely used, the b + value derived from the newer magnitude difference method has been shown to be more robust to short-term aftershock incompleteness (van der Elst, 2021).This method is also well suited for our relative magnitude measurements that are designed to quantify magnitude differences rather than the absolute magnitude.To estimate the uncertainty of our measurements we perform bootstrap resampling with 2,500 iterations.
A singular b-value can give valuable information about the characteristics of a certain region.However, the temporal and spatial variations in MFD and b-value are of potentially more importance when analyzing the characteristics and evolution of a specific sequence.Thus, we investigate the spatial and temporal variations in b-value.It is advantageous to use a window with a constant number of events rather than a constant time period so that the uncertainty of the b-value due to sample size does not change throughout the analysis.We evaluate the temporal variations in b-value using a moving window of 400 events.For our analysis of spatial variations, we calculate b-value at each hypocentral location using the closest 400 events as long as the farthest event does not exceed 5 km.

Relative Magnitude Method Benchmark Test
We first evaluate the quality of magnitude estimates produced by the relative magnitude method by benchmarking the results of our method with coda-envelope derived magnitudes for events in Oklahoma and southern Kansas (Shelly et al., 2021).Note, these events are not part of the Prague sequence and are a separate data set used for methodological benchmarking.Coda-envelope magnitudes are considered a stable method of estimating earthquake magnitude due to the superior averaging properties of coda waves compared to direct phase waves (Mayeda & Walter, 1996).We apply the relative magnitude method using five stations (O2.GORE, GS.KAN08, GS.KAN13, GS.KAN14, and GS.OK032) previously utilized by Shelly et al. (2021) and recalculate 683 out of 760 events using a cross-correlation coefficient of 0.55 (Figure 4).To compare our results with the cataloged coda-derived moment magnitudes (M coda ), we use a scaling constant of ⅔ and calibrate them by adding a constant equal to the mean of cataloged M coda (2.658) because we consider these magnitudes to be reliable.
We find that the differences between our calibrated relative magnitudes and M coda are less than 0.5 for 98.9% of the events used for benchmarking with a maximum magnitude difference of 0.66.Additionally, we compare our absolute relative magnitudes from this data set to 337 events with M L cataloged by the Oklahoma Geological Survey (OGS) again, we find very good agreement where 80.7% of the relative magnitudes have a difference <0.5 from the OGS magnitudes with a maximum magnitude difference of 0.97.If we compare M L and M w from the Skoumal and Cochran catalogs for events with M > 1.75 (following the Shelly benchmark events), 20.8% of those events have a magnitude difference greater than 0.5 with a maximum magnitude difference of 1.57.These results show that the relative magnitude method can produce magnitude estimates that are in good agreement with widely accepted and robust forms of magnitude calculations and has the potential to improve upon the magnitude estimation for the Prague sequence, though we note that the coda-envelope catalog does not include any M < ∼1.5 events where the majority of discrepancy between the Cochran and Skoumal catalogs take place.

Prague Sequence Magnitude Recalculation
After applying SNR and CC threshold tests to all waveforms and waveform pairs, we can recalculate 81% of the events in the combined catalog, all of which occurred between March of 2010 and March of 2012.As mentioned previously, the reader should note that while we do provide an absolute number for our relative magnitudes (M rel ), these are not scaled to the "standard" magnitude scale.Instead, the least squares method ensures that the mean of the relative magnitudes is equal to zero.Thus, M rel is not intended to be equal to M w or M L and our results focus on the shape of the MFD and b-value rather than the absolute value of the magnitudes.In our new catalog of recalculated M rel we observe many events that form a peak in the MFD between M rel = ∼−1-0.This is unrelated to the aforementioned clustering problem because the peak is not centered on M rel = 0. Therefore, we suspect that these small magnitude event detections may represent events where the waveform amplitudes are contaminated by noise.Furthermore, these events are likely below the true M C so we elect to remove them from further analysis.
We determine an overall b-value using both the maximum-likelihood method (b GR ) (Aki, 1965) and the b + method (van der Elst, 2021) for all three catalogs (Skoumal, Cochran, and relative magnitude catalogs) and find that b + is slightly higher than b GR for all three cases (Figures 5a, 5c, and 5e).The difference between b GR and b + for the same catalog is likely due to the dependence of b GR on M C .We can observe curvature in the MFDs of the Cochran catalog as well as our own relative magnitude catalog which is visible in Figure 5 and quantitatively shown by b-value stability curves in Figure S3 in Supporting Information S1.Curvature makes estimation of the b-value difficult because we assume that the MFD follows a power law distribution and attempt to fit a linear trend to gradually curved data.
One explanation for gradual curvature in the relative magnitude MFD is biased sampling of events in the relative magnitude process since ∼19% of events are excluded due to not meeting the SNR or CC thresholds.Additionally, our high cross-correlation threshold may exclude potential false-positive detections in the template-matched catalogs leading to an apparent lack of small magnitude events.Biased sampling is potentially concerning for MFDs because the b-value will be affected if incompleteness is unequally distributed across the full magnitude range.However, our results show that the MFD of the entire catalog is not significantly different from the MFD of the catalog with only the recalculated events (Text S5 and Figure S6 in Supporting Information S1), and there is a similar amount of curvature in the Cochran and relative magnitude MFDs (Figure S3 in Supporting Information S1), so it is unlikely that the curvature is caused entirely by sampling bias in this case.Another possible explanation is that the gradual curvature of the MFD is due to short-term aftershock incompleteness (STAI).In a catalog of aftershocks, large events may mask the signal of a smaller successive event, causing an artificially low number of small event detections.This in turn may lead to errors in calculating M C and underestimation of b-value.However, a method that only utilizes event pairs where the first earthquake is smaller than the successive event, as the b + method does, will be more robust to the effects of STAI.Therefore, we choose to use b + for further analysis of the spatial and temporal variations in b-value.For each spatial and temporal window, we choose a minimum magnitude difference of 0.2 following the procedure in van der Elst et al., 2021 though we do recognize that the choice of minimum magnitude difference may influence the results of b + when there is a deficiency of small magnitude events.We discuss this in the supplementary information (Text S4 in Supporting Information S1).

Temporal Variations in b +
A single estimate of b-value can give an idea about the historical magnitude distribution for a region but the temporal variations in the MFD and b-value will give more information about how a specific sequence evolves over time.Thus, we evaluate the temporal variations in b + for the Prague sequence using a window of 400 events that moves for every new event.We plot the b + as a function of the date for the final event in each window (Figure 6).

Spatial Variations in b-Value
We present the spatial variations in b + for the Prague sequence aftershocks to determine which locations may have been at a higher risk for relatively large events (Figure 7).The areas in the southwest extent of the aftershock sequence experienced lower b-values than the locations in the northeast where the M 5.7 mainshock and the M 4.8 foreshock had previously occurred.Overall, b + decreases with distance from the mainshock rupture area.When considering depth (Figure 7, bottom and right) we note that events in the basement rock below ∼5 km depth have low b-value similar to the events that occur on the farthest southwest section of the Meeker-Prague fault segment.For events in the upper sedimentary layers b + is consistent with distance to the mainshock rather than depth.
We also note a cluster of 328 events that occur southeast of the main faults that were activated during the Prague sequence (Figure 7, "Southeast Cluster").These events occurred from 1 day after the mainshock to 14 February 2012 and only one other event occurred in this location before the Prague Sequence.While the b + for this cluster is not calculated in Figure 7 due to having less than 400 events within 5 km of each earthquake in the cluster, we calculate the b-value of this group of aftershocks separately and find b + = 0.85 ± 0.09 which is comparable to the value closest to the mainshock location.

Relative Magnitudes
When benchmarking the relative magnitude method with the M coda results of Shelly et al. (2021) and OGS local magnitudes, we find that the vast majority of M rel estimates are less than 0.5 magnitude units different from both M coda and M L,OGS , suggesting that a relative magnitude recalculation is a viable way to quantify the relationships between earthquake magnitudes.This method is similar to the estimation of magnitude in template-matching studies, but the main difference is that we require a higher cross correlation and thus, more similar waveforms to constrain amplitude ratio than what is routinely used in template matching.Furthermore, we pair each event with every other event whereas template matching only measures the detected event in reference to its own template.
The greatest benefit of the relative magnitude method is that it relies only on cross correlation between event pairs measured at the same station.This eliminates the need to determine path or station corrections that require certain assumptions of prior geologic knowledge which may not always be readily available or accurate.We are also able to compute reliable relative magnitudes for a small number of stations since rupture directivity and radiation pattern do not need to be averaged out over many station measurements.Relative magnitudes can be measured from as little as a single station, though we recommend multiple stations at different azimuths to maximize the number of high cross-correlation pairs.However, the relative magnitude method has limitations related to the necessary sample size and computational cost.An optimal data set for this method would be a large or high-quality data set with closely located events as well as multiple stations with sufficient signal-to-noise ratios and azimuthal coverage.In the central and eastern United States, this kind of high-quality data set is not always achievable due to sparse network coverage, though this can be achieved with template-matched catalogs.Conversely, if the data set is very large (>∼20,000 events), significant allocation of computational resources is necessary to process the amplitude ratio calculations and matrix inversions.The large number of densely located events in the Prague sequence are a good trial data set for this method, but there have been difficulties regarding station coverage.The addition of rapid deployment stations allows for extensive data collection and spatial resolution, but they only record the period between the first large foreshock (M 4.8, 5 November 2011) and March of 2012, making magnitude recalculation difficult both before and after this time period.
Another potential drawback for the relative magnitude method is the exclusion of events that have either low SNR or poor correlation.Low SNR may be more likely to affect small or distant events that are at or below the noise threshold of a particular station.Poor cross correlation may affect waveform pairs of any magnitude but may particularly impact moderate and large events due to increased waveform complexity and a lack of similar magnitude events (Shelly, 2016).However, the exclusion of events due to the CC threshold may be advantageous in some ways as it would eliminate potential false detections from template-matched catalogs where the detection threshold is often lower than what is required in this study.
In our case ∼19% of events are excluded due to a combination of low SNR and poor correlation over the entire magnitude range.For the Prague catalog we investigate the effects of the missing events on b-value using the cataloged magnitudes (Text S5 and Figure S6 in Supporting Information S1) and find that the MFD is similar between the entire catalog of events and the recalculated events.However, many events are still excluded, and we add a note of caution for future research that attempts to interpret spatial and temporal variations in b-value (b GR or b + ) of a relative magnitude catalog if the incompleteness is not continuously uniform across space and/or time.
Future work to improve the relative magnitude method may include focusing on specific segments of the waveform as opposed to using the entire waveform to determine amplitude ratios.Previous studies (Mayeda et al., 2003;Mayeda & Walter, 1996) note that coda waves tend to average over fine scale path structure and represent source properties better than direct phases so the use of coda waves to determine relative amplitudes may allow for better amplitude ratio accuracy for pairs with lower cross correlation.Additionally, the chosen frequency band (0.5-20 Hz) may lead to amplitude saturation for the largest events in the sequence and subsequently, a possible underestimation of magnitude difference for event pairs that include the largest events.However, we note that most events in this sequence are of generally small magnitude (<3) and are likely not affected by amplitude saturation.For a catalog with a larger magnitude range, we would recommend further exploration of the effects of measuring amplitude at different frequency bands.

b + for the Full Catalog
b-Values reported for other induced sequences show a very wide range of values (0.6-2.89).Very high b-values are usually associated with hydraulic fracturing in which new fractures are created in rock with low tensile strength, and low b-values are associated with wastewater disposal and enhanced geothermal systems where previously existing tectonic faults are reactivated (Goebel et al., 2016;Maxwell et al., 2009;Mousavi et al., 2017, and references therein).
Overall, the b + values calculated from our relative magnitudes are very low compared to many other induced and tectonic sequences, but our results are consistent with McMahon et al.'s (2017) results for the Prague sequence and the fact that b-values in intraplate cratonic regions are generally less than one (Okal & Sweet, 2007).In addition, prior studies have identified low b-values in the range of 0.6-0.9 during periods of injection activity (Goebel et al., 2016;Huang & Beroza, 2015).The role of specific injection parameters on the b-value for the Prague sequences was not investigated in this study but is a potential avenue for investigation of the Prague sequence.
It is possible that our choice of minimum magnitude difference (0.2) has caused a slight underestimation of b + (discussed in supplements, Text S4 and Figures S3-S5 in Supporting Information S1); however, our results are within a range of reasonable b-values for a sequence that, although initiated by pore-pressure perturbations, is generally thought to have released preexisting tectonic stress.

Temporal and Spatial Variations in b +
Using our new catalog of relative magnitudes, we identify temporal and spatial variation in b + within the range of ∼0.6-0.8 (with standard deviation in a range of ∼0.035-0.075)for the time period immediately after the M 4.8 foreshock until ∼5 months after the Prague mainshock (when we are no longer able to calculate b-value due to data availability).Because of a lack of station coverage, we are not able to establish a robust background b-value or identify any temporal variations in b + that may occur before the initiation of the Prague sequence.After the mainshock we observe a rise in b + similar to prior studies of temporal Gutenberg-Richter b-value variations (Gulia & Weimer, 2019;Gulia et al., 2015;Tormann et al., 2015;van der Elst, 2021;Walter et al., 2017;Weimer et al., 2002).In fact, the post-mainshock temporal variations are very similar to those described by Gulia et al. (2018) who suggest that b + increases immediately after a mainshock until reaching a peak value approximately 1-2 months later.In our relative magnitude catalog b + reaches peaks on 10 November and on 30 December 2011, 10 days and almost 2 months, respectively, after the Prague mainshock.In comparison with b + calculated from the original Skoumal and Cochran catalogs (Figure S7 in Supporting Information S1), all three catalog's b + values reach the minimum sometime between the M 5.7 mainshock and the M 4.8 aftershock, and experience a general increase over the next ∼10 days.During this time period, our b + values are similar but less oscillatory than the Cochran catalog.After ∼10 days, the Cochran catalog shows a drop in b + which is not observed by our catalog or the Skoumal catalog.In addition to the temporal variations of b + , we observe a decrease in b + with distance from the area where the M 4.8 foreshock and M 5.7 mainshock ruptured where b + is the highest.This decrease in b + as a function of distance from the mainshock is also in line with Gulia et al.'s (2018) observations of tectonic earthquake sequences.
There is some debate about whether the size of induced earthquakes is related solely to the total volume of injected fluid (McGarr, 2014), or whether induced earthquakes are initially triggered by injection but ultimately release tectonic stress already existing on the fault.The Prague sequence is suggested to follow the latter (Keranen et al., 2013;Sumy et al., 2014), though the total moment released during the Prague sequence is still within the limit predicted by McGarr and Barbour (2017)'s relation of seismic moment and injected volume.However, whether an induced earthquake sequence follows a specific pattern of b-value variation may be related to both the pressure perturbations caused by injection and the preexisting tectonic stresses.The low b + values in the southwest portion and high b + in the vicinity of the mainshock could be due to missing detections of small events during the aftershock sequence, but we mitigate the effects of STAI as much as possible by recalculating magnitudes and using the b + method.If we consider b-value (b + in this case) as a stress proxy where b + is inversely correlated with stress, then the distribution of b + indicates that the tectonic stress in the northeast part of the fault was relaxed by the M 4.8 foreshock and the M 5.7 mainshock while stress perturbations still existed on the southwest portion of the Meeker-Prague fault, leading to a higher proportion of relatively large events, though a further analysis of Coulomb stress would be required to confirm this.In addition, b-values may also be influenced by the fault heterogeneity.Long fault segments that are able to host larger ruptures may produce more large events, and hence lower b-values as opposed to the smaller, compartmentalized fault segments along the main Wilzetta fault system in the eastern part of the sequence.McMahon et al. (2017) also finds much higher b-value for shallow events (1.05) than deep events that occur in the crystalline basement (0.52).We observe low b + for many of the deeper events.However, in our results the depth distribution of high versus low b + appears to be more related to the location of the mainshock and large foreshock, rather than to depth, although we observe a jump in b + near 5-6 km depth.
The smaller cluster of events that occurred southeast of the Wilzetta/Meeker-Prague fault system is also of interest.Almost all of these events occur approximately 1 day to 2 months after the mainshock with ∼two-thirds of the events occurring within the first 10 days after the mainshock.This temporal correlation with the mainshock suggests that stress transfer from the main sequence may have triggered this small cluster, likely on a small, unmapped fault segment nearby.The b-value for this cluster is less than one (0.85 ± 0.09), and despite the overall low b + for the sequence as a whole, this cluster is one of the highest among the spatial groups.It is possible that the size of the fault segment may be limiting the size and distribution of earthquakes hosted on that fault.

Implications for Induced Seismicity and Seismic Hazard
The ability to identify statistical indicators of heightened probability of earthquake rupture that can be monitored in near real-time would be incredibly important for induced earthquakes, as it could provide an avenue for adjusting injection activity to mitigate seismicity.Since we are unable to establish background b-value for the Prague region, we cannot quantify a percent change of b-value after the mainshock as prior studies have done (Gulia & Weimer, 2019;Gulia et al., 2018).However, our results indicate that the Prague sequence, despite being initially triggered by pore-pressure perturbations caused by fluid injection, shows distributions of b-value variations that are similar to tectonic sequences.Thus, we suggest that the same assumptions regarding the use of b-value variations as earthquake precursors may also apply to induced sequences where the main driver of seismicity is relaxation of tectonic stress initiated by fluid injection.Due to the variability of the effects of injection on earthquake activity, we suggest that more induced sequences in a variety of stress conditions and fault structures should be investigated further and compared with tectonic sequences in similar geologic settings, if possible, to understand whether injection activity correlates with short-term variations in b-value.
The regions of lowest b-value during the aftershock sequence can be considered more hazardous in terms of earthquake size distribution, but it is also important to distinguish between precursory changes that may indicate a "large" rupture, the maximum magnitude that may be produced, and the actual ground shaking and damage caused by an earthquake.In this study, we show that the Prague sequence follows similar b-value variations as tectonic sequences and posit that the spatial and temporal variations in b-value may allow us to understand where stress has been relaxed on the fault system and where stress concentrations may persist ∼5 months after the mainshock.However, our findings alone cannot lead to any assumptions about the maximum magnitude of a potential future rupture or the amount of damage that may be caused.

Conclusions
As earthquake rates in the central and eastern United States remain high mainly due to wastewater disposal injection and hydraulic fracturing, improved knowledge of how these sequences evolve is necessary to mitigate damage caused by induced earthquakes.In our study we attempt to further understand the 2011 Prague, Oklahoma sequence through the magnitude-frequency distribution and spatiotemporal variations in the b-values of the aftershocks.We first recalculate earthquake magnitude using a relative amplitudes-based method that allows us to take advantage of two high-resolution template-matched catalogs with different magnitude scales.We recalculate a majority of the events that occurred between March 2010 and March 2012, and find that the relative magnitudes combined with a new b-value estimation method that utilizes successive magnitude differences provides reliable estimates of b-value.
Using this method, we examine the temporal and spatial variations of the b-value and identify a ∼5 month period of low b-values (in reference to the globally expected b-value of 1) occurring during and after the Prague mainshock, indicating a heightened probability of relatively large events particularly during the aftershock sequence in reference to the global average b-value of one.We also identify specific subregions of the Meeker-Prague area that exhibit very low b-values where tectonic stresses may not have been sufficiently relaxed by the three M > 4 earthquakes and are more likely to host earthquakes in the future.While b-values are limited to only providing information on the frequency distribution and probability of events of a given magnitude as opposed to predicting damage potential, these new relative magnitudes and spatiotemporal variations in the MFD are important steps to understanding the statistics and our hazard assessment capabilities of induced seismicity.

Figure 2 .
Figure2.M L versus M w for events recorded by the Skoumal and Cochran catalogs respectively.Dotted gray lines represent a margin of ±0.5 magnitude units while the solid gray line represents a 1-to-1 relationship between M L and M w .The inset plot shows the distribution of differences between the magnitude types.

Figure 1 .
Figure 1.Earthquakes that occurred in the Meeker-Prague region between March of 2010 and March of 2012 with the M 4.8 foreshock, M 5.7 mainshock, and M 4.8 aftershock highlighted as yellow stars.Open circles represent earthquakes with size scaled to magnitude.Green squares represent injection wells scaled by cumulative volume injected until the end of 2010(Sumy et al., 2017).Focal mechanisms for the three largest events are retrieved from the global CMT catalog(Dziewonski et al., 1981;Ekström et al., 2012) and fault lines are retrieved from the Oklahoma Fault Database(Holland, 2015).

Figure 4 .
Figure 4. Left: Coda-envelope derived magnitudes (M coda ) for Oklahoma and southern Kansas compared with uncalibrated relative magnitudes (light purple) and relative magnitudes calibrated by the mean of cataloged M coda (dark purple).Right: OGS local magnitudes compared with uncalibrated relative magnitude (light blue) and relative magnitudes calibrated by the mean of cataloged M L, OGS (dark blue).

Figure 5 .
Figure 5. Left column (a, c, and e): Noncumulative MFDs (light colors) and MDFDs (dark colors) for the Skoumal catalog (turquoise/teal), Cochran catalog (pink/purple), and the relative magnitude catalog (orange/red) for events from March 2010 to March 2012.M c is expressed by gray dashed lines.Right column (b, d, and f): Magnitude as a function of event number for each catalog.Light gray portions of the catalog in (e, f) represent recalculated events that are excluded from further analysis.

Figure 6 .
Figure 6.b + as a function of event date for the Relative Magnitude Catalog.Dotted lines represent the event times of the M 4.8 foreshock, M 5.7 mainshock, and M 4.8 aftershock.Inset: b + as a function of event date for the first 10 days of the sequence.Gray shaded region for both the main and inset plots represents ± one standard deviation.

Figure 7 .
Figure 7. Map of b + at each aftershock hypocenter using the closest 400 events.Gray markers are events for which b + is not calculated because there are less than 400 events within 5 km of the hypocenter.Bottom and right figures express spatial b + variations with depth.While this map only calculates b + for aftershocks, we have included the locations of the M 4.8 foreshock, M 5.7 mainshock, and M 4.8 aftershock as black stars for reference.