The 2010 Haiti earthquake revisited

(1)

The 2010 Haiti earthquake revisited

An acoustic intensity map from remote atmospheric infrasound observations

Shani-Kadmiel, Shahar; Averbuch, Gil; Smets, Pieter; Assink, Jelle; Evers, Läslo

DOI

10.1016/j.epsl.2021.116795

Publication date

2021

Document Version

Final published version

Published in

Earth and Planetary Science Letters

Citation (APA)

Shani-Kadmiel, S., Averbuch, G., Smets, P., Assink, J., & Evers, L. (2021). The 2010 Haiti earthquake

revisited: An acoustic intensity map from remote atmospheric infrasound observations. Earth and Planetary

Science Letters, 560, 1-11. [116795]. https://doi.org/10.1016/j.epsl.2021.116795

Important note

To cite this publication, please use the final published version (if applicable).

Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

(2)

Contents lists available atScienceDirect

Earth

and

Planetary

Science

Letters

www.elsevier.com/locate/epsl

The

2010

Haiti

earthquake

revisited:

An

acoustic

intensity

map

from

remote

atmospheric

infrasound

observations

Shahar Shani-Kadmiel

a,b,

∗

,

Gil Averbuch

a,b

,

Pieter Smets

a,b

,

Jelle Assink

b

,

Läslo Evers

b,a

a_Faculty_of_Civil_Engineering_and_Geosciences,_Department_of_Geoscience_and_Engineering,_Delft_University_of_Technology,_Delft,_the_Netherlands b_R&D_Seismology_and_Acoustics,_Royal_Netherlands_{Meteorological}_Institute_(KNMI),_De_Bilt,_the_Netherlands

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

Articlehistory:

Received21July2020

Receivedinrevisedform23January2021 Accepted29January2021 Availableonlinexxxx Editor:R.Bendick Keywords: ShakeMap infrasound seismo-acoustics beamforming numericalmodeling

In the days following the January 12, 2010 Mw 7 Haiti earthquake the shaking intensity near the epicenter was overestimated and the spatial extent of the potentially damaging shaking was underestimated. Thiswas dueto the lack ofseismometers in the near-source regionat the time of theearthquake. Besides seismicwaves,earthquakesgenerateinfrasound, i.e.,inaudibleacoustic waves inthe atmosphere. Here weshow that infrasound signals, detectedatdistant ground-based stations, canbeusedtogenerateamapoftheacousticintensity,whichisproportionaltotheshakingintensity. Thisisdemonstratedwithinfrasoundfromthe2010HaitiearthquakedetectedinBermuda,over1700 kmaway.Wavefrontparametersareretrievedinabeamformingprocessandarebackprojectedtomap themeasuredacousticintensitytothesourceregion.Thebackprojectionprocessaccountsforhorizontal advection effects dueto winds and inherent uncertainties with regard to the time of detectionand the back azimuth resolution. Furthermore, we resolve the ground motion polarity in the epicentral regionand usesynthetics generated by anextended infrasound source model to supportthisresult. InfrasoundmeasurementsareconductedgloballyfortheveriﬁcationoftheComprehensive Nuclear-Test-BanTreatyand althoughthe network wasdesignedto provideglobalcoveragefornuclear explosions inthe atmosphere, it is showninthis paper thatthere is alsoglobal coverageforthe estimation of acousticshaking intensity.Inthisstudy, welaythegroundworkthatcanpotentiallymake infrasound-basedShakeMaps ausefultool alongside conventional ShakeMapsand avaluable tool forearthquake disastermitigationinsparselymonitoredregions.

1. Introduction

TheJanuary12,2010Mw7Haitiearthquakeisoneofthemost devastating earthquakedisasterinrecenthistory.This earthquake was precededbyseveralhistoricalearthquakesthat haveseverely damaged the city of Port-au-Prince throughout history (Bilham,

2010). However, no research-quality seismic stations were oper-atinginHaitipriortothisearthquake(Houghetal.,2010;Douilly etal.,2013).

To rapidly assess the potential impact and the necessary re-sponsemeasuresfollowinganearthquake,theU.S.Geological Sur-vey(USGS)launchedthePromptAssessmentofGlobalEarthquake for Response(PAGER) in 2007 (Jaiswal et al., 2011). A key com-ponentofthePAGERsystemistheShakeMap,whichindicatesthe

*

Correspondingauthorat:FacultyofCivilEngineeringandGeosciences, Depart-ment of Geoscience and Engineering,Delft Universityof Technology,Delft, The Netherlands.

E-mailaddress:s.shanikadmiel@tudelft.nl(S. Shani-Kadmiel).

earthquake impact through the distribution of shaking intensity. Theseverity ofgroundshakinggenerallydecreases withdistance from the epicenter, however, near-surface geology, topography, andthesource radiationpattern, contribute tolocalvariations in groundshakingintensity.Whereasearthquakesource characteris-tics,e.g.,locationandmagnitude,canberapidlydeterminedusing distantseismicstations,rapidgenerationofanaccurateShakeMap requiresgroundmotionmeasurements fromstationsinthe near-sourceregion(Earleetal.,2009).Duetothelackofseismometers during the 2010 Haitiearthquake disaster, the initial ShakeMaps overestimatedtheshakingintensityneartheepicenterand under-estimatedthe spatial extent of the potentially damaging shaking (Fig.1).

MuchscientiﬁcworkhasbeencarriedoutinHaitifollowingthe 2010earthquake disaster.These effortsincluded constructingthe initialseismichazardmaps(Frankeletal.,2010),deployingseismic stationstorecordandpreciselylocateaftershocks,gaininginsight intothefaultingmechanismsandgroundmotioncharacteristicsin the region (Houghet al., 2010; Mercier de Lépinay etal., 2011), estimatingpeakgroundaccelerationfromrigidbodydisplacement

https://doi.org/10.1016/j.epsl.2021.116795

(3)

Fig. 1. ShakeMaps.(a)TheShakeMapestimatedbytheUSGSapproximatelyonedayaftertheearthquake.(b)AUSGSupdatedShakeMapcompiledonJanuary27,2017. (Source:USGS).(Forinterpretationofthecolorsintheﬁgure(s),thereaderisreferredtothewebversionofthisarticle.)

(Hough et al., 2012), and damage assessment usingsatellite and aerial imagery aswell assurveying damage on theground (Cor-baneetal.,2011;Francisetal.,2016).Basedon theseefforts,the USGSupdatedtheinitialShakeMapestimatestobetterexplainthe surveyeddamage.The updatedShakeMap inFig.1b,compiledby theUSGSinJanuary 2017,is muchmoredetailedthantheinitial ShakeMap that did not incorporateany groundmotion measure-mentsnorreporteddamage.

Inadditiontoseismicwaves,earthquakesandunderground nu-clearexplosionsgeneratelow-frequencyacousticwavesinthe at-mosphere,knownasinfrasound. Bolt(1964) describedexceptional atmosphericwavesrecordedbyamicrobarographinBerkeley, Cal-ifornia, after the Mw 9.2 Alaskan earthquake ofMarch 27,1964, over3000kmaway.Multiplemechanismsplayarole inthe gen-erationofinfrasoundwavesfromasubsurfacesource.

Arrowsmith et al. (2010) discussed seismo-acoustic coupling mechanisms forearthquakes aswell asvariousother naturaland human-made sources. Intuitively, the mechanical disturbance of the ground-atmosphere interface compresses (and decompresses) the atmosphereabove and generates an acoustic pressure wave. Physically, boundaryconditions require thecontinuity on normal stressandverticaldisplacementalongtheinterface.Thusthe con-tribution of the mechanical coupling mechanism to the acoustic pressureperturbationinﬂatregions,isattributedtoseismicwaves thatcontainaverticalcomponent(e.g.,P,SV,andRayleigh). How-ever,steeptopographycanfacilitatecouplingofadditionalphases, evenonesthatarehorizontallypolarizedasshownbyGreenetal. (2009).

