Assessing Earthquake Rates and b -value given Spatiotemporal Variation in Catalog Completeness: Application to Atlantic Canada

Spatiotemporal variations in the magnitude of completeness M c make it challenging to confidently assess seismic hazard or even to simply compare earthquake rates between regions. In this study, we introduce new techniques to correct for heterogeneous M c in a treatment of the eastern and Atlantic Canada earthquake catalog (1985–2023). We first introduce new methodology to predict M c ( x, t ) based on the distribution of seismometers. Second, we introduce a modified maximum-likelihood estimator (MLE) for b (the b -value) that accounts for spatiotemporal M c variation, allowing the inclusion of more earthquakes. Third, we compute the ratio of detected/predicted M> 1 earthquakes as a function of M c and apply it to create a calibrated M> 1 event-rate map. The resulting map has advantages over a moment-rate map, which is effectively sensitive only to the very largest earthquakes in the dataset. The new MLE results in a modestly more precise b when applied to the Charlevoix Seismic Zone, and a substantial increase in precision when applied to the full Atlantic Canada region. It may prove useful in future hazard assessments, particularly of regions with highly heterogeneous M c and relatively sparse catalogs.


Introduction
The rate of earthquake occurrence in a given region is generally reported as either an event rate, a moment rate, or, most formally, with a Gutenberg-Richter (GR) model (Ishimoto and Iida, 1939;Gutenberg and Richter, 1944).An event rate is the simplest way to communicate earthquake density, it is the number of earthquakes per unit time, generally considering only those above some threshold magnitude.A moment rate sums the seismic moment of all earthquakes in the region, and is effectively only sensitive to the largest earthquakes in the region.GR models express the number of earthquakes as a function of magnitude N (M ) through a log-linear relation log 10 N (M ) = a − bM , where the constant a describes the overall abundance of earthquakes and b (the b-value) describes the relative abundance of small earthquakes to large ones, and typically b ≈ 1.This is a central equation to probabilistic seismic hazard assessments, and reliable estimates of a and b are therefore critical to seismic hazard analysis.
The magnitude-frequency distribution (MFD) in actual earthquake catalogs always deviates from the strict log-linear model.They are often characterized by a double-truncated GR model, which has an upper magnitude limit M max based on the maximum fault size in the region, as well as a magnitude of completeness M c , below which there will be fewer earthquakes detected than predicted because of our limited ability to detect them.M c can be affected by noise conditions as well as geological factors that affect seismic attenuation, but it primarily depends on the distribution of seismometers in the region.High M c increases uncertainty in hazard assessments and complicates even event-rate estimates.Knowledge of background rates of naturally-occurring earthquakes is critical to the responsible management of any activity that can pose risk of induced seismicity, such as hydraulic fracturing and wastewater injection in the oil and gas industry, as well as geological carbon storage (Schultz et al., 2020;Cheng et al., 2023).This paper introduces new methodologies to calculate b and compare event rates across regions that account for spatiotemporal M c variations.It will focus on eastern and Atlantic Canada, which generally has low to moderate levels of seismicity and seismometer coverage, but does contain several prominent onshore and offshore seismic zones.

