The 2010 Haiti earthquake revisited: an acoustic intensity map from remote atmospheric infrasound observations

In the days following the January 12, 2010 M w 7 Haiti earthquake the shaking intensity near the epicenter was overestimated and the spatial extent of the potentially damaging shaking was underestimated. This was due to the lack of seismometers in the near-source region at the time of the earthquake. Besides seismic waves, earthquakes generate infrasound, i.e., inaudible acoustic waves in the atmosphere. Here we show that infrasound signals, detected at distant ground-based stations, can be used to generate a map of the acoustic intensity, which is proportional to the shaking intensity. This is demonstrated with infrasound from the 2010 Haiti earthquake detected in Bermuda, over 1700 km away. Wavefront parameters are retrieved in a beamforming process and are backprojected to map the measured acoustic intensity to the source region. The backprojection process accounts for horizontal advection effects due to winds and inherent uncertainties with regard to the time of detection and the back azimuth resolution. Furthermore, we resolve the ground motion polarity in the epicentral region and use synthetics generated by an extended infrasound source model to support this result. Infrasound measurements are conducted globally for the veriﬁcation of the Comprehensive Nuclear-Test- Ban Treaty and although the network was designed to provide global coverage for nuclear explosions in the atmosphere, it is shown in this paper that there is also global coverage for the estimation of acoustic shaking intensity. In this study, we lay the groundwork that can potentially make infrasound-based ShakeMaps a useful tool alongside conventional ShakeMaps and a valuable tool for earthquake disaster mitigation in sparsely monitored regions.

In the days following the January 12, 2010 M W 7 Haiti earthquake the shaking intensity near the epicenter was overestimated and the spatial extent of the potentially damaging shaking was underestimated (a). This was due to the lack of seismometers in the near-source region at the time of the earthquake.
Post-earthquake scientific efforts included surveying damage on the ground, formulating GMPEs using seismic stations rapidly deployed to record aftershocks, and more complex forward propagation models. Based on these efforts the USGS updated the initial ShakeMap to better explain the observations.
The updated ShakeMap (b), compiled by the USGS in January 2017, is much more detailed than the initial ShakeMap that did not incorporate any ground motion measurements nor reported damage. sured at IS51 to source patches in the epicentral region (Fig. 5a). In addition, we are able to detect the ground motion polarity of the compressional and dilatational quadrants around the epicenter (Fig. 5b). This observation is supported by a simulation of an extended infrasound source using the Rayleigh integral (Fig. 6).

Data acquisition and beamforming results
IS51 is a four-element array located on Bermuda island, situated 1738 km from the epicenter in Haiti (Fig. 2a). The array is equipped with MB2005 microbarometric sensors that have a flat frequency response between 0.01 and 27 Hz. A rosette wind-noise designed to include back azimuth and apparent velocity values of interest. The back azimuth values range between 180 • and 260 • and are spaced by 1 • . This range is selected to avoid detection of microbarom sources in the Atlantic. The apparent velocity values range between 0.28 km/s and 6 km/s. Between 280 and 450 m/s (the infrasonic signal range), the values are evenly spaced by 5 m/s, and between 450 m/s and 6 km/s (the seismic signal range), the spacing increases logarithmically from 16 to 200 m/s. This yields 8829 slowness vectors that are evaluated at each of the 26,700 processing time-windows. We used only three out of the four elements due to anomalous high noise levels at H1 (see Fig. S1 in Supplementary material), likely due to a malfunctioning wind-noise reduction system.
Characteristic wavefront parameters are extracted in the beamforming process, namely, the direction of arrival back azimuth (BAZ), the speed of horizontal propagation over the array apparent velocity (AV), and the signal coherency in terms of SNR ( Fig. 3c-e, respectively). Generally, epicentral infrasound signals are empirically characterized by a celerity (epicentral distance divided by the total travel-time) range of 0.34 to 0.31 km/s for stratospheric propagation and 0.31 to 0.28 km/s for thermospheric propagation (Evers and Haak, 2007). However, coherent signals (SNR > 0.7) are only detected between ∼5500 and ∼6100 seconds (celerity range of 0.32 to 0.28 km/s), mostly corresponding to the thermospheric celerity range. The right column in Fig. 3 focuses on detections that correspond to infrasound signals from the epicentral region. Only detections that fit our selection criteria (200 • < BAZ < 220 • , 320 < AV < 400 m/s, and SNR > 0.7) are shown. These are later used in the backprojection process. Inclination angles measured up from the horizontal are calculated using the relation φ = arccos(c/AV), where c=342 m/s is the local speed of sound (assuming zero wind) at IS51 and AV is the observed apparent velocity for each detection. This yields observed inclination angles from 7 • to 25 • .

