Local Explosion Detection and Infrasound Localization by Reverse Time Migration Using 3-D Finite-Difference Wave Propagation

Infrasound data are routinely used to detect and locate volcanic and other explosions, using both arrays and single sensor networks. However, at local distances (<15 km) topography often complicates acoustic propagation, resulting in inaccurate acoustic travel times leading to biased source locations when assuming straight-line propagation. Here we present a new method, termed Reverse Time Migration-Finite-Difference Time Domain (RTM-FDTD), that integrates numerical modeling into the standard RTM back-projection process. Travel time information is computed across the entire potential source grid via FDTD modeling to incorporate the effects of topography. The waveforms are then back-projected and stacked at each grid point, with the stack maximum corresponding to the likely source. We apply our method to three volcanoes with different network configurations, source-receiver distances, and topography. At Yasur Volcano, Vanuatu, RTM-FDTD locates explosions within ∼20 m of the source and differentiates between multiple vents. RTM-FDTD produces a more accurate location for the two Yasur subcraters than standard RTM and doubles the number of detected events. At Sakurajima Volcano, Japan, RTM-FDTD locates the source within 50 m of the active vent despite notable topographic blocking. The RTM-FDTD location is similar to that from the Time Reversal Mirror method, but is more computationally efficient. Lastly, at Shishaldin Volcano, Alaska, RTM and RTM-FDTD both produce realistic source locations (<50 m) for ground-coupled airwaves recorded on a four-station seismic network. We show that RTM is an effective method to detect and locate infrasonic sources across a variety of scenarios, and by integrating numerical modeling, RTM-FDTD produces more accurate source locations and increases the detection capability.


INTRODUCTION
Numerous methods have been applied to detect, locate, and characterize volcanic and anthropogenic explosions using lowfrequency acoustic waves (infrasound). Long-range (defined here as >15 km distance) infrasound localization techniques typically use intersecting back-azimuths from multiple infrasound arrays (termed Direction of Arrival (DOA) or triangulation) (e.g. Brachet et al., 2010;Landès et al., 2012;Matoza et al., 2017), with some methods incorporating Bayesian likelihood functions for back-azimuth and arrival times (e.g. Modrak et al., 2010). Networks of distributed single infrasound or seismic sensors are also used to locate infrasonic sources using back-projection, gridsearch techniques such as Reverse Time Migration (RTM) (Walker et al., 2010;Sanderson et al., 2020). Large networks can also be divided into smaller arrays to increase waveform coherency, and then processed to detect and locate events (de Groot-Hedlin and Hedlin, 2015). At closer ranges (local, defined here as <15 km), DOA methods using arrays have been applied for years (Szuberla et al., 2009), but networks of single sensors surrounding a source are becoming more common. In this scenario the propagating infrasonic wave may not be planar across the network, as it typically is for infrasound arrays, therefore Time Difference of Arrival (TDOA) information must instead be derived using waveform similarity between stations. Numerous TDOA methods employing hyperbolic travel time curves (e.g. Szuberla et al., 2009;Rowell et al., 2014), RTM via semblance (e.g. Ripepe and Marchetti, 2002;Cannata et al., 2009;Jones and Johnson, 2011), and station-pair double difference (Fee et al., 2016) have been successful in locating infrasonic events at these ranges. However, for regions with topography or other obstructions that block the direct path of sound, such as is often the case at volcanoes and in mountainous regions, recent studies (e.g., McKee et al., 2014) have found that these traditional techniques fail to accurately locate the source. In these situations numerical simulations can model the complex wave propagation and be integrated into a location framework to provide an accurate infrasound source location Lees, 2014, 2015). Related work has similarly shown that numerical simulations are needed to accurately determine other source parameters (e.g. mass flux, explosive yield) under complex propagation conditions (e.g. Lacanna and Ripepe, 2020). A heterogeneous and moving atmosphere (winds) may also affect the location accuracy, although these factors appear to produce minor changes in cases explored thus far (Kim and Lees, 2014;McKee et al., 2014). Kim and Lees (2014) introduced a 3-D Finite-Difference Time Domain Algorithm (termed infraFDTD3D) that can accurately model the infrasonic wavefield in regions of complex topography in a windless and homogeneous atmosphere, and can be run on a desktop computer utilizing GPU parallelization. The authors applied the Time Reversal Mirror (TRM) algorithm at Sakurajima (Kim and Lees, 2014) and Santiaguito  volcanoes and derived high-resolution, accurate source locations of explosive events despite topographic obstructions and complex acoustic propagation conditions. Their TRM method 1) reverses the recorded waveforms in time, 2) applies a distance-corrected amplitude, 3) backpropagates the reversed waveforms over the search grid using infraFDTD3D, and 4) finds the maximum convergence of the back-propagated waves at each time step. The likely source location corresponds to the location of peak convergence. Although it provides high-resolution source localization and can account for complex wave propagation, TRM is computationally expensive and a new waveform simulation needs to be conducted for each event.
Here we present a new technique, termed Reverse Time Migration-Finite-Difference Time Domain (RTM-FDTD), which integrates numerical modeling via FDTD into the traditional RTM back-projection process. Realistic travel time information is computed from each station across the entire potential source grid from 3-D Green's functions derived via infraFDTD3D. We then back-project the waveforms using this travel time information and stack them at each grid point. The maximum stack location corresponds to the likely source. This single travel time simulation can then be used for additional events, since we assume a static atmosphere where the travel time is not changing. The stack function can also be used as a "detector" for coherent infrasound events recorded across the network. We apply our method to detect and locate explosions in three different settings: Yasur Volcano, Vanuatu; Sakurajima Volcano, Japan; and Shishaldin Volcano, Alaska. These scenarios provide variable network geometries, sourcereceiver distances, complex wave propagation, and multiple closely spaced source locations to evaluate the performance of RTM-FDTD. The Shishaldin example uses ground-coupled airwaves (GCAs) recorded on nearby seismic stations. We restrict our ranges to less than 15 km distance. We also compare the RTM-FDTD results with those from more traditional methods using direct straight-line propagation (RTM), as well as the TRM method at Sakurajima. We find that both RTM and RTM-FDTD are able to locate infrasound sources at local distances to within tens of meters of the actual source when propagation conditions are simple. However, when complex topography blocks propagation paths, the poorly resolved location may be improved via numerical modeling. Accurate travel times not only improve source location estimates but also increase the number of detected events by properly aligning the waveforms. Because eruption source processes may occur within complex and rapidly evolving vent systems, accurate locations and well-resolved eruption catalogs can improve characterization of volcanic processes that may be used to inform volcano monitoring for hazards. Our open-source Python RTM code is available on GitHub at https://github.com/uafgeotools/rtm and can be applied to a range of scales and network geometries. The travel times derived via FDTD must be provided separately.

