A UK M_w catalogue derived from coda envelopes

SUMMARY

The United Kingdom (UK) experiences low-to-moderate levels of seismicity; only 12 onshore earthquakes with local magnitude (M_L) ≥ 4.0 have been recorded in the past 20 yr. It is therefore difficult to estimate moment magnitude (M_w) using conventional moment tensor inversion for the majority of UK seismicity, resulting in limited reliable estimates of M_w. To address this, we calibrated coda envelopes at 16 broad-band seismic stations distributed across the UK to produce an M_w catalogue for 100 events with M_w ≥ 2.13 that occurred since 2006. This was achieved using the open-source Coda Calibration Tool, which requires independent source parameter estimates for calibration. For 13 UK events between 2006 and 2022, we used spectral modelling to estimate apparent stress (0.32–1.74 MPa) and moment tensor inversion to estimate M_w (3.35–4.52). These independent source parameters formed a subset of the inputs into the final calibration, which used seismic data from 33 events with coda-derived values of 2.57|$\le $| M_w|$\le $|4.49. The resultant coda calibration parameters were applied to 67 further events (M_w ≥ 2.13). The coda envelopes exhibit slow seismic coda decay across the UK, with significant energy up to 20 Hz, consistent with other regions of low tectonic activity. This M_w catalogue, and the application of the calibration to future UK seismic events, will be useful for both assessing seismic hazard and event characterization.

Coda waves, Earthquake source observations, Waveform inversion

1 INTRODUCTION

There are numerous magnitude scales for measuring seismic event size. Richter (1935) developed the first magnitude scale that measures the maximum-recorded signal amplitude, without distinguishing between seismic phases, for observations at seismic stations within 600 km of the earthquake. This magnitude scale, commonly referred to as the local magnitude (M_L), was developed for observations of earthquakes in Southern California but continues to be used in local and regional seismic catalogues globally. The M_L formulation underestimates earthquake magnitudes for distant (>600 km) earthquakes (where S-wave attenuation is not fully captured), deep earthquakes (where surface waves are smaller) and for large earthquakes (it does not take into account the duration of the shaking). Additionally, M_L saturates at large magnitudes due to the amplitude measurements being made across a limited frequency. Moment magnitude (M_w), derived from seismic moment (M₀), is regarded as the most reliable magnitude scale, since it does not saturate at large magnitudes (Hanks & Kanamori 1979), and it can be directly related to properties of the earthquake fault (e.g. Aki 1966).

Robust estimation of M_w is vital in seismic monitoring for multiple reasons. First, it is needed for creating accurate seismic catalogues, which are used in seismic hazard studies to highlight areas of increased seismicity. Additionally, M_w estimates are used for characterizing seismic sources. For example, event characterization is required to discriminate between earthquakes and underground explosions, and M_w values can be used to estimate explosion yields (e.g. Pasyanos & Chiang 2022). Accurate estimates of M_w are also important in seismic monitoring of subsurface resource operations, including mining, hydraulic fracturing, and geothermal energy. Finally, attempts at linking narrow-band magnitudes such as M_L to M_w are often associated with significant bias (Shelly et al. 2022) so an M_w methodology that can be applied for both small and large events is preferable.

Conventionally, estimates of seismic source parameters are calculated from long-period moment tensor inversions derived from waveform fitting using 1-D velocity models (e.g. Petersen et al. 2021). Accurately estimating M_w from moment tensor inversions is dependent upon knowledge of the subsurface seismic velocity model. Inaccuracies in the seismic velocity model result in inaccurate synthetic Green's functions and therefore an unreliable M_w estimate. This is a significant disadvantage of the waveform modelling approach since the subsurface velocity structure is frequently poorly known and likely to vary, meaning 3-D seismic velocity models are required for different areas and source–station pairs. For smaller magnitude events (M_w < 3.5), this approach is unstable due to low signal-to-noise ratios (SNRs) in the 10–20 second period band. In the UK, there are limited reliable M_w estimates for seismic events below M_w 3.5. A lack of an accurate knowledge of the UK 2-D and 3-D crustal structure, required for intermediate (100–15 s) and short-period (2.0–0.2 s) waveform modelling, and the large computational expense in deriving complex Green's functions at shorter periods makes this technique unfeasible. This is important since low magnitude M_w estimates are known to vary between studies, meaning it is difficult to compare event M_w estimates if they are obtained via different techniques and with differing velocity models. For low magnitude events, which often have low SNRs at intermediate periods, the effects of 2-D and 3-D structure are more pronounced in the signal passband. Therefore, identifying methods to estimate magnitudes for such events remains an active area of research (e.g. Laurendeau et al. 2022).

Modelling high-frequency seismic amplitude spectra is another approach for calculating source parameters (e.g. Holt et al. 2021), such as apparent stress and M_w. These estimates are often influenced by heterogeneous path and site variations, and require significant corrections that often require large amounts of data for path and site calibration (e.g. Oth et al. 2011). As a result, creating a homogenous set of M_w estimates using either long-period moment tensors or modelling seismic amplitude spectra in the UK is unsuitable for large catalogues that include events as small as M_w 2.5.

However, seismic coda display minimal sensitivity to 3-D path heterogeneities at the low frequencies required to constrain M_w (Mayeda & Malagnini 2010). Determining source model parameters such as seismic moment, corner frequency and apparent stress, using seismic S-wave coda can be robust and stable, even when considering low magnitude events (e.g. Mayeda et al. 2003). Seismic coda are generated by the scattered wavefield, and therefore have low sensitivity to the source radiation pattern and directivity, with the result that source spectra, and hence seismic moment estimates from coda modelling, are often much more stable and robust than those derived from traditional direct wave approaches (Mayeda et al. 2003). We use the Java-based coda calibration tool (CCT), which follows the empirical method outlined by Mayeda et al. (2003) and offers an efficient method for calibrating and processing envelopes of seismic coda (Barno 2017). The coda calibration process results in estimates of earthquake source parameters, including M_w and apparent stress, which are up to four times less variable than those derived from conventional direct wave considerations (Mayeda et al. 2020; Morasca et al. 2022 ). This stability arises because coda envelopes are less affected by the radiation pattern and source directivity than direct waves (Aki & Chouet 1975; Mayeda & Malagnini 2010) since scattering acts to average out a significant amount of the spatial heterogeneity. Coda methods are also applied across a wide passband, such that low magnitude events that are rich in high frequencies can also be analysed using this method; hence the applicability of coda methods for the analysis of low magnitude events.

The UK lies on the northwest part of the Eurasian tectonic plate, with the nearest plate boundary (the Mid-Atlantic ridge) over 1500 km away. The continental crust was formed in either the Precambrian or lower Palaeozoic age (Bluck et al. 1992) by the amalgamation of several fault bounded terranes. Tectonic processes in the Caledonian (460 Ma) and Variscan (290 Ma) orogenies produced the current terrane configuration (Woodcock & Strachan 2012). As a result of its location, the UK is usually characterized as having low-moderate levels of seismicity (Musson 2007), with ML^BGS [local magnitude calculated by the British Geological Survey (BGS)] values seldom exceeding 5.0 (Fig. 1). Observations dating back several hundred years suggest that people regularly feel earthquakes in the UK, but they rarely cause significant damage. Only 12 onshore earthquakes with ML^BGS > 4.0 have been recorded in the past 20 yr and Musson (2007) estimates a 100-yr return period for an ML^BGS 5.5 event.

Figure 1.

(a) Events processed in this study (coloured circles), overlying regional seismicity 2006–2022 (grey circles). This event data set includes validation (blue) and reference (purple) events with independent source parameter estimates, additional calibration events (orange) and CCT-processed events (yellow). Event symbol sizes are scaled by the event ML^BGS estimate. Event locations are from the BGS earthquake bulletin (e.g. Galloway 2020). (b) Distribution of stations (triangles) used throughout our study. Stations used for coda calibration are orange triangles; whereas stations used only for independent source parameter estimates (via moment tensor inversion and spectral modelling) are blue triangles.

The spatial distribution of seismicity varies significantly across the UK (Fig. 1). The majority is found on the western flank of the UK in NW Scotland, across northern England, and in Wales, while eastern Scotland and southeast England are relatively aseismic. These earthquakes are generally considered to be due to the reactivation of existing faults by present-day deformation. However, surface ruptures are rarely reported for any UK earthquakes (modern or historical) making it difficult to associate earthquakes with specific faults.