Seismo-acoustic coupling and propagation from Haiti to Bermuda
Waveguides in the atmosphere facilitate long-range infrasound propagation and its detection at ground-based stations (Waxler and Assink, 2019). These waveguides are defined by temperature these spatial uncertainties to have a Gaussian distribution around the resolved location for each detection, resulting in a detection patch.
The acoustic intensity I associated with each detection window i at time t d with start time t 0 = t d − T w /2 and end time t 1 = t d + T w /2 is calculated as I i = t1 t0 (p 2 )dt, where T w is the processing time-window length, p is the pressure, and dt is the temporal resolution of the signal (Fig. 3f). As a detection window in time is mapped to a detection patch in space, this value is used to estimate the acoustic intensity associated with the part of the signal that was contributed by that patch. The region outlined in the acoustic intensity map (Fig. 5a), covers a much larger radiated acoustic pressure field, p(r r r, ω), from a planar vibrating surface can be computed by the Rayleigh integral (Wapenaar and Berkhout, 1989;Green et al., 2009): where r r r = (x, y, z) is a receiver location, r r r 0 is a location on the vibrating surface, ω is angular frequency, k is the medium wavenumber, v(r r r 0 , ω) is the complex spectral component of the surface vertical velocity, and the kernel, exp(−ik|r r r − r r r 0 |)/|r r r − r r r 0 |, is the Wavefront parameters are retrieved in a beamforming process and are backprojected to map the measured acoustic intensity to the source region. The backprojection process accounts for horizontal advection effects due to winds and inherent uncertainties with regard to the time of detection and the back azimuth resolution. Furthermore, we resolve the ground motion polarity in the epicentral region and use synthetics generated by an extended infrasound source model to support this result.
We lay the groundwork that can potentially make infrasound-based ShakeMaps a useful tool alongside conventional ShakeMaps and a valuable tool for earthquake disaster mitigation in sparsely monitored regions.

ARRAY PROCESSING
Characteristic wavefront parameters are extracted in the beamforming process, namely, the direction of arrival back azimuth (BAZ), the speed of horizontal propagation over the array apparent velocity (AV), and the signal coherency in terms of SNR. Generally, epicentral infrasound signals are empirically characterized by a celerity (epicentral distance divided by the total travel-time) range of 0.34 to 0.31 km/s for stratospheric propagation and 0.31 to 0.28 km/s for thermospheric propagation. However, coherent signals (SNR > 0.7) are only detected between ∼5500 and ∼6100 seconds (celerity range of 0.32 to 0.28 km/s), mostly corresponding to the thermospheric celerity range. The right column focuses on detections that correspond to infrasound signals from the epicentral region. Only detections that fit our selection criteria (200°< BAZ < 220°, 320 < AV < 400 m/s, and SNR > 0.7) are shown. sured at IS51 to source patches in the epicentral region (Fig. 5a). In addition, we are able to detect the ground motion polarity of the compressional and dilatational quadrants around the epicenter (Fig. 5b). This observation is supported by a simulation of an extended infrasound source using the Rayleigh integral (Fig. 6).

