Source Separation and Medium Change of Contained Chemical Explosions from Coda Wave Interferometry

The Source Physics Experiment (SPE) is a multiphase experiment to better understand explosion physics and, thereby, improve explosion monitoring (Snelson et al., 2013). The first part, Phase I, examined the explosion source in a hard‐rock medium and took place in a granite outcrop of the Climax Stock section (Townsend et al., 2011) of the Nevada National Security Site (Fig. 1, inset). Events from this series of chemical explosions in saturated granite are usually referred to simply as the SPE events (i.e., SPE‐1, SPE‐2, etc.).

All shots were recorded on a local network of short‐period seismometers (Fig. 1), and we employ these recordings in the analysis. Figure 2 shows recordings of the SPE Phase I chemical explosions from this network. The great difference in the waveforms is primarily due to differences in path heterogeneity, which is difficult to model. However, differential measurements of sources recorded at the same site do not suffer from these complications and can get at near‐source properties more directly (Poupinet et al., 1996). Poupinet et al. (2008) gives a good review of differential measurement techniques. One of the methods—coda wave interferometry (CWI)—has shown promise in finding very small changes in source and medium properties. Snieder and Hagerty (2004) measured changes in volcanic tremor location, and Grêt et al. (2006) measured medium velocity perturbation due to an imposed stress change using CWI. We will use CWI to infer small changes in the near‐source environment of the SPE explosions due to damage and test the ability to retrieve known source location offsets.

It is important to quantify the damage component of the seismic wavefield to better understand the effects of source, path, and site. Johnson and Sammis (2001) found that the damage radius is up to 10 times the post‐shot cavity radius and is a separate source of seismic radiation that must be accounted for. Patton and Taylor (2011) attempt to quantify the contribution of late‐time damage to seismic moment and find that it is significant for certain emplacement conditions. Modeling of the 1993 “tuna‐can” shaped Non‐Proliferation Experiment explosion by Stevens and O’Brien (2012) shows a somewhat spherical damage region, and Martin et al. (2011) report a velocity change of 20% in the near‐source region due to damage in granite for an experiment in Vermont that could affect shear‐wave generation (Stroujkova et al., 2012). The use of CWI with the well‐studied SPE explosions offers a unique opportunity to assess its use in other experiments and at other nuclear test sites. This observational study will focus on applying the CWI method to SPE data, and future work will focus on detailed numerical and physical modeling of how such damage can occur.

The first chemical explosion, SPE‐1, intended as a calibration shot, consisted of 87.9 kg trinitrotoluene (TNT)‐equivalent (1 kg TNT‐equivalent is 4184 kJ) of sensitized heavy ammonium nitrate fuel oil (SHANFO) loaded into a right circular cylinder canister with a diameter and height of, approximately, 0.5 m and lowered to a source centroid depth of 55.1 m. The next two shots, SPE‐2 and SPE‐3, were approximately the same size, and consisted of 997 and 905 kg TNT‐equivalent of SHANFO, respectively, loaded into right circular cylinder canisters with a diameter of 0.8 m, and height of 3 and 2.7 m, respectively, and lowered to source centroid depths of 45.7 and 45.8 m, respectively. These shots were intended to examine effects of rock damage on explosion signals.

A variety of instruments were deployed in five lines extending radially away from the shot point. In this study, we focus on the closest surface instruments, which are short‐period (4.5 Hz corner) GS11D vertical‐component geophones. The geophones were deployed at 100 m intervals and installed on concrete pads set into the surface. Geophones at a distance less than 200 m from SPE‐2 and SPE‐3 exceeded the sensor measurement range (i.e., “clipped”), so we do not include them in the analysis (indicated by × in Fig. 1). In addition, recordings of SPE‐1 at L4‐06, SPE‐2 at L2‐09, L4‐07, and L5‐03, and SPE‐3 at L2‐01, L2‐03, L3‐03, L4‐07, L5‐09, and L5‐10 had problems (e.g., dropouts, sensor noise, and altered sensor locations) and are not employed in the analysis. These traces are zeroed out in Figure 2.

CWI is a powerful technique for detecting small changes in source and structure parameters from repeated explosions (Snieder et al., 2002). The improved detection capability is due to exploitation of the Earth as a multiply scattering medium so that it can be used as an interferometer. Small changes in the scattering medium can be detected as small time shifts in the coda of seismic waves that sample the medium. The theory presented here follows Robinson et al. (2007).

