Remotely imaging seismic ground shaking via large-N infrasound beamforming

Seismic ground motion creates low-frequency atmospheric sound (infrasound) that is detectable at remote sensor arrays. However, earthquake infrasound signal analysis is complicated by interference between multiple waves arriving at sensors simultaneously, reducing the accuracy and detail of ground motion detection. Here we show that individual waves in complicated wavefields can be resolved by recording infrasound on large-N arrays and processing with CLEAN beamforming. Examining both a local (M3.5, purely tropospheric infrasound propagation) and regional earthquake (M6.5, upper-atmospheric returns), we detect infrasound from tens of km away and up to several hundred km away respectively. Source regions span arcs of approximately 90°, indicating that although detection bias does occur (most likely from atmospheric winds) the recorded infrasound sources are widely dispersed and not simply epicentral. Infrasound-based remote detection of ground motion over wide areas can complement point measurements by seismometers and spur innovations in earthquake research and real-time hazard monitoring.


Introduction
Ground shaking exerts normal tractions on the atmosphere, exciting infrasound waves; due to the strong wave speed contrast between the ground and air, this seismic-coupled infrasound propagates at steep incidence angles to the ground surface (fig.1).Consequently, earthquake infrasound detections are typically either 1) "primary" or "local" infrasound that converts from seismic to infrasound at the recording site and propagates upward with low slowness, or 2) diffracted "secondary" infrasound that converts from seismic to infrasound at topographic features and propagates to recording sites subhorizontally (i.e., with high slowness) 1.,2. .A single infrasound array can record secondary infrasound from nearby or distant ground shaking, making secondary infrasound a promising low-latency ground motion monitoring and research tool.Infrasound waves that travel up to tens of km ("proximal" secondary infrasound 3. ) follow simple, nearly linear paths; however, waves traveling longer distances are influenced by refractions in the troposphere 4. or upper atmosphere 5. , which depend on time-variable atmospheric wind and temperature structures.Although atmospheric models can be used to predict infrasound arrival patterns, unresolvable short-time-scale atmospheric variability affects their accuracy 4. .Proximal infrasound, confined to the troposphere, experiences milder atmospheric path effects, but effects from topographic obstruction and diffraction can be significant 6. .The potential value of secondary infrasound in ground motion monitoring motivates efforts to backproject it to its sources 3.,7. .Key steps in the backprojection process include 1) calculating the backazimuth of secondary infrasound, 2) estimating seismic/infrasound travel times from the earthquake hypocenter to candidate conversion points to the infrasound array, and 3) identifying seismic-to-infrasound conversion points considering the total travel time and infrasound backazimuth.The backprojection results can be complicated by ambiguity in the type (and speed) of seismic waves and potential multipathing.

Slowness component labels describe horizontally-diffracted secondary infrasound (3 s/km horizontal slowness) and vertical primary infrasound (near-zero horizontal slowness).
The first step of the backprojection process-finding the direction toward infrasound sources-requires beamforming analyses of data from infrasound arrays (multiple sensors separated by approximately tens of meters).Beamforming with infrasound arrays is a common method for locating diverse phenomena, e.g.volcanic explosions, nuclear tests, ocean swells, lahars, and avalanches 8.-12. .Usual beamforming methods identify a single best-fit propagation direction per time interval analyzed, which is acceptable when the wavefields being analyzed have one dominant signal and minor background noise.However, secondary earthquake infrasound arrives from many directions simultaneously and cannot be fully described using single-solution beamforming.We present two methodological improvements that enable beamforming of simultaneous arrivals from multiple directions, enabling scientists to more fully exploit information on infrasound sources and atmospheric structure contained in the wavefield.
First, we use the CLEAN beamforming algorithm, which identifies distinct wavefield components by iteratively deconvolving the array response from the slowness spectrum (fig.2).
CLEAN has been applied to various physical problems, e.g., microseism seismic signals 13. and microbarom infrasound 14. , but has seen little use for studying rapidly changing wavefields like secondary earthquake infrasound.
Second, we improve the resolution of CLEAN beamforming by recording data with an unprecedented "large-N" infrasound array (meaning many more sensors than the 3 required to beamform a single wave).Large-N arrays benefit from having a compact array response (enabling distinction of signals arriving from different directions, fig. 3) while avoiding issues with spatial aliasing and loss of coherence between sensors 15. , with ancillary benefits of noise reduction and redundancy in case of sensor failure.
The utility of large-N infrasound arrays, likely to include enhanced detection of weak sources and changes in atmospheric structure, remains largely untested.However, the related field of seismology has recently been transformed by the ability to distinguish weak signals and multiple sources through the use of large-N seismometer arrays, leading to revolutionary improvements in seismic imaging at local 16.,17.and continental 18. scales, with benefits scalable to array size up to at least thousands of sensors.

