Determination of Local Magnitude for Induced Earthquakes in the Western Canada Sedimentary Basin: An Update

Introduction

With the increased seismicity rate in the Western Canada Sedimentary Basin (WCSB), accurate determination of local magnitude (M_L) using data from regional seismograph networks is a vital task for induced seismicity monitoring and mitigation of seismic hazard.  Specific mitigation strategies have been put in place by the regulators when the induced events exceed a predefined M_L threshold (Kao et al., 2016, 2018).  Regulations in both British Columbia and Alberta impose suspension orders to temporarily halt injection operations in case of an induced event with M_L 4.0 or larger. 

Currently, Natural Resources Canada (NRCan) follows the Richter (1935) scale, which was originally designed for southern California, to determine the M_L of seismic events in western Canada.  One key difference between the Richter’s calculation of M_L and the NRCan’s is that the former measures the maximum amplitude from horizontal components, whereas the later uses the vertical component.  Moreover, the original Richter method uses the largest recorded amplitude (on a Wood-Anderson, WA, seismometer) regardless of seismic phases, whereas NRCan uses the window encompassing the S phase.

Yenier (2017; referred to hereafter as Y17) observed an increase in the M_L values reported by NRCan with distance (Figure 1) and concluded that the distance correction terms in the original Richter scale may not be suitable for WCSB.  Using data from horizontal components, Y17 also obtained distance correction terms to account for the attenuation effects across WCSB.  In a similar study, Babaie Mahani and Kao (2019; referred to hereafter as BMK19) estimated the distance correction terms (i.e. proper adjustment to the original Richter M_L table) when the vertical component is used in the calculation.

In this report, we provide an update on the distance correction terms of BMK19 using newly available data.  We also provide comparisons between the result obtained by the updated model and those of Y17 and BMK19.  Moreover, using the dataset from a local seismograph array established by one of the operators in the region, we analyze the performance of each model and verify their consistency in the determination of M_L for induced earthquakes within WCSB.

Dataset

In BMK19, the analysis was done with the vertical WA amplitudes taken from the comprehensive catalogue of local and regional earthquakes in WCSB compiled by Visser et al. (2017) for the period of 2014-2016 (inclusive).  The dataset includes 18918 amplitude values recorded at epicentral distances between 0.4 and 600 km from 4916 earthquakes.  Taking advantage of the recent increase in the number of local and regional seismograph stations, Visser et al. (2020) provided another comprehensive dataset for the period of 2017-2018 (inclusive).  The new dataset includes 63883 vertical WA amplitudes from 10693 earthquakes within the southern Montney play in northeast British Columbia recorded at epicentral distances between 0.1 and 600 km.  The new dataset provides an opportunity to review the correction factors obtained by BMK19 and to improve the accuracy of these factors in light of a bigger dataset.  The combined dataset (i.e. Visser et al., 2017 and 2020) includes a total of 82801 vertical WA amplitudes at epicentral distances between 0.1 and 600 km from 15609 earthquakes.  The other dataset used in this study is provided by a local operator, consisting of peak amplitude measurements from a local seismograph array in west central Alberta.  This dataset includes 19770 WA amplitude values (vertical and horizontal components) from 948 earthquakes.  Note that the private dataset is only used for the purpose of comparison of the results.  Figure 2 shows the number of data samples as a function of epicentral distance.

Methodology

Richter (1935) devised the first systematic measure of the size of an earthquake using a small dataset from southern California.  The most important input parameter in the magnitude calculation is the maximum zero-to-peak amplitude measured from the horizontal components of ground displacement, regardless of the wave type, recorded on WA torsion seismometers.  Since then, this method has been widely adopted around the globe with various modifications for different tectonic regions (e.g. Hutton and Boore, 1987; Savage and Anderson, 1995; Uhrhammer et al., 2011; Ottemoller and Sargeant, 2013).  To preserve the original sensor type, WA amplitudes are synthetically obtained through deconvolution of the recording sensor’s instrument response (non-WA sensors) from the recorded waveform and convolution with the WA type instrument response.

The Richter magnitude scale is defined as

where A is the recorded amplitude, in mm, on a WA seismometer, A0 is the distance correction term for the recorded amplitude, and S is a correction term for each station that accounts for any systematic amplitude difference such as local site effect.  Richter (1935) determined the S terms for each station by taking the average of the magnitude differences between the values obtained for individual stations and the mean value for the event.  A0 at each distance was obtained by defining the zero magnitude (ignoring the S term) at a distance of 100 km (A = 0.001 mm).  Here, we obtain A0 for WCSB through regression of the WA amplitudes, as compiled by Visser et al. (2017, 2020), against a given model that accounts for the attenuation of ground motion amplitudes in the WCSB.

The model for the distance correction can be defined as (Rezapour and Rezaei, 2011)

where Rhypo is hypocentral distance, and n and k are the corrections for geometrical spreading and anelastic attenuation, respectively.  Equation (2) is written in such a way that for the reference distance of 100 km, the corresponding amplitude is 0.001 mm (-log(A0) = 3) to maintain the consistency with the original definition of M_L.

Equations (1) and (2) can be combined and re-arranged to give

where ε is the residuals in log units.  We obtain n, k, M_L, and S by constructing the linear inverse equation of d = Gm (Menke, 1984) where d is observed data (log(A) + 3.0) and m is the model parameter vector (n, k, M_L, and S) to be determined.  G is the kernel matrix, which relates d to m.  In solving for the model parameters, the correction term for each station (S) is obtained simultaneously.  Here, the correction terms for stations are computed relative to the average site condition by constraining the S terms to attain zero when averaged over all stations.