The crustal strain field and its relation to the observed distribution of UK seismicity is unclear, with suggestions that it is related to glacio-isostatic rebound (e.g. Main et al. 1999), as well as upwelling of hot material from the original Iceland plume head, causing uplift of the crust and driving tensional loading stresses (e.g. Arrowsmith et al. 2005). The difficulty in obtaining stable M_w estimates via waveform modelling in an area of low magnitude and low-to-moderate seismicity levels, such as the UK, emphasizes the need for the coda calibration approach.

The BGS operates a nationwide network of seismograph stations in the UK and produces an annual bulletin of earthquakes located in the UK mainland and coastal waters. The BGS bulletin contains locations, local magnitudes (ML^BGS) and phase data but does not include estimates of M_w. Instead, M_w estimates for UK earthquakes are currently only available on an ad hoc basis as part of individual studies (e.g. Edwards et al. 2008; Sargeant & Ottemöller 2009; Ottemöller & Sargeant 2010). Here, we use CCT to create a new homogeneous M_w catalogue for the UK. To calibrate CCT-derived M_w estimates, we first calculate independent M_w and apparent stress estimates using moment tensor inversions and spectral modelling. The CCT method reduces the influence of site and path effects and allows us to create an M_w catalogue that includes low magnitude events. Once the region is calibrated, this approach allows M_w values for new events to be rapidly estimated using the derived calibration parameters, therefore increasing the number of events that are included in our M_w catalogue.

Seismic coda analysis for M_w has been undertaken in various tectonic settings, including high heat flow regions such as the western US (e.g. Mayeda & Walter 1996), cratonic regions like eastern Canada (e.g. Bent et al. 2021), the complex faulting region of central and northern Italy (Morasca et al. 2008; Morasca et al. 2022), the Arabian Peninsula (e.g. Chiang et al. 2021), and western India (e.g. Malagnini et al. 2006). Our study thus enables comparisons to be made between UK seismic coda envelopes and those in different tectonic regimes with distinctly different seismic attenuation (Q), heat flow, and tectonic history.

2 DATA

We use three-component broad-band (100 Hz sampling rate) seismic data from the BGS local distance seismic networks (Fig. 1b). The majority of the data we use are from the BGS seismic network and from UKNET sites run by the Atomic Weapons Establishment (AWE) Blacknest, which are available to download from Observatories and Research Facilities for European Seismology (ORFEUS). Data from stations in Ireland (Irish National Seismic Network), France (CEA/DASE Seismic Network), Germany (German Regional Seismic Network), Luxembourg (Luxembourg Seismic Network) and the Netherlands (Netherlands Seismic and Acoustic Network) were also used and data were downloaded from the Incorporated Research Institutions for Seismology (IRIS).