Data and regional seismicity
This study examines the Canadian National Earthquake Database (CNED, see section "Data and code availability") from 1985-2023.We selected this timespan for two reasons: i) 1985 is the first year for which the catalog is publicly available, and ii) it is challenging to determine what seismometers were operational and used in routine earthquake detection at any given time before approximately the mid-1990s.Seismometer locations and operating dates were downloaded from the IRIS database, and supplemented with an annual publication on the network that ended in the late 1980s (Munro et al., 1988).Note that these sources indicate when stations became active, but not periods of time they were nonoperational.A map of epicenters and a magnitudetime plot are shown in Figure 1, and Figure 2a shows the location of all stations that may have contributed to the catalog.
The catalog contains 13612 earthquakes, with a variety of magnitude types.We follow the procedure of Halchuk et al. (2015) to convert all magnitudes to M W ; the method is partly based on the m N -M W relation computed by Bent (2011), but does not follow it directly.We note that the linear scaling suggested by Bent (2011) would directly affect b estimates; more details on the treatment of magnitudes is found in Appendix 1.The magnitude-time plot (Figure 1b) shows a drastic change around 1995, when the seismograph network and detection routines underwent several major changes (Bent, 2011).It also has notably few events in the range of approximately 0 < M W < 1 (the range is higher pre-1995 than post-1995), whereas there are lobes of relatively abundant events both above and below.The lowermagnitude lobe consists mainly of onshore events reported with only a local magnitude M L , for which the conversion to M W may be more dubious.Proper scaling of onshore M L to M W may warrant future investigation, but for this study we opt to simply ignore M W < 0 for subsequent analyses, removing the vast majority of these onshore M L events.
Following the approach taken in Canadian Seismic Hazard models (Kolaj et al., 2020), we do not attempt to decluster the earthquake catalog.Declustering is often performed as part of probabilistic seismic hazard analysis to remove foreshocks and aftershocks, such that the remaining earthquakes can be considered a Poisson process in time (Gerstenberger et al., 2020).However, declustering techniques require arbitrary thresholds to define what constitutes a foreshock or aftershock, and can cause unintended bias of b and hazard level estimates (Gerstenberger et al., 2020;Mizrahi et al., 2021).We also assume, given the vast geographic scale and generally low event rates and high M c of this dataset, that temporal variation of M c due to short-term aftershock incompleteness (Stallone and Falcone, 2021; van der Elst, 2021) is not a significant issue.Map boundaries were chosen in order to include the Charlevoix Seismic Zone (CSZ, Lamontagne et al., 2003a;Yu et al., 2016 ) and Lower St. Lawrence Seismic Zone (LSZ, Lamontagne et al., 2003b;Plourde and Nedimović, 2021) in the west, as well as the less-studied seismic zones at the Laurentian slope and fan in the south (Adams and Basham, 1989;Bent, 1995), and in the the Labrador Sea to the north (Bent and Hasegawa, 1992;Bent and Voss, 2022).These areas are labeled in the earthquake density map in Figure 2b.The map includes all M W >1 earthquakes and is significantly affected by spatial M c variations.The CSZ and LSZ stand out as the most earthquake-dense regions on the map, with another prominent area in northern New Brunswick around the epicenter of the 1982 Miramichi M 5.7 event (Wetmiller et al., 1984).
Different trends emerge if we map moment-density of the same earthquake catalog (Figure 2c).The CSZ and LSZ are still prominent, but they have much more similar moment-rates to the offshore seismic zones than in the event-rate map.Spatial M c variation biases moment-rates much less than event-rates, but the trade-off is that moment-rate is controlled almost entirely by the largest, infrequent events, so it is more highly affected by the limited catalog duration.The highest mapped moment-rate results from the largest event in the catalog, the M W 6.3 Ungava Bay earthquake of December 1989 (Bent, 1994), which falls in the northwestern corner of the map.The anomaly just east of it is a M W 5.7 that occurred in March 1989.Another prominent anomaly lies just north of the CSZ, and results from the second-largest earthquake in the dataset, the 1988 Saguenay M W 5.9 (Haddon, 1995).

Earthquake density mapping in Atlantic Canada
We assume that the MFD in Atlantic Canada follows a GR model allowing us to predict the distribution of undetected small earthquakes based on the distribution of larger earthquakes.This requires i) knowledge of M c as a function of both space and time, i.e.M c (x, t), ii) a method for estimating b given spatiotemporal variations in M c , and iii) a function to predict the ratio of undetected to detected earthquakes for a given M c (x, t).
New methodologies to address these three issues are presented in the following subsections, along with results from their application to Atlantic Canada.