Fig. 2. Demonstration of the iterative CLEAN method. In each iteration, we identify the maximum power of the slowness spectrum and its slowness coordinates (sx,sy). We define the component to remove as the array response centered at that slowness coordinate with power equal to the spectrum's max power times a factor phi (we use phi = 0.1 when analyzing earthquakes in this paper, but use a larger value in this figure for demonstration clarity).
Although the component to subtract and the results are shown as slowness maps, the actual subtraction is performed in the cross-spectral domain.Finally, we add a delta function with the corresponding slowness coordinates and power to the clean spectrum.This permits the iterative removal of the most powerful signal from the slowness spectrum and accounts for it as a discrete wave in the clean spectrum instead.

Fig. 3. Effect of large-N arrays in conventional beamforming (not CLEAN). (A) Array geometries used in this and other figures: the PARK 22-element large-N array and a threeelement subset. (B) The three-element array's response has a broad main lobe and close sidelobes, compared to (C) the compact response of the large-N array. Dashed circles in (B)-(C)
are wavenumbers corresponding to infrasound with frequency of 2 Hz (inner) and 10 Hz (outer).

(D) The slowness spectrum of secondary earthquake infrasound calculated using the threeelement array is blurry and less detailed than (E), the slowness spectrum of the same data calculated using the large-N array. Dashed circles in (D)-(E) indicate 3 s/km (the slowness of sound).
Despite its benefits, the large-N strategy (in this paper, N=17-22 sensors) is difficult because of the lifetime costs and installation/maintenance effort for instrumentation.Passive-source seismology overcame similar challenges through the use of low-cost, low-effort "nodal" seismometers 16. .In infrasound, the Gem infrasound logger 19.fills a similar niche with its low cost, low power needs, and rapid cable-free deployment; its use in our large-N array made installation and maintenance of this site feasible.

Results
We examine two earthquakes recorded during low-noise nighttime conditions at the PARK infrasound array (fig.3), which we installed to monitor local aftershocks following the M6.5 2020-03-31 Stanley, Idaho earthquake 20.,21. .The first earthquake we study is a M3.5 aftershock 22. that occurred beneath PARK and produced more than one minute of correlated infrasound, beginning with primary infrasound and followed by proximal secondary infrasound from many backazimuths.The second earthquake we study is a M6.5 regional earthquake 23.("Monte Cristo") that occurred 718 km away with a backazimuth of -160°.The regional earthquake includes tens of minutes of infrasound, including primary infrasound, proximal secondary infrasound, and distant secondary infrasound.Of the 22 sensors originally deployed in PARK, 20 recorded the local aftershock and 17 recorded the regional earthquake.

Local M L 3.5 Event
With the local earthquake example (0.6 km S, 7.5 km deep, ~1 km uncertainty via relocation analysis 22. ), primary infrasound from incident seismic waves is detected in the seconds after the earthquake onset, and is followed by tens of seconds of secondary infrasound (fig.4).Results from the CLEAN analysis are plotted as image colors and show initial low-slowness primary infrasound followed by a longer phase dominated by high-slowness secondary infrasound.
Secondary infrasound waves come from multiple directions that converge over time toward the northwest.Standard beamforming analyses only identify the slowness with maximum semblance or power for each time window; such results are plotted over the CLEAN results as dots.
Considerably more details can be observed in the CLEAN results than with standard beamforming.
Regional M L 6.5 Monte Cristo Event With the regional earthquake example 23.(720 km, -160° backazimuth), we observe a lower signal-to-noise ratio than the local event, with pervasive background noise before and during the event (e.g., the 90°-180° quadrant, and around -120°, fig.5C).We suppressed the noise by calculating the pre-event power distribution at each azimuth and plotting only the signals that significantly exceed background noise (fig.5D).Because all arrivals travel at minimum the full distance from the hypocenter to the array, converting from seismic to infrasound waves en route, and because seismic waves travel much faster than infrasound, we expect waves that converted near the array to arrive earliest, waves that converted more distantly to arrive at intermediate times, and waves converted near the epicenter to arrive later still.