DATA
We use data from three volcanoes with different network geometries, number of stations, data types, topography, and scales.

Yasur, Vanuatu
Yasur Volcano, Vanuatu is a persistently active volcano and prodigious infrasound source. We use data collected over five days in July-August 2016 from a six-station infrasound network surrounding the active vents ( Figure 1A). At the time of recording, Yasur had two main explosive vents (termed Vent A and C) within the main 400 m wide crater, with small explosions occurring every ∼1-5 min and large explosions more sporadically from both vents. Vents A and C are located in subcraters and are separated by ∼145 m. Based on visual evidence, we assume each vent location lies at the bottom of each subcrater and ∼150 m below the main crater rim. A less active subvent (Vent B) is located near Vent A. We note that the vent locations are estimated from digital elevation models (DEMs) and therefore contain some inherent uncertainty. A thick plume and other observational challenges restricted our ability to accurately measure the vent locations during the experiment (see Iezzi et al., (2019) for a discussion of the challenges with DEM creation and vent locations), and there was some qualitative evidence that the vents migrated up to tens of meters during the experiment. In this manuscript we analyze 10 s of data from a single event from both Vent A and C, along with 12 h of typical activity.
Source-receiver distances range from ∼250-800 m for our network, with four stations along the edge of the main crater (YIF2-YIF5) and two along the flanks (YIF1, YIF6) ( Figure 1A). Each infrasound station consisted of a single Chaparral Physics Model 60-UHP infrasound sensor (flat response 0.03-245 Hz) connected to a DiGOS DATA-CUBE 3 digitizer. An infrasound sensor tethered to an aerostat was also deployed above the crater rim (Jolly et al., 2017;Iezzi et al., 2019). However, we do not use this sensor here due to larger uncertainties in the sensor location. Azimuthal coverage of the network is good with the largest azimuthal gap being 80°. The four crater stations have line-of-sight views to at least one of the active vents while the two flank stations have topographic obstructions along the propagation paths. We use a highresolution (∼2 m) DEM of Yasur topography constructed from a UAV survey and satellite images (Iezzi et al., 2019) for the infraFDTD3D modeling and to constrain the vent locations. See Marchetti et al., (2013), Jolly et al., (2017), Iezzi et al., (2019), and Fitzgerald et al., (2020) for more details on Yasur explosions and the deployment during this time period.