Estimating M c (x, t)
In practice, it is challenging to estimate M c (x, t) precisely.M c is typically estimated from GR plots either by visual inspection or using any number of algorithms (several popular ones are reviewed by Woessner and Wiemer, 2005), but that requires many earthquakes (typically hundreds or more).We resolve this by using a predictive M c (x, t) model based on the distribution of seismometers (Mignan et al., 2011), which uses the distance to the n th closest seismometer d n = d n (x, t).Mignan et al. (2011) produces empirical powerlaw models of the form , where c i are constants, based on a dense earthquake catalog from Taiwan.They produce models for distance to the 3 rd -, 4 th -, and 5 th -nearest seismometer, each with similarly close fits.Here we form a similar model for our dataset, but to partially reduce the effect of temporarily nonoperational stations we define a (admittedly arbitrary) weighted station distance metric as a weighted sum of the distances to the 4 th -6 th nearest stations: (1) d(x, t) = 0.70d 4 + 0.25d 5 + 0.05d 6 .
We evaluate d(x, t) for each earthquake in the dataset using the distribution of seismometers that were active when it occurred.
Earthquake density in our Atlantic Canada catalog is orders of magnitude lower than the Taiwan catalog used by Mignan et al. (2011)-it has about one tenth the number of earthquakes in a region 20 times larger.As such, there are few areas with enough earthquakes in a small radius (e.g.50 km) to reliably estimate a, b, or M c .We therefore need an alternative way to fit the powerlaw model, and take the following approach.We first compute d(x, t), as defined above, for each earthquake given its origin time and epicentre.We then sort the earthquakes by their d(x, t) and bin into groups of 300 with 50% overlap, resulting in 78 groups with maximum d(x, t), or d max , of 21 to 1310 km.For each group of 300 earthquakes, we apply the method of Ogata and Katsura (1993) to estimate M c .The method assumes the number detected earthquakes of a given magnitude N (M ) depends on the actual number N 0 (M ) and a "thinning" function q: (2) The thinning function is assumed to be a cumulative normal distribution function: where µ is the magnitude at which 50% are detected and σ describes the width of the thinning function; a low σ indicates a steep falloff in detection below M c , whereas a high σ indicates a more gradual falloff.We use the nonlinear optimization toolbox of MATLAB ® to find the set of b, σ, and µ that maximizes the loglikelihood function defined by Ogata and Katsura (1993) (see their Equation 8or Si and Jiang, 2019 for details).M c (d max ) is then taken to be µ + 2.4σ, i.e. the magnitude where we expect 99% of earthquakes are detected.We provide the optimizer limits on b, selecting 0.85 ≤ b ≤ 1.05, as we found the output b to vary dramatically otherwise.Also, given that we threshold our catalog at M W ≥0, we set 0 as the lower integral bound in Eq. 3, rather than −∞.
Results for three example distance bins with d max = 23, 192, and 684 km are shown in Figure 3a-c, and the overall results shown in Figure 3d.Note that we attempted to estimate uncertainties of each M c by repeating the Ogata and Katsura (1993) method in 200 bootstrap iterations (and these are shown in Figure 3).However, we find that systematic trends in the M c (d) data are more relevant to model fits rather than "noise" that is characterized by the bootstrap confidence intervals, so we opt not to use these uncertainties in the modelfitting process.The thinning width σ(d) generally covaries with M c (d), although it has a prominent peak at d<100 km, where there is little range between M c and the cutoff magnitude of M W 0, so σ is less well constrained; this results in an overall σ-M c correlation coefficient of 0.60.
The resulting M c (d) is poorly fit by a power law due to a bend to unexpectedly low M c in the d range of ∼100-500 km, and we therefore omit that distance range to compute the power-law model shown in Predicted yearly earthquake density based on the CNED catalog and the magnitude-of-completeness analysis of this study.All maps were first computed on a coarse grid (∼15 km spacing), then converted to a finer grid and smoothed with a 2D Gaussian filter.
far greater freedom than the power-law, and as a result fits the data much more closely.For the remaining sections we use the smooth-increasing M c (d) model.However, we repeat the analyses using the power-law model and plot the results in Supplemental Figures S1-S4, which demonstrate that the choice has little impact on our overall conclusions.