The ground-atmosphere and ocean-atmosphere interfaces are typicallyconsideredtobereﬂectiveboundariestotheseismicand

acousticwaveﬁelds.However,atlowfrequencies,thisboundary be-comestransparenttopartofthewaveﬁeld(Godin,2008,2011)and infrasonicsignalscanbegeneratedbyanother,lessintuitive mech-anism,evanescent wave couplingasobserved followingthe 2004 Mw8.1Macquarieridgeearthquake(Eversetal.,2014).

Seismo-acousticcoupledsignalscanoriginatefromthe epicen-ter(MutschlecnerandWhitaker,2005),aswellasfromsecondary sourcesofinfrasound away fromtheepicentralregion (LePichon etal., 2006; Greenet al., 2009; Marchetti etal., 2016).Given an origintime andlocation,seismo-acoustic signalscanbe backpro-jectedtolocatethepointontheEarth’ssurfaceatwhichcoupling has occurred. Backprojections of infrasonic signals from earth-quakes(Shani-Kadmieletal.,2018)andnuclearexplosions(Assink etal.,2018)provideinsightintotheseismo-acousticcoupling pro-cessandcanmapthedistributionoftheepicentralandsecondary sourcesofinfrasound. FollowingtheTohokuearthquakeofMarch 11,2011,andtheAmatriceearthquakeofAugust24,2016,ground motionsinferredfrombackprojectionwere showntobe in agree-ment withthe measured peak surface pressure in the epicentral region(Walkeretal.,2013;Hernandezetal.,2018).

AlthoughseismicsignalswerenotmeasuredinHaitiduringthe 2010 earthquake, infrasound was generated over the region and was detected by an International Monitoring System (IMS) array IS51onBermuda island,1738kmaway (Fig. 2). Inwhat follows, weﬁrstusearrayprocessingtoresolveplane-waveparametersof theinfrasonicsignalsdetectedinBermuda(Fig.3).Next,we exam-inethegenerationandpropagationofinfrasoundfromthe earth-quakesourceinthesubsurface(Fig.4).We thenusetheresolved wavefrontparameters andatmosphericpropagation conditionsto backproject the detections and map the acoustic intensity

(4)

mea-Fig. 2. Overviewmapshowingthe epicenter(star)and infrasoundarrayIS51on Bermudaisland(triangle).ThegreatcirclepathconnectingIS51withtheepicenter measures1738kmlongwiththeoreticalbackazimuthof209◦.Theboundsofthe backazimuthrangeusedforbeamformingareindicatedingray.Thearray conﬁgu-rationofIS51isshownininsetframe.

suredatIS51 tosourcepatches inthe epicentralregion (Fig.5a). In addition, we are able todetect the groundmotion polarity of the compressional anddilatational quadrants around the epicen-ter(Fig. 5b). Thisobservationis supported by asimulation ofan extendedinfrasoundsourceusingtheRayleighintegral (Fig.6). 2. Data acquisition and beamforming results

IS51 is a four-elementarray located on Bermudaisland, situ-ated 1738 km fromthe epicenter in Haiti(Fig. 2a). The array is equipped withMB2005 microbarometric sensors that havea flat frequencyresponsebetween0.01and27Hz.Arosettewind-noise reduction system is used to reduce wind noise over the infra-sonic frequencyband by spatially averaging the pressure field in the vicinity of each infrasound sensor (Hedlin andRaspet, 2003; Gabrielson, 2011). The pressure field is continuously recorded at a rate of 20 samples-per-second, and the waveform data are detrended, tapered, and band-pass filtered before time-domain beamforming(Fig. S1inSupplementarymaterial). Asecond-order Butterworth band-passfilterbetween0.45and2Hzis chosento reduceinterferencefromlow-frequencysignalsinthemicrobarom bandaswellashigherfrequencywind-noise.

A time-domain beamforming technique (Melton and Bailey,

1957) is used for the detection of coherent infrasound and the estimationof plane-waveparameters.It isadelay-and-sum tech-niquethat enhancescoherentsignalsandsuppressesthe incoher-entbackgroundnoise. Consequently,coherentsignalswith ampli-tudes belowthe background noise levels can be detected. A de-tection is based on the evaluation of a Fisher ratio (Fr), which corresponds to thesignal-to-noise ratio(SNR):Fr

=

N

·

SNR2

+

1, where N is the number of array elements (Melton and Bailey,

1957).

The waveform data are processedin time-windows Tw of 30

seconds with 99% overlap between successive windows. A large overlap yields detections witha high temporal resolution and is beneﬁcial for the backprojection algorithm. The samples are de-layed and summed over a horizontal slowness grid. The grid is

designedtoincludebackazimuthandapparentvelocity valuesof interest. The back azimuthvalues range between180◦ and 260◦ and are spaced by 1◦. This range is selected to avoid detection ofmicrobaromsources intheAtlantic.The apparent velocity val-uesrangebetween 0.28km/sand 6km/s. Between280and450 m/s (the infrasonic signal range), the values are evenly spaced by 5 m/s, andbetween 450 m/sand 6 km/s (the seismic signal range),the spacingincreaseslogarithmically from16to 200m/s. This yields 8829 slowness vectors that are evaluated at each of the 26,700 processingtime-windows. We usedonly three out of thefourelements dueto anomaloushighnoise levelsat H1(see Fig. S1inSupplementarymaterial), likelyduetoamalfunctioning wind-noisereductionsystem.

Characteristicwavefrontparametersareextractedinthe beam-forming process, namely, the direction of arrival back azimuth (BAZ),thespeedofhorizontalpropagationoverthearrayapparent velocity (AV),andthesignalcoherencyintermsofSNR(Fig. 3c-e, respectively). Generally, epicentral infrasound signals are empir-ically characterized by a celerity (epicentral distance divided by thetotaltravel-time) rangeof0.34to0.31km/sforstratospheric propagationand0.31to0.28km/s forthermosphericpropagation (EversandHaak,2007).However,coherentsignals(SNR > 0.7)are onlydetected between ∼5500and∼6100 seconds(celerityrange of0.32to0.28km/s), mostlycorresponding tothethermospheric celerity range. The right column in Fig. 3 focuses on detections that correspond to infrasound signalsfromthe epicentral region. Onlydetectionsthatﬁtourselectioncriteria(200◦ < BAZ < 220◦, 320 < AV < 400 m/s, and SNR > 0.7) are shown. These are later used in the backprojectionprocess. Inclination angles mea-sured up from the horizontal are calculated using the relation

φ

=

arccos

(

c

/

AV

)

, where c=342 m/s is the local speed of sound (assumingzerowind)atIS51andAVistheobservedapparent ve-locity for each detection. This yields observed inclination angles from7◦to25◦.

3. Seismo-acoustic coupling and propagation from Haiti to Bermuda

Waveguidesintheatmospherefacilitatelong-rangeinfrasound propagation and its detection at ground-based stations (Waxler andAssink,2019). Thesewaveguides are defined bytemperature andwindgradientsthatrefractpartoftheacousticwavefieldback toward the ground. The typical atmospheric waveguides, classi-fiedbythelayersoftheEarth’satmosphere,arethetropospheric, stratospheric, andthermospheric waveguides. In the troposphere and stratosphere, the adiabatic speed of sound is usually insuf-ficient to form waveguides. Therefore, wind in the direction of propagationisessentialtofacilitateground-to-groundpropagation. Thisdependenceonwind limitssoundpropagationinspecific di-rections.