When the source and receiver are separated, $K$ must be calculated numerically to give $⟨τ(t)⟩$ by equation (5). This technique is especially suited for assessing changes in velocity due to damage from the SPE Phase I explosions SPE‐2 and SPE‐3, where direct effects at near‐regional distances are difficult to detect for these approximately one tonne shots. Results using this technique will be presented in a later section.

We employ equation (2) to estimate the separation $δ$ between SPE‐1 and SPE‐2. As can be seen in Figure 2, SPE‐2 has a greater source–time duration than SPE‐1, due to its larger yield, which effectively low‐passes the recordings. In order to appropriately compare the recordings, we low‐pass the SPE‐1 recordings below the source corner frequency at 40 Hz with a six‐pole causal Butterworth filter. Before measuring the cross correlation, we upsample the 500 Hz sample rate to a 1000 Hz sample rate (⁠$Δt=1ms$⁠). Alternatively, one could measure the shifts from the slope of the phase of the cross spectra. We find that the dominant period is about 40 ms, so we use $tw=200ms$⁠.

Figure 4 shows the source separation calculation at L4‐05 for SPE‐1 and SPE‐2. The $δ$ estimation becomes more stable after, approximately, 5 s, where the multiple scattering assumption is valid. The average $δ$ after this time is 9.1 m, whereas, the actual is 9.4 m (SPE‐1 centroid depth of 55.1 m − SPE‐2 centroid depth of 45.7 m). The late‐time estimations are less stable due to the increased effect of noise on the measurements. When we employ estimates from all the stations, the expected displacement is 9.2 m.

In order to make a statement on medium change, we forward model possible scenarios and compare their predictions with observations. We model the perturbed region as a sphere, which is expected for such a deeply buried shot, with slowness perturbation that decreases to zero at the boundary. The diffusion constant $D=65×106m2/s$ used to calculate the sensitivity kernel is from a study by Wesley (1965) at the Nevada Test Site, as reported by Dainty and Toksöz (1981). Figure 5a shows a cross section of the scattering kernel $K$ for the SPE‐2 to L4‐05 source–receiver geometry where the source depth is 45.7 m and the epicentral distance is 500 m. Figure 5c shows the three localized velocity perturbation scenarios proposed as models for spherical damage. Each scenario has a slowness perturbation of 25% at the center and decreases to zero at various spherical radii. These models are then used to make predictions of $⟨τ(t)⟩$ at 500 m for comparison with observations at L4‐05 in Figure 6. The predicted shift times for a damaged sphere where the velocity perturbation $δv/v$ decreases from a maximum of 25% at the center to zero at a radius of 40 m (scenario 1 in Fig. 5c) are almost an order of magnitude greater than observed. However, the predicted shift times for a damaged sphere with radius less than 10 m (scenario 3 in Fig. 5c) are less than a few record samples and are difficult to definitively observe; therefore, this scenario cannot be rejected and serves as an upper bound. Alternative scenarios that are still bounded by the observation are balanced variations of this scenario. Balanced in the sense that a stronger, smaller perturbation or a weaker, larger one will predict shift times that remain bounded by the observations. For example, a maximum perturbation of 40% with a radius of 5 m or a maximum perturbation of 10% with a radius of 25 m are also scenarios that cannot be rejected. The feasibility of these scenarios will be discussed in the next section.

The constraint from the CWI analysis agrees with the ground truth from geological assessment. Townsend et al. (2012) report damage effects out to a radius of almost 10 m, where there are distributed microfractures out to 3 m and discrete fractures from 3 to 5 m. In addition, Schwering et al. (2019) perform a cross‐borehole change detection analysis and invert differential arrival times collected from a pair of 40 m distant boreholes for compressional velocity changes in the medium between the boreholes. The inverted image shows an approximately 25% velocity perturbation centered near the collocated SPE‐2 and SPE‐3 site that tapers out to zero over a 40 m wide by 20 m tall ellipsoidal region. Subsequent work by Hoots et al. (2020) shows that such an image could be produced by a spherical change to the medium that is not constrained by the imaging geometry. The simple spherical model presented here suggests that such a low‐velocity region would be detected in the shift times of the SPE Phase I network.