Data acquisition and beamforming results
IS51 is a four-element array located on Bermuda island, situated 1738 km from the epicenter in Haiti (Fig. 2a). The array is equipped with MB2005 microbarometric sensors that have a flat frequency response between 0.01 and 27 Hz. A rosette wind-noise reduction system is used to reduce wind noise over the infrasonic frequency band by spatially averaging the pressure field in the vicinity of each infrasound sensor (Hedlin and Raspet, 2003; Gabrielson, 2011). The pressure field is continuously recorded at designed to include back azimuth and apparent velocity values of interest. The back azimuth values range between 180 • and 260 • and are spaced by 1 • . This range is selected to avoid detection of microbarom sources in the Atlantic. The apparent velocity values range between 0.28 km/s and 6 km/s. Between 280 and 450 m/s (the infrasonic signal range), the values are evenly spaced by 5 m/s, and between 450 m/s and 6 km/s (the seismic signal range), the spacing increases logarithmically from 16 to 200 m/s. This yields 8829 slowness vectors that are evaluated at each of the 26,700 processing time-windows. We used only three out of the four elements due to anomalous high noise levels at H1 (see Fig. S1 in Supplementary material), likely due to a malfunctioning wind-noise reduction system.
Characteristic wavefront parameters are extracted in the beamforming process, namely, the direction of arrival back azimuth (BAZ), the speed of horizontal propagation over the array apparent velocity (AV), and the signal coherency in terms of SNR (Fig. 3c-e,  respectively). Generally, epicentral infrasound signals are empirically characterized by a celerity (epicentral distance divided by the total travel-time) range of 0.34 to 0.31 km/s for stratospheric propagation and 0.31 to 0.28 km/s for thermospheric propagation (Evers and Haak, 2007). However, coherent signals (SNR > 0.7) are only detected between ∼5500 and ∼6100 seconds (celerity range of 0.32 to 0.28 km/s), mostly corresponding to the thermospheric celerity range. The right column in Fig. 3 focuses on detections that correspond to infrasound signals from the epicentral region. Only detections that fit our selection criteria (200 • < BAZ < 220 • , 320 < AV < 400 m/s, and SNR > 0.7) are shown. These are later used in the backprojection process. Inclination angles measured up from the horizontal are calculated using the relation φ = arccos(c/AV), where c=342 m/s is the local speed of sound (assuming zero wind) at IS51 and AV is the observed apparent velocity for each detection. This yields observed inclination angles from 7 • to 25 • .

Seismo-acoustic coupling and propagation from Haiti to Bermuda
Waveguides in the atmosphere facilitate long-range infrasound propagation and its detection at ground-based stations (Waxler and Assink, 2019). These waveguides are defined by temperature and wind gradients that refract part of the acoustic wavefield back toward the ground. The typical atmospheric waveguides, classified by the layers of the Earth's atmosphere, are the tropospheric, stratospheric, and thermospheric waveguides. In the troposphere

BACKPROJECTIONS
The severity of ground shaking generally decreases with distance from the epicenter, however, nearsurface geology, topography, and the source radiation pattern, contribute to local variations in ground shaking intensity. This variability is captured by the coupled acoustic pressure field over the disturbed region. The radiated signals from the different sub-patches in the near-source region remain ordered in time throughout the propagation from the epicentral region to Bermuda, and the variability of acoustic pressure perturbations from one location to another can be retrieved.
The acoustic intensity I (integral of the square of the pressure) associated with each detection window is mapped to the detection patch that best corresponds to the associated travel-time and back-azimuth.
Similarly, the pressure-time integral S in each detection window, which yields a positive or negative overall sum, indicates whether the detection patch mostly moved upward or downward. these spatial uncertainties to have a Gaussian distribution around the resolved location for each detection, resulting in a detection patch.

S. Shani-Kadmiel, G. Averbuch, P. Smets et al. Earth and Planetary Science Letters 560 (2021) 116795
The acoustic intensity I associated with each detection window i at time t d with start time t 0 = t d − T w /2 and end time t 1 = t d + T w /2 is calculated as I i = t1 t0 (p 2 )dt, where T w is the processing time-window length, p is the pressure, and dt is the temporal resolution of the signal (Fig. 3f). As a detection window in time is mapped to a detection patch in space, this value is used to estimate the acoustic intensity associated with the part of the signal that was contributed by that patch. The region outlined in the acoustic intensity map (Fig. 5a), covers a much larger radiated acoustic pressure field, p(r r r, ω), from a planar vibrating where r r r = (x, y, z) is a receiver location, r r r 0 is a location on the vibrating surface, ω is angular frequency, k is the medium wavenumber, v(r r r 0 , ω) is the complex spectral component of the surface vertical velocity, and the kernel, exp(−ik|r r r − r r r 0 |)/|r r r − r r r 0 |, is the