In the troposphere, propagation is typically not efficient over long-ranges(distancesofmorethan100km)duetotherelatively limited size of the jet stream. In the stratosphere, advantageous conditionsbythemuchlargercircumpolarvortexsupportefficient propagationoververylong-ranges(WaxlerandAssink,2019) (dis-tancesexceeding2000km),though,stratosphericwindsare signif-icantlyvaryingona(sub)seasonal scale.Inthethermosphere,the temperaturegradientaloneissufficienttorefractinfrasoundtothe ground regardless of propagation direction. However, long-range propagationislessefficientcomparedtoastratosphericwaveguide duetothehigherattenuationinthethermosphere(Sutherlandand Bass,2004).

Ina horizontallylayered atmosphere, theﬁrst-order effectsof temperatureand horizontalwind in the directionof propagation areapproximatedbytheeffectivespeedofsound ceff(Godin,2002). Following Snell’s law, the horizontal componentof the speed of

(5)

Fig. 3. Beamformingresults.Theleftcolumnsummarizesallbeamformingresults.(a)Thespectrogramofthebestbeam(grayscaleindicatingpowerisshowntotheright of (d)).(b)Thebest beam.Subsequentframesshowthefollowing wavefrontparametersas retrievedineachtime-window:(c)Back azimuth(BAZ;truebackazimuth indicatedbydottedgrayline),(d)Apparentvelocity(AV;localspeedofsoundindicatedbydottedgrayline),and(e)Fisherratio(F-ratio).Verticalgraylinesindicatethe expectedtimerangeofepicentralinfrasound.Theredboxindicatestheenlargedrangeplottedintherightcolumn.Therightcolumnfocusesontheepicentralinfrasound detectionsthatareusedinthebackprojectionalgorithm.(f)AcousticintensityI usedtoscalethedetectionsinthebackprojectionalgorithmforFig.1c.(g)Normalizedsum ofthepressureusedtoscalethedetectionsinthebackprojectionalgorithmforFig.1d.Thebestbeam,low-passﬁlteredto1/30Hzandnormalized,isalsoplotted.(h-j) sameas(c-e).X-symbolsinframe(i)correspondtothetravel-timeandapparentvelocityofcalculatedeigenraysplottedinFig.4.ThecolorscaleindicatestheSNRofthe detection.TraveltimeinsecondssinceorigintimeandUTCtimeareindicatedonthebottomaxisandcelerity(averagehorizontalpropagationvelocity)isindicatedonthe topaxis.

propagation(referredtoasAVinsection2)remainsconstantwhile crossing betweenlayers.Therefore,a negativegradient intheceff proﬁlecausesthewave torefractupward,anda positivegradient causesthewave torefractdownward.Atthereturnheight ofthe wave, the AVof the wave andthe ceff inthe medium are equal.

The ceff profile,thereforedescribes, to first-order,the partofthe acousticwavefieldthatcangettrappedintheatmospheric waveg-uidesandthuscontributetolong-rangepropagation.Furthermore, return heightsare thus constrained to the parts ofthe ceff pro-filethatexceedtheinitialvelocityattheground(ceff ratio > 1)in

(6)

Fig. 4. Couplingandpropagationmodeling.(a)10evenlyspacedeffectivesoundspeedprofilesalongthegreatcirclepathconnectingthe epicenterinHaititoIS51on Bermuda.Topaxisshownforscaleandindicatesceffoftheleftmostprofile.Theclimatologyprofilesareindicatedbythinblacklines,theERA5ECMWFprofilesareindicated

inthickgraylines.(b,left)Averagedeffectivespeedofsoundproﬁle(valuesonthetopaxis)andseismicvelocityproﬁle(valuesonthebottomaxis)usedintheFFPandceff

raytracer.Verticaldashedlineindicatesceffatgroundlevelandregionsofpotentialreturnheightforaceffratio > 1relativetothepreviouswaveguideareoutlinedinred.

Horizontaldottedlinesindicateboundariesbetweentroposphere,stratosphere,mesosphere,andthermosphere.Forcomparison,theERA5ECMWFprofilesarealsoaveraged andplottedasthickgrayline.(b,right)Verticalsectionshowingacousticintensitytransmissionloss(TL)alongthepropagationpathfromasubsurfacesourceinHaitito IS51onBermudaisland.Thesourcesareindicatedbystarswithsizecorrespondingtorelativesourcemagnitude.IS51isindicatedbyatriangle.Eigenraysconnectingthe sourceregionandIS51calculatedusingthesameceffprofileareoverlaid.(c)VerticalcrosssectionofeffectivesoundspeedcalculatedfromERA5ECMWFspecifications

showingraysback-propagatedalongthetheoreticalbackazimuth±15◦.Solidlinesindicaterayscorrespondingtotheobservedrangeofinclinationangles(7◦- 25◦)and dashedlinesindicateraysoutsidethatrange.

thecaseofa singlewaveguideorthatofthewaveguidebelowin thecaseofmultiplewaveguides.Dependingonthedirectionofthe windwithrespecttothedirectionofpropagation,multiple waveg-uidescanco-exist.

Thepartofthewaveﬁeldconsistingofhighapparentvelocities, whichcorrespondtosteeperanglesofpropagation,donotbecome trappedintheatmosphericwaveguidesandarenotreturnedtothe ground.Inparticularcases,thesesignalscanbedetectedby satel-litesinlowearthorbits.Forinstance,at270kmaltitude,theGOCE missiondetecteddensityvariationsandverticaldisplacements cre-ated by post-seismic infrasound from the Tohoku earthquake of March11,2011(Garciaetal.,2013).

For the analysis of infrasound propagation conditions at the time ofthe earthquake,empirical atmosphericmodels, knownas climatologies, fortemperature (MSIS-00(Picone etal., 2002))and horizontalwind (HWM14(Drob etal., 2015))are comparedwith actual conditionsprovidedbythehigh-resolution ERA5reanalysis from the European Centre for Medium-Range Weather Forecasts (ECMWF) (Fig. 4a). The resulting ceff proﬁles from these, differ

in two waveguides: (1) In the troposphere, a weak tropospheric ductispresentintheECMWFproﬁles,whichisunresolvedinthe climatology. However, as mentioned above, this waveguide does not sustain long-range propagation andacts as a ﬁlter that pre-vents lower apparent velocities from reaching higher ducts. (2) Thestratosphericwaveguideisweakerthanintheclimatologyand weaker than the tropospheric waveguide, potentially creating an elevatedduct.ThiseffectstrengthensoverBermudaandtherefore makes itunlikely todetect stratosphericreturns accordingtothe atmosphericconditionsofECMWF.

AFastField Program(FFP)(Averbuch et al., 2020a) isutilized to simulatelong-rangeinfrasound propagation froma subsurface source inHaiti to IS51on Bermuda island. The FFPis a coupled seismo-acoustic solver that provides an exact solution for wave propagationinlayered media in thefrequency-wavenumber (f-k) domain. Within each layer, an exact solution of the wave equa-tiondescribesthe propagating waves,andthe layers arecoupled bytheboundaryconditions,i.e.,continuityofstressand displace-ment.Thisformofsolutionsupportsthesimulationofwave

(7)

prop-agation through stark discontinuities in the density andvelocity proﬁles likethe ground-atmosphereandocean-atmosphere inter-faces(Averbuchetal.,2020a,b).

Seismo-acousticcouplingisgovernedbythreepartsofthe seis-mic waveﬁeld spectrum. (1) Homogeneous body waves, (2) In-homogeneous (evanescent) body waves, and (3) Surface waves. However, onlytheparts ofthewaveﬁeldthat consistofapparent velocitiesthatcangettrappedintheatmosphericwaveguides con-tributetolong-rangepropagationandareofinterestinour model-ing.Moreover,thesewavesundergoenhancedsoundtransmission acrosstheground-atmosphereinterface,generatinginfrasonic sig-nalswithlargeramplitudesthanexpected(Godin,2011).