Sakurajima, Japan
Sakurajima Volcano, Japan typically produces ash-and ballistic-rich vulcanian explosions from either the Showa or Minamidake craters (Yokoo et al., 2013). These explosions generate high-amplitude infrasound signals that have been well-characterized by past studies (e.g. Yokoo et al., 2013;. We use data from a five-station infrasound network deployed around Sakurajima between 18 and July 27, 2013 ( Figure 1B) as part of a community infrasound experiment . Explosions during the deployment emanated from the bottom of the approximately 100 m deep by 300 m wide Showa crater. We analyze a single high signal-to-noise (SNR) explosion from Sakurajima and compare it to the TRM location from Kim and Lees (2014), who found that all 10 explosions analyzed during this period had similar source locations.
Source-receiver ranges for the Sakurajima dataset are between 2.4 and 6.3 km and the largest azimuthal gap is ∼95°. Sound propagation paths from the bottom of Showa crater to two stations (HAR and SVO) have notable topographic obstructions . Three stations (SVO, ARI, KOM) consist of National Center for Physical Acoustics infrasound sensors while the other two (HAR, KUR) have Hyperion IFS-5201 infrasound sensors. The two sensor types have similar instrumentation with built-in digitizers sampling at 500 Hz, flat responses between 0.02 and 250 Hz, and clip at ∼1,000 Pa. The DEM we use has a 5 m resolution around the Showa crater region and 10 m resolution elsewhere. See Fee et al., (2014), McKee et al., (2014), and Kim and Lees (2014), among others, for more information on this dataset and related activity and analysis.

Shishaldin, Alaska
Shishaldin Volcano is a large conical stratovolcano on Unimak Island, Alaska. The summit elevation is 2,857 m and the edifice extends nearly 16 km in diameter. It is one of the most active volcanoes in the Aleutian Arc with over 24 eruptions since 1775 (Miller et al., 1998). The most recent eruption began in July 2019 and produced a variety of eruptive activity, including numerous Strombolian explosions, debris flows, and weak subPlinian eruptions. Eruptions typically emanate from a summit crater that varies in size but is roughly 250 m wide and 100 m deep in our DEM.
The Alaska Volcano Observatory (AVO) operates a threestation seismo-acoustic network on Shishaldin, with two additional seismic stations on nearby Isanotski Volcano ( Figure 1C). We focus on a short time period of the most recent eruption. In October 2019 Strombolian explosions produced ground-coupled airwaves (GCAs) visible on the four nearest seismic stations at distances of ∼6.4-14.5 km from the summit crater. Although these explosions were also recorded on the two colocated local infrasound stations on Shishaldin (SSLN and SSBA), we focus on 60 s of data from October 22 where GCAs were recorded on the seismic stations SSLN, SSBA, SSLS, and ISNN. These four stations permit an evaluation of the effectiveness of acoustic source localization of GCAs using RTM over a relatively large (16 × 16 km) region. The largest azimuthal gap between stations is 148°. For the topography, we use a 5 m resolution DEM derived via Interferometric Synthetic Aperture Radar (InSAR) (Data.gov, 2017).

METHODS
In this section we discuss the RTM and RTM-FDTD methods employed in this manuscript, as well as the parameters used for each scenario. RTM-FDTD follows the same procedure as RTM, but involves travel times derived from FDTD numerical modeling as opposed to travel times derived from assuming a straight-line, unobstructed path. The back-projection can be evaluated through stacking of the waveforms or evaluation of semblance. Processing parameters used in this manuscript are listed in Table 1.

RTM
Reverse Time Migration (RTM) is a common method to locate seismic and acoustic sources (Ishii et al., 2007;Walker et al., 2010;Sanderson et al., 2020) and consists of back-projecting, or migrating, recorded waveforms over a grid of trial source locations assuming an apparent propagation velocity. The preprocessed waveforms (e.g., filtered, enveloped) back-projected to each potential source location are then stacked, with waveform energy adding constructively at the likely source location assuming accurate propagation times. Our implementation of RTM follows Walker et al., (2010) and Sanderson et al., (2020). However, those studies focused on relatively large source-receiver distances and grid sizes (hundreds to thousands of kilometers), while our study analyzes acoustic signals recorded over distances of a few to tens of kilometers.

Grid Construction
Prior to data processing, we construct a grid of trial source locations. The grid encompasses the likely source location and is centered on the volcanic vent. All scenarios have notable topography that may affect the acoustic wave propagation, so we construct a grid from the aforementioned high-resolution DEMs. Grid sizes for the Yasur, Sakurajima, and Shishaldin scenarios are 1.4 × 1.4 km, 8.0 × 8.0 km, and 16.0 × 16.0 km, respectively ( Figure 1). Accurate resolution of a potential source at each grid node requires that the grid spacing is smaller than the distance an acoustic wave travels for each waveform sample (Walker et al., 2010). We downsample the data to 80, 40, and 40 Hz for Yasur, Sakurajima, and Shishaldin, respectively, to speed-up the calculations. For an acoustic velocity of 340 m/s, 40 and 80 Hz sampling rates correspond to a maximum grid resolution of 8.5 and 4.25 m, respectively. For Yasur, Sakurajima, and Shishaldin we set the grid node spacing to 4, 6, and 20 m, respectively. The extent of the Shishaldin grid requires a larger grid size to ensure reasonable computational speeds. The DEM is resampled over the grid in UTM coordinates using cubic spline interpolation to produce an (x, y, z) grid of potential source locations. Sources are assumed "clamped" to the surface for this analysis. Figure 1 shows the grids (black dashed box) and network layouts for each scenario.

Data Processing
Appropriate data processing prior to stacking is key for successful application of RTM (e.g. Sanderson et al., 2020), especially when source-receiver distances are large and/or waveforms may not be coherent. First, data are read in for the analysis window length, plus additional data to pad for travel time removal across the entire grid. We read in enough data for the furthest station from the grid center to record a signal in the window of interest from the furthest grid point. We then band-pass filter each waveform, u, to isolate the signal of interest. Filter bands ( Table 1) were chosen to maximize the SNR for each scenario. Next we downsample to the desired sample rate to speed up the grid search process. Infrasound signal coherency decreases as a function of range and frequency (Green, 2015), and GCAs may not be coherent across a network (Fee et al., 2016). Therefore, standard multi-channel cross-correlation techniques, such as semblance, may not be effective. For this reason we construct the waveform envelope via the analytical signal of the Hilbert transform. For longer source-receiver distances and at higher frequencies, smoothing of the signal may also be necessary to increase stacking power (Walker et al., 2010). For the Shishaldin scenario we apply a 0.5 s Hann smoothing window to the waveform envelope. Lastly, for short-duration (<60 s) analysis windows, we normalize each waveform to reduce the effect of amplitude loss as a function of distance. Longer data segments with numerous events in the timeseries and varying SNR may require an automatic gain control (AGC) to regularize the envelopes. Here we analyze a continuous 12-h time section for Yasur. See Walker et al., (2010) and Sanderson et al., (2020) for more details and discussion on RTM waveform processing, and Table 1 for the parameters used here.

Back-Projection and Stacking
The processed waveforms are then back-projected to every grid node (potential source) to produce a set of computed station to grid node travel times. This RTM travel time, Δt ijl , is computed assuming straight-line propagation from the 3-D grid node surface (x, y, z) to the station (n) at a fixed apparent velocity, or celerity (range divided by time). The range is the Euclidean (slant) distance between the two points. This is a common assumption in local explosion detection and localization (Szuberla et al., 2009). We choose a fixed celerity for each of our scenarios based on average temperatures in the region: 343.5 m/s (20°C) for Yasur (Iezzi et al., 2019), 349.3 m/s (30°C) for Sakurajima (McKee et al., 2014), and 334.0 m/s (4°C estimate) for Shishaldin. The waveforms are then corrected for travel time and are stacked, or summed together, and then divided by the number of stations to form a stack function A: where h is the processed waveform (envelope), N is the number of stations, and i, j, k, l are the indices associated with x and y coordinates, time t, and station n, respectively. The z coordinate is not specified here since all grid nodes are clamped to the DEM surface.
The maximum of the stack function represents the most likely source location and the maximum stack value for our scenarios is normalized to 1. We perform this normalization to provide a more direct comparison between scenarios with different numbers of stations and the semblance calculation. The stack function can also be evaluated with more complex detector functions based on relative SNR within a longer timeseries (e.g. Walker et al., 2010). Since we primarily focus on single events in this study and comparison of algorithms, we just analyze the stack function. The results can be visualized by plotting an intuitive map-view "slice" of the stack function, A ijk , at each time step.
Stacking via semblance, a measure of normalized multichannel coherence, has also been effective at source localization, particularly at local distances. We apply semblance during the stacking process for comparison and evaluation of our algorithms and scenarios. Semblance is calculated by first computing the beam, u, for the entire network at each grid point for a time window of interest: where u is the processed waveform. The beam is defined as the mean of time-shifted waveforms in an array or network. Semblance can then be calculated by taking the average beam power divided by the average station power (Neidell and Taner, 1971): where M is the number of samples in the window. Semblance is then defined over the interval [0, 1]. We acknowledge that other definitions of semblance exist (e.g. Ripepe and Marchetti, 2002;Johnson et al., 2011), although they all rely on the similar principle of multi-station coherence. We use a 5 s time window with 50% overlap for semblance processing in all of our scenarios. We define an "event" or detection to be the maximum of the stack function when A or S exceeds 0.6 for the time window of interest.

RTM-FDTD
The RTM-FDTD method follows the same process outlined above, but implements travel times computed by solving the acoustic wave equations using the Finite-Difference Time Domain (FDTD) algorithm of Kim and Lees (2014). We compute the travel times Δt ijl by inserting an impulsive source and then propagating it out over the 3-D topographic grid ( Figure 1). Rather than run a simulation from each one of the (tens of) thousands of grid nodes to every station, which would be very computationally expensive, we instead run N simulations with a source at each station that propagates out across the grid. Each grid node is then a synthetic receiver. We then use the principle of reciprocity (Morse et al., 1970) to assume the travel time from grid node (receiver) to the station (source) is the same as the station to grid node. This reduces the computational burden considerably. The source-time function is a Blackman-Harris (Gaussian) function and the travel time is derived from the time difference between the peak in the source-time function and peak in the synthetic receiver waveform. The source of the acoustic wavefield is represented by the mass flow rate, which Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 620813 5 corresponds to pressure changes in air following the finitedifference formulation of Ostashev et al., (2005). We use a smoothed step function as the time history of mass flow rate, which generates pressure waveforms corresponding to a Blackman-Harris function. Then, the travel time is derived from the time difference between the peaks in the pressure waveforms at the source and synthetic receivers. The interface between the ground and atmosphere is assumed to be a rigid boundary, perfectly reflecting incident acoustic waves without absorption. The corner frequency of the simulation source-time function corresponds to at least 15 grid nodes per smallest wavelength to ensure numerical stability ( Table 1). Each FDTD simulation grid encompasses the associated RTM grid (Section 3.1.1) and celerities are the same as in the RTM backprojection (Section 3.1.3). Note we do not consider or correct for changes in amplitude due to geometric spreading or diffraction Lacanna and Ripepe, 2020), as this adds additional complexity and assumptions not needed for accurate source location in our scenarios.
Our FDTD computation assumes a 1-D, static atmosphere with a fixed velocity. This assumption has been shown to be valid in similar previous studies (e.g. Lacanna and Ripepe, 2013;Kim and Lees, 2014;Fee et al., 2017;Iezzi et al., 2019), where topography has a dominant effect on the wave travel time and structure at local distances. Kim and Lees (2014) show how changes in the celerity due to wind or temperature variability did not have an appreciable change in location accuracy at Sakurajima using the same dataset we analyze in Section 4.2. We note that the scenarios considered here all had relatively low winds (<10 m/s), and that the 1-D static atmosphere may not be valid in other scenarios. We note that future work adding a wind advection term to the FDTD calculation may violate the reciprocity assumption.
The travel time difference between the FDTD and straight-line calculations can be considerable when topographic barriers exist. Figure 2A shows a snapshot of the FDTD simulation originating from Yasur station YIF6 at 3.0 s. Red indicates a positive pressure of the acoustic wave along the ground surface while blue indicates a negative rarefaction. The spherically radiating wave is clearly distorted around the crater area. Figure 2B shows the difference between the FDTD computed and straight-line propagation travel times at each grid node with the source at station YIF6. The straight-line propagation underestimates the travel time by up to 0.42 s at the bottom of the crater, despite the distance only being a few hundred meters from the source, an error of ∼14%. The travel times on the far side of the crater (southeast) are also notably underestimated by the straight-line assumption. Similar travel time underestimation exists for other Yasur and Sakurajima stations. Note the FDTD travel times are always equal to or greater than the straight-line calculation, as the interaction of acoustic waves with topography creates a time delay.