Signals from low-moderate magnitude events (M_w < 4.0), such as those in the UK, are rarely observed at distances beyond 800 km. The moment tensor inversions and spectral modelling used data from stations within 800 km of the source. Using seismic stations at this distance range also ensured we had sufficient azimuthal coverage for reliable moment tensor solutions. For coda calibration, data from only UK stations were used (13 BGS and three UK-Net stations were used (Fig. 1b), to avoid source–station paths with significantly different crustal properties. These UK stations were selected based on their longevity, reliability and spatial distribution within the UK.

3 INDEPENDENT EARTHQUAKE SOURCE CONSTRAINTS: ESTIMATING M_w AND APPARENT STRESS

The coda calibration process (Fig. 2) requires independent reference earthquake source estimates (M_w and apparent stress) as constraints to determine each station's frequency-dependent site term correction. We refer to these events as reference events. Validation events, which also have independent M_w estimates, are used to evaluate the performance of the coda calibration. Only limited M_w estimates for UK seismicity are reported and are generally only available on an ad hoc basis as part of individual studies (e.g. Sargeant & Ottemöller 2009; Ottemöller & Sargeant 2010). As such, the methods used to estimate M_w differs and the estimated M_w can vary for the same event. Therefore, the independent estimates of M_w used in this study were obtained by moment tensor inversion (Section 3.1). For the high frequency site terms, independent estimates of apparent stress are used in the coda calibration and are preferable to assuming, a priori, a source scaling for the region (e.g. Morasca et al. 2022). Apparent stress estimates for UK seismicity are limited. Therefore, spectral modelling (Section 3.2) was used to estimate independent apparent stress and M_w for the coda calibration process.

Figure 2.

Workflow for this study.

3.1. Moment tensor inversion to determine M_w

M_w estimates were calculated using MTtime, a Python package for time-domain moment tensor inversion (Chiang & USDOE National Nuclear Security Administration 2020) that estimates the deviatoric focal mechanism using the approach outlined by Minson & Dreger (2008). We follow the same procedure as described in previous studies (e.g. Chiang et al. 2016, 2021). Observed waveform data are instrument corrected to displacement and then bandpass filtered (with the highpass corner between 0.02 and 0.05 Hz and the lowpass corner between 0.06 and 0.13 Hz, depending on the frequency content observed for each event).

Green's functions (GFs) are calculated for each source–receiver pair in the same frequency range as the observed data, using the frequency–wavenumber integration method (Wang & Hermann 1980) in Computer Programs in Seismology (Herrmann 2013). The best-fitting moment tensor solution is computed by inverting the time-domain waveform data using a least-squares minimization approach (Minson & Dreger 2008). This technique measures the fit between the observed and synthetic data over the length of the trace by the variance reduction (VR):

$$\begin{eqnarray} \mathrm{ VR} = \left( {1 - \frac{{\mathop \sum \nolimits_n {{w}_n}{{{\left( {{{{\boldsymbol{d}}}_{\boldsymbol{n}}} - {{{\boldsymbol{s}}}_{\boldsymbol{n}}}} \right)}}^2}}}{{\mathop \sum \nolimits_n {{w}_n}{{{\boldsymbol{d}}}_{\boldsymbol{n}}}^2}}} \right) \times 100, \end{eqnarray}$$

(1)

where |${\boldsymbol{d}}$| is the data, |${\boldsymbol{s}}$| is the synthetic waveforms, and w is the inverse distance weighting factor at each station n (Chiang et al. 2016).

This method estimates the source depth by performing the inversion for a range of depths, with the preferred solution being that with the largest VR. The generation of GFs requires a velocity model. For this study, we use the ak135 velocity model (Kennett et al. 1995), however we tested a range of alternative regional and global velocity models, including various BGS velocity models (e.g. Lownet and general UK, Mid Wales, Borders; Galloway 2020), Preliminary Reference Earth Model (PREM; Dziewonski & Anderson 1981), and a Western US model (Herrmann et al. 2011). The ak135 velocity model produced the largest VR for the 13 events for which we calculated moment tensors. To account for uncertainties in the hypocenter parameters and the velocity model, small time-shifts between the observed and synthetic seismograms are allowed (Pasyanos & Chiang 2022). These time-shifts are limited to half the minimum filtered period to avoid cycle skipping.

3.2. Moment tensor solutions

We were able to calculate stable deviatoric moment tensor solutions for 13 UK events that are distributed across the region of interest (Fig. 3, red circles), exhibit a range of magnitudes (ML^BGS between 3.5 and 5.2), and have high-quality data leading to high values of VR. The resulting moment tensor solutions have M_w values between 3.35 and 4.52 (see Table 1 for full moment tensor inversion solutions).

Figure 3.

Moment tensors estimated for UK earthquakes from this study (red) and Baptie (2010; blue). Numbers correspond to events in table 1 of Baptie (2010), and letters correspond to events in Table 1 (this study), with events in common highlighted with a dashed box. Circle size for each event is scaled by ML^BGS value.

Table 1.

Moment tensor inversion results for earthquakes in our study. Event origin time, latitude, and longitude are from the BGS earthquake bulletin (e.g. Galloway 2020).

Event	Origin time (YR-MM-DY Time)	MT depth (km)	DC/CLVD (%)	MLBGS	Mw	Mrr	Mtt	Mff	Mrt	Mrf	Mtf	Exponent (dyne-cm)	Locality
A	2006-12-26 10:40:04	6	91/9	3.5	3.41	−2.12	9.64	−7.52	−1.87	−3.71	−9.94	20	Dumfries
B	2007-04-28 07:18:11	2	85/15	4.3	4.03	−5.52	8.60	−3.08	3.56	−6.44	−4.16	21	Folkestone
C	2008-02-27 00:56:47	20	93/7	5.2	4.52	3.18	−3.36	0.18	2.86	−1.72	−4.22	22	Market Rasen
D	2009-04-28 10:22:09	6	82/18	3.7	3.50	−1.67	−15.06	16.73	0.14	−1.02	−8.55	20	Morecambe
E	2011-07-14 06:59:10	6	96/4	3.9	3.65	5.13	−0.92	−4.21	6.32	−12.62	−28.38	20	English Channel
F	2013-05-29 03:16:28	11	99/1	3.8	3.59	−5.95	18.93	−12.98	−4.15	−10.38	−16.88	20	Lleyn Peninsula
G	2014-02-20 13:21:30	8	81/19	4.1	3.63	3.37	12.03	−15.40	2.43	−7.61	−23.15	20	Bristol Channel
H	2015-05-22 01:52:17	13	86/14	4.2	3.67	1.95	−1.02	−0.93	−1.44	−2.10	−1.28	21	Ramsgate
I	2017-08-04 14:43:38	9	99/1	4.0	3.67	0.03	−28.88	28.84	−3.09	−0.84	17.55	20	Moidart
J	2018-02-17 14:31:07	15	86/14	4.6	4.03	−8.49	−81.18	89.67	0.50	−3.67	−75.75	20	Cwmllynfell
K	2018-06-09 22:14:28	10	97/3	3.8	3.35	3.47	−11.09	7.62	0.86	−5.32	−2.07	20	Grimsby
L	2020-09-08 08:45:28	3	86/14	3.5	3.47	−1.12	−5.17	6.29	−1.35	0.22	−14.56	20	Leighton Buzzard
M	2022-05-30 14:36:58	9	88/12	3.8	3.45	−0.88	−6.31	7.19	0.25	2.25	−13.58	20	Wem

An example deviatoric moment tensor solution from an earthquake in Cwmllynfell, South Wales (event J in Table 1) is shown in Fig. 4. This solution demonstrates a good fit between the observed (black) and synthetic (blue) waveforms, and the VR for this preferred solution is 89 per cent. In addition, this example shows the lack of sensitivity to source depth, where VR appears almost constant between 10–20 km. This feature of the inversion arises since we are predominantly fitting long-period surface waves. Nevertheless, our estimated depth is consistent with that estimated by the BGS (to within ∼5 km). The Cwmllynfell earthquake is an example of a strike-slip event, which is dominated by Love waves. It is well known that Love waves add information to the moment tensor components but are not sensitive to depth (Ben-Menahem & Harkrider 1964), whereas Rayleigh waves include information about the earthquake mechanism and are sensitive to the source depth. The approach of computing solutions for a range of depths provides insight as to how the inferred double couple (DC) component of the focal mechanism changes with source depth (Fig. 4b). For event J, the DC component is small at the shallowest and deepest depths, and large (86 per cent) at the depth of the preferred solution. For shear faulting, the DC component should comprise most of the focal mechanism (e.g. Vavryčuk 2015). Therefore, the large DC component in our preferred solution is a good indication of a reliable moment tensor solution.

Figure 4.

Moment tensor inversion results for event J calculated using MTtime. (a) Summary of the best-fitting inversion results and waveform fits at each station. Blue triangles on the lower hemisphere stereographic projection show station azimuth. From left to right, the values provided for each station are |$\Delta $| (distance in km), |$\theta $| (azimuth), time shift in samples (TS) and variance reduction (VR). (b) Moment tensor solution and M_w as a function of source depth. The colour of the moment tensor reflects the percentage DC (double couple) component of each solution. The vertical black line indicates the best-fitting solution.

The azimuthal station coverage for each event influences the stability of the moment tensor solution, and the non-DC component of the moment tensor solution. Petersen et al. (2021) calculate moment tensors for three earthquakes (M_w 3.5, 3.9 and 4.1) recorded by the dense AlpArray Seismic Network. This network provides excellent azimuthal coverage for moderate to large earthquakes in this region. Petersen et al. (2021) take advantage of the large number of stations in this network to investigate the stability of the moment tensor inversion to gaps in the azimuthal station distribution. They observe a clear increasing non-DC component trend with decreasing azimuthal coverage for the M_w 3.9 event; with a 180° station azimuthal coverage the non-DC component is only 10 per cent, but this increases to 40 per cent when the station azimuthal coverage is only 45°. The moment tensor calculated for the Cwmllynfell, South Wales earthquake (Fig. 4) had an azimuthal coverage of 270|$^\circ $|⁠. This azimuthal coverage is typical for the other events in our data set, excluding a few events in northern England and Scotland, which only have azimuthal coverage's of 180|$^\circ $|⁠. By only selecting events that have good station azimuthal coverage (>180|$^\circ $|⁠) for each of our moment tensor calculations, ensures that the non-DC component should remain low (Petersen et al. 2021).

We compare focal mechanisms estimated in this study to those found by Baptie (2010; Fig. 3). Baptie (2010) focal mechanisms are based on P-wave first motions and are therefore constrained to be pure DC solutions. Collectively, these mechanisms demonstrate a range of faulting regimes, even for closely separated earthquakes (e.g. events 5 and 17, separated by 25 km). For the three events in common between this study and Baptie (2010), the focal mechanisms are comparable.

3.3. Spectral modelling to determine apparent stress and M_w

We use SpecMod, a Python-based toolbox for processing and modelling seismic spectra (Edwards et al. 2010), to provide independent estimates of both apparent stress and M_w (allowing comparisons with estimated M_w values from moment tensor inversion). Spectral modelling for larger magnitude events (M_w > 3.5) is carried out to obtain stress drop estimates, such that the coda calibration process has input information regarding the high frequency site terms. The CCT-derived M_w estimates are not sensitive to the high frequency fall-off of the spectra (which is described by the stress drop), meaning that the stress drop estimates have minimal effect on the CCT M_w estimates. However, the high frequency terms of the coda calibration allow comparisons to be made between coda decay in different tectonic regimes, and therefore it is important that these terms are calibrated.

We follow a similar procedure to other studies (e.g. Holt et al. 2021). The horizontal components of the observed waveforms are instrument corrected to velocity and rotated to radial and transverse components. Signal and noise windows are extracted from the transverse component, which is the preferred component for isolating the S-wave signal. The signal window starts at 95 per cent of the elapsed time between the P and S arrivals and ends 20–30 s later. The noise window ends 2 s before the P arrival and is of equal length to the signal window. The signal and noise spectra are computed using the multitaper approach (Prieto et al. 2009). Finally, for each velocity spectrum we only fit the model between frequencies where the spectral SNR is greater than 3. If the spectral SNR is greater than 3 across the entire observed signal spectra, the lower frequency limit is set at 0.2 Hz (where the low-frequency surface wave appears to dominate) and the upper frequency limit is set at 40 Hz. The velocity spectra, |$A( f )$| are modelled using the Brune (1970) model using the following equation derived by Holt et al. (2021) (modified from Edwards et al. 2010):

$$\begin{eqnarray} {{\log }_{10}}A\left( f \right) = {{\log }_{10}}2\pi f + {{\log }_{10}}{{\Omega }_0} - {{\log }_{10}}\left[ {1 + {{{\left( {\frac{f}{{{{f}_\mathrm{ c}}}}} \right)}}^2}} \right] - \frac{{\pi f{{t}^*}}}{{\ln 10}}, \\ \end{eqnarray}$$

(2)

where f is frequency (Hz), |${{\Omega }_0}$| is the long-period spectral level at the source (m/s), |${{f}_{\mathrm{c}}}$| is corner frequency (Hz), and |${{t}^*}$| is an attenuation parameter. We constrain |${{\Omega }_0}$|⁠, |${{f}_{\mathrm{c}}}$| and |${{t}^*}$| using Powell's minimization technique (Powell 1964). Following Holt et al. (2021), we perform this inversion in two stages. First, we let all parameters be free and compute the best-fitting |${{\Omega }_0}$|⁠, |${{f}_{\mathrm{c}}}$| and |${{t}^*}$| at each station. We then calculate the average weighted |${{f}_{\mathrm{c}}}$| for the event, where each station's weight equates to the inverse of the hypocentral distance to the event. Second, we repeat the inversion with |${{f}_{\mathrm{c}}}$| fixed to this average value. Performing the inversion in this way helps resolve the trade-off between |${{f}_{\mathrm{c}}}$| and |${{t}^*}$| (Holt et al. 2021).

Before calculating |${{M}_0}$| from |${{\Omega }_0}$|⁠, we use |${{f}_{\mathrm{c}}}$| to calculate the rupture radius, r (assuming a circular fault; Madariaga 1976):

$$\begin{eqnarray} r = \frac{{k{{V}_{\mathrm{s}}}}}{{{{f}_{\mathrm{c}}}}}, \end{eqnarray}$$

(3)

where k is a constant related to the radiation pattern and |${{V}_{\mathrm{s}}}$| is the shear-wave velocity (m s⁻¹). For our study, we use |$k = $| 0.21 (Klinger & Werner 2022) and |${{V}_{\mathrm{s}}} = $| 3460 ms⁻¹ (velocity in ak135 to a depth of 20 km, the maximum source depth in our study). Madariaga (1976) demonstrated that f_c is expected to vary with azimuth around the focal sphere, relative to the nodal plane, by a factor of ∼2.5 for S waves. This effect is averaged out by assigning k as 0.21 for S waves. We then use r with |${{\Omega }_0}$| to calculate |${{M}_0}$|⁠:

$$\begin{eqnarray} {{M}_0} = \frac{{4\pi {{\Omega }_0}V_{\mathrm{s}}^3\rho r}}{{F\theta G\left( R \right)}}, \end{eqnarray}$$

(4)

where |$\rho $| is density (kg m⁻³), F is the free surface amplification factor, |$G( R )$| describes the geometric spreading as a function of hypocentral distance |$( R )$| and |$\theta $| is an empirically derived parameter related to the S-wave radiation pattern, averaged over the full focal sphere. We use |$\rho = $| 2800 kg m⁻³, |$F = $| 2 and |$\theta = $| 0.55 (Boore & Boatwright 1984). Note, we assume that site amplification is negligible. The site corrections in this area fluctuate between low values of amplification and deamplification (Edwards et al. 2008), meaning the resultant averaged source estimates are minimally affected by this assumption. For the geometric spreading correction, we follow Herrmann & Kijko (1983):

$$\begin{eqnarray} G\left( R \right) = \left\{ {\begin{array}{@{}*{1}{c}@{}} {{{R}^{ - 1}},\ R < 100\ \mathrm{ km}}\\ {\ {{{\left( {100 \times R} \right)}}^{ - \frac{1}{2}}},\ R \ge 100\ \mathrm{ km}} \end{array}} \right.\mbox{.} \end{eqnarray}$$

(5)

The geometric spreading correction has previously been used by Ottemöller & Sargeant (2010) for UK seismicity. This enables us to calculate M_w (Hanks & Kanamori 1979):

$$\begin{eqnarray} {{M}_\mathrm{w}} = \frac{2}{3}{{\log }_{10}}{{M}_0} - 6.03. \end{eqnarray}$$

(6)

Stress drop, |$\sigma $|⁠, is calculated using:

$$\begin{eqnarray} \sigma = \frac{{7{{M}_0}}}{{16{{r}^3}}}. \end{eqnarray}$$

(7)

And apparent stress, |${{\tau }_\mathrm{a}}$| (Singh & Ordaz 1994), is given by

$$\begin{eqnarray} {{\tau }_\mathrm{a}} = \frac{\sigma }{{4.3}}. \end{eqnarray}$$

(8)

For each event, this methodology results in estimates of M_w and apparent stress at every station. The mean network M_w and apparent stress for each event are determined.

3.4. Spectral modelling results

We estimate stable M_w and apparent stress for 13 UK events (Table 2): spectra-derived M_w values range from 3.40 to 4.45, and apparent stress ranges from 0.32 to 1.74 MPa. An example of our spectral modelling for an earthquake in Cwmllynfell, South Wales (event J) recorded at GB.KESW is shown in Fig. 5. This station is 317 km from the earthquake source and the amplitude spectra is modelled between 0.2–30 Hz (where SNR > 3). The goodness of fit between the observed and synthetic spectra is quantified using the reduced chi-square statistic (⁠|$\chi _{red}^2$|⁠), which suggests reasonable agreement for this station (⁠|$\chi _{red}^2 = $| 1.04). Any stations that demonstrate a poor fit (⁠|$\chi _{red}^2$| >> 1) or a narrow estimation bandwidth (<10 Hz) are removed from the analysis to improve robustness. In total, data from a further nine stations are used for this event, to estimate an M_w of 3.95 and apparent stress of 0.78 MPa. We determine M_w and apparent stress uncertainties by calculating the standard deviation of the population of single station estimates. For this event, the standard deviations of M_w and apparent stress are 0.11 and 0.31 MPa, respectively.

Figure 5.

Spectral analysis of event J using SpecMod. (a) Broad-band unfiltered, transverse component velocity recorded at GB.KESW showing manual Pg and Sg arrival picks and noise (black) and signal (blue) windows. The spectra of these windows are then shown in (b), where the dotted red lines represent the bandwidth (SNR ≥ 3) over which the signal is modelled. (c) Modelled signal spectrum. The signal spectrum is blue over the modelled bandwidth (0.2–30 Hz) and grey outside of this bandwidth. Likewise, the computed model (green) is solid within this bandwidth and dashed outside of it.

Table 2.

Spectral analysis results from SpecMod.

Event ID	|${{f}_{\mathrm{c}}}$| (Hz)	Rupture radius (m)	M0 (N m)	Mw	Apparent stress (MPa)
A	1.88	386.49	1.84E + 14	3.48	0.32
B	1.19	610.59	1.30E + 15	4.04	0.58
C	1.08	672.78	5.21E + 15	4.45	1.74
D	2.03	357.93	2.60E + 14	3.58	0.58
E	1.95	372.62	4.61E + 14	3.75	0.91
F	2.10	346.00	2.87E + 14	3.61	0.70
G	1.74	417.59	4.20E + 14	3.72	0.59
H	2.20	330.27	3.78E + 14	3.69	1.07
I	2.14	339.53	3.34E + 14	3.65	0.87
J	1.47	494.29	9.26E + 14	3.95	0.78
K	2.88	252.29	1.58E + 14	3.44	1.00
L	2.22	327.30	2.28E + 14	3.54	0.66
M	2.86	254.06	1.40E + 14	3.40	0.87

Despite few published apparent stress estimates for events in the UK, there is one event that we can compare to the results of this study. For the 2007 April 28 Folkestone (event B) earthquake, we calculate an apparent stress of 0.6 MPa with standard deviation of 0.5 MPa. For the same event, Ottemöller & Sargeant (2010) estimate apparent stress (converting their stress drop values to apparent stress using equation 8) of 0.7|$\pm $|0.8 MPa. Even though the two studies use different methods and data sets the apparent stress estimate for these this event is comparable.

We were able to use waveforms from the 2008 Market Rasen earthquake and one of it's aftershocks to estimate apparent stress independently through the coda spectral ratio method outlined by Mayeda et al. (2007) and applied to numerous global crustal data sets (Malagnini et al. 2014). The coda spectral ratio method provides very stable, averaged apparent stress estimates that are free of path, site, and source mechanism effects. The SpecMod apparent stress value was lower than the preliminary spectral ratios (Jorge Roman-Nieves personal communication, 2024) that suggest a higher corner frequency, closer to ∼4–5 Hz. Therefore the CCT derived apparent stresses may well be a lower bound based upon preliminary spectral ratios. However, the CCT estimated Mw's will remain unchanged because Mw is computed in a frequency band that is unaffected by the choice of our reference events' apparent stress. The reference M_w's set the low-frequency site terms, whereas the reference apparent stresses set the higher frequency site terms. Therefore, t* and whatever path terms exist at longer periods, are completely accounted for by the use of these independent M_w's.

4 CODA CALIBRATION

4.1. Coda calibration methodology to determine M_w

To estimate M_w, we use the Coda Calibration Tool (CCT; Barno 2017) and follow the method described in Mayeda et al. (2003), which has been widely used in other studies (e.g. Chiang et al. 2021; Morasca et al. 2022; Shelly et al. 2022). Horizontal component velocity seismograms are filtered into 15 narrow passbands between 0.05 and 19 Hz (0.05–0.1 Hz, 0.1–0.2 Hz, 0.2–0.3 Hz, 0.3–0.5 Hz, 0.5–0.7 Hz, 0.7–1.0 Hz, 1.0–1.5 Hz, 1.5–2.0 Hz, 2.0–3.0 Hz, 3.0–4.0 Hz, 4.0–6.0 Hz, 6.0–8.0 Hz, 8.0–10.0 Hz, 10.0–15.0 Hz, 15.0–19.0 Hz). For each frequency band, |${{f}_\mathrm{b}}$|⁠, the two bandpass filtered horizontal components are stacked to produce a single envelope (e.g. Fig. 6a). A model envelope, |$E( {{{f}_\mathrm{b}},x,t} )$|⁠, is then fit to the data using the appropriate propagation distance, |$x( {km} )$|⁠, and propagation time, |$t( s )$|⁠, using:

$$\begin{eqnarray} E\left( {{{f}_\mathrm{b}},x,t} \right) &=& W\left( {{{f}_\mathrm{b}}} \right) \times S\left( {{{f}_\mathrm{b}}} \right) \times T\left( {{{f}_\mathrm{b}}} \right) \times P\left( {{{f}_\mathrm{b}},x} \right) \\&&\times H\left( {t - \frac{r}{{v\left( {{{f}_\mathrm{b}},x} \right)}}} \right) \times {{\left( {t - \frac{r}{{v\left( {{{f}_\mathrm{b}},x} \right)}}} \right)}^{ - \gamma \left( {f,r} \right)}}\exp \\&&\left[ {\beta \left( {{{f}_\mathrm{b}},x} \right) \times \left( {t - \frac{r}{{v\left( {{{f}_\mathrm{b}},x} \right)}}} \right)} \right], \end{eqnarray}$$

(9)

where W is the S-wave source amplitude, S is the frequency-dependent site effect, T is the S-to-coda transfer function, P is the path correction, H is the Heaviside step function, v is the propagation velocity (km s⁻¹) corresponding to the coda envelope peak amplitude, and |$\gamma $| and |$\beta $| are parameters associated with the coda shape (e.g. Mayeda et al. 2003; Morasca et al. 2022; Shelly et al. 2022). The S-to-coda transfer function (T) as well as frequency-dependent site response terms (S) remove the effective t* using independent reference earthquake source estimates (M_w and apparent stress). The coda methodology depends on correction terms in narrow bands, that accounts for path, S-to-coda transfer function, and site effect. The spectral level and moment-rate spectral shape, which is then converted to radiated seismic energy, is then free of all path effects. This has the added benefit of avoiding the issue of the frequency dependence of t*, which can be significant (e.g. Yoshimoto et al. 1993). The reader is directed to Mayeda et al. (2003) for an in-depth description of these parameters. The coda is fit across a specified time length, between the coda envelope peak amplitude and an end limit. The end limit is located immediately before the coda amplitude returns to the pre-arrival noise amplitude, provided there are no aftershocks (i.e. ‘f-markers’ in Fig. 6a). For UK earthquakes, the coda length decreases for higher frequency passbands (Fig. 6a), we typically fit between 200 to 400 s of coda. This allows us to derive moment rate spectra by fitting the modelled moments of each frequency passband to create the spectrum (e.g. Fig. 6b).

Figure 6.

(a) Coda envelope fits for four narrow-band frequencies at station GB.CWF. The f-marker/red bars represent the end of envelope fits, where the synthetics (green) transition from solid to dashed. (b) The resultant moment rate spectra used to calculate M_w for event C using the CCT. The symbols denote coda amplitudes at different stations, which are then averaged (black line) and modelled using the CCT to compute M_w (dashed red line). This model fit is compared to a reference, REF, (dashed black line), representing values of M_w and apparent stress from moment tensor inversion and spectral analysis, respectively. Confidence interval, UQ1 and UQ2, represent values that lie within the standard error and twice the standard error, respectively.

Figure 7.

Comparison of magnitude estimates. Moment tensor inversion (MTtime) M_w estimates versus M_w estimates from spectral analysis (SpecMod) (a) and CCT (b). In both (a) and (b), dotted lines represent values |$\pm $|0.1 magnitude units from the 1:1 line. Data labels correspond to event IDs in Table 2. M_w estimates calculated using CCT versus ML^BGS (c) and ML^P (d). Linear regressions for our data and the data in Green et al. (2020) are computed using the orthogonal distance regression method and shown in black and green, respectively. The dotted line is the 1:1 line.

We perform CCT calibration iteratively (see Fig. 2) by incorporating progressively lower magnitudes (ML^BGS > 3.0) as additional calibration events.This ensures that the initial calibration is generated from the largest magnitude and most reliable events and allows us to assess the robustness of the calibration (by analyzing how this changes with the inclusion of additional events). This iterative approach also allows quality control of the f-marker positions and removal of envelopes with significant noise, therefore the calibration parameters become increasingly robust. Initially this calibration procedure uses just the reference and validation events (Fig. 1 and Table 3) before including progressively more calibration events with each iteration. These calibration events (Fig. 1 and Table 3), which are vital when constraining path terms, were geographically well distributed around the UK, and specifically chosen to improve ray-path coverage across the UK and therefore providing well-calibrated path terms. Once source, site, transfer function, and path effects are corrected for, the envelope amplitudes are still dimensionless. Following Mayeda et al. (2003), the non-dimensional amplitudes are converted to seismic moment, M₀, using previously determined reference events. These events have independently estimated M_w and apparent stress values from our moment tensor inversions and spectral modelling, and are used to constrain the frequency-dependent site terms (e.g. Morasca et al. 2022; Shelly et al. 2022). We use the validation events, which also have independent estimates of M_w, to validate the reliability of the calibration. The peak envelope amplitudes, representing M₀, across all passbands are combined to give a resultant source spectrum, which is fit using a Brune (1970) source model to determine M_w and apparent stress (Fig. 6b). The derived calibration parameters can then be used to directly measure M_w and apparent stress for different events within the calibrated region and potentially for future events in routine processing.

Table 3.

All events processed using the CCT (Coda Calibration Tool). Event origin time, latitude and longitude from BGS earthquake bulletin. Event type includes Validation (V), Reference (R), Calibration (C), CCT processed (P), Quarry (Q), Newdigate (N), and hydrofrac (H) .

Origin time (Year-Month-Day Time)	Lat (⁠|$^\circ $|⁠)	Lon (⁠|$^\circ $|⁠)	MLBGS	MLP	Mw	Apparent stress (MPa)	Event type
2006–12-26 10:40:04	55.09	-3.64	3.5	3.5	3.38	0.16	C/V
2007–04-28 07:18:11	51.08	1.17	4.3	4.3	4.09	0.25	C/V
2008–02-27 00:56:47	53.40	-0.33	5.2	4.8	4.49	2.69	C/R
2008–05-28 20:09:08	54.69	-2.95	2.5	2.5	2.61	0.06	P
2008–09-30 20:46:32	58.08	-3.18	2.5	2.5	2.75	0.02	P
2008–10-10 04:28:39	56.83	-5.53	3.5	3.5	3.30	0.16	C
2008–10-26 18:06:25	52.20	-2.63	3.5	3.4	3.14	0.60	P
2008–11-03 09:53:42	56.37	-5.52	2.5	2.6	2.35	0.15	P
2009–03-03 14:35:55	51.12	1.18	3.0	3.1	2.93	0.15	P
2009–04-11 11:39:07	53.70	-0.25	3.0	2.9	2.84	0.32	P
2009–04-28 10:22:09	54.16	-2.99	3.7	3.7	3.48	0.16	C/R
2009–06-05 19:42:07	51.62	-3.65	2.9	2.9	2.67	0.14	P
2009–09-15 00:30:22	53.74	1.12	3.3	3.1	3.34	0.02	P
2009–09-30 13:04:56	53.70	-0.69	2.9	2.8	2.56	0.11	P
2009–10-06 05:05:55	51.77	-4.10	2.5	2.6	2.44	0.31	P
2010–07-30 23:39:39	51.67	-2.42	2.7	2.6	2.33	0.19	P
2010–11-13 09:37:00	53.22	-1.81	2.2	2.2	2.57	0.01	P/Q
2010–12-03 09:53:00	56.61	-5.53	2.2	2.1	2.31	0.02	P/Q
2010–12-04 01:53:25	53.98	0.87	2.6	2.7	2.69	0.01	P
2010–12-21 22:59:12	54.39	-3.15	3.5	3.5	3.26	0.34	C
2011–01-03 21:03:28	54.17	-1.65	3.6	3.1	3.03	0.23	C
2011–01-23 06:02:49	56.82	-5.78	3.5	3.4	3.26	0.30	C
2011–01-27 06:04:13	50.00	-1.91	2.5	2.4	2.46	0.08	P
2011–05-17 14:14:00	56.60	-5.60	2.2	2.3	2.77	0.01	P/Q
2011–06-23 13:43:38	50.57	-3.73	2.7	2.6	2.62	0.09	P
2011–07-14 06:59:10	50.12	-0.74	3.9	3.6	3.57	0.40	C/V
2011–08-21 08:37:23	56.85	-5.67	2.9	2.9	2.79	0.09	P
2011–09-08 19:02:51	56.59	-5.64	4.5	4.1	3.90	1.13	P
2012–01-18 18:33:29	49.63	-4.92	3.5	3.3	3.32	0.07	P
2012–02-20 05:35:48	55.78	-6.35	2.6	2.6	2.55	0.05	P
2012–02-26 22:31:16	54.65	-0.84	2.9	2.5	2.74	0.01	P
2012–02-29 09:14:26	55.78	-6.34	2.8	2.8	2.76	0.03	P
2012–03-04 23:23:53	56.24	-4.77	2.8	2.9	2.67	0.10	P
2012–06-01 12:16:46	52.41	-2.62	2.8	2.6	2.53	0.18	P
2012–10-31 15:59:19	55.97	-5.94	2.7	3.0	2.63	0.21	P
2012–12-14 23:03:03	50.96	-0.82	2.9	2.9	2.75	0.19	C
2013–01-18 05:20:44	52.80	-1.25	2.9	2.7	2.57	0.13	C
2013–02-27 23:57:01	52.90	-1.04	2.5	2.3	2.31	0.05	P
2013–05-15 17:43:48	57.67	-5.58	2.8	2.9	2.80	0.11	P
2013–05-18 19:18:02	56.78	-5.71	2.9	3.0	2.79	0.12	P
2013–05-29 03:16:28	52.88	-4.72	3.8	3.9	3.59	0.52	C/V
2013–06-26 22:28:01	52.88	-4.72	2.7	2.9	2.86	0.41	P
2013–06-26 22:28:29	52.88	-4.70	2.5	2.9	2.70	0.10	P
2013–07-16 04:04:00	57.72	-5.72	2.8	2.8	2.65	0.13	P
2013–08-25 05:37:48	53.86	-3.38	2.5	2.5	2.66	0.03	P
2013–08-25 09:58:36	53.88	-3.40	3.3	3.3	3.37	0.04	C
2013–08-27 10:06:06	56.64	-4.37	2.8	2.8	2.45	0.16	P
2013–08-31 06:36:11	53.89	-3.40	2.6	2.7	2.72	0.04	P
2014–02-20 13:21:30	51.36	-4.16	4.1	3.7	3.63	0.30	C/V
2014–04-17 06:07:14	52.72	-0.73	3.2	3.0	2.82	0.26	C
2014–04-18 06:50:51	52.72	-0.73	3.5	3.4	3.21	0.15	C
2014–06-18 08:44:39	53.40	-1.38	2.8	2.7	2.55	0.21	P
2014–06-20 16:01:48	55.79	-6.38	2.5	2.6	2.41	0.14	P
2014–07-03 18:36:08	56.82	-5.26	2.9	3.0	2.84	0.16	P
2014–07-11 11:54:32	49.15	-2.41	4.3	3.9	3.85	0.45	P
2014–07-23 16:26:41	49.10	-2.40	3.3	3.1	3.04	0.21	P
2014–07-25 19:05:25	53.72	1.39	2.7	2.6	2.67	0.02	P
2014–10-28 19:16:54	53.06	-1.19	2.6	2.7	2.80	0.03	P
2014–11-14 11:23:58	55.10	-3.65	2.6	2.9	2.58	0.17	P
2015–01-27 18:30:17	51.07	-1.30	2.9	2.8	2.97	0.02	C
2015–01-28 22:25:53	52.73	-0.72	3.8	3.6	3.38	0.22	C
2015–02-25 22:55:07	49.09	-2.38	3.0	2.9	2.86	0.18	P
2015–05-22 01:52:17	51.30	1.44	4.2	4.1	3.69	0.77	C/R
2015–05-26 15:41:03	53.12	-4.36	3.0	3.0	2.73	0.20	C
2015–05-30 19:20:12	54.33	-1.86	2.6	2.3	2.19	0.13	P
2015–09-22 21:40:11	52.70	-0.72	2.8	2.6	2.50	0.15	P
2017–03-03 09:28:11	52.91	-2.15	2.6	2.6	2.44	0.14	P
2017–08-04 14:45:34	56.80	-5.87	3.4	3.5	3.23	0.39	C
2017–11-01 20:59:22	55.88	-5.43	2.6	2.8	2.70	0.08	P
2018–02-17 14:31:07	51.77	-3.83	4.6	4.6	4.03	1.17	C/R
2018–02-28 07:33:51	54.64	-3.34	3.4	3.2	2.99	0.14	C
2018–03-09 08:14:08	52.76	-3.60	2.7	2.7	2.50	0.27	P
2018–04-01 11:10:58	51.16	-0.26	2.6	2.6	3.24	0.01	P/N
2018–04-29 18:19:00	55.89	-5.58	2.6	2.7	2.64	0.04	P
2018–05-01 06:15:59	55.89	-5.58	2.5	2.7	2.48	0.12	P
2018–06-09 22:14:28	53.65	-0.05	3.8	3.7	3.39	0.28	C/V
2018–06-27 12:28:24	51.18	-0.24	2.6	2.6	3.42	0.01	P
2018–07-05 10:53:24	51.15	-0.27	3.0	3.0	3.52	0.01	P/N
2018–08-02 17:42:52	52.61	0.45	2.7	2.5	2.33	0.41	P
2018–08-28 04:59:21	54.42	-0.88	2.8	2.4	2.26	0.24	P
2018–09-15 18:39:09	54.57	-1.64	3.1	2.8	2.64	0.15	C
2019–02-27 03:42:21	51.16	-0.25	3.1	3.2	3.47	0.01	C
2019–04-08 01:49:27	52.84	-2.60	2.5	2.6	2.58	0.07	P
2019–08-26 07:30:47	53.79	-2.96	2.9	2.6	2.78	0.02	C/H
2019–12-05 22:49:18	51.12	-2.98	3.2	3.2	2.84	0.17	C
2020–01-23 05:57:00	54.59	-1.31	3.1	2.9	2.74	0.09	C
2020–08-21 17:41:52	57.18	-5.77	2.6	2.6	2.49	0.20	P
2020–09-08 08:45:28	51.93	-0.74	3.5	3.6	3.51	0.07	C/V
2020–09-22 08:32:15	51.92	-0.68	3.0	3.0	3.00	0.06	C
2020–10-04 18:43:42	56.27	-3.75	2.5	2.7	2.71	0.03	P
2020–10-21 02:49:34	53.35	-0.57	2.5	2.3	2.13	0.13	P
2021–05-15 15:02:47	58.55	-4.72	2.5	2.6	2.62	0.07	P
2021–05-16 15:16:49	51.60	-2.80	2.9	3.0	2.78	0.13	P
2021–07-26 10:01:21	55.81	-3.20	2.5	2.6	2.45	0.22	P
2021–11-16 01:44:07	56.03	-5.54	3.4	3.4	3.14	0.10	C
2021–12-11 17:32:58	57.85	-5.11	2.5	2.5	2.38	0.14	P
2021–12-11 17:37:53	57.85	-5.11	2.5	2.5	2.43	0.11	P
2022–02-21 22:59:29	52.54	-1.94	2.8	2.8	2.53	0.23	P
2022–05-30 14:36:57	52.80	-2.62	3.8	3.7	3.38	0.57	C/V
2022–06-01 01:51:22	53.80	-2.24	2.7	2.7	2.47	0.12	P

4.2. Coda calibration results

Our calibration data set includes 33 geographically separated UK events (Fig. 1 and Table 3) with ML^BGS ranging between 2.9 and 5.2. Of these 33 events, eight are validation events and four are reference events. Event I was not used as a validation or reference event, due to a significant aftershock present in the waveform data that could influence the calibration. We limit the number of reference events to four to avoid path bias resulting from too many paths from a particular region (Holt et al. 2021). The resultant M_w and apparent stress values of our calibration data set range from 2.57 to 4.49 and 0.01 to 2.69 MPa. The remaining 21 events in the calibration data set are calibration events used to improve the path terms.

Coda envelopes and moment rate spectra for validation event C, the 2008 February 27 Market Rasen earthquake, are shown in Fig. 6. These envelopes exhibit low noise, and the synthetic envelopes compare well to the observed. The subsequent moment rate spectra agree with the validation spectra, and the difference in estimated M_w between the coda estimate and validation estimate is 0.03. The CCT quantifies uncertainty by calculating confidence intervals for both M_w and apparent stress estimates based on the misfit of the inversion. The confidence interval limits are labelled UQ1 and UQ2, representing values that lie within the standard error and twice the standard error, respectively. In Fig. 6(b), UQ1 and UQ2 follow the model fit closely, suggesting the uncertainty for this event is low. When compared to model fits in the calibration data set, the mean M_w UQ2 lower and upper values are [-0.04, +0.04] and the mean apparent stress UQ2 limits are [−0.10, +0.12] MPa.

Once the CCT calibration is finalized, we use the resultant calibration parameters (i.e. |$W( {{{f}_\mathrm{b}}} )$|⁠, |$S( {{{f}_\mathrm{b}}} )$|⁠, and |$P( {{{f}_\mathrm{b}},x} )$| in eq. 9) to calculate M_w and apparent stress for an additional 64 UK events with ML^BGS|$\ge $| 2.5. These coda-derived M_w and apparent stress estimates range from 2.13 to 3.90 and 0.01 to 1.17 MPa. Considering a similar magnitude range for UK seismicity, our CCT-derived apparent stress estimates are comparable to those reported by Edwards et al. (2008). In addition to natural seismicity we also calculate CCT-derived M_w and apparent stress for quarry events, a hydrofracture event and events from the 2018–2019 Newdigate swarm (Hicks et al. 2019). We analyse three ML^BGS 2.2 quarry events located in NW Scotland and in central England (Fig. 1). CCT determines M_w values of 2.31, 2.57 and 2.77, and apparent stresses of 0.02, 0.01 and 0.01 MPa (Table 3) for these events. However, signals for these events have low SNRs at frequencies less than 2 Hz and are only observed at limited, close by seismic stations, therefore hindering the ability to model the spectra. Furthermore, it has been shown (e.g. Murphy et al. 2009) that shallow earthquakes and explosions produce near-surface scattering that results in a peak in the low-frequency spectra, resulting in overestimated M_w. This peak is observed in the spectra for the three quarry events, and therefore suggests our M_w estimates are likely to be too large.

In 2018, two horizontal wells were drilled, to depths of ∼ 2 km, into the gas-bearing Bowland Shale at Preston New Road, UK. Hydraulic fracturing operations at the second of these wells took place between 15–21 and the 23rd August 2019 (Kettlety et al. 2020). At 07:30 (UTC) on 26th August 2019 an ML^BGS 2.9 earthquake occurred near the Preston New Road site, the largest of 135 induced events recorded in this region in 2019 (Galloway 2020). Using CCT we determine a moment magnitude for the ML^BGS 2.9 hydro fracture event on the 26th August 2019. This is likely to be an over estimate of the magnitude due to the event having a depth shallower than 5 km, which produces peaks in the low-frequency spectra.

We also calculate moment magnitudes for three events of the 2018–2019 Newdigate earthquake swarm (Hicks et al., 2019). Similar to the quarry and hydro fracture events the CCT moment magnitudes for the Newdigate swarm is likely to be too high due to peaks in the low-frequency spectra caused by the shallow depths of the events. This earthquake swarm occurred within the aseismic Weald Basin, in southeast England. These shallow earthquakes (<3 km depth) received significant public and media attention due to their proximity (∼3 km) to operating oilfields in the area and the possibility that they were induced by this subsurface industrial activity. After detailed analysis of seismic observations offered by the installed temporary seismic network and reported operational data provided by the operators, Hicks et al. (2019) suggest these events are unlikely induced from near by industrial activity and are unique examples of very shallow, natural earthquakes within a sedimentary basin. However, this conclusion relies on reported operational data provided by the operators.

5 DISCUSSION

For a subset of 12 events (Sections 3.2–3.4) we have determined M_w using three independent methods: moment tensor inversion, spectral analysis, and coda calibration, allowing us to compare the results of each method. Additionally, this study is the first to analyse seismic coda for UK seismicity, providing a unique opportunity to compare coda observed for UK seismicity with coda observed in different tectonic settings.

5.1. Magnitude comparisons

Moment tensor estimated M_w compares well to spectral analysis estimated M_w (Fig. 7). The mean difference between calculated M_w is 0.06, with a maximum difference of 0.10. The good agreement between M_w estimates using two independent methods suggests they are reliable. Additionally, we compare our coda calibration M_w estimates to those determined from moment tensor inversion (Fig. 7b). These M_w estimates are also in good agreement, with a mean difference of 0.03 and maximum of 0.08, indicating that the CCT has been effective in calibrating moment magnitudes across the UK.

We investigate the relationship between M_w with two local magnitude scales: ML^BGS (a local magnitude scale based on S-wave amplitudes) and ML^P (a local magnitude scale based on P-wave amplitudes). Using M_w estimates from spectral modelling (Sargeant & Ottemöller 2009 and Ottemöller & Sargeant 2010), Green et al. (2020) quantified the relationship between M_w with both ML^BGS and ML^P for UK events. Here we examine whether CCT M_w estimates produce similar relationships with ML^BGS and ML^P to those reported by Green et al. (2020). Calculating M_w using spectral modelling is reliant on several parameter choices (e.g. |${{V}_{\mathrm{s}}}$|⁠, k, |$\rho $| and |$\theta $| in eqs 3–6). The previously reported M_w estimates for the UK (Sargeant & Ottemöller 2009 and Ottemöller & Sargeant 2010) use different values for these parameters to those used in this study. Therefore, for consistency between the two M_w data sets, these previously reported M_w estimates (Sargeant & Ottemöller 2009; Ottemöller & Sargeant 2010) are recalculated using the same parameter values as defined in this study and termed |${{M}_{W,SO}}$|⁠. Using orthogonal distance regressions, where the quoted uncertainties are |$\pm $| 1 standard error, the relationships between |${{M}_{W,SO}}$| with ML^BGS and ML^P are:

$$\begin{eqnarray} {{M}_{\mathrm{ w,SO}}} = \left( {0.71 \pm 0.05} \right)\mathrm{M}{{\mathrm{L}}^{\mathrm{BGS}}} + \left( {0.77 \pm 0.16} \right),2.7 \le \mathrm{M}{{\mathrm{L}}^{\mathrm{BGS}}} \le 5.2 \end{eqnarray}$$

(10)

$$\begin{eqnarray} {{M}_{\mathrm{ w,SO}}} = \left( {0.91 \pm 0.05} \right)\mathrm{M}{{\mathrm{L}}^{\mathrm{P}}} + \left( {0.21 \pm 0.15} \right),2.4 \le \mathrm{M}{{\mathrm{L}}^{\mathrm{P}}} \le 4.8, \\ \end{eqnarray}$$

(11)

ML^P is estimated for all events in the CCT data set using the methodology of Green et al. (2020); the half peak-to-trough displacement values of the P-wave packet are measured, corrected for both propagation path and site effects, and resultant station magnitudes averaged across the recording network. The relationship between CCT M_w with ML^BGS and ML^P is calculated using the same regressions as before:

$$\begin{eqnarray} {{M}_\mathrm{w}} = \left( {0.79 \pm 0.04} \right)\mathrm{M}{{\mathrm{L}}^{\mathrm{BGS}}} + \left( {0.50 \pm 0.13} \right),2.2 \le \mathrm{M}{{\mathrm{L}}^{\mathrm{BGS}}} \le 5.2, \end{eqnarray}$$

(12)

$$\begin{eqnarray} {{M}_\mathrm{w}} = \left( {0.88 \pm 0.04} \right)\mathrm{M}{{\mathrm{L}}^{\mathrm{P}}} + \left( {0.26 \pm 0.12} \right),2.1 \le \mathrm{M}{{\mathrm{L}}^{\mathrm{P}}} \le 4.8, \\ \end{eqnarray}$$

(13)

Comparisons of Equation 10 with Equation 12, and Equation 11 with Equation 13 indicate that our results, using updated |${{M}_{W,SO}}$|values and a significantly expanded data set (from 46 events to 100 events), are comparable with those of Green et al. (2020) (Figs 7c and d). However, our observations indicate that the gradient of the M_w: ML^P relationship is not consistent with a direct scaling between ML^P and M_w (i.e. the gradient is not 1). Therefore this does not support the tentative hypothesis of Green et al. (2020) that ML^P scales with M_w due to the P-wave magnitude scale being less influenced than ML^BGS by estimation errors in along-path attenuation. In addition, by extending our results to lower magnitudes than Green et al. (2020) we observe more clearly that ML^BGS estimates are greater than M_w estimates at low magnitude, as is expected due to anelastic attenuation modifying the recorded spectra for smaller earthquakes (Deichmann et al. 2017; Dost et al. 2018).

5.2. Scaling of apparent stress and M_w

We analyse the relationship between our CCT estimates of moment magnitude and apparent stress (Fig. 8). The CCT apparent stress is calculated by integrating the squared velocity spectrum, where the frequency range is extrapolated to zero and infinity at either end. In Fig. 8, we include only events for which the extrapolated region is less than 45 per cent, meaning at least 55 per cent of the energy is observed. In all cases, this ensures the corner frequency is within the observed frequency range. Additionally, we categorize the UK events by location and origin time, allowing us to isolate the natural earthquakes from those related to anthropogenic sources (e.g. quarry and hydrofracture events). This allows us to focus only on sources where a Brune source model is applicable and thus provide more reliable apparent stress estimates. For comparison, we also include CCT results from Eastern Canada (Bent et al. 2022) and Central Italy (Morasca et al. 2022 ).

Figure 8

(a) Relationship between CCT moment magnitude and apparent stress estimates. Results include UK events analysed in this study (blue circles are onshore and English Channel earthquakes and orange symbols are North Sea or non-earthquake sources), events from Eastern Canada (black diamonds) and Central Italy (grey squares). (b) Map showing BGS locations for UK events, where symbol size is scaled by CCT M_w.

Fig. 8 suggests two principal conclusions. First, events with larger M_w result in larger values of apparent stress. This relationship is most apparent for the UK and Central Italy CCT results. This conclusion is tentative however, since the trend of M_w with apparent stress has previously been explained by frequency limitations for lower magnitude events, where lower magnitude events are deficient in low frequencies and therefore appear to have lower apparent stress due to observational limitations. Secondly, the CCT apparent stress estimates from different tectonic regimes form distinct groups. The apparent stress values from events in Eastern Canada are approximately an order of magnitude greater than events of comparable moment magnitude in Central Italy. The UK events lie between Eastern Canada and Central Italy and are generally closer to the Eastern Canada events.

5.3. Comparisons of coda in other tectonic regimes

We compare our coda envelopes to those from events in previously well-studied regions, including the Western US, Eastern Canada, and Italy (Fig. 9). Each event is recorded by a broad-band station approximately 100 km from a source of similar magnitude; events in the UK, Western US, and Italy are estimated as M_w ∼4.5, whilst the Eastern Canada event has an M_w estimate of 4.97 (Table 4).

Figure 9.

Coda envelopes from different regions for a range of frequencies. Representative events from different areas (see Table 4) are of comparable magnitude and recorded by broad-band stations all approximately 100 km from each event.

Table 4.

Event and station information for the coda shown in Fig. 9. M_w estimates and event information are from the Saint Louis University Moment Tensor Determinations and the Northern California Earthquake Data Center Catalog (Western US event only).

Region	Origin time (Year-Month-Day Time)	Lat (⁠|$^\circ $|⁠)	Lon (⁠|$^\circ $|⁠)	Mw	Station code	Station Lat (⁠|$^\circ $|⁠)	Station Lon (⁠|$^\circ $|⁠)	Event-station distance (km)
Eastern Canada	2002-04-20 10:50:48	44.53	−73.73	4.97	MNT	45.50	−73.62	108
Western US	2006-08-03 03:08:13	38.35	−122.57	4.50	JRSC	37.40	−122.24	110
UK	2008-02-27 00:56:47	53.40	−0.33	4.7	CWF	52.74	−1.31	98
Italy	2016-09-03 01:34:13	42.77	13.13	4.3	CAFI	43.33	11.97	113

At low frequencies (0.05–0.10 Hz), the coda envelopes are similar between different regions, and the steepness of the coda decay is comparable for each example event. However, as frequency increases, the shape of the coda envelopes become increasingly distinct and appear to divide into two separate categories: (a) the UK and Eastern Canada, and (b) Italy and Western US, where the decay of coda envelopes at high frequencies (>8 Hz) are much faster in Italy and Western US compared to the UK and Eastern Canada.

The differences in envelope shape can be explained by regional variations in attenuation (e.g. Mayeda et al. 1992). Coda Q describes the decay rate of narrow-band local coda under the assumption of the single scattering model, first defined by Aki (1969) and subsequently used to characterize tectonic variation in a various regions (e.g. Havskov et al. 2016). Coda generation is due to 3-D scattering of seismic waves, which primarily results from subsurface impedance (Sato et al. 2012). Variation in Q has been related to a variety of factors, including water content in the crust, heat flow, the degree of tectonic activity, and crustal age (e.g. Hauksson & Shearer 2006). Italy and Western US are associated with high heat flow and frequent tectonic activity, resulting in low Q (e.g. Phillips et al. 1988). In contrast, the UK and Eastern Canada have low heat flow and low levels of tectonic activity, resulting in high Q (e.g. Woodgold 1990). The effect of these factors on coda decay increases with frequency, as seen in Fig. 9, with differences in envelope shape being most apparent at 15–19 Hz. This is because lower frequency waves sample deeper parts of the lithosphere, which are more homogeneous than the shallowest parts of the lithosphere that are sampled by high frequency waves. These regional differences in coda are consistent with previous local magnitude attenuation parameters calculated for the UK, which fall between the rapid amplitude reductions with range observed in Southern California and the slower amplitude decay observed in intra-plate areas such as Norway and the North-eastern US (Ottemöller & Sargeant 2013). Low stress drop events are often located in areas of low Q, less competent rock, and higher heat flow (e.g. Zoback 1991). For example, in the Western US, a region of high heat flow, stress drop values range between 0.1 and 4.7 MPa for events of M_w|$\le $|4.5 (Mayeda et al. 1996). In contrast, the UK is an area of high stress drop (and apparent stress), with values of stress drop up to approximately 12 MPa for a similar magnitude range (e.g. Edwards et al. 2008). These significant differences in coda decay between regions highlights that region-specific coda calibration is necessary; hence the derived calibration parameters in this study will be useful for future studies of seismic coda in the UK.

6 CONCLUSION

M_w’s have been calculated for 100 UK seismic events post 2006 (with 2.13|$\le $|M_w|$\le $|4.49) using coda measurements. The data analysis, including the calibration of path and site effects were undertaken using the open-source software CCT (Barno 2017) based upon the methodology of Mayeda et al. (2003).

Our analysis calibrates coda at 16 broad-band UK seismic stations, which were chosen for their wide geographical coverage and long deployment history. To achieve the calibration of path and site effects we use 33 events distributed across the UK. To calibrate source size effects, M_w values for a subset of 12 events were independently estimated using moment tensor inversion. We apply the calibrated magnitude estimation technique to a further 67 events. The CCT derived Mw values are within ± 0.1 magnitude units of these estimated for the validation events, providing confidence that the calibration is robust.

Our extended M_w catalogue for the UK allows us to compare the estimated M_w values with local magnitude scale estimates (ML^BGS, ML^P) at lower magnitudes than in previous studies (e.g. Green et al. 2020). This study confirms that M_L is greater than M_w for small magnitude UK earthquakes, which is expected from both theory and observations (e.g. Deichmann 2017). UK earthquake coda decay rates are comparable to those in other intraplate (e.g. Eastern Canada) regions rather than tectonically active regions (e.g. Western US, Italy).

The UK coda derived M_w catalogue has the advantage over previous studies that it has been generated using a common approach. The derived coda calibration parameters will allow rapid M_w estimation for future UK events to lower magnitude values (i.e. M_w > 2.2) than is currently feasible with moment tensor inversion techniques.

ACKNOWLEDGEMENTS

We would like to thank the BGS for providing the data and the CCT working group for their support and feedback throughout the duration of the study. We would also like to thank AWE for supporting and facilitating the completion of this project, with a specific mention to their support of Charlie Peach. The authors thank Jorge Roman-Nieves for his spectral ratio calculations for the 2008 Market Rasen earthquake. The authors thank Katherine Whidden and one anonymous reviewer for their comments and suggestions on this article.

DATA AVAILABILITY

Seismic data used in this study are available from ORFEUS (https://www.orfeus-eu.org) and IRIS (http://ds.iris.edu/ds/nodes/dmc/). MTtime (v1.1.0) was used for moment tensor inversion (https://github.com/LLNL/mttime). SpecMod (v0.1.1) was used for spectral modelling (https://github.com/uofuseismo/SpecMod). The Coda Calibration Tool (v1.0.18.3) was used to calibrate the UK and create the M_w catalogue (https://github.com/LLNL/coda-calibration-tool). Event information for international earthquakes (Table 4) are from the Saint Louis University Moment Tensor Determinations (last accessed 2022 August 23) (https://www.eas.slu.edu/eqc/eqcmt.html) and the Northern California Earthquake Data Center Catalog (last accessed 2022 August 23) (https://www.ncedc.org). Event information for UK earthquakes (Tables 1–3) comes from the BGS seismic bulletins (last accessed 2022 August 23) (https://www.earthquakes.bgs.ac.uk/publications/bulletins/bulletins_list.htm).