The computational domainconsistsof coupledsolid-ﬂuid me-dia:Inthesubsurface,seismicvelocitiesVP andVS correspondto

thestructure oftheuppercrustbelowHaiti(Posseeetal., 2019). Inthe atmosphere,tenevenly spacedclimatologyprofiles are av-eragedalongthegreatcirclepathfromHaititoBermuda(Fig.4b, left)todefinetheceff profileusedintheFFPsimulations(Fig.4b, right).

Sourceslipinversionsfromteleseismicdatashowthatslipwas distributed betweenthe hypocenter, at a depth of13 km,and 2 km depth, with the largest slip occurring at a depth of 10 km (Hayes et al., 2010). To account for this distribution, pressure-wave propagationisinitiated by ﬁveevenlyspaced pointsources withmagnitudedecreasingasafunctionofthesquarerootofthe vertical distance from10 kmto 2 kmdepth (Fig. 4b, right). The simulatedfrequencyis0.5Hz,andphasevelocitiesarebound be-tween280m/sand450m/s,correspondingtotheinfrasonicsignal range in the beamforming process. This also corresponds to the effective speed of sound range that can become trapped in the atmospheric waveguides. The pressure ﬁeld resulting from each source is stacked and the acoustic intensity (I) is calculated as I

=

p2

_/(

_ρ

_·

_c

e f f

)

,wherep isthepressure,and

ρ

isthedensity.The

transmissionloss(TL)iscalculatedasTL

= −

10log10

(

I

/

I0

)

,where I0 isthe acousticintensityjust above groundover thesource. In thisstudy weassumealossless propagationmediuminwhich at-tenuationisonlyduetogeometricalspreading.

Fig. 4b (right)showstheinfrasound propagationfromsources inthesubsurfaceunderHaititoBermudaintermsofTL.For sim-plicity, the values in the subsurface are masked. Return heights correspond to the parts of the ceff proﬁle that are outlined in Fig. 4b (left), where the ceff ratio > 1 relative to the ground or thepreviouswaveguide.

In addition, the averaged ceff profile is used to calculate ray paths in a one-dimensional (1D) layered windless atmosphere. EigenraysconnectingtheepicenterandIS51areshowninFig.4b. Eigen rays 1,2,and3fit theobservationsin termsofthe travel-timewithvaluesof5621,5588,and5985seconds,respectively,but onlyray number3iswithintheobserved inclinationanglerange (see section2andFig.3i).Recallingthattheaveragedceff profile, whichisbasedonclimatologydata,differsfromtheECMWF pro-files, itisunlikely that eigenrays1and2wouldhavepropagated to Bermuda; they would either get trapped in the tropospheric waveguideanddecayorbecometrappedintheelevatedductover Bermuda.Ray4isalsowithintheobservedinclinationanglerange butrays4,5,and6havetheoreticaltravel-timevaluesbeyondthe observedrangeofcoherentdetections(6322,6715,and6136 sec-onds,respectively).

To model propagation paths in further detail and investigate the inclinationrange ofpossiblerays incomingto IS51,we prop-agated rays using an in-house developed ray tracing algorithm, cast in spherical coordinates, that accounts for the full effect of thethree-dimensional(3D)inhomogeneouswindandtemperature ﬁelds (see section 2.2 in Smets (2018)). For the forward propa-gation (see rays on the left side of Fig. 4c), we used the high-resolution ERA5 ECMWF atmospheric speciﬁcations provided for

22:00UTC (∼7 minutes after the earthquake’s origintime). Rays arelaunched attheepicentertowardIS51overan azimuthrange of

±

15◦ ofthetheoreticalazimuthat0.5◦ azimuthal spacingand atinclination angles of0◦ to 40◦ every 0.5◦ above the horizon-tal.Forthebackpropagation(seeraysontherightsideofFig.4c), we used the atmospheric speciﬁcations provided for 23:00 UTC. Rays are launched at the central coordinates of IS51 toward the epicenterover an azimuthrange of

±

15◦ of thetheoretical back azimuth at 0.5◦ azimuthal spacing and at inclination angles of 0◦ to 40◦ every 0.5◦ above the horizontal. In order to facilitate backwardpropagation,thehorizontalcomponentsofthewindare reversed.

Intheforwardpropagationcase,rayspropagatesteeplyupward through the troposphere and stratosphere and exit through the topboundaryduetotheweaktroposphericductandevenweaker stratospheric duct over Haiti. In the back propagation case, rays withinclinationanglesbelow14.5◦ remain trappedinthe tropo-sphereuntil thewaveguide isweak enough toescape. Rays with inclinationanglesabove14.5◦ haveahighenoughapparent veloc-itytopropagate throughthetroposphere, butare not trappedby thestratosphere andcontinue upward. In eithercase, the strato-sphericwaveguideintheERA5ECMWFspeciﬁcationsistooweak to refract these rays back toward the ground. It is therefore as-sumedhereinthat infrasoundpropagationfromHaititoBermuda wasfacilitatedbyathermosphericwaveguide.

4. From detections in time to acoustic intensity map

The severity of groundshaking generally decreases with dis-tancefromtheepicenter,however,near-surface geology, topogra-phy,andthesourceradiationpattern,contributetolocalvariations ingroundshakingintensity.Thisvariabilityiscapturedbythe cou-pled acousticpressure ﬁeld over thedisturbed region (Walkeret al.,2013).Giventhatthewaveismostsensitivetotheatmosphere atthe returnheight (Assink etal., 2019) andthat the epicentral infrasoundpropagationwas facilitatedby athermospheric waveg-uide,itisvalidtoassumethattheatmosphereisconstant through-out the duration ofthe detected signal (∼600seconds inFig. 3). Thus the radiated signals from the different sub-patches in the near-sourceregion remain ordered intime throughout the prop-agationfromtheepicentralregiontoBermuda,andthevariability ofacousticpressureperturbationsfromonelocationtoanothercan beretrieved.

Thebackprojectionalgorithmmakesuseofthetravel-timeand backazimuthassociatedwitheachdetectiontime-windowtomap detections in time to their point of origin on the Earth’s sur-face.Building onthe backprojectionalgorithm outlined in Shani-Kadmiel et al. (2018), corrections for advection effects are in-corporated. Each cell in the grid-search domain is prescribed a backazimuth andinfrasound propagationvelocitycorrection val-ues,whicharederivedfromthecross-trackandalong-trackwinds. Furthermore, uncertainties with regard to the time of detection andbackazimutharetreated;thetimeofadetectionisassumed to be at the center of each processing time-window however, a detectioncanappearanywherewithin thedetectionwindow.This uncertaintytranslatestoaspatialerrorinthealong-trackdirection thatamountsto

±

Tw

·

ci

/

2 (∼5kmfor30second-long processing

windows Tw andinfrasound celerity ci of ∼0.3 km/s).

Addition-ally, uncertaintyin the detection’s back azimuth due to the res-olution of the beamforming (1◦), translates to a spatial error in thecross-track direction.Incontrast tothedetectiontime uncer-tainty that translates to a constant along-trackspatial errorsize, thesizeofthespatialerrorinthecross-trackdirectionduetothe backazimuthuncertaintygrowswithreceiver-source (backprojec-tion)distance.Forexample,attheepicenter,1738kmaway from IS51,this error amounts to ∼

±

15 km.In thisstudy, we assume

(8)

Fig. 5. (a)AcousticintensitymapfrombackprojectionofinfrasounddetectionsshowingtheacousticintensityI measuredatIS51mappedontothesourceregion.TheII

=5contourlinesfromtheinitial(brokenline)andupdated(solidline)USGSShakeMaps(Fig.1)areoverlaidingreen.(b)Sourceradiationpatternfrombackprojectionof infrasounddetections,redindicatesupwardmotion,blueindicatesdownwardmotion.Thebeachballrepresentationofthemomenttensor(NettlesandHjörleifsdóttir,2010) andnodalplanesareoverlaid.ThedirectiontoIS51onBermudaislandisindicatedbyanarrow(∼45◦tothenodalplanes).