RESULTS
In this section we compare RTM and RTM-FDTD event detection and source location results for our three scenarios. We compare detections and source localization from RTM-FDTD using waveform stacking, which we term RTM-FDTD (sum), and semblance, termed RTM-FDTD (semblance), vs. standard RTM assuming straight-line propagation and stacking of waveforms, termed RTM (sum).

Yasur Volcano
We first apply RTM and RTM-FDTD to one explosion originating from each vent at Yasur Volcano. We find that although all algorithms are able to locate and distinguish events between the two vents, RTM-FDTD produces the most accurate location, irrespective of using envelope stacking or windowed-semblance. Figure 3A shows the filtered (gray) and processed (black) waveforms for 10 s surrounding a Vent A explosion recorded on the six station network ( Figure 1A). The explosion has a clear impulsive onset and coda lasting Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 620813 6 ∼7 s. The SNR is high at most stations, although notably lower for station YIF6 on the crater flank. Figure 3B) shows a slice of the RTM-FDTD (sum) stack near the vent (700 × 700 m region) for the time of the peak in the RTM-FDTD stack function (red line, Figure 3C). The peak RTM-FDTD (sum) stack value of 0.96 is close to theoretical maximum of 1, demonstrating that the maximum of all processed waveforms stack coherently. This peak also provides the time of the explosion signal maximum amplitude back-projected to the source, which is likely very close to the explosion onset for this explosion. The RTM-FDTD (sum) source location ( Figure 3B, red star) is only 24 m to the southwest from the assumed vent location. RTM-FDTD (semblance) run in 5 s overlapping windows provides a similar location: 25 m to the southwest ( Figure 3B, gray square) and peak value of 0.81. The RTM (sum) location assuming straight-line propagation is further away, ∼53 m to the southeast ( Figure 3B, blue diamond), but still reasonably close to the assumed vent. The maximum RTM (sum) stack is shown as a blue line in Figure 3C, with a peak amplitude of 0.91.
The stack functions in Figure 3C for RTM-FDTD (sum) and RTM (sum) have similar shapes and amplitudes. However, the RTM-FDTD (sum) stack in red has a higher peak at ∼02:17:48.5 UTC and secondary peak about a second later, indicating improved alignment from more accurate travel times. This second peak, at ∼02:17:49.5 UTC, is much more clear in the RTM-FDTD maximum stack amplitude in Figure 3C than the envelope stacks in the bottom panels of Figure 3A, suggesting that the latter may be a secondary source. The travel time reduced, processed envelopes for the maximum stack location in Figure 3B are shown in the bottom two panels of Figure 3A. The superior alignment is apparent in the red RTM-FDTD (sum) waveforms. The blue RTM stacks are also delayed with respect to the red RTM-FDTD stacks. This is related to the longer travel times computed by the FDTD modeling ( Figure 2) that have been subtracted from the stacks. The timing difference between the RTM-FDTD and RTM stacks agrees with the travel time differences of up to ∼0.4 s in Figure 2B. The noisier YIF6 envelope is also clear as an outlier in both panels. Note that Figure 3C shows the maximum stack amplitude at each time sample, while the bottom two panels of Figure 3A show the normalized envelopes for the maximum stack location in Figure 3B. RTM-FDTD also provides high-resolution source localization for a Yasur Vent C explosion. Figure 4 shows the waveforms, map of the peak stack slice, and peak stack functions for RTM-FDTD (sum) and RTM (sum). This event also has an impulsive infrasound signal but it is shorter in duration (∼3 s). The explosion has a clear impulsive onset producing a peak of ∼25 Pa on the crater rim stations, which is clearly visible in the waveform envelopes. The maximum RTM-FDTD (sum) stack amplitude is 1, i.e. the theoretical maximum. The RTM-FDTD (sum) location ( Figure 4B, red star) is 20 m to the northeast of the assumed vent location, while the RTM-FDTD (semblance) location ( Figure 4B, gray square) is only 14 m to the northeast. The maximum RTM-FDTD (semblance) value is 0.83. The RTM (sum) location is further away at 59 m to the southeast with a maximum stack amplitude of 0.94. Figure 4C both have a clear peak and similar trend, with the RTM-FDTD (sum) stack having a higher peak and a few additional increases in energy, indicating improved alignment. The overall shape of the stack function generally resembles the waveform envelopes, although some artifacts exist in the grid stacking for both this and the Vent A explosion.