In modeling the attenuation, an important factor is the shape or functional form of geometrical spreading, which is characterized by n.  Linear, bi-linear, and tri-linear functional forms have been used to model the attenuation of ground motion from geometrical spreading (e.g. Babaie Mahani and Atkinson, 2012).  Following BMK19, we use a bi-linear function with the transition distance at 85 km as

to model the effect of geometrical spreading (n), where b₁ and b₂ are the coefficients to be determined.

In theory, the M_L value obtained from the horizontal component of ground motions should be consistent with that from the vertical component, as long as the correction terms are properly calibrated.  Historically NRCan chose to use the vertical-component peak amplitude in the M_L calculation for two reasons.  First, single-component (vertical), short-period sensors were installed in many seismograph stations of the Canadian National Seismograph Network (CNSN).  Therefore, using vertical-component data would maximize the number of available stations in the calculation, which is especially important for regions with very sparse station coverage.  We note, however, that all CNSN stations are now equipped with 3-component broadband sensors.  The second reason is that the level of background noise on the vertical component is generally lower than that on the horizontal components.  To maximize our data resolution, we follow NRCan’s convention of using vertical-component measurements in BMK19 and this study.

Results

Using Equations (3) and (4), model parameters b₁, b₂, k, M_L, and S are obtained through the least-square inversion in MATLAB.  Here, we only use events with at least 8 observations and stations that recorded at least 8 events (3575 events with 33994 vertical WA amplitudes at 56 stations).  Figure 3 shows the results of regression by plotting the residuals versus hypocentral distance and magnitude.  There is a clear over-prediction at hypocentral distances < 2 km.  Although the epicentral distance of the recorded events can be as low as 0.1 km, the focal depths of most events are more than 2 km.  Therefore, samples with hypocentral distances < 2 km may correspond to very shallow events in the earthquake catalogues or mislocated events due to insufficient data constraint.  For this reason, the validity of our model would be for hypocentral distances of > 2 km.  Beyond 2 km there are no trends in the residuals with distance.  Similar results can be seen for the relationship between residual and magnitude, as no trend can be recognized across the entire magnitude range.  The lack of an apparent trend means that the inversion was successful assuming the functional forms for the attenuation and geometric spreading. 

Our updated distance correction term, -log(A₀), for WCSB is

Figure 4 (top panel) shows the comparison of the -log(A₀) term from Y17, BMK19, and the updated model with respect to hypocentral distance. All data samples are normalized by the average amplitude in the distance bin of 90-110 km such that the distance correction term is 3 at 100 km. The three models are similar for distances around 100 km and above. At closer distances, the updated model shows small improvement compared to BMK19, whereas Y17 gives significantly smaller values for -log(A₀) (and therefore smaller M_L). This can also be seen in Figure 4 (bottom panel) which shows, for each event i, the difference between the M_L value derived from the peak amplitude at each station j (M_L,ij) and the event’s final M_L value (M_L,i, defined as the median of all M_L,ij). The observed under-prediction by Y17 at hypocentral distances < 10 km is mainly due to the difference between vertical and horizontal amplitudes. As mentioned before, whilst Y17 use the horizontal amplitudes in the M_L calculation, BMK19 and the updated model in this study both use the vertical component of ground motion. The lower amplitudes of shear waves on the horizontal components compared to the vertical component at short distances can explain the discrepancy between the models. In Table 1, the correction term for each station (S) is compiled.

The Private Dataset

Figure 5 shows the comparison between M_L values as reported in the private dataset and those obtained using Y17, BMK19, and the updated model. In the top row, we use all of the available WA amplitudes, including both vertical and horizontal amplitudes, whereas in the middle and bottom rows, only the horizontal and vertical components are used, respectively. It is clear in Figure 5 that, for the top and middle row, M_L values obtained by Y17 match the reported values in the private dataset very well, while BMK19 and the updated model over-estimate for all magnitudes. On the other hand, BMK19 and the updated model match the values reported in the private dataset very well when the M_L calculation is based on measurements from the vertical component despite the small number of available samples. As expected, Y17 under-predicts the M_L value for all magnitudes in this case.

Conclusion

Using a comprehensive earthquake dataset for the period of 2014-2018 (inclusive), we provide an update on the distance correction term in the M_L calculation of local earthquakes in WCSB. A total of 33994 vertical Wood-Anderson amplitudes from 3575 earthquakes recorded at 56 local and regional seismograph stations are used in the regression. Our revised formula can better characterize the attenuation effect associated with the unique geological and tectonic configuration of WCSB. As a result, it provides higher consistency between M_L values calculated at individual stations and the final event M_L over a wide hypocentral distance range from 2 to 600 km, and therefore is a reasonable alternative to the original Richter formula for regulatory purposes in northeastern BC and western Alberta.

Acknowledgements

Waveform data used in this study are obtained from the Canadian Hazard Information Service and the Data Management Center of the Incorporated Research Institutions for Seismology (IRIS). The private dataset is provided to NRCan under a data sharing agreement. This research is partially supported by NRCan’s Energy Innovation Program, Geological Survey of Canada Environmental Geoscience Program, and a research grant from Geoscience BC. NRCan contribution 20200022.