Infrasound arrives from many directions initially, whereas the latest-arriving infrasound comes strictly from the northwest. Repeating the analysis using only three stations yields consistent but less-precise results for both (D) horizontal slowness and (E) backazimuth.
In contrast to the local event's minute of detectable signals, waves from the regional event arrive in a complicated pattern spanning over 40 minutes (fig.5).Like the local event, we interpret the first signals as primary infrasound.Body wave arrivals modeled using the 1-D IASP91 Earth model range from 95-214 s, in approximate agreement with the primary infrasound arrival times of 120-280 s.Secondary infrasound follows the primary infrasound, and due to its early arrival time (150-500 s after earthquake), we infer it originated locally (within tens of km of the array).During an ensuing gap in secondary infrasound from the regional earthquake, two local signals are detected: primary and secondary infrasound from a small local earthquake between 670-710 s, and an inferred vehicle at 910-935 s (whose rapidly changing backazimuth agrees with highway speeds on nearby ID-21).Secondary infrasound from the regional earthquake arrives from the south and east between 900-2200 s; its arrival times indicate that it originated far from both the array and the earthquake epicenter.Finally, infrasound arrives from the epicenter's backazimuth at times similar to the 2470 s arrival time predicted using the infraGA atmospheric ray-tracing model.A subsequent signal from 2600-2800 s with variable horizontal slowness and rapidly changing backazimuth is attributed to an aircraft.The full array is required to see most of these features (fig.5D); pervasive noise obscures them with the 3element array (fig.5G).

(E)
With the reduction from 17 to 3 sensors, primary infrasound from the regional earthquake is still clear but the local earthquake is not.Horizontal slowness is poorly resolved in general; unlike in (B), 3 sensors fail to detect the sustained horizontally-propagating infrasound.

(F) The 90°-180° backazimuth range is still dominant here, but the -130° to -110° range is not visible. The features seen in (B) are less visible. (G)
The features that were obvious in (D) are unresolvable.

Sources of secondary infrasound
Backprojected possible secondary infrasound sources correspond partly to topographic features (fig.6).For the local event (beneath PARK), secondary infrasound generation strictly occurs on the near side of mountain crests, suggesting that topographic obstruction prevents infrasound radiated on the far side of crests from reaching the infrasound array (fig.6a).In the regional event, non-local secondary infrasound is first observed from the ranges just south of the Snake River Plain and not from the nearer lowlands (fig.6b).In both events, mountains and basins both appear within source regions, possibly due to uncertainty from seismic waves propagating at multiple speeds.Source regions do not strictly coincide with mountain ranges in either event.We credit the detection of a wide range of secondary infrasound sources to the superior resolving capabilities of CLEAN beamforming with a large-N array.

(B) Map of the western US
showing secondary infrasound in the regional earthquake identified in fig.5d between 1000-2500 s after the earthquake (local secondary infrasound is excluded because of the very high uncertainty in seismic travel times).Secondary infrasound observed at PARK originates mostly from the southwest, south, and southeast, and the source region follows some mountain ranges.
In both earthquakes, secondary infrasound overwhelmingly arrives from approximately 90°wide sectors that coincide with the dominant background noise backazimuths (fig.4c, 5c-d), suggesting that both noise and secondary infrasound may be ducted by favorable atmospheric conditions.Importantly, the direction toward the epicenter is not strongly related to the backazimuths of secondary infrasound: the local earthquake epicenter approximately coincides with the array, and the regional earthquake epicenter is on the edge of the detected source region.