thesespatial uncertaintiesto havea Gaussiandistributionaround the resolved location for each detection, resulting in a detection patch.

The acousticintensity I associated with each detection win-dow i at time td with start time t0

=

td

−

Tw

/

2 and end time

t1

=

td

+

Tw

/

2 is calculated as Ii

=

t1 t0

(

p

2

₎

_dt, _where _T

w is the

processing time-window length, p is the pressure, anddt is the temporal resolution of the signal (Fig. 3f). As a detection win-dow intime is mapped to adetection patchin space, thisvalue isusedtoestimatetheacousticintensityassociatedwiththepart ofthe signalthat was contributedby thatpatch. The region out-linedintheacousticintensitymap(Fig.5a), coversamuchlarger extent than the initial USGS ShakeMap estimates (Fig. 1a). For comparison, the region where shaking intensity was suﬃciently strong to cause damage is outlined by the instrumentalintensity II = 5contour line, extracted fromthe updated USGS ShakeMap

estimates (Fig. 1b). The same contour linefrom the initial USGS ShakeMap(Fig. 1a)outlines amuchsmallerregion.Inbothcases, The USGS ShakeMaps, aswell as ouracoustic intensitymap, are models.And eventhoughouracousticintensitymap andthe up-dated ShakeMap do not perfectly overlap, they both provide a better representation of the geographical extent of the observed damage.

5. Source radiation map

Similarly,thepressure-timeintegral (S)ineachdetection win-dow iscalculated as Si

=

t1

t0

(

p

)

dt and yields a positiveor

neg-ative overall pressure sum (Fig. 3g); This indicates whether the detection patch mostly moved upward or downward. It is noted that thesevaluesfollowthe lowpassedbestbeamtrace (Fig. 3g). The resolved radiation pattern in Fig. 5b outlines three distinct regions. The two blue patches reﬂect an overall downward mo-tion andcoincide withthe dilatationalquadrantsdepictedby the earthquakemomenttensor(NettlesandHjörleifsdóttir,2010).The red patchinthemiddlereﬂectsan overallupward motionandis smearedacrossthecompressionalquadrants.

Tobetterunderstandthisresult,aconceptualmodelissetupto simulatethe acousticpressure ﬁeld froman extendedinfrasound source,asillustrated inFig.6a.Ina homogeneoushalf-space, the

radiated acoustic pressure ﬁeld, p

(

rrr

,

ω

)

, from a planar vibrating surfacecan becomputedby theRayleighintegral (Wapenaar and Berkhout,1989;Greenetal.,2009):

p

(

rrr

,

ω

)

=

ik

ρ

c 2

π

S0 e−ik|rrr−rrr0|

|

rrr

−

rrr0

|

v

(

rrr0

,

ω

)

dS0

,

(1)

whererrr

= (

x

,

y

,

z

)

isareceiverlocation,rrr0isalocationonthe vi-bratingsurface,

ω

isangularfrequency,k isthemedium wavenum-ber, v

(

rrr0

,

ω

)

is the complex spectral component of the surface vertical velocity, and the kernel, exp

(

−

ik

|

rrr

−

rrr0

|)/|

rrr

−

rrr0

|

, is the Green’sfunctionforahomogeneousmedium.

InthisstudyweuseaGaussianpulseSTF

v

(

rrr0

,

t

)

=

ω

0

√

2

π

e −ω2 0(t−t0)2/2

_,

₍₂₎

witha central frequency

ω

0 = 0.1 Hz to drive the source veloc-ity,whichwhenintegrated,describesthepermanentdisplacement (see Section 2 ofthe Supplementary material). Using the Fourier transform, the spectral components of the STF are represented as: v

(

rrr0

,

ω

)

=

∞

−∞ v

(

rrr0

,

t

)

e−iωtdt

,

(3) and the waveforms at the receiver are obtained by the inverse transform p

(

rrr

,

t

)

=

1 2

π

∞

−∞ p

(

rrr

,

ω

)

eiωtd

ω

.

(4)

Toevaluate theintegral, the surface isdiscretized into individual pistonswitha400m radius(examplecodeisprovided Section2 oftheSupplementarymaterial).

Pistons in the top-left and bottom-right quadrants (ﬁrst and third)are prescribed a positive (upward) STF, andpistonsin the top-rightandbottom-left quadrants (secondand fourth)are pre-scribedanegative(downward)STF(Fig.6a).Theactivationofeach

(9)

Fig. 6. Extendedsourcemodelingandbackprojection.(a)Extendedsourcesetupwithfourquadrants:redindicatesupwardmotion,blueindicatesdownwardmotion.White contourlinesindicateisochronsofpistonactivationtimeinseconds.(b)Sourceradiationpatterninferredfrombackprojectionofsyntheticinfrasoundsignalsusingone arraylocated150kmawayfromtheepicenter,45◦tothenodalplanesinthedirectionindicatedbythearrow.(c)Superpositionofsourceradiationpatternsinferredfrom backprojectionofsyntheticinfrasoundsignalsusingtwoarrayslocated150kmawayfromoriginat45◦tothenodalplanesinthedirectionindicatedbythearrows. piston is offset in time to mimic a radially propagating seismic

wave withamoveoutvelocity of3km/s fromthe simulated epi-centeratthecenteroftheextendedsource.

Syntheticwaveformsarecalculatedforafour-elementarrayin thefar-ﬁeld,150kmawayfromtheepicenterat45◦ tothenodal planestomimictheorientation ofIS51withrespecttothenodal planesoftheHaitiearthquake.Wavefrontparametersareextracted in a beamforming process, as described in Data acquisition and beamformingresults,andthenusedinthebackprojectionprocess. Thepressure-time integral S ineachdetectionwindowisusedto inferthegroundmotionpolarityofeachdetectionpatch.

As before, three distinct patches are resolved (Fig. 6b) with themiddlepatchsmearedacrosscompressionalquadrantsoneand three.Thesource-receiverdistanceinthesyntheticmodelis three-times the width of the extended infrasound source and roughly ﬁve-timesthewidthoftheobservedinfrasoundsourceinthecase oftheHaitiearthquake(Fig.5b).Recallingthatthesizeofthe spa-tialerrorinthecross-trackdirectiongrowswithdistancefromthe receiver, it is clear why these quadrants become merged in the backprojection process. Werepeat the process to simulate wave-forms at a second array at -45◦ to the nodal planes. Since the travel-timeandback azimuthdetected atone array are indepen-dentofthosedetectedatanotherarray,thebackprojectionprocess canbecarriedoutindependentlyperarrayandtheresultscanbe combinedinseveralways(Hernandezetal., 2018). Here,because ofthesimplicityofthesourceandthefactthateachpistonmoves eitherup ordown withnofurtheroscillations, wecan stackand averagethetwosourceradiationmapsandresolvethepolarityin eachofthefourquadrants(Fig.6c).

6. Discussion and conclusions

Following the2010Haiti earthquake,infrasound was detected at an IMS station on Bermuda island. Although a glance at the waveforms in Fig. S1(Supplementary material)mightnot be en-couraging, we wereable toextract coherentsignalsthat arewell withinthelevelofthebackgroundnoiseandassociatethesewith theepicentralregioninHaitiusingwavefrontparametersextracted in the beamforming process (Fig. 3). Two seismo-acoustic cou-plingmechanismsandtheir contributiontoinfrasoundgeneration over thesource regionare discussed: (1)our FFPsimulations

in-dicate that energy from sources in the subsurface, evanescently coupledtotheatmosphereandgeneratedinfrasound(Fig. 4), and (2) extended infrasound source modeling demonstrates that the mechanicalcouplingprocess,thatis,theperturbationofthe acous-ticpressureﬁeldinresponsetoshakingoftheground-atmosphere interface,preservesthegroundmotionpolaritythatisdetermined bythefaultingmechanism(Fig.6).

Our analysis of the atmospheric propagation conditions in-cluded climatological profiles as well as high-resolution ERA5 ECMWFatmosphericspecifications.From thisanalysisitfollowed that propagation of infrasound was most likely facilitated by a thermosphericwaveguide(Fig.4).Thisisfurthersupportedbythe celerityrangeandthecleanwavetrain(single-ductwaveform,no multipathing signature)ofthe observed infrasound detections.In fact,ifinsteadofcalculatingthecelerity withrespecttothe hori-zontaldistancetotheepicenterwe correctforthehorizontal dis-tancetothedetectionpatch,thecelerityvaluesremainbelow0.3 km/s.Infrasound propagationmodelingthroughthethermosphere remains a challenge as wind and temperature specifications are, forthemostpart,limitedtoclimatologicalaverages. Suchmodels donotyet accuratelydescribe thisregion oftheatmospherethat is characterized by large diurnal variations in wind and temper-aturedue to (non-linear interactions between) atmospherictides and(breaking)gravitywaves(Drobetal.,2015;Blancetal.,2017; Assinketal.,2019).Inaddition,non-linearpropagationeffects, ab-sorptionofinfrasoundaswell astheinterplaybetweenthesetwo aspects are significant in the upper atmosphere (Lonzaga et al.,

2015)andnotwellunderstood.

Backprojections of infrasound signals have been shown to be incorrelationwithearthquakegroundmotions priortothiswork butforshorter propagationranges and onlyfor stratospheric in-frasound(Marchettietal.,2016;Walkeretal.,2013;Hernandezet al.,2018). Forthefirsttime, infrasoundthat haspropagated over 1700kmisusedtooutline theregionwhereshakingintensityis sufficientlylargetoleadtodamage(Fig.5a).Inadditiontoshaking onland,ourbackprojectionilluminatessourcesofinfrasoundover theoceanin theepicentralregion.Infrasound fromthesesources evanescently coupled from shaking of the ocean floor. Further-more,groundmotionpolarityinthe epicentralregion isresolved (Fig.5b). Thisprovides an unprecedentedinsight intoearthquake

(10)

Fig. 7. AcousticintensitymappotentialoftheIMS.Thegreengradientshadingindicatescoverageofthecurrentlyinstalledinfrasoundstations(fullcircles).Greenshading correspondstodistanceoutto2000km.Contourlinesspaced30minutesapartindicatetravel-timetotheneareststation,calculatedonthebasisofthermospheric propa-gation.Thelightredshadingindicatesregionsthatwillbecoveredbyplannedstations(emptycircles).Darkredregionscorrespondtolandmassthatwillremainuncovered bytheIMSforthisapplication.

sourcecharacterizationbasedoninfrasounddetectionsat ground-basedstationsinthefar-ﬁeld.

In this study, we demonstrate the potential of remote infra-sounddetectionsformappingtheacousticintensityoveran earth-quakesource region.Wedefer derivationofabsoluteground mo-tions in terms of peak ground velocity or acceleration to future studies that should: (1) account for the propagation term in a moreprecisemanner,and(2)comparemeasuredgroundmotions inthenear-sourceregionwithderivedgroundmotionsfrom back-projections.Suchderivations canthen beincorporatedasanother informationlayerintheShakeMapgenerationprocess.

Infrasound technology is part of the IMS, a global network forthe veriﬁcationoftheComprehensive Nuclear-Test-BanTreaty (CTBT)(Dahlman et al., 2009). Currently, 52out of 60IMS infra-sound arrays are operational andstreaming data in real-time to the International Data Center in Vienna. Considering a 2000 km radiusaroundeachstation,thetotallandmasscoveragebyatleast onestationis90%withanadditional7%thatwillbecoveredonce the network is complete (Fig. 7). Large regions are covered by more than one station.Other infrasound stations,known to pro-vide additional coverage, can be incorporated. These include the USArray-Transportable-Array(TA)stationsequippedwitha micro-barometer alongsideeach seismometersince2007, stationsofthe ARISE(Blanc etal., 2017) (AtmospheredynamicsResearch InfraS-tructureinEurope) initiative,andother nationallyoperated infra-soundstations(Pilgeretal.,2018).

Nippress andGreen (2019) have shownthat infrasound prop-agation ina thermospheric waveguideis effectiveout to a range of 2000kmwithcelerity values around0.28 km/s.The expected infrasound travel-time over 2000 km assuming a thermospheric waveguide is approximately 2 hours under typical propagation conditions.Thismeansthatanacousticintensitymapcanbe pro-ducedfasterthan other methodssuch asdamageanalysison the groundorfromaerialandsatellite imagery.Thiscan bedone for earthquakesalmostanywhereonlandorclosetoshore.The tech-niquespresentedhere,together withthecoverageextent,makeit plausibletouseinfrasoundasaglobalearthquakedisaster mitiga-tiontechniquefortheﬁrsttime.

CRediT authorship contribution statement

Shahar Shani-Kadmiel: Conceptualization, Methodology,Project administration, Software, Visualization, Writing – original draft. Gil Averbuch: Methodology, Software, Writing – review & edit-ing. Pieter

Smets: Methodology,

Software,Visualization,Writing – review & editing. Jelle

Assink: Methodology,

Software, Writing – review & editing. Läslo

Evers: Funding

acquisition, Supervision, Writing–review&editing.

Declaration of competing interest

Theauthorsdeclarethattheyhavenoknowncompeting ﬁnan-cialinterestsorpersonalrelationshipsthatcouldhaveappearedto inﬂuencetheworkreportedinthispaper.

Acknowledgements

The infrasound data used for this articleare provided by the ComprehensiveNuclear-Test-BanTreatyOrganization(CTBTO) and are available to member states. Data can be requested for aca-demicpurposesviathevirtualDataExplorationCenter(vDEC).The CTBTO and InternationalMonitoring System(IMS) station opera-torsarethanked forthehigh-qualitydataandproducts.The high resolutionERA5ClimateReanalysisdataisfreelyavailable viathe CopernicusClimateChangeService(C3S)ClimateDataStore(CDS). Figures in this article (1-6) are made with Matplotlib (Hunter,

2007) and (7) The Generic Mapping Tools (Wessel et al., 2019). Wethankreviewer SusanHoughandan anonymousreviewerfor theirthoroughandconstructivereviewsofthemanuscript.Wealso thank Kees Wapenaar, Auke Barnhoorn, and Thomas Reinsch for theirpre-submissionreviewanddiscussions.

Funding: The contributions by Shahar Shani-Kadmiel, Pieter Smets,andLäsloEversarefundedthroughaVIDIprojectfromThe DutchResearchCouncil(NWO),Project864.14.005.GilAverbuchis fundedthroughtheMarieCurieAction WAVESfromtheEuropean UnionwithinH2020,GrantNumber641943.

(11)

Appendix A. Supplementary material

Supplementarymaterialrelatedtothisarticlecanbefound on-lineathttps://doi.org/10.1016/j.epsl.2021.116795.

References

Arrowsmith,S.J.,Johnson,J.B., Drob,D.P.,Hedlin,M.A.H.,2010.The seismoacous-ticwaveﬁeld:anewparadigminstudyinggeophysicalphenomena.Rev. Geo-phys. 48,23.https://doi.org/10.1029/2010RG000335.

Assink, J.D., Averbuch, G., Shani-Kadmiel, S., Smets,P.S.M., Evers, L.G., 2018. A seismo-acousticanalysisofthe2017NorthKoreannucleartest.Seismol. Res. Lett. 86,2025–2033.https://doi.org/10.1785/0220180137.

Assink, J., Smets,P.,Marcillo, O., Weemstra,C., Lalande, J.-M.,Waxler,R., Evers, L., 2019. Advances ininfrasonic remote sensing methods.In: LePichon, A., Blanc,E.,Hauchecorne,A.(Eds.),InfrasoundMonitoringforAtmosphericStudies. SpringerInternationalPublishing,Cham,pp. 605–632.

Averbuch,G., Assink,J.D.,Evers, L.G.,2020a. Long-rangeatmosphericinfrasound propagation from subsurface sources. J. Acoust. Soc. Am. 147, 1264–1274.

https://doi.org/10.1121/10.0000792.

Averbuch,G., Waxler, R.M.,Smets, P.S.M., Evers,L.G., 2020b. Probabilistic inver-sionfor submergedsourcedepthandstrengthfrominfrasound observations. J.Acoust.Soc.Am. 147,1066–1077.https://doi.org/10.1121/10.0000695. Bilham, R., 2010. Invisible faults under shaky ground. Nat.Geosci. 3, 743–745.

https://doi.org/10.1038/ngeo1000.

Blanc,E.,Ceranna,L.,Hauchecorne,A.,Charlton-Perez,A.,Marchetti,E.,Evers,L.G., Kvaerna,T.,Lastovicka,J.,Eliasson,L.,Crosby,N.B.,Blanc-Benon,P.,LePichon, A.,Brachet,N.,Pilger,C.,Keckhut,P.,Assink,J.D.,Smets,P.S.M.,Lee,C.F.,Kero, J.,Sindelarova,T.,Kämpfer,N.,Rüfenacht,R.,Farges,T.,Millet,C.,Näsholm,S.P., Gibbons,S.J.,Espy,P.J.,Hibbins,R.E.,Heinrich,P.,Ripepe,M.,Khaykin,S.,Mze, N.,Chum,J.,2017.Towardanimprovedrepresentationofmiddleatmospheric dynamicsthankstotheARISEproject.Surv.Geophys.,D2S04.https://doi.org/10. 1007/s10712-017-9444-0.

Bolt,B.A.,1964.Seismic airwavesfromthe great1964Alaskanearthquake. Na-ture 202,1095–1096.https://doi.org/10.1038/2021095a0.

Corbane,C.,Saito,K.,Dell’Oro,L.,Bjorgo,E.,Gill,S.P.,EmmanuelPiard,B.,Huyck, C.K.,Kemper,T.,Lemoine,G.,Spence,R.J.,Shankar,R.,Senegas,O.,Ghesquiere, F., Lallemant, D., Evans, G.B., Gartley, R.A., Toro, J., Ghosh, S., Svekla, W.D., Adams,B.J., Eguchi,R.T.,2011. Acomprehensiveanalysisofbuildingdamage in the 12 January2010 Mw7 Haiti earthquake usinghigh-resolution satel-liteand aerialimagery.Photogramm.Eng.RemoteSens. 77,997–1009.https:// doi.org/10.14358/PERS.77.10.0997.

Dahlman,O.,Mykkeltveit,S.,Haak,H.,2009.NuclearTestBan:ConvertingPolitical VisionstoReality.

Douilly, R.,Haase,J.S.,Ellsworth,W.L.,Bouin,M.-P.,Calais,E.,Symithe,S.J., Arm-bruster, J.G.,deLepinay,B.M., Deschamps,A.,Mildor, S.-L.,Meremonte, M.E., Hough,S.E.,2013.Crustalstructureandfaultgeometryofthe2010Haiti earth-quakefromtemporaryseismometerdeployments.Bull.Seismol.Soc.Am. 103, 2305–2325.https://doi.org/10.1785/0120120303.

Drob,D.P.,Emmert,J.T.,Meriwether,J.W.,Makela,J.J.,Doornbos,E.,Conde,M., Her-nandez,G.,Noto,J.,Zawdie,K.A.,McDonald,S.E.,Huba,J.D.,Klenzing,J.H.,2015. Anupdatetothehorizontalwindmodel(HWM):thequiettimethermosphere: empiricalmodelofthermosphericwinds.EarthSpaceSci. 2,301–319.https:// doi.org/10.1002/2014EA000089.

Earle, P.S., Wald, D.J., Jaiswal, K.S., Allen, T.I., Hearne, M.G., Marano, K.D., Ho-tovec,A.J.,Fee,J.,2009.Promptassessmentofglobalearthquakesforresponse (PAGER):asystemforrapidlydeterminingthe impactofearthquakes world-wide.

Evers,L.G.,Brown,D.,Heaney,K.D.,Assink,J.D.,Smets,P.S.M., Snellen,M., 2014. Evanescent wavecoupling ina geophysicalsystem:airborne acousticsignals from the Mw 8.1 Macquarie Ridge earthquake. Geophys. Res. Lett. https:// doi.org/10.1002/2013GL058801.

Evers, L.G., Haak, H.W., 2007. Infrasonic forerunners: exceptionally fast acoustic phases.Geophys.Res.Lett. 34,1–5.https://doi.org/10.1029/2007GL029353. Francis, G., Stuart, G., Galen, E., Ross, G., Joaquin, B.A., Toro, Paul A., John,

B., Lawrence, C., Ron, E., Michael, E., Shubharoop, G., Matt, H., Zenghui, H., Charles, H., Melisa, H., Walter, S., Enrica, V., Su, W., 2016. Post-disaster building damage assessment using satellite and aerial imagery interpretation, fieldverification and modeling techniques. 2010 Haiti earth-quake - Final report. A World Bank, GFDRR and ImageCat report. http:// learningfromearthquakes.org/2010-01-12-haiti/11-resources/121-post -disaster-building-damage-assessment-using-satellite-and-aerial-imagery-interpretation -field-verification-and-modeling-techniques.(Accessed October2020).

Frankel,A.,Harmsen,S., Mueller,C.,Calais,E.,Haase,J.,2010.Documentationfor initialseismichazardmapsforHaiti.

Gabrielson, T.,2011.Insitu calibrationofatmospheric-infrasoundsensors includ-ingtheeffectsofwind-noise-reductionpipesystems.J.Acoust.Soc.Am. 130, 1154–1163.https://doi.org/10.1121/1.3613925.

Garcia,R.F.,Bruinsma,S.,Lognonné,P.,Doornbos,E.,Cachoux,F.,2013.GOCE:the ﬁrstseismometerinorbitaroundtheEarth.Geophys.Res.Lett. 40,1015–1020.

https://doi.org/10.1002/grl.50205.

Godin,O.A.,2002.Aneffectivequiescentmediumforsoundpropagatingthrough aninhomogeneous,movingﬂuid.J.Acoust.Soc.Am. 112,1269–1275.https:// doi.org/10.1121/1.1504853.

Godin,O.A.,2008.Soundtransmissionthroughwater–airinterfaces:newinsights into an old problem. Contemp. Phys. 49, 105–123. https://doi.org/10.1080/ 00107510802090415.

Godin,O.A.,2011.Low-frequencysoundtransmissionthroughagas–solidinterface. J.Acoust.Soc.Am. 129,EL45–EL51.https://doi.org/10.1121/1.3535578. Green,D.N.,Guilbert,J.,LePichon,A.,Sebe,O.,Bowers,D.,2009.Modelling

ground-to-aircouplingfortheshallowML4.3Folkestone,UnitedKingdom,earthquake of28April2007.Bull.Seismol.Soc.Am. 99,2541–2551.https://doi.org/10.1785/ 0120080236.

Hayes,G.P.,Briggs, R.W.,Sladen,A.,Fielding,E.J.,Prentice, C.,Hudnut,K.,Mann, P.,Taylor,F.W.,Crone,A.J.,Gold,R.,Ito,T.,Simons,M.,2010.Complexrupture duringthe12January2010Haitiearthquake.Nat.Geosci. 3,800–805.https:// doi.org/10.1038/ngeo977.

Hedlin,M.A.H., Raspet,R.,2003.Infrasonic wind-noisereductionbybarriersand spatialﬁlters.J.Acoust.Soc.Am.

Hernandez,B.,LePichon,A.,Vergoz,J.,Herry,P.,Ceranna,L.,Pilger,C.,Marchetti, E.,Ripepe,M.,Bossu,R.,2018.Estimatingtheground-motiondistributionofthe 2016Mw6.2Amatrice,Italy,earthquakeusingremoteinfrasoundobservations. Seismol.Res.Lett.https://doi.org/10.1785/0220180103.

Hough,S.E.,Altidor,J.R.,Anglade,D.,Given,D.,Janvier,M.G.,Maharrey,J.Z., Mere-monte,M.,Mildor,B.S.-L.,Prepetit,C.,Yong,A.,2010.Localizeddamagecaused bytopographic ampliﬁcation during the 2010 M 7.0 Haiti earthquake. Nat. Geosci. 3,778–782.https://doi.org/10.1038/ngeo988.

Hough,S.E.,Taniguchi,T.,Altidor,J.-R.,2012.Estimationofpeakgroundacceleration fromhorizontalrigidbodydisplacement:acasestudyinPort-au-Prince,Haiti. Bull.Seismol.Soc.Am. 102,2704–2713.https://doi.org/10.1785/0120120047. Hunter,J.D.,2007. Matplotlib: a2D graphicsenvironment. Comput. Sci.Eng. 9,

90–95.https://doi.org/10.1109/MCSE.2007.55.

Jaiswal,K.,Wald,D.,Earle,P.,Porter,K.,Hearne,M.,2011.Earthquakecasualty mod-elswithin the USGSprompt assessmentof globalearthquakesfor response (PAGER)system.In:Spence,R.,So,E.,Scawthorn,C.(Eds.),HumanCasualties inEarthquakes. In:AdvancesinNaturaland TechnologicalHazardsResearch, vol. 29.Springer,Dordrecht,pp. 83–94.

LePichon,A.,Mialle,P.,Guilbert,J.,Vergoz,J.,2006.Multistationinfrasonic observa-tionsoftheChileanearthquakeof2005June13.Geophys.J.Int. 167,838–844.

https://doi.org/10.1111/j.1365-246X.2006.03190.x.

Lonzaga,J.B.,Waxler,R.M.,Assink,J.D.,Talmadge,C.L.,2015.Modellingwaveformsof infrasoundarrivalsfromimpulsivesourcesusingweaklynon-linearraytheory. Geophys.J.Int. 200,1347–1361.https://doi.org/10.1093/gji/ggu479.

Marchetti,E.,Lacanna,G.,LePichon,A.,Piccinini,D.,Ripepe,M.,2016.Evidenceof largeinfrasonicradiationinducedbyearthquakeinteractionwithalluvial sedi-ments.Seismol.Res.Lett. 87,678–684.https://doi.org/10.1785/0220150223.

Melton,B.S.,Bailey,L.F.,1957.Multiplesignalcorrelators.Geophysics 22,565–588.

MercierdeLépinay,B.,Deschamps,A.,Klingelhoefer,F.,Mazabraud,Y.,Delouis,B., Clouard,V.,Hello,Y.,Crozon,J.,Marcaillou,B.,Graindorge,D.,Vallée,M., Per-rot,J.,Bouin,M.-P.,Saurel,J.-M.,Charvis,P.,St-Louis,M.,2011.The2010Haiti earthquake:acomplexfaultpatternconstrainedbyseismologicandtectonic ob-servations.Geophys.Res.Lett. 38.https://doi.org/10.1029/2011GL049799. Mutschlecner,J.P.,Whitaker,R.W.,2005.Infrasoundfromearthquakes.J.Geophys.

Res.,Atmos. 110.https://doi.org/10.1029/2004JD005067.

Nettles,M., Hjörleifsdóttir, V.,2010.Earthquake sourceparametersfor the2010 JanuaryHaiti main shock and aftershocksequence: Haiti earthquake source parameters.Geophys.J.Int. 183,375–380.https://doi.org/10.1111/j.1365-246X. 2010.04732.x.

Nippress,A.,Green,D.,2019.Globalempiricalmodelsforinfrasonicsignalcelerity, backazimuthanddurationfromgroundtruthdata.In:AGUFallMeeting.San Francisco,California,6- 13December.

Picone,J.M.,Hedin,A.E.,Drob,D.P.,Aikin,A.C.,2002.NRLMSISE-00empiricalmodel oftheatmosphere:statisticalcomparisonsandscientiﬁcissues:TECHNIQUES.J. Geophys.Res. 107,SIA15-1.https://doi.org/10.1029/2002JA009430.

Pilger,C.,Ceranna,L.,Ross,J.O.,Vergoz,J.,LePichon,A.,Brachet,N.,Blanc,E.,Kero, J.,Liszka,L.,Gibbons, S.,Kvaerna,T.,Näsholm, S.P.,Marchetti,E.,Ripepe,M., Smets,P.,Evers,L.,Ghica,D.,Ionescu,C.,Sindelarova,T.,BenHorin,Y.,Mialle,P., 2018.TheEuropeaninfrasoundbulletin.PureAppl.Geophys. 175,3619–3638.

https://doi.org/10.1007/s00024-018-1900-3.

Possee,D.,Keir,D.,Harmon,N.,Rychert,C.,Rolandone,F.,Leroy,S.,Corbeau,J., Stu-art,G.,Calais,E.,Illsley-Kemp,F.,Boisson,D.,Momplaisir,R.,Prépetit,C.,2019. ThetectonicsandactivefaultingofHaitifromseismicityandtomography. Tec-tonics 38,1138–1155.https://doi.org/10.1029/2018TC005364.

Shani-Kadmiel, S., Assink, J.D., Smets, P.S.M., Evers, L.G., 2018. Seismoacoustic coupledsignals from earthquakesinCentral Italy:epicentral and secondary sourcesofinfrasound.Geophys.Res.Lett. 45,427–435.https://doi.org/10.1002/ 2017GL076125.

(12)

Smets, P., 2018. Infrasound and the dynamical stratosphere: a new application for operational weather and climateprediction.https://doi.org/10.4233/uuid: 517f8597-9c24-4d01-83ed-0f430353e905.

Sutherland,L.C.,Bass,H.E.,2004.Atmosphericabsorptionintheatmosphereupto 160km.J.Acoust.Soc.Am. 115,1012–1032.https://doi.org/10.1121/1.1631937. Walker,K.T.,LePichon,A.,Kim,T.S.,deGroot-Hedlin,C.,Che,I.Y.,Garcés,M.,2013. Ananalysisofgroundshakingandtransmissionlossfrominfrasoundgenerated bythe2011Tohokuearthquake.J.Geophys.Res.,Atmos. 118,831–851.https:// doi.org/10.1002/2013JD020187.

Wapenaar,C.,Berkhout,A.(Eds.),1989.ElasticWaveFieldExtrapolation,vol.2. El-sevier.

Waxler,R.,Assink,J.,2019.Propagationmodelingthroughrealisticatmosphereand benchmarking.In:LePichon,A.,Blanc,E.,Hauchecorne,A.(Eds.),Infrasound Monitoringfor AtmosphericStudies.Springer InternationalPublishing, Cham, pp. 509–549.

Wessel, P., Luis, J.F., Uieda, L., Scharroo, R., Wobbe, F., Smith, W.H.F., Tian, D., 2019.The genericmappingtools version6.Geochem. Geophys.Geosyst. 20, 5556–5564.https://doi.org/10.1029/2019GC008515.