The stack functions in
We also apply the RTM algorithms to multiple hours of data and find that RTM-FDTD detects more events and provides a more stable, accurate location. Figure 5 displays 12 h of filtered waveforms for all six Yasur stations on 28-29 July 2016. Note the numerous high-amplitude explosions apparent throughout the time period, particularly the first 4 h, where the crater rim stations recorded signals up to 100 Pa. Station YIF6 has a lower SNR and increased wind noise is apparent, especially after 28 July 22:00 UTC. This station is further from the vent and positioned along the crater flank ( Figure 1A). Figure 6 compares the RTM and RTM-FDTD results for this time period using windowed semblance. Figures 6A,C) show a 2-D histogram of the event locations per grid point during this period for RTM and RTM-FDTD, respectively, along with the Vent A and C locations and 95% confidence ellipses. The maximum semblance values in each overlapping 5 s window are shown in Figures 6B,D). RTM-FDTD has a consistently higher maximum semblance value and the maximum semblance values are generally higher the first 4 h and then decrease, consistent with the higher amplitude explosions and the decrease in SNR visible in Figure 5. RTM-FDTD detects over twice as many events during this period as RTM: 1,589 compared to 746. The RTM-FDTD event locations are also closely clustered near Vents A and C ( Figure 6C), with a single grid point for Vent A and C containing 480 and 402 events, respectively. The RTM events are distributed more broadly within the crater region ( Figure 6A), with the peak grid points for Vent A and C having only 93 and 25 events, respectively. The 95% confidence ellipses provide a measure of the location uncertainty for this scenario and are much smaller for RTM-FDTD. For RTM, the Vent A ellipse has a major half-axis of 122 m and minor half-axis of 55 m, while for RTM-FDTD they are about 25% of these values at 28 m and 15 m, respectively. For Vent C, the RTM major half-axis is 68 m and minor halfaxis is 55 m, while for RTM-FDTD they are about a third that of RTM at 21 m and 17 m. Both algorithms locate very few events outside the crater region.