Effectiveness of advanced array analysis
For both the local (fig.4) and regional earthquakes (fig.5), the large-N array provides significantly more detailed information about earthquake and ambient noise infrasound wavefields than the three-element sub-array, regardless of whether single-solution or CLEAN beamforming is used (fig.4).In particular, the large-N array provides excellent horizontal slowness resolution, indicated by peak slowness matching sonic slowness (3 s/km) for most time windows during the aftershock's secondary infrasound period starting at 5 s (fig.4b).By contrast, estimates of horizontal slowness using the three-element array have significant scatter during the secondary infrasound period (fig.4d), resulting in many signals with lower slowness estimates being misclassified as primary infrasound (seismic arrivals) and omitted from the secondary infrasound image (fig.4e).The poor horizontal slowness resolution of the threeelement array is even more pronounced for the regional event (fig.5b), probably due to its low signal-to-noise ratio (fig.5a), which complicates the separation of primary infrasound.Even a compact source like the aircraft, obvious with the large-N array, is nebulous with the 3-element sub-array.Additionally, the process of identifying and removing likely background noise is less effective with the sub-array than with the large-N array (fig.5d, 5g).Finally, large-N arrays have the significant benefit of redundancy; the loss of 2-5 sensors in a 22-element array had little effect on its resolving capabilities, whereas 3-element arrays lose their ability to calculate a wave's backazimuth with the loss of a single sensor.
The benefits of CLEAN vary between our two case studies.In the high signal-to-noise ratio local earthquake, the traditional beams (dots in fig.4) yield an incomplete representation of the wavefield, but do capture major features revealed by CLEAN (fig.4b-c).However, in the low signal-to-noise ratio regional earthquake, CLEAN is required for a step that we found essentialcalculating the background noise distribution at each azimuth and plotting only the signals that significantly exceed background noise (fig.5c-d).
By design, CLEAN approximates a wavefield as the convolution of the array response with a finite number of delta functions in slowness space.Consequently, it is inherently prone to misrepresenting diffuse sounds as a combination of discrete waves 13. .When studying diffuse wavefields like secondary infrasound, this weakness can be mitigated by using large-N arrays because their compact array response leads to diffuse wavefields being approximated as many closely-spaced points instead of a few sparse points.
Large-N arrays incur costs including funds for instrument purchasing, field effort and site footprint, and computational expense.These can be mitigated by choice of instrumentation (using low-cost, rapid-deploy, low-maintenance devices) and analysis settings (e.g., number and size of analysis time windows).Researchers must assess their constraints and analytical needs when determining array size.
The cost of CLEAN is mainly configuration effort (running tests to optimize tuning parameters) and runtime.For example, analyzing the 85-second local earthquake on a modern laptop, traditional beamforming (the array_processing() function in obspy 24. ), took 0.6 seconds (3 elements) and 16 seconds (20 elements), whereas the CLEAN analysis took 58 seconds (3 elements) and 3368 seconds (20 elements).Computational expense of advanced methods may restrict their use in real-time applications, or when large amounts of data must be processed.

Advancing acoustic earthquake monitoring
We distinguish primary infrasound from secondary infrasound solely on the basis of horizontal slowness.Co-deploying a seismic array alongside the infrasound array would reduce the potential for ambiguity between steeply incident secondary infrasound (with unusually low horizontal slowness) and primary infrasound coupled to slow seismic waves (with unusually high horizontal slowness).Additionally, it would elucidate air-to-ground and ground-to-air coupling (common in primary earthquake infrasound and for high-amplitude infrasound waves 25. ); this will reduce ambiguity in secondary infrasound backprojection by determining conditions that are favorable to seismic-infrasound conversions.We also note the value of monitoring ground shaking remotely on other planets, where installing large networks of seismic sensors is costly (e.g., Mars) or impossible (e.g., Venus).Ground-based or balloon-borne 26.large-N acoustic arrays may accelerate the study of tectonics and planetary structure elsewhere in the solar system.

Fig. 6 .
Fig. 6.Maps of inferred secondary infrasound sources.Because seismic waves propagate at a