Moment tensor event identification for collapses

SUMMARY

We introduce a seismic identification method for collapse events using moment tensors (MTs). We start by computing full (six-element) MT solutions for 43 identified collapse events from around the world, and statistically characterizing the population on the MT hypersphere. We then test a large data set of over 1000 full MTs for the western U.S. against the distribution of collapses using a MT-based identification method similarly as used for testing explosions. Known collapses and explosions are readily identified, along with other anomalous events in the Geysers and central California coast. Misidentification rates are determined for various screening angles with optimal misidentification rates between earthquakes and collapses on the order of 3 per cent. The method is demonstrated to be very effective at identifying non-earthquake sources with a 97–98 per cent accuracy. It is likely to be transportable to other regions, and can be used for event identification anywhere full MT solutions are routinely calculated.

INTRODUCTION

Moment tensor (MT) solutions describe the point-source force couples of seismic sources. As such, they are ideally suited to characterize the type of event, be it an earthquake, explosion or the collapse of an underground cavity. Although MTs are now routinely determined for many events at the global and regional scales, improvements in methodology (e.g. Minson & Dreger 2008) have enabled the calculation of full, six-element MT solutions, which do not make the assumptions often made about non-double couple mechanisms and volume change. The additional degrees-of-freedom, however, require improved inversions, including better azimuthal coverage and more accurate velocity models.

For most earthquakes, the full MTs are generally still primarily double couple solutions that represent slip on a fault plane, with the solutions providing relevant source information, including the earthquake size (seismic moment), orientation, and sense of slip (normal, reverse, strike-slip and oblique). For non-double couple sources, full MTs are necessary in order to capture physical mechanism beyond simple planar slip, such as crack openings/closings, or the volume changes caused by explosions and implosions (e.g. Kanamori et al. 1993; Shuler et al. 2013a, b; Rodríguez-Cardozo et al. 2021).

U.S. Government reports (NAS 2002; NRC 2012) have raised the prospect of using mining operations and large chemical explosions to mask an underground nuclear explosion, a scenario known as mine masking. The location of many observed collapses at known nuclear test sites further increases the chance that they may be confused with explosions, as recently demonstrated by a collapse following the 2017 declared North Korean nuclear test (e.g. Walter et al. 2018). Like explosions, collapses are relatively shallow compared to earthquakes and lead to shallow event hypocentres and exhibit similar signatures of shallow events, such as spectral peaking, an observed low-frequency bump in the source spectra of explosions and shallow-source earthquakes thought to be related to Rg-to-S scattering (Murphy et al. 2009).

There is an expectation that source mechanisms for mining events termed rock-bursts, strain-bursts or coal-bumps have explosive isotropic MT components since rocks expand during failure. However, seismic analysis of mining events from inverted in-mine seismic records show significant implosive MT components (e.g. Julia et al. 2009). Malovichko (2020) demonstrated through comparisons between inversions and simulations that the near-source excavation needs to be included in the Green's functions in order to correctly resolve these deformational forces. Malovichko (2020) also states that not taking the excavation into account with complex GFs is appropriate when the excavation is considered a part of the overall deformation for the source mechanism given that the wavelengths of simulated or inverted seismic waves exceed the size or volume of deformation. A mining rock-burst by itself, involves the pre-existing mining cavity since the cavity changes the local stress field. Whereas the rock-burst, as precursor or simultaneous subevent with the collapse, directly involves the pre-existing mining cavity and the collapse mechanism dominates at lower frequency when it has a larger dimension than the rock-burst.

A regional long-period MT inversion is a point source method that does not solve for a complex source mechanism, using wavelengths greater than 10 km, much longer than the largest cavity dimension. A regional MT solution also does not resolve the details of the source mechanism including those described in Hasegawa et al. (1989) and Malovichko (2020) as it uses lower-frequency (<0.2 Hz) seismic waves recorded by sensitive broadband sensors at far-field distances (>100 km) as opposed to near-field (1–10 Hz) waves from geophones inside the tunnels. As long as this limitation in regional MT analysis is recognized, then the method still lends itself to far-field monitoring for source identification.