Estimating b given spatiotemporal M c variations
The standard and most-accepted way to estimate b is the Aki-Utsu maximum-likelihood estimator (MLE) (Utsu, 1965;Aki, 1965): where M is the mean magnitude of the catalog, including only earthquakes with M ≥M c .There have been several techniques introduced to allow time-varying M c in the MLE (e.g.Weichert, 1980;Kijko and Smit, 2012;van der Elst, 2021).Taroni (2021) presented a convenient modification of the MLE that does not require the evaluation of b for subcatalogs.Ignoring two minor corrections, their MLE can be written: where M cE is M c (t) evaluated for each earthquake.In this study, we further generalize their method by allowing spatial variation of M c in addition to temporal variation, i.e. we consider M cE = M c (x, t), evaluated using d(x, t) for each earthquake.Derivation of Equation 5 is shown in Appendix 2, beginning from the magnitude probability density function of Aki (1965) and including an extension to incorporate a maximum magnitude (Page, 1968).The notion of considering the completeness level for each earthquake is not only useful in applying the MLE, but more generally we can examine the MFD using the relative magnitude M * = M − M cE .
Figure 4 displays regular M W and M * GR plots for both the CSZ and the full Atlantic Canada region, as defined by the map area in Figures 1 and 2. In both cases, M * -derived b are lower than the M W -derived estimate, but not significantly so according to the 95% confidence limits from bootstrapping.Confidence limits are highly dependent on the number of events included, and thus the M c chosen.We therefore plot b vs. M c (or b * vs. M * c ) for each GR plot (Figure 4c,d,g,h).These plots demonstrate that b * is more stable over M * c than their equivalent M W -derived estimates.Note that here we are ignoring potential variations of b to get average results over broad areas and times, despite observing varying b in the previous section.

Correcting density maps for undetected earthquakes
Here we apply our new MLE (Eq.5) in order to estimate r = r(M c ), defined as the ratio (total M W ≥1 earthquakes)/(recorded M W ≥1 earthquakes) expected for a given M c (x, t).If the thinning parameter σ was consistent in space and time, we could estimate the function r(M c ) using only the MFD in Figure 4d.However, because we noted in Section 3.1 that σ is not constant, we expect more accurate results if we estimate r(M c ) from multiple MFDs, formed using narrower ranges of M c (or, equivalently, narrower ranges of d).We therefore consider the following procedure for a series of M c spanning 1.0 to 5.2, using increments of 0.1: 1. Select all earthquakes with M cE ≤ max(M c , 2) to form the M * MFD.
2. Estimate b using the MLE described in Eq. 5 and Appendix 2.
3. Fit a thinning function q(M * |σ, µ) on the MFD, with b constrained to the MLE value.
4. Compute the ratio r(M c ) as: where the integral limits define the range between M W = 1 and the maximum magnitude in the MFD.
The resulting r(M c ) function is shown in Figure 5; note that it approaches a log-linear relationship with slope of ∼1, which is expected as b ≈ 1.0 in Figure 4f,h.As an alternative model that does not depend on fitting a thinning function, we can simply extrapolate the GR model fit in step 2 to predict the total number events This ratio is plotted as the green curve in Figure 5 and produces similar results.
We opt not to fit a best-fit curve and instead directly interpolate the r(M c ) results to compute r E = r(M cE ) for each earthquake.The ratios r E represent the predicted number of M W ≥1 earthquakes each (recorded) earthquake represents, and can be used to make the calibrated event-rate map shown in Figure 2d.This is in practice very similar to how the scalar moment of each earthquake is used to make the moment-rate map.Note that we cannot simply multiply the grid-cell values from Figure 2b by a ratio like r(M c ) because M c varies in time, as well as space.Finally, the resulting map shows that offshore seismic zones have comparable earthquake densities to the CSZ and LSZ, and it is much smoother than the moment-rate map because it is not dominated by the infrequent, largest earthquakes.