Sakurajima Volcano
We apply the RTM and RTM-FDTD algorithms to a high SNR explosion from Sakurajima Volcano. We find that although all algorithm variations provide a reasonable source location estimate within 171 m of the assumed vent, RTM-FDTD provides the most accurate source location estimate (47 m). Figure 7A) shows the filtered and processed waveforms for an explosion from Showa crater on July 18, 2013, while Figure 7B) shows the various RTM locations associated with the maximum stack in Figure 7C), displayed over a 1,400 × 1,400 m region centered on the assumed vent location. The explosion is typical of Sakurajima events, with a relatively high-amplitude, impulsive onset clearly visible on all five stations, followed by a variable coda lasting 5-10 s ( Figure 7A). The waveform envelopes look mostly similar although there is a dual pressure increase for KUR due to crater reflections . All RTM variations locate the source to the east or southeast of the expected source location (green circle in Figure 7B) along the edge of or outside of Showa crater. RTM-FDTD produces the most accurate source location: RTM-FDTD (sum) 81 m southeast (red star) and RTM-FDTD (semblance) 47 m east (gray square) ( Figure 7B), with a maximum stack amplitude of 0.99 and semblance value of 0.95. The RTM (sum) result is further away, at 171 m to the southeast (blue diamond in Figure 7B), and lies outside Showa Crater. The RTM (sum) stack maximum is the same as RTM-FDTD (sum) at 0.99. Sakurajima RTM and RTM-FDTD locations compare favorably to the high-resolution TRM locations of Kim and Lees (2014). TRM locates this event 55 m to the southeast of the assumed vent, in a very similar location as RTM-FDTD (semblance) and RTM-FDTD (sum). Kim and Lees (2014) located 10 events in this dataset using two different celerities and found very similar locations, suggesting that the offset is not related to a changing atmosphere or source location. We note that McKee et al., (2014) used straight-line propagation assumptions and a 3-D search grid not clamped to the surface for their semblance-based analysis of this dataset, and found the average source location to be offset horizontally approximately 130 m to the northeast. Notably, their average location's 3-D slant offset from the source was 370 m, and over two hundred meters below the surface. They note poor vertical resolution due to a roughly 2-D network all located below the source. Our RTM (sum) location assuming straight-line propagation is closer to the assumed vent than McKee et al., (2014), most likely because we force our source to be located on the DEM rather than stacking over a fully 3-D grid. Waveform processing also differs slightly.
We also located nine other explosive events from Sakurajima, previously examined by Kim and Lees (2014), using the same processing parameters as for the July 18, 2013 event noted above. We find the RTM and RTM-FDTD location results to be consistent with the example shown in Figure 7. The mean locations for these nine events are offset 117 m to the southeast for RTM and 34 m to the east for RTM-FDTD, comparable to the 171 m and 47 m offsets, respectively, in Figure 7B.