Although MTs are often considered an attribute of large events, the regional MT inversion approach used here has been used in near real-time by regional seismic networks over the last 30 yr (e.g. Pasyanos et al. 1996; Clinton et al. 2006; Herrmann et al. 2011) for magnitudes ranging between 3 and 6. Furthermore, using shorter periods, full MTs have been calculated for mine events as small as magnitude one (Caputa et al. 2021), small (M_w < 2.5) volcanic events (Alvizuri et al. 2018), and explosions with M_w < 0.5 and yields as low as 100 kg TNT equivalent (Pasyanos & Chiang 2022).

Here, we have assembled a database of 43 full MT solutions for known collapses from around the world. We use this data set to characterize the population of collapse events, against which we can compare other events. We then test a data set of over 1000 full MTs for the western United States against the populations of collapses and explosions where we search for and find anomalous events. The end result is an effective and well-understood ability to characterize events as collapses which, along with analogous identification methods for explosions, can be used for the routine event identification of primary seismic source types.

METHOD

As a general descriptor, seismic MTs have the capacity to provide detailed event characteristics. In previous work, Ford et al. (2009a) reliably identified isotropic events using full MT inversions. Cesca et al. (2013) established a procedure for analysing natural and induced events in Europe by full MT inversion and decomposition. More recently, Ford et al. (2020) introduced a statistical method for screening on the MT hypersphere. A sample of explosions was used to characterize the population, and screening is achieved by calculating the angle on the hypersphere between the explosion population mean and any new events. A significant test of the robustness of this method was performed using the full MT database of 130 nuclear and 10 chemical explosions at the Nevada National Security Site (NNSS) of Pasyanos & Chiang (2022). We assessed the performance of the discriminant by calculating the misidentification rates of both explosions and earthquakes as a function of the screening angle, which is a distance on MT hypersphere. Where the two misidentification rates (of earthquakes wrongly characterized as explosions, and explosions wrongly characterized as earthquakes) were equal, the misidentification rate was about 5 per cent, showing it to be an effective discriminant between these populations.

In this study, we perform a similar analysis for collapses by performing full MT solutions for dozens of events, statistically characterizing the population, testing it against a large data set, and assessing its performance. We have assembled a collection of 43 collapses from around the world (Fig. 1), including eastern (Pennsylvania, Virginia) and western (Nevada, Utah, Wyoming) United States, Europe (Germany, Poland, Sweden), Russia (Kola Peninsula, Urals, Siberia), South Africa, China, North Korea and Australia. Some of the events, such as those in Nevada and North Korea, are collapses following nuclear explosions whereas others result from mining activities (e.g. Pechmann et al. 1995; Walter et al. 2018).

Information including origin time, location and description of events were taken from multiple sources including the Human-Induced Earthquake (HIQuake) database (e.g. Wilson et al. 2017; Foulger et al. 2018), internet news reports, and mine operators for mine collapses, and the document NV-209 (DOE 2015) and Springer et al. (2002) for collapses associated with U.S. nuclear tests.

For each of these events, we perform full MT inversions using the MTINV code (Ichinose et al. 2003). The time-domain regional MT inversion method solves for the MT using local to regional distance (<∼1000 km) long-period (5–100 s) body and surface waves. We use the formulation from Jost & Herrmann (1989) for the deviatoric MT and formulation from Herrmann & Hutchensen (1993) and Minson & Dreger (2008) for the full, six-element inversion. Estimating point source MTs from long-period filtered three-component (vertical, radial, transverse) waveforms are ideal because the wavelengths are significantly larger than the asperities of the fault rupture by 10–100 times. The wavelengths are also larger than the crustal heterogeneities, only propagating for a few wavelength cycles between the source and receiver allowing for use of calibrated 1-D earth models.

The displacement Green's functions are computed using the frequency–wavenumber reflectivity method (e.g. Wang & Herrmann 1980; Saikia 1994; Zeng & Anderson 1995). The linearized inversion method is written here as the forward problem,

where $\overrightarrow{u}$ is the ground motion vector, $\mathit{G}$ is the Green's function (GF) matrix and