Damage (or a crushed zone) around an explosion‐created cavity in a high‐strength material like granite extends 2–5 cavity radii around the source (Rodionov et al., 1971; Office of Technology Assessment, 1989) with a radial fracture zone out to about 4–15 cavity radii (Adushkin and Spivak, 2015). Rougier and Patton (2015) show that an assumed cavity with a radius of, approximately, 1 m is consistent with the near‐field observations of reduced displacement potential that are very similar for both SPE‐2 and SPE‐3. If we take the upper limits of damage radii, this equates to about a 5 m radius of damage (crushed zone) and a 15 m radius of radial fractures. Rougier and Patton (2015) also estimate a “bulking” region of reduced density within about 6.5 m radius that fits their model for reduced displacement potential. These scenarios are consistent with the analysis presented here. In addition, Stroujkova et al. (2015) found similar scaled damage radii in their study of chemical explosions in granite.

There are several reasons why the elastic scattering theory adopted here may not be applicable to the realistic deformation induced by large chemical explosions. The damage region does not behave linearly, so scatterers within this region may not contribute in accordance with the averaging theory used to derive the shift time dependence on source separation and medium change. A possible explanation for the small shift times between SPE‐2 and SPE‐3 is that, although, the nonlinear region changed quite a bit, the scattering medium outside this region and in the elastic regime did not change at all. However, Robinson et al. (2011) found success determining earthquake source separations of tens to hundreds of meters for magnitude 1–3 events from stations up to 26 km away, which certainly would have produced a nonlinear zone around the fault‐slip region. And the success of applying the theory to find a fairly accurate assessment of source separation between SPE‐1 and SPE‐2 would argue that the theory may be applicable, such that the scattering medium is perturbed enough to prove the scattering theory valid.

Finally, the reader is referred to the outstanding work by Robinson et al. (2011, 2013). Not only are their aforementioned applied results relevant to this work, but their theoretical and numerical investigations of CWI also allow for determination of the criteria (scattering regime, signal to noise, and sensor bandwidth) for extending CWI to a regional context with source separations of hundreds of meters for magnitude 3–4 events from stations a few hundred meters away. Such a scenario exists for many nuclear test sites, and the use of CWI with double‐difference and phase arrival techniques could improve regional location for nuclear explosion monitoring.

The well‐controlled experimental conditions provided by the SPE allow us to test the capability of CWI to measure source medium damage and source separation. The damage caused by the, approximately, one tonne TNT‐equivalent SPE‐2 did not significantly perturb the direct wave arrival time or coda waves such that a travel‐time perturbation was observed using CWI. We could, however, place an upper bound on the velocity perturbation caused by any damage that is constrained to a spheroidal region smaller than 10 m radius with velocity perturbation less than 25%. This constraint is consistent with observations obtained from a drillback to the source. CWI analysis between the, approximately, 100 kg TNT‐equivalent SPE‐1 and, approximately, 1 tonne TNT‐equivalent SPE‐2 estimated the source separation between those events to be, approximately, 6–18 m, with an expected separation near 9.2 m. The actual source separation is 9.4 m. This analysis shows that CWI is a powerful technique to probe small changes in the source and medium. Future work should focus on the combined effects of damage and source separation as would be observed in a more complicated scenario.

Surface seismic data from the Source Physics Experiment Phase I are provided to the Incorporated Research Institutions for Seismology (IRIS) using network identifier SN (doi: 10.7914/SN/SN). Full data and reports are also part of Assembled Datasets at IRIS (http://ds.iris.edu/SeismiQuery/assembled.phtml, last accessed May 2018 and search dataset named “Source Physics Experiment”). Some plots were made using the Generic Mapping Tools v. 4.2.2 (www.soest.hawaii.edu/gmt, last accessed September 2002). Seismic analysis was done with Seismic Analysis Code v. 7.10.5 (https://ds.iris.edu/ds/nodes/dmc/software/downloads/SAC/102-0/, last accessed October 2011).

The authors acknowledge there are no conflicts of interest recorded.

The authors are grateful for helpful comments by Associate Editor Steven J. Gibbons and two anonymous reviewers that greatly improved the article. The Source Physics Experiment (SPE) would not have been possible without the support of many people from several organizations. The authors wish to express their gratitude to the National Nuclear Security Administration, Defense Nuclear Nonproliferation Research and Development, and the SPE working group—a multiinstitutional and interdisciplinary group of scientists and engineers. The authors thank the Survey Department of Mission Support and Test Services, LLC for the ground‐surface coordinates. This research was performed under the auspices of the U.S. Department of Energy by the Lawrence Livermore National Laboratory under Contract Number DE‐AC52‐07NA27344.