Shishaldin Volcano
We apply RTM and RTM-FDTD to a series of explosions from Shishaldin Volcano, which produced GCAs on the local seismic network. Figure 8 shows the waveforms, location comparison, and maximum of the stack for 60 s of data for at least two GCAs produced by explosions. The GCAs are most clear on the three Shishaldin stations closest to the volcano (SSLN, SSBA, and SSLS), and less so on ISNN ( Figure 8A). RTM-FDTD (sum) and RTM (sum) both locate the source close to the assumed vent: 40 m east for RTM-FDTD (sum) and 45 m north for RTM (sum) ( Figure 8B), which is quite close considering the large sourcereceiver ranges. Both algorithms have a peak stack of around 0.89. Note also that there are multiple peaks in the maximum stack amplitude in Figure 8C), likely due to multiple explosions occurring during this analysis period. Two explosions are clearly visible in the Shishaldin station waveforms in Figure 8A), especially on stations SSLN, SSBA, and SSLS. Figure 8B) only displays the location for the second event with an approximate source time of 14:51:42. RTM-FDTD (semblance) does not produce a realistic detection for this time window. The maximum semblance values are low (<0.6) and the location of the maximum semblance is over 8 km from the source and is not visible in the zoomed-in view of the summit in Figure 8B). We explored other data processing parameters for use with windowed semblance, including longer smoothing windows, but all returned low semblance values and poor locations. We note the peak in the stacks is only modestly above the background, likely due to the low SNR of the signal on ISNN and the relatively large grid. However, the RTM-FDTD (sum) and RTM (sum) source locations are both quite close to the assumed vent location (<50 m). We also note a lack of notable topography between the source and receivers, which is consistent with similar results between RTM and RTM-FDTD.

DISCUSSION
Our results demonstrate that RTM and RTM-FDTD are effective at detecting and locating infrasonic sources recorded on both seismic and infrasound stations at distances <15 km. In situations where topography blocks the direct propagation from source to receiver, such as at Yasur and Sakurajima volcanoes, RTM-FDTD provides more accurate travel times-that increase the number of events detected-and a more accurate source location estimate. RTM-FDTD clearly reduces the overall uncertainty and variance in source location estimates, as shown at Yasur in Figure 6. Highfrequency GCA waveforms are not coherent between seismic stations for our Shishaldin scenario, suggesting that caution should be used applying semblance to GCAs. Waveform envelope stacking was still effective in detecting and locating the explosive source of GCAs. Proper tuning of the waveform processing parameters is key in this situation.
In addition to just basic explosion detection and localization, the location accuracy demonstrated here by RTM-FDTD can FIGURE 8 | Shishaldin explosion waveforms and location comparison. Figure layout is the same as in Figure 3, except waveforms are filtered between 5 and 15 Hz and these are seismic data. RTM (sum) and RTM-FDTD (sum) produce similar results and locate the source to within 50 m of the base of the crater, while RTM-FDTD (semblance) does not provide a realistic location due to poor coherence.
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 620813 provide insight into eruption dynamics, such as differentiation of vents within a crater, monitoring changes in vent location, distinguishing flank vs. summit eruptions, and related phenomena. Rockfalls and others mass flows can likely be similarly detected and located (e.g. Johnson and Palma, 2015). The signals from these sources would likely have more extended durations but should still be locatable assuming waveform coherence. Non-point-sources (e.g. fissures) may radiate sound from multiple locations, which in principle could be located and distinguished using RTM. The algorithms presented here can already permit improvements in volcano monitoring and eruption detection and characterization of both infrasonic events and GCAs. RTM, RTM-FDTD, and TRM are all computationally limiting at some level. The grid-search process in RTM and RTM-FDTD is slow for large grids and high sample rates. A standard laptop computer can construct a grid, process the data, and run the gridsearch in near-real-time, although the grid spacing would have to be relaxed slightly from our scenarios. Testing of larger grid sizes did not show an appreciable decrease in location accuracy for our scenarios, suggesting that the grid node size justification based on sampling frequency in Section 3.1.1 is conservative and that moderately larger grid node sizes would provide a good trade-off between computational speed and accuracy. We note that future implementations of the grid-search process could implement multi-core processing for notable speed improvements. The infraFDTD3D calculations for RTM-FDTD and TRM are also relatively computationally expensive, even when run using GPU parallelization. They typically take 1-2 days for our scenarios and computational resources (three NVIDIA Tesla K80 GPU cards). The advantage of RTM-FDTD vs. TRM lies in the fact that the travel times are only computed once for each station, compared to a separate simulation for each time-reversed waveform in TRM. We note that TRM may be preferred when the waveforms are considerably distorted during propagation, which could cause the envelope stacking or semblance to be ineffective. This may occur in regions of complex topography, where a strong reflection modifies the timing of the peak pressure signal, or at high frequencies (i.e. >10 Hz) where diffraction will be more prominent. However, this was not the case for the scenarios and frequencies explored here, even at Sakurajima where strong reflections likely occur around Showa crater. A coarser RTM grid will also substantially decrease the computation time for the grid search and would be preferred for near-real-time detection and localization, and standard RTM may be effective enough for these applications.
Waveform distortion between stations should be considered for each scenario. Our envelope stacking focuses on the peak of the stack function, which is typically the first impulse in an explosion waveform (e.g. Figure 7A). This appears to mitigate the effects of waveform distortion at both Sakurajima and Shishaldin. Care should be taken in using semblance for GCAs and in scenarios when waveform distortion is considerable. Higher frequencies, topographic reflections, ground coupling, and large sourcereceiver distances will all decrease coherence between waveforms, suggesting that envelope stacking may be more effective in these scenarios. Similar to work by Fee et al., (2016), our analysis gives confidence that RTM and RTM-FDTD can be used to locate GCAs from explosions on local seismic networks.
We suggest a number of avenues for future work on local infrasonic source location and detection. We note that other methods could be used to compute accurate acoustic travel times, and infraFDTD3D is only one accurate implementation with moderate computational cost. Simpler, more efficient methods for travel time calculations that still incorporate topography should be explored. For example, one could solve the Eikonal equation by the finite-difference method to compute the travel times, such as is commonly done in seismology (e.g. Vidale, 1988). Such a formulation for infrasonic waves encountering steep topography should be validated against infraFDTD3D and other accepted methods. If valid, this method would likely provide substantial decreases in computational speed. 3-D and timedependent atmospheric structure do not appear to limit the accuracy and effectiveness of the algorithms discussed here; in recent work, McKee et al., (2014) and Kim and Lees (2014) explored the effect of wind and different celerities and found only minor changes in source location estimates. However, these effects will likely be important in scenarios that have stronger winds and larger source-receiver distances (e.g. winds >15 m/s or sourcereceiver distances greater than 20 km). Variations in celerity could also be considered in the grid search process, as demonstrated on larger grids (Walker et al., 2010;Walker et al., 2011;Sanderson et al., 2020). Near-real-time applications of RTM should also be explored, such as for monitoring networks at volcano observatories. We note the RTM algorithms can be applied to infrasound networks, GCAs on seismometers, or a combination. Corrections for transmission loss, such as those demonstrated in TRM by , could also be integrated into RTM and related back-projection techniques. We did not explore this concept here as the source locations were sufficiently accurate. The infraFDTD3D algorithm could be used to estimate the pressure loss as a function of range, which could then provide a relative weight for each station.