are the six MT elements that are formulated to be linearly related to the GFs. These GFs are computed in ground motion units of displacement (or occasionally velocity) for a source–receiver distance and depth using a 1-D layered earth velocity models depending on tectonic region. These regions include the WUS model (Ritsema & Lay 1995) for the western United States, CUS model (Herrmann et al. 2011) for the central and eastern United States, SAFRICA0 (Qiu et al. 1996) for South Africa, MDJ2 (Ford et al. 2009b) for the Korean Peninsula and finally ak135 model (Kennett et al. 1995) for tectonically stable regions globally.

The low-frequency (f) corners of the filter passband can be different for each station and event to optimize the signal-to-noise depending on noise level, source magnitude, and the response of the broadband sensors. Typically, low-f corners range between 0.01 and 0.05 Hz. The high-f corner of the passband depends on the quality of the earth model and source-to-receiver distance, where calibrated 1-D models at local distance allow for low-f corner ∼0.2 Hz, while default models only fit up to f ∼ 0.05 Hz. MTs can be estimated using a single station (e.g. Kim & Kraeva 1999); however, at least two to three stations optimally distributed around the focal sphere leads to better constrained solutions (e.g. Ichinose et al. 2003), see e-supplement for the number of stations used for each event. We finally iterate the MT inversion over a range of depths (from 1 to 33 km, every 1 km) and origin-times (shift up to 10 s) to solve for the best-fitting solution based on the highest variance reduction. In all cases, the sources are found to be very shallow and origin time shifts are generally found to be small.

Fig. 2 is an example of waveform modelling and MT inversion for a mine collapse in Russia in 2013. The waveform fits have high (79 per cent) variance reduction, and the resulting solution is 32 per cent double couple (reverse mechanism), 21 per cent compensated linear vector dipole (CLVD) and 47 per cent (negative) isotropic. As is typical for collapse events (e.g. Chiang et al. 2018), the diagonal MT elements (M_xx, M_yy, M_zz) are large and negative, indicating volume reduction. This MT sits at lune coordinates of 50°S and 11°E, plotting close to the position for a closing crack. A table of MT solutions for all of the collapses is provided in Table 1.

Ford et al. (2020) investigated source type screening using the seismic MT and its description on the hypersphere, specifically the 5-sphere 𝕊⁵ of 6-degree unit vectors representing the normalized symmetric MT (Tape & Tape 2015). Following their approach we assume the data set of collapse events is a random sample from the population of collapse seismic MTs and use the sample to parametrize a unimodal distribution on the hypersphere, namely the p-dimensional Langevin (von Mises–Fisher) distribution M_p(μ,κ) which can be described by a location parameter μ and a concentration parameter κ. The explosion population was previously estimated as ${\hat{\mathit{\mu }}}^T$ = (0.450, 0.524, 0.713, 0.0272, 0.0245, ­−0.112) and $\hat{\kappa }$= 73.7. After verifying that the sample is well described by the Langevin distribution we find the maximum likelihood estimates for the collapse population based on the 43 collapses here is ${\hat{\mathit{\mu }}}^T$ = (­−0.333, −­0.344, −­0.873, 0.0663, −­0.0683, −­0.0111) and $\hat{\kappa }$= 64.8. The collapse population is, in some ways, the opposite of the explosion population. The first three terms, which correspond to the diagonal components of the MT, are large and of opposite polarity, with the third term, which corresponds to the M_zz component the largest value. The three other terms have significantly smaller amplitudes. The smaller κ value indicates that the explosions have a slightly tighter distribution than the collapses.

RESULTS

We have plotted the 43 collapse events, along with 140 explosions and ∼1000 presumed earthquakes, on the fundamental lune of the MT eigensphere (Tape & Tape 2012; Fig. 3). Events range in magnitude from M_w 3.4 to 5.5 for the collapses, M_w 0.4 to 5.6 for the explosions and M_w 2.1 to 7.0 for the earthquakes. Collapses associated with underground nuclear explosions are plotted black. Although based on a small number of events, these particular events don't appear to be fundamentally different than the rest of the population. For explosions, we use the MT database for explosions from Pasyanos & Chiang (2022). For the earthquake population, we use the same full MT solution data set that was used for event identification in Pasyanos & Chiang (2022). This includes the data sets of Dreger et al. (2000), Minson & Dreger (2008), Ford et al. (2009a), Boyd et al. (2015) and additional MTs calculated at LLNL. Events are distributed across the region, but are concentrated in California. The resulting collection is a data set of over 1000 full MT solutions in the western United States.