RAYLEIGH INTEGRAL
To better understand this result, a conceptual model is set up to simulate the acoustic pressure field from an extended infrasound source, as illustrated in (a). Pistons in the top-left and bottom-right quadrants (first and third) are prescribed a positive (upward) STF, and pistons in the top-right and bottom-left quadrants (second and fourth) are prescribed a negative (downward) STF. The activation of each piston is offset in time to mimic a radially propagating seismic wave with a moveout velocity of 3 km/s from the simulated epicenter at the center of the extended source.
Synthetic waveforms are calculated for a four-element array in the far-field, 150 km away from the epicenter at 45°to the nodal planes to mimic the orientation of IS51 with respect to the nodal planes of the Haiti earthquake. Wavefront parameters are extracted in a beamforming process, and then used in the backprojection process. The pressure-time integral S in each detection window is used to infer the ground motion polarity of each detection patch.
Extended source modeling and backprojection. (a) Extended source setup with four quadrants: red indicates upward motion, blue indicates downward motion. White contour lines indicate isochrons of piston activation time in seconds. (b) Source radiation pattern inferred from backprojection of synthetic infrasound signals using one array located 150 km away from the epicenter, 45°to the nodal planes in the direction indicated by the arrow. (c) Superposition of source radiation patterns inferred from backprojection of synthetic infrasound signals using two arrays located 150 km away from origin at 45°to the nodal planes in the direction indicated by the arrows.
S. Shani-Kadmiel, G. Averbuch, P. Smets et al. Earth and Planetary Science Letters 560 (2021)  piston is offset in time to mimic a radially propagating seismic wave with a moveout velocity of 3 km/s from the simulated epicenter at the center of the extended source. Synthetic waveforms are calculated for a four-element array in the far-field, 150 km away from the epicenter at 45 • to the nodal planes to mimic the orientation of IS51 with respect to the nodal planes of the Haiti earthquake. Wavefront parameters are extracted in a beamforming process, as described in Data acquisition and beamforming results, and then used in the backprojection process. The pressure-time integral S in each detection window is used to infer the ground motion polarity of each detection patch.
As before, three distinct patches are resolved (Fig. 6b) with the middle patch smeared across compressional quadrants one and dicate that energy from sources in the subsurface, evanescent coupled to the atmosphere and generated infrasound (Fig. 4), an (2) extended infrasound source modeling demonstrates that th mechanical coupling process, that is, the perturbation of the acous tic pressure field in response to shaking of the ground-atmospher interface, preserves the ground motion polarity that is determine by the faulting mechanism (Fig. 6).
Our analysis of the atmospheric propagation conditions in cluded climatological profiles as well as high-resolution ERA ECMWF atmospheric specifications. From this analysis it followe that propagation of infrasound was most likely facilitated by thermospheric waveguide (Fig. 4). This is further supported by th

DISCUSSION AND CONCLUSIONS
• Backprojections of infrasound signals have been shown to be in correlation with earthquake ground motions prior to this work but for shorter propagation ranges and only for stratospheric infrasound (Marchetti et al., 2016;Walker et al., 2013;Hernandez et al., 2018). For the first time, infrasound that has propagated over 1700 km is used to outline the region where shaking intensity is sufficiently large to lead to damage. • We demonstrate the potential of remote infrasound detections for mapping the acoustic intensity over an earthquake source region. • The expected infrasound travel-time over 2000 km assuming a thermospheric waveguide is approximately 2 hours under typical propagation conditions. This means that an acoustic intensity map can be produced faster than other methods such as damage analysis on the ground or from aerial and satellite imagery. This can be done for earthquakes almost anywhere on land or close to shore. The techniques presented here, together with the coverage extent, make it plausible to use infrasound as a global earthquake disaster mitigation technique for the first time.
Acoustic intensity map potential of the IMS. The green gradient shading indicates coverage of the currently installed infrasound stations (full circles). Green shading corresponds to distance out to 2000 km. Contour lines spaced 30 minutes apart indicate travel-time to the nearest station, calculated on the basis of thermospheric propagation. The light red shading indicates regions that will be covered by planned stations (empty circles). Dark red regions correspond to landmass that will remain uncovered by the IMS for this application. Fig. 7. Acoustic intensity map potential of the IMS. The green gradient shading indicates coverage of the currently installed infrasound stations (full circles). G corresponds to distance out to 2000 km. Contour lines spaced 30 minutes apart indicate travel-time to the nearest station, calculated on the basis of thermos gation. The light red shading indicates regions that will be covered by planned stations (empty circles). Dark red regions correspond to landmass that will rem by the IMS for this application.

S. Shani-Kadmiel, G. Averbuch, P. Smets et al. Earth and Planetary Science Letters 560
source characterization based on infrasound detections at groundbased stations in the far-field. In this study, we demonstrate the potential of remote infrasound detections for mapping the acoustic intensity over an earthquake source region. We defer derivation of absolute ground mo-