With the increased seismicity rate in the Western Canada Sedimentary Basin (WCSB), accurate determination of local magnitude (ML) 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 ML 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 ML 4.0 or larger.
Currently, Natural Resources Canada (NRCan) follows the Richter (1935) scale, which was originally designed for southern California, to determine the ML of seismic events in western Canada. One key difference between the Richter’s calculation of ML 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 ML 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 ML 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 ML for induced earthquakes within WCSB.
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.
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 ML.
Equations (1) and (2) can be combined and re-arranged to give
where ε is the residuals in log units. We obtain n, k, ML, 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, ML, 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 b1 and b2 are the coefficients to be determined.
In theory, the ML 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 ML 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.
Using Equations (3) and (4), model parameters b1, b2, k, ML, 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(A0), for WCSB is
Figure 4 (top panel) shows the comparison of the -log(A0) 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(A0) (and therefore smaller ML). This can also be seen in Figure 4 (bottom panel) which shows, for each event i, the difference between the ML value derived from the peak amplitude at each station j (ML,ij) and the event’s final ML value (ML,i, defined as the median of all ML,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 ML 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 ML 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, ML 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 ML calculation is based on measurements from the vertical component despite the small number of available samples. As expected, Y17 under-predicts the ML value for all magnitudes in this case.
|Table 1. Correction term (S in Equation 1) for each seismograph station.|
Using a comprehensive earthquake dataset for the period of 2014-2018 (inclusive), we provide an update on the distance correction term in the ML 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 ML values calculated at individual stations and the final event ML 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.
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.
About the Author(s)
Alireza Babaie Mahani (Ali Mahani) Ph.D., P.Geo. is the director of Mahan Geophysical Consulting Inc. (www.mahangeo.com). He has more than 15 years of experience in earthquake ground-motion analysis and seismic hazard assessment. Since 2015, he has been the project seismologist for the Induced Seismicity Monitoring Project in British Columbia working with Geoscience BC, BC Oil and Gas Commission, Geological Survey of Canada, and Canadian Association of Petroleum Producers. He has performed several research projects in northeast British Columbia including evaluation of the monitoring network for accurate earthquake detection, the basis for local magnitude determination, construction of velocity models for the Montney unconventional resource play, analysis of focal mechanism and orientation of stress in the region, and ground motion analysis of large induced earthquakes.
Honn Kao Ph.D., P.Geo. obtained his BSc in Geophysics from the National Central University, Taiwan, in 1985, and his MSc and PhD in Geophysics from the University of Illinois at Urbana-Champaign in 1991 and 1993, respectively. He worked at the Institute of Earth Sciences (IES), Academia Sinica, Taiwan, as an Assistant Research Fellow (1993-1996), eventually was promoted to Full Research Fellow in 2000. He joined the Geological Survey of Canada in 2002. He has been an Adjunct Professor at the School of Earth and Ocean Sciences, University of Victoria since 2006. He is also an IES Adjunct Research Fellow since 2019. Currently, he is the leader of the Induced Seismicity Research Project of the Department of Natural Resources of the Government of Canada.
Babaie Mahani, A. and G. M. Atkinson (2012). Evaluation of functional forms for the attenuation of small-to-moderate-earthquake response spectral amplitudes in north America, Bulletin of Seismological Society of America 102, 2714-2726.
Babaie Mahani, A. and H. Kao (2019). Accurate determination of local magnitude for earthquakes in the western Canadian sedimentary basin, Seismological Research Letters 90, 203-211.
Hutton, L. K., and D. M. Boore (1987). The ML scale in southern California, Bulletin of Seismological Society of America 77, 2074-2094.
Kao, H., D. W. Eaton, G. M. Atkinson, S. Maxwell, and A. Babaie Mahani (2016). Technical meeting on the traffic light protocols (TLP) for induced seismicity: summary and recommendations, Geological Survey of Canada open file 8075, 20 pp.
Kao., H., R. Visser, B. Smith, and S. Venables (2018). Performance assessment of the induced seismicity traffic light protocol (IS-TLP) for northeastern British Columbia and western Alberta, The Leading Edge, 37(2), 117-126, doi: 10.1190/tle37020117.1.
Menke, W. (1984). Geophysical Data Analysis: Discrete Inverse Theory, Academic Press, San Diego, California, 289 pp.
Ottemoller, L. and S. Sargeant (2013). A local magnitude scale for the United Kingdom, Bulletin of Seismological Society of America 103, 2884-2893.
Rezapour, M. and R. Rezaei (2011). Empirical distance attenuation and the local magnitude scale for northwest Iran, Bull. Seismol. Soc. Am. 101, 3020-3031.
Richter, C. F. (1935). An instrumental earthquake magnitude scale, Bulletin of Seismological Society of America 25, 1-31.
Savage, M. K., and J. G. Anderson (1995). A local-magnitude scale for the western Great Basin-eastern Sierra Nevada from synthetic Wood-Anderson seismograms, Bulletin of Seismological Society of America 85, 1236-1243.
Uhrhammer, R. A., M. Hellweg, K. Hutton, P. Lombard, A. W. Walters, E. Hauksson, and D. Oppenheimer (2011). California integrated seismic network (CISN) local magnitude determination in California and vicinity, Bulletin of Seismological Society of America 101, 2685-2693.
Visser, R., B. Smith, H. Kao, A. Babaie Mahani, J. Hutchinson, and J. E. McKay (2017). A comprehensive earthquake catalogue for northeastern British Columbia and western Alberta, 2014-2016, Geological Survey of Canada Open File 8335, 28 pp.
Visser, R., H. Kao, B. Smith, C. Goerzen, B. Kontou, R.M.H. Dokht, J. Hutchinson, F. Tan, and A. Babaie Mahani (2020). A comprehensive earthquake catalogue for the Fort St. John-Dawson Creek region, British Columbia, 2017-2018, Geological Survey of Canada Open File, in press.
Yenier, E. (2017). A local magnitude relation for earthquakes in the western Canadian sedimentary basin, Bulletin of Seismological Society of America 107, 1421-1431.