We test these events against both the estimated collapse and explosion populations, which are shown on the maps in Fig. 4. The seismic MT screening statistic is defined as the angle δ from the population mean vector $\hat{\mathit{\mu }}$ to a newly estimated unit MT $\mathit{x}$

It is important to note that the angle is calculated on the hypersphere. Since the hypersphere cannot be rendered in a 2-D plot, the fundamental lune (Tape & Tape 2012) shown in Fig. 3, which is a particular reduced-dimension representation, is simply shown for convenience. For every event in our study, the screening angle is calculated against both the collapse and explosion population. Fig. 4 plots the smaller of the two angles on the respective colour scales. As a result, explosion-like events plot as redder symbols, collapse-like events plot as blacker symbols and events that are the least explosion- and collapse-like as the purest green. Most events on the map are bright green, indicating that the event is inconsistent with both collapse and explosion mechanisms, and is hence very ‘earthquake-like’. Exceptions are for known collapses (squares) in Nevada, Utah and Wyoming and underground explosions (stars).

In Fig. 4, we also focus on three subregions where many of the non-earthquake populations are clustered. The first is at the Nevada National Security Site (NNSS) in southern Nevada (Fig. 4b) where the United States conducted numerous underground nuclear explosions up until 1992, as well as a number of chemical explosions in more recent years. At NNSS, we see red symbols associated with explosions in Yucca Flat in the east and mesas to the northwest, consistent with previous observations. In addition, we see earthquakes from ambient seismicity in the region, as well as several collapses associated with underground tests. The second is at The Geysers (Fig. 4c), the world's largest geothermal field, located north of San Francisco Bay. Here, we see that the earthquake-like ambient seismicity of the region is punctuated by explosion-like earthquakes in the immediate vicinity of the Geysers, presumably as a result of crack openings from geothermal production (Boyd et al. 2015).

Perhaps the most unexpected find was the collapse-like events in the central California coast located between San Francisco and Los Angeles (Fig. 4d). The events in the region are associated with the 2003 M_w 6.6 earthquake near San Simeon, CA. One possibility, of course, is that the large non-double couple components are simply noise caused by any number of reasons including poor velocity model, noisy waveforms, poor coverage, origin mislocations, etc. It would be unusual, however, for this to be so localized and we have no reason to suspect that the velocity models are poor for this region in particular. Analyses of this well-studied earthquake sequence indicate that it was a shallow, reverse-motion event with long (>20 km) along-strike extended slip (e.g. Rolandone et al. 2006). Perhaps complex faulting contributes to having large non-double couple components (e.g. Shuler et al. 2013a,b; Contreras-Arratia & Neuberg 2020; Rodríguez-Cardozo et al. 2021). Other analysis indicated that the sequence ruptured within the weak, water-saturated Franciscan complex (Hauksson et al. 2004). Perhaps the shallow thrust mechanism closed near-surface cracks, resulting in uncompensated closures that could be the cause of the observed significant increase in streamflow (Wang et al. 2004). Detailed analysis of the earthquake sequence, however, is beyond the scope of this paper.

We calculate the misidentification rates against both explosions and collapses in Fig. 5. In Fig. 5(a), we are testing events against an explosion screening. The explosion curve shown in red are Type II (that is, false negative) errors, whereas the earthquake (green) and collapse (black) curves are Type I (false positive) errors. When the screening angle is set to 0, then no events are categorized as explosions since no events satisfy δ < A. This results in a 100 per cent misidentification rate for the explosions and a 0 per cent misidentification rate for the other two populations. As the tested screening angle increases, then increasing number of explosions are properly categorized as explosions (true positive), but earthquakes and collapses may be improperly categorized as explosions (false positives). In Fig. 5(b), we are testing events against a collapse screening (δ < B). In this case, the collapse curve shown in black are now Type II errors, whereas the earthquake (green) and explosion (red) curves both represent Type I errors.