CONCLUSION
RTM-FDTD can detect, locate, and distinguish multiple explosions from multiple sources. The algorithm is particularly effective at detecting and locating explosions that may otherwise be challenging due to complex acoustic propagation. Similar to previous work, RTM assuming straight-line propagation works reasonably well at local distances (<15 km) assuming sufficient station coverage and direct propagation paths. However, by integrating travel times computed by solving the acoustic wave equations using numerical modeling, we are able to achieve higher location accuracy and stack amplitudes with RTM-FDTD. This improvement is pronounced in regions with considerable topographic barriers, such as-for our example-deployment configurations at Yasur and Sakurajima Volcanoes, which both have eruptive vents located near the bottom of sizable craters. RTM-FDTD located events within tens of meters of the actual vents in these scenarios, despite the infrasound networks being spaced up to 6 km from the source. The more accurate travel times provided by Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 620813 infraFDTD3D improved the alignment of the waveforms during stacking and increased the number of detected events. For a 12 h period at Yasur, we detected roughly twice as many events with RTM-FDTD than RTM. Location accuracy was also increased, with error ellipses 25-30% smaller for RTM-FDTD compared to RTM. Both RTM-FDTD and RTM were also able to detect and locate ground-coupled airwaves on seismometers over a large (16 × 16 km) region at Shishaldin Volcano. RTM-FDTD produces comparable results to the high-resolution Time Reversal Mirror algorithm, yet only has to be run once for each station rather than for each event. This results in notable computational savings. Our open-source RTM code provides a straightforward means to detect and locate infrasonic sources from natural and anthropogenic sources over a wide variety of ranges and sources. Our results suggest that although simple acoustic propagation assumptions may permit a reasonable source location and alignment of waveforms, numerical modeling will increase the location accuracy and detection capability, especially when topography is present in volcanic and mountainous terrain. Future work should explore the effect of dynamic and heterogeneous atmospheres across larger ranges, investigate real-time applications, and test more efficient methods for computing accurate travel times.

DATA AVAILABILITY STATEMENT
The datasets analyzed for this study can be found at the Incorporated Research Institutions for Seismology Data Management Center (https://ds.iris.edu/ds/nodes/dmc/) under network codes 3E (Yasur), YO (Sakurajima), and AV (Shishaldin). The RTM code is available at https://github.com/ uafgeotools/rtm. RTM-FDTD travel times for each scenario are also provided as Supplementary Material. The infraFDTD3D source code is available upon request to KK.

AUTHOR CONTRIBUTIONS
DF performed the analysis and wrote the manuscript, with input from all co-authors. DF and LT wrote the RTM code, based partially on code by RS, while KK wrote the infraFDTD3D code. DF, AI, RM, and AJ helped collect the Yasur data; and DF, RM, and colleagues collected the Sakurajima data. All authors contributed to early data interpretations and provided motivation for this study.
FUNDING DF, LT, AI, JJL, and MH acknowledge support from the Alaska Volcano Observatory through the U.S. Geological Survey Volcano Hazards Program. DF also acknowledges support from NSF Grant EAR-1614323. AI also acknowledges support from NSF Grant EAR-1952392. RS and RM acknowledge NSF grants EAR-1614855 and EAR-1847736. ADJ was supported by the New Zealand Ministry of Business Innovation and Employment (MBIE) Strategic Science Investment Fund (SSIF) and the GNS Science Capability Fund.