Discussion and Conclusions
Although we opted to use the smooth-increasing M c (d) model rather than the power-law fit (which was highly sensitive to the particular data range included), we are not suggesting that a power-law is inappropriate for the region.The Mignan et al. (2011) M c (d) power-law fitting results show significant scatter for individual M c estimates, but they converge to a best-fit model because they have sufficient data to average M c from many MFDs   per distance bin.It is therefore unsurprising our M c (d) data show substantial scatter, but this does not fully explain why our M c (d) seems to have systematic deviations for a power-law, and remains constant over a range of ∼50 < d < 200 km.It could be that the true M c (d) follows a power-law more closely and that our M c (d) estimates for this distance range are underestimated due to non-log-linear effects in their MFDs, although we have no hypothesis to suggest why this should be the case over an extensive range of d(x, t).We demonstrate with Supplemental Figures S1-S4 that this issue does not substantially affect our following results or conclusions, but nevertheless this topic may warrant further investigation in future studies.The assessment of M c (x, t) more generally is discussed further near the end of this section.
Upon visual inspection, the GR plot produced by Equation 5 for the CSZ is quite similar to the raw M W plot (Figure 4a,b), but the b vs. M c plots show a modest improvement in the stability of b when using M * (Figure 4c,d).M * provides a much more dramatic improvement for the full Atlantic Canada catalog, as it produces much smaller confidence intervals on b, more stable b vs. M c , and a lower thinning width σ (Figure 4e-h).The lower σ (0.39 for M * vs. 0.65 for M W ) suggests a more angular MFD, rather than the gradual curvature that is caused by spatiotemporal heterogeneity in M c (Mignan, 2012).The full Atlantic Canada catalog is an extreme case, where M c varies greatly in both space and time, and there may be interesting b variation within the sample.Nevertheless, these observations suggest that M *derived b are more reliable than traditional estimates and that future studies, even those in areas where earthquakes and seismometers are more abundant, should consider spatiotemporal M c variation when estimating b.
The predicted event-rate map (Figure 2d), as well as the associated methodology introduced here, are a robust way to compare earthquake rates between regions with different levels of seismometer coverage.It may be preferable to moment rate maps because it is sensitive to all earthquakes, instead of (effectively) only very large ones.With regard to the regulation of geological fluid-injection activities and induced-earthquake risk, the method is no replacement for local seismic monitoring efforts, but it may provide the best-available baseline earthquake rate estimates.Due to variable b and non-log-linear effects in MFDs, event rates will not always correlate perfectly with hazard.We also note that, for the purposes of hazard assessment, any map of recorded earthquakes is only useful to the extent that previous earthquakes locations can help predict the sites of future large earthquakes; our analysis does nothing to account for the possibility that areas of elevated intraplate seismicity today are extended aftershock sequences (Basham and Adams, 1983; Toda and   Ogata and Katsura (1993) model fits in panels a, b, e, and f, were computed with b and M c fixed to be consistent with the MLE fit and resulted in thinning widths σ of 0.39, 0.39, 0.65, and 0.39, respectively.Figure 5 The predicted ratio of total/detected M W ≥1 earthquakes r, as a function of M c , computed using Equation 6 for 1.0≤M c ≥5.2 using increments of 0.1.The ratio approaches the log-linear relationship log 10 r = 1.00M c − 1.92 (computed using least-squares over the range 3.9≥ M c ≤5.2).The green curve is an alternative r(M c ) model computed by extrapolating the GR model (without the thinning function), as described in the text.
Stein, 2018).Finally, we must acknowledge some of the limitations of our M c (x, t) analysis.In addition to uncertainty in station metadata and the "noise" caused by periods where a seismometer was nonoperational that we have not accounted for, treating all seismometers equally is also a severe limitation.Noise levels, and signal-tonoise ratio, vary between seismometers for many reasons; geology at the site, instrument type, and proximity to anthropogenic or ocean-wave noise sources all being important factors.Instrumentation quality and noise levels also generally improve throughout the study period.Although we do not expect it to be a major factor in this dataset, short-term aftershock incompleteness (Stallone and Falcone, 2021) can also cause M c variation in time.Schorlemmer and Woessner (2008) use a full phase pick catalog to empirically determine the likelihood of an earthquake being picked at a particular seismometer as a function of magnitude and distance; then, taking these functions at all seismometers, they equate M c (x, t) to a threshold probability of the earthquake being detected at four or more seismometers.Mahani et al. (2016) measure ambient noise levels at each seismometer and compare them with theoretical earthquake amplitudes in order to spatially map M c .We expect that incorporating either system in place of, or in combination with, our simplified d-M c relation would further improve the predicted event rate and b estimates.