As seen in Fig. 3, there are generally fewer presumed earthquakes near the lower right half of the lune near the negative (closing) crack and negative (implosive) isotropic positions. As a result, the misidentification rate is much lower for earthquakes versus collapses than versus explosions for the same screening angle. The misidentification rates of the two populations cross at a screening angle of around 60° where the rate is ∼3 per cent. In comparison, for explosion screening, we find a misidentification rate of 5 per cent at a screening angle of 40°. On the assumption that Geysers events are fundamentally different than other tectonic sources, we have also tested the event screening by removing events contained in the boundary plotted on Fig. 4(c). The results are shown by the dashed green lines in Fig. 5. By removing these crack opening events, the optimal screening angle increases to 50°, where the misidentification rate is about 2 per cent. The misidentification rates for collapses are essentially unchanged.

In our analysis, we make the statistical assumption made in other identification methods (e.g. m_b:M_S, high-frequency P/S ratios) that the dataset is a random sample from the population and that error in the individual MT estimates is smaller than the sample variance. Although the MT inversion provides formal uncertainties, these are generally considered to be underestimates of the true uncertainty because they do not take into account inaccuracies in the velocity model, centroid location or other factors. Although several approaches have been proposed to address these shortcomings (e.g. Mustać et al., 2018), for purposes of event identification, however, the uncertainties do not factor in, since event ID utilizes the maximum likelihood MT solution. MT estimation error can be quite large when there are a small number of noisy station measurements and the velocity model used to estimate the Green's functions for the inversion are uncertain. In these cases the measurement error may be greater than spread in the distribution and therefore should not be used in any screening or identification analysis.

CONCLUSIONS

MT solutions are a common way of describing a seismic source. Relative to first motion focal mechanisms, MTs provide more source information for events, including seismic moment (and derived moment magnitude) and non-double couple components. Similarly, relative to constrained MTs, full MT solutions provide additional valuable information on event source type and, if an explosion, on yield (Pasyanos 2022). Some recent studies have suggested that, for earthquakes, non-double-couple components are artefacts of the inversion (e.g. Rösler & Stein 2022). The significance of non-double components appears to be true for non-earthquake sources like collapses and explosions, but may also be for some earthquakes as well, such as those observed in the Geysers' geothermal region, as was found in previous studies (e.g. Guilhem et al. 2014).

This study demonstrated that full MT solutions can reliably identifying non-earthquake seismic sources. We have tested this on a large dataset of solutions in the western United States where we identify known collapses and explosions, as well as some unusual earthquakes concentrated in the Geysers and central California coast. Populations for collapses and explosions are well separated on the hypersphere, allowing for effective discrimination with high (97–98 per cent) accuracy rates. Rates of identifying explosion and collapse sources can be further improved by increasing the screening angle, but at the cost of higher misidentification of earthquakes.

Although we would expect this method to be transportable outside of the western United States, in order to rigorously test this, we need to apply it to other regions. At this time, however, there are limited regions where full MT solutions are routinely calculated. As this type of analysis becomes more commonplace, we can apply the method to other populations of collapses.

Acknowledgement

We thank Doug Dreger for a helpful discussion about the San Simeon earthquake, and Bill Walter for a critical read of the manuscript. We thank reviewers Felix Rodríguez-Cardozo, Brent Delbridge and editor Carl Tape. This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344, and is document LLNL-JRNL-843094. This Ground-based Nuclear Detonation Detection (GNDD) research was funded by the U.S. National Nuclear Security Administration, Office of Defense Nuclear Nonproliferation Research and Development (NNSA DNN R&D).

Author contribution statement: MEP led the study, performed the event screening and wrote the first draft of the paper. GAI calculated all of the MT solutions for collapses, and contributed to the paper. SRF ran the statistical analysis on the collapses and contributed to the paper.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

DATA AVAILABILITY

The MTINV code is available at SourceForge (https://sourceforge.net/projects/mtinv/). Data products for this study were accessed through The Human-Induced Earthquake Database (HiQuake, https://inducedearthquakes.org). [Last accessed 16 Nov. 2022]. Waveform fits for all of the events along with network and station information are provided at doi:10.5281/zenodo.7435275.

MT solutions were performed using data available at the IRIS Data Management Center (https://www.iris.edu), which is part of the Earthscope Consortium (http://www.earthscope.org). All other data used in this paper came from published sources listed in the references.

References