Figure 1
Figure 1 (a) Map of M W ≥ 1 earthquakes, with dots coloured and sized by magnitude.A 600 m bathymetric contour, corresponding roughly to the continental shelf edge, is plotted in blue.(b) Scatter plot of earthquake magnitudes in time (blue dots).The grey markers and line mark the mean magnitude by year.The red markers and line indicate the total number of earthquakes per year, corresponding to the second y-axis on the right side.

Figure 2
Figure 2 (a) Time-averaged weighted station-distance metric d over the study period of 1985-2023.Contours indicate distance in kilometres.Red triangles indicate locations of seismometers active for ≥3 years of the study period.Canadian provinces/regions and the United States (U.S.) are labeled.(b) Uncorrected yearly M >1 earthquake density (N km −2 y −1 ) from the CNED catalog.(c) Moment density (J km −2 y −1 ) for the same catalog.(d) The main result of this study (Section 3):Predicted yearly earthquake density based on the CNED catalog and the magnitude-of-completeness analysis of this study.All maps were first computed on a coarse grid (∼15 km spacing), then converted to a finer grid and smoothed with a 2D Gaussian filter.

Figure 3
Figure3(a-c) The MFD and resultingOgata and Katsura (1993) model fit for three distance bins of 300 earthquakes.95% confidence intervals from bootstrapping are reported for each model parameter.(d) Overall M c (d) results (blue dots) with a power-law fit (red curve) and a best-fit smooth, increasing model (black dashed curve).Cyan and grey error bars indicate the 50% and 95% confidence intervals from bootstrapping, respectively.The green curve indicates the (4 th -nearest station) model ofMignan et al. (2011), M c = 5.96d 0.0803 4 − 5.80.Orange triangles indicate the corresponding σ for each distance bin, as indicated on the right-hand side y-axis.

Figure 4
Figure 4 (a) Standard GR plot for the CSZ.Red and blue markers show the non-cumulative and cumulative MFDs, respectively.MLE-derived b is printed with 95% confidence intervals from 500 bootstrap iterations, as well as the standard error σ b , which depends directly on the number of events with M ≥M c (N ): σ b = b/ √ N .The M c used is shown with the dashed black line.(b) CSZ GR plot using M * to account for variable M c (x, t).(c) b vs. M c for the CSZ, using raw M W , cyan and grey error bars show the 50% and 95% intervals from 500 bootstrap iterations.(d) b vs. M * c for the CSZ, using M * , (e-h) like a-d except for the full map region of Figure 1.TheOgata and Katsura (1993) model fits in panels a, b, e, and f, were computed with b and M c fixed to be consistent with the MLE fit and resulted in thinning widths σ of 0.39, 0.39, 0.65, and 0.39, respectively.