Evaluating the temporal capability of empirical climatologies for rapid long-range volcanic infrasound propagation estimates using a multidecadal data set of persistent Vanuatu volcanic eruptions

SUMMARY

Powerful infrasound (acoustic waves $<$20 Hz) can be produced by explosive volcanic eruptions. The long-range propagation capability, over hundreds to thousands of kilometers, of atmospheric infrasound motivates the development of regional or even global scale volcano-infrasound monitoring systems. Infrasound propagation paths are subject to spatiotemporal atmospheric dynamics, which lead to deviations in the direction-of-arrival (back-azimuth) observed at sensor arrays and contribute to source location uncertainty. Here, we further investigate the utility of empirical climatologies combined with 3-D ray tracing for providing first-order estimates of infrasound propagation paths and back-azimuth deviation corrections. The intended application is in scenarios requiring rapid or pre-computed infrasound propagation calculations, such as for a volcano-infrasound monitoring system. Empirical climatologies are global observationally based function fitting models of the atmosphere, representing robust predictors of the bulk diurnal to seasonal atmospheric variability. Infrasound propagation characteristics have previously been shown to have strong seasonal and diurnal components. At the International Monitoring System infrasound station IS22, New Caledonia, quasi-continuous multiyear infrasound array detections show oscillating azimuthal variations for arrivals from volcanoes in Vanuatu, including Yasur (⁠$\sim$400 km range), Ambrym (⁠$\sim$670 km range) and Lopevi (⁠$\sim$650 km range). We perform 3-D ray tracing to model infrasound propagation from the Ambrym and Yasur volcano locations to IS22 every six hours (00:00, 06:00, 12:00 and 18:00 UTC) for every day of 2004 and 2019 for Ambrym and Yasur, respectively and evaluate the results as compared to the multiyear observations. We assess a variety of models and parametrizations, including both empirical climatologies and hybrid descriptions; range-independent and range-dependent atmospheric discretizations; and unperturbed and perturbed range-independent empirical climatologies. The hybrid atmospheric descriptions are composed of fifth generation reanalysis descriptions (ERA 5) from the European Centre for Medium-Range Weather Forecasts below $\sim$80 km altitude combined with empirical climatologies above. We propose and employ simple parametric perturbations to the empirical climatologies, which are designed to enhance the stratospheric duct and compensate for missing gravity wave perturbations not included in the climatologies, and thereby better match observations. We build year-long back-azimuth deviation interpolations from the simulations and compare them with three different multiyear array detection data sets from IS22 covering from 2003 up to 2022. Through a systematic comparison, we find that the range-independent empirical climatologies can capture bulk azimuth deviation variability and could thus be useful for rapid infrasound propagation calculation scenarios, particularly during favourable sustained propagation ducting conditions. We show that the hybrid models better describe infrasound propagation during periods of weak stratospheric ducting and during transient atmospheric changes such as stratospheric wind reversals. Overall, our results support the notion that climatologies, if perturbed to compensate for missing gravity wave structure, can improve rapid low-latency and pre-computed infrasound source discrimination and location procedures.

1 INTRODUCTION

Explosive volcanic eruptions are formidable sources of infrasound (acoustic waves $<$20 Hz), perhaps even the most powerful kind of infrasound recorded on Earth (Matoza et al. 2022b; Vergoz et al. 2022). Infrasound from a typical Volcanic Explosivity Index (VEI) 4 eruption (e.g. Matoza et al. 2011a, 2018) can propagate hundreds to thousands of kilometers due to low intrinsic attenuation (Sutherland & Bass 2004) and the formation of atmospheric ducts. These waveguides result from the combination of vertical gradients in horizontal winds and temperature, primarily at stratospheric and thermospheric altitudes (e.g. Le Pichon et al. 2009b; Dabrowa et al. 2011; Matoza et al. 2011a, 2018, 2019, 2022b; Waxler & Assink 2019; Vergoz et al. 2022;).

The International Monitoring System (IMS) of the Comprehensive Nuclear-Test-Ban Treaty Organization (CTBTO) includes a global infrasound network that currently has 53 certified of 60 planned infrasound arrays (Christie & Campus 2010; Marty 2019). The IMS, designed to detect explosions with a yield greater than 1 kT TNT equivalent anywhere on Earth (Le Pichon et al. 2009a; Green & Bowers 2010; Marty 2019), has recorded numerous eruptions (Campus & Christie 2009; Dabrowa et al. 2011; Matoza et al. 2019). For example, 68 of 110 volcanic events worldwide (⁠$\sim$61 per cent) were detected by at least one infrasound station from 2002 to 2009 (Dabrowa et al. 2011), while up to 41 IMS stations were recording (Hedlin et al. 2002; Christie & Campus 2010; Hupe et al. 2022).

Previous studies have shown that we can use several infrasound arrays to detect, locate and track the chronology of distant volcanic eruptions, and provide eruption source parameters for ash transport and dispersal models (e.g. Matoza et al. 2007, 2011a, b, 2017, 2018; Garcés et al. 2008; Fee et al. 2010; Dabrowa et al. 2011; Green et al. 2013; Marchetti et al. 2013; Caudron et al. 2015; Taisne et al. 2019; Perttu et al. 2020). Although most practical algorithms require detection by at least two or three infrasound arrays (stations) to locate the source (e.g. Evers & Haak 2005; Arrowsmith et al. 2008; Modrak et al. 2010; Morton & Arrowsmith 2014; Matoza et al. 2011a, 2017), single-station detection (i.e. detection by one array) of volcanic infrasound signals can still provide valuable information on eruption occurrence (e.g. Fee et al. 2010; De Angelis et al. 2012; Marchetti et al. 2019; Matoza et al. 2019; Ortiz et al. 2020; De Negri et al. 2022; Morelli et al. 2022; Gheri et al. 2023). Long-range infrasound also provides observational constraints on the activity of very remote volcanoes, such as those in the South Sandwich Islands, that can be challenging to observe with satellite remote sensing and for which local instrumentation is difficult to maintain (De Negri et al. 2022).

The dynamic nature of the atmosphere presents a major challenge for infrasound monitoring. Spatiotemporal variability at temporal scales including hourly, daily and seasonal, and at spatial scales from the microscale to synoptic scale, is partially captured by various different atmospheric specification products (Garcés et al. 1998; Drob et al. 2010; Smets et al. 2016; Schwaiger et al. 2019). One of the most important atmospheric phenomena influencing infrasound propagation are the vector wind fields, which can create infrasound waveguides (along-path winds), but also impose cross-path winds (cross-winds) leading to deviations from the great-circle path connecting the true source and receiver locations (back-azimuth deviations), and upwinds resulting in shadow zones (e.g. Evers & Haak 2005; Matoza et al. 2011a; Green et al. 2012; De Negri & Matoza 2023). The back-azimuth deviations can be up to $\sim {10}^{\circ }$ (e.g. Le Pichon et al. 2005b; Antier et al. 2007; Assink et al. 2014a; Morelli et al. 2022), resulting in source mislocations of up to hundreds of kilometers at remote distances (⁠$>$250 km) if using cross-bearing (i.e. triangulation) methods (Matoza et al. 2011b, 2017, 2018; De Negri & Matoza 2023).

Morton & Arrowsmith (2014) showed that by combining 3-D ray tracing techniques with realistic atmospheric profiles, lookup tables for celerity (ratio of range over total traveltime) and back-azimuth priors can be produced to describe global seasonal effects that could improve infrasound detection, location and association. A similar approach was used by Drob et al. (2010), who proposed combining the $\tau -p$ equations  (Garcés et al. 1998) and hybrid ground-to-space (G2S) atmospheric descriptions of Drob et al. (2003), to account for hourly, daily and seasonal changes for real time detection, location and characterization automated algorithms at the International Data Center of the CTBTO.

Building on the previous ideas, Matoza et al. (2018) and De Negri & Matoza (2023) showed that the mislocation of known volcanic sources could be reduced by an order of magnitude by applying corrections to observed back-azimuths before locating the source with a brute-force, grid-search and cross-bearings method (Matoza et al. 2017). De Negri & Matoza (2023) compared the results of using different atmospheric models in this approach, and found that statistical models, also known as empirical climatologies [i.e. the Horizontal Wind Model HWM14 (Drob et al. 2015) and the MSIS2.0 model (Emmert et al. 2020)], can be used to generate back-azimuth corrections that perform comparably well to those calculated using the more advanced and accurate hybrid models [e.g. the European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis v5 (ERA 5) for the first $\sim$80 km of the atmosphere merged with the empirical climatologies for higher altitudes]. However, their study was limited to two eruption case studies lasting only several days, a small time window embedded in much larger temporal changes.

It is clear that, in terms of temporal and spatial accuracy, the best available descriptors of the lower atmosphere (⁠$<$80 km) are the numerical weather prediction (NWP) models (e.g. ECMWF). This is the main reason why the use of empirical climatologies is usually overlooked for infrasound propagation. However, there is an unavoidable trade-off between accuracy and computation resources that scales substantially with the model size. Thus, we usually need to choose either a static atmosphere for a large model or a small model with better accuracy but much more demanding calculations (and a good internet connection). We argue that there is a middle ground where the empirical climatologies and spherical range-independent ray tracing can play a crucial role for first-order global characterization of infrasound propagation that has not been exploited for volcanic near-real time location and characterization. At present, global infrasound propagation results remain accessible only as reference summaries and there are no easy operational tables that could be directly implemented at volcanic observatories (e.g. lookup tables). Additionally, infrasound propagation can involve the thermosphere, which is not included in current NWP products, and a hybrid model fusing the NWP in the lower atmosphere with an empirical climatology for the upper atmosphere is required, but at greater computational complexity (Drob et al. 2003; Schwaiger et al. 2019). Infrasound arrivals that refract back to the surface from thermospheric heights (⁠$\sim$110 km) have been observed (e.g. Le Pichon et al. 2005b; Green et al. 2011; De Negri et al. 2022) and are commonly modelled with ray tracing techniques. The empirical climatologies map altitudes up to the exobase (⁠$\sim$600 km), allowing thermospheric returns to be simulated within the propagation models.

The methodology introduced by De Negri & Matoza (2023) uses readily accessible numerical tools and offer a reproducible workflow for calculation of lookup back-azimuth tables for arbitrary source–receiver locations, considering infrasound arrivals up to thermospheric altitudes. The rapid, first-order approach also reduces the computation burden that more accurate models require. Here, we evaluate the accuracy of their methodology with observations from the IMS infrasound array IS22, New Caledonia, which has recorded persistent infrasound from volcanoes in the Vanuatu volcanic archipelago (⁠$\sim$400–700 km) during the last 20 yr (2003–2022) (e.g. Le Pichon et al. 2005a, 2006; Antier et al. 2007; Assink et al. 2014a; Morelli et al. 2022). We calculate year-long infrasound back-azimuth estimates for arrivals at IS22 using both empirical and hybrid atmospheric specifications, and compare these results to the directions of Yasur and Ambrym volcanoes, which have exhibited persistent activity over the past 20 yr of IS22 recording. We show that for long-range propagation models with year-long timescales, the simpler, first-order approach based on the empirical climatologies and range-independent (spherically layered) atmosphere provides comparable results with the hybrid models that use ECMWF specifications for the lower atmosphere. As such, the calculated lookup tables can be robust descriptors of seasonal changes and provide useful information for global monitoring and characterization of infrasound sources. Operationally, the first-order models could then be improved upon when needed by more accurate descriptions (NWP), or when resources allow (computation power, time, etc.).

In Section 2, we introduce the area of study, highlight the most infrasonically active volcanoes in the region, and present the data sets we use to compare with our models. In Section 3, we introduce our methodology to estimate year-long infrasound back-azimuth deviations, and the different atmospheric and ray tracing propagation schemes we consider. In Section 4, we summarize the annual back-azimuth prediction results, compare them with the multiyear data set, and measure the accuracy of each model on a year-long time window. In Section 5, we discuss the results and the limitations of our approach, and pinpoint some interesting results. Finally, in Section 6 we review our findings and suggest the next steps forward.

2 THE VANUATU ARCHIPELAGO VOLCANIC ARC

Along the South Vanuatu Trench in the South Pacific Ocean, more than 80 islands stretch northeastward to form the Vanuatu archipelago volcanic arc (Fig. 1; Eissen et al. 1991; Simkin & Siebert 1994). The active volcanoes of this archipelago have a rare diversity of eruptive styles which range from effusive to highly explosive, hence becoming a valuable laboratory for volcanologists to better understand eruption dynamics (Beaumais et al. 2016; Coppola et al. 2016; Meier et al. 2016; Vergniolle & Métrich 2016).

As shown in Coppola et al. (2016), from 2000 to 2015 Yasur (Tanna island), Lopevi (Lopevi island), Ambrym (Ambrym island), Lombenben (Ambae island, herein referred as Ambae) and Mont Garet (Gaua island, herein referred as Gaua; see also Park et al. 2021) volcanoes were the most active volcanoes of Vanuatu archipelago. Although most local multiparameter studies in the region target Yasur due to its essentially continuous explosive activity (Marchetti et al. 2013; Meier et al. 2016; Spina et al. 2016; Maher et al. 2022; ; Matoza et al. 2022a).

2.1 IS22: multidecadal long-range volcanic infrasound records

At $\sim$400 km southwest of the Vanuatu archipelago, the IMS infrasound station IS22 (Grande Terre island, New Caledonia) has been recording signals from Vanuatu volcanoes since 2003 (Fig. 1). The station IS22 (or I22FR) is composed of four MB2000 microbarometers separated by $\sim$1 km (Fig. 1, bottom right panel). Each microbarometer has a flat frequency response (within 3 dB) between 0.01 and 8 Hz, with noise levels less than 2 mPa between 0.02 and 4 Hz, and a sampling frequency of 20 Hz (Marcillo et al. 2012).

Coherent infrasound signals propagating across IS22 can be extracted from non-coherent infrasound with the progressive multichannel correlation (PMCC) time-domain method (Cansi 1995). The PMCC method can obtain the coherent signal’s mean frequency, trace velocity, amplitude and azimuth remarkably well under low signal-to-noise ratios (⁠$<$0.5) (Cansi 1995; Matoza et al. 2013). Once a passing wave is identified, the PMCC groups coherent waveforms with similar plane-wavefront properties across sensors into what we call a detection (Cansi & Klinger 1997).

Among many natural and human-made sources detected by IS22 (e.g. microbaroms, surf, mining explosions, etc.; see Fig. 2), the persistent eruptive activity of Yasur volcano (19.532${}^{\circ }$S, 169.447${}^{\circ }$E; 361 m summit; distance $\sim$400 km), and the combined eruptive sequences of Ambrym volcano (16.25${}^{\circ }$S, 168.12${}^{\circ }$E; 1334 m summit; distance $\sim$700 km), and Lopevi volcano (16.507${}^{\circ }$S, 168.346${}^{\circ }$E; 1413 m summit; distance $\sim$700 km) consistently appear during the last $\sim$20 yr (2003–2022; see Fig. 2). These incessant detections have provided an important test-bed to compare and evaluate the accuracy of Numerical Weather Prediction atmospheric models (NRL-G2S and ECMWF), and stochastic models such as the Horizonal Wind Model (HWM) and the Mass Spectrometer Incoherent radar based model (MSIS) (Drob et al. 2003; Le Pichon et al. 2005a, b, 2006; Antier et al. 2007; Assink et al. 2014a). They have also been used to develop passive atmospheric inversion techniques targeting high-altitude winds (Le Pichon et al. 2005a, b; Lalande et al. 2012), and to propose improvements for the IMS infrasound network and test new techniques for remote detection of volcanic eruptions (Le Pichon et al. 2009b, 2010 ; Tailpied et al. 2017; Morelli et al. 2022).

2.2 Open access infrasound data products

In this study, we use open access infrasound data products developed by Hupe et al. (2022), which represent summaries of array processing detection lists rather than the raw lists. These data products started as the result of synthesizing 2003 to 2020 IMS infrasound detections produced from PMCC array processing with a 26-band logarithmic-spaced parametrization, and since then have been updated yearly. The open access data products (see Hupe et al. 2021) use a quality value to gather the most robust PMCC detections and classify them into four specific frequency ranges, which target the most common natural and human-made phenomena: 0.02–0.07 Hz “very low frequency” (e.g. mountain associated wave infrasound), 0.15–0.35 Hz “low-frequency microbarom”, 0.45–0.65 Hz “high-frequency microbarom”, and 1–3 Hz “high-frequency” (e.g. volcanic eruptions and surf).

Here, we will use the “high-frequency” (1–3 Hz) open access products (herein, “PMCC products”) for IS22 covering from years 2003 to 2023 (Fig. 2; see also Fig. S1 in the Supporting Information), as most of the infrasound at IS22 from Lopevi/Ambrym and Yasur volcanoes centres around 1 and 2 Hz, respectively (Le Pichon et al. 2005a; Antier et al. 2007; Assink et al. 2014a; Morelli et al. 2022).

For comparison and validation of these data summary products, we compare with full array detection lists using 15-band (2003–2017) (Matoza et al. 2013; Ceranna et al. 2019) and 26-band (2003–2023) (Hupe et al. 2022) PMCC parametrizations in the Supporting Information (see Figs S4–S8 and S11–S13 in the Supporting Information). Unlike the summaries by Hupe et al. (2022), these detection lists are the direct output of PMCC array processing of the waveforms recorded at the station IS22.

2.3 Chronological and “folded” observations

Volcanic, surf and possibly mining infrasound (“anthropogenic”), are the three main types of “high-frequency” (1–3 Hz) detections at IS22 from 2003 to 2022 (see labels in Fig. 2a). Non-volcanic infrasound falls at much higher back-azimuths than Yasur (⁠$\sim {43}^{\circ }$⁠) and the “Ambrym group” (i.e. Ambrym/Lopevi/Gaua/Ambae; $\sim$5–14${}^{\circ }$⁠). The distinct quasi-continuous infrasound around the back-azimuth of Yasur roughly resembles the shape of “bottom half sinusoids” during summer months (Fig. 2b), with periods of one year. The detections around the back-azimuth of Ambrym appear more sporadically, resembling the slopes of a sinusoid in years with more detections (e.g. 2004–2008 in Fig. 2c).

The periods of volcanic activity catalogued by the Vanuatu Meteorology and Geo-hazards Department (Figs 2b and c in grey vertical lines and shaded areas; also by relative RSAM from a couple of available stations near Yasur volcano in Fig. 2b) serve as valuable complementary observations (e.g. McKee et al. 2021a, b; Matoza et al. 2022a). Their observations indicate that overlapping of infrasound sources among the Ambrym group seems feasible. In particular, during the second half of 2010 (Fig. 2c), a large group of detections appears around the azimuth of Gaua volcano, coincident with a long eruptive period that occurred from September 2009 to July 2010 (Bani et al. 2016).

Sorting the detections by day of year (i.e. “folding” into a year; Figs 2d–f), shows that most of the 2003–2022 IS22 infrasound detections come from the directions of Yasur and the Ambrym group volcanoes. The abundant detections that chronologically resemble “bottom half sinusoid” curves from the back-azimuth of Yasur (Figs 2a and b) are sorted into a “capped arc” now (Fig. 2e). This feature shows up between October and March (Fig. 2e; days $\sim$280 to 365, and 1 to $\sim$100), with minimum values by the end of December to the start of January, and transition periods from March to April (negative to positive azimuth deviations) and October to November (positive to negative). Relatively large groups of detections also appear sporadically between May and August, decoupled from the seasonal deviation trend (i.e. the “capped” part of this arc), with azimuth deviations ranging from approximately $-4$ to $+6^{\circ }$ (Fig. 2e). The “folded” Ambrym group (Figs 2d and f) shows detections resembling the sides of a “main arc” around Ambrym/Lopevi. In this case, the features are less clear, as the numbers of detections are much lower than Yasur’s case. However, Fig. 2(c) does suggest sinusoidal back-azimuth variation during periods of higher number of detections (e.g. 2004–2008). Similarly, we consider it feasible that the traces of a “secondary arc” exist centred near Ambae/Gaua back-azimuths (approximately $-3$ to $-5^{\circ }$ below the main arc). The hypothesis of a secondary arc stems from the fact that, if Ambae and/or Gaua produced detectable infrasound during shorter periods of time (e.g. in the middle of 2010), these detections should be equally influenced by the travel path effects that make Ambrym/Lopevi detections have a sinusoidal back-azimuth deviation through the year. For both cases, most of their detections appear from March to April and August to October, while decoupled groups of detections appear from June to August, similar to the case of Yasur volcano (Fig. 2e).

3 ANNUAL LONG-RANGE PROPAGATION MODELS

3.1 Methodology

We attempt to model the bulk seasonal to daily effects to first-order (e.g. the “capped arc” and “main arc” features; Fig. 2e), which could provide valuable insights and potentially seasonal predictions of the azimuth deviation from each source. To do so, we use the methodology introduced by De Negri & Matoza (2023) to simulate the observed detections at IS22 from the locations of Yasur and Ambrym volcanoes (see the “ARCADE” open-source repository link in Data Availability) in year-long time windows.

The ARCADE procedure (De Negri & Matoza 2023) is designed to automatically look for ground intercepts (arrivals) at pre-defined receiver locations (stations) coming from an arbitrary number of source locations (volcanoes), at specified dates and times. By default, it couples empirical climatologies (HWM14/MSIS2.0) (Drob et al. 2015; Emmert et al. 2020) and models infrasound propagation with 3-D ray tracing (infraGA) (Blom & Waxler 2012, 2017; Blom 2020). The atmospheric models can also be discretized as range-independent (layered spherical) or range-dependent (3-D grid for the area), allowing different levels of detail in the specifications. The modularity of this approach permits using different atmospheric specifications that we will test in terms of accuracy.

3.1.1 Hybrid model

The empirical climatologies are stochastic models based on interpolations of historical data collections, which are improved over time as new global data sets are gathered. This means that they are robust atmospheric descriptions, but smooth out smaller spatiotemporal changes that occur in the lower parts of the atmosphere (troposphere/lower stratosphere), as the quasi-biennal oscillation (QBO; Anstey & Shepherd 2014; Kidston et al. 2015) and sudden stratospheric warming (SSW) events (Drob et al. 2010; Assink et al. 2014a; Smets et al. 2016; Drob 2019). To try to account for this issue, ARCADE (De Negri & Matoza 2023) does include the option of using a hybrid model that merges the ECMWF ERA 5, or ERA 5 (for simplicity) descriptions (horizontal winds, pressure and temperature) below $\sim$78 km altitude, with pure (raw) empirical climatologies at higher altitudes. The ERA 5 reanalysis models, available from 1959 to the latest processed forecast model, describe the lower atmosphere in greater detail, accounting for specific year and date–time variations with a minimum time step of one hour and a minimum spatial step of 31 km (Bell et al. 2021).

The hybrid models will be the only type able to distinguish each year, as the empirical climatologies are interpolations of multiyear data into one characteristic year. These models will therefore be more accurate in representing transient non-seasonal changes, in detriment of representing repetitive (seasonal) phenomena in a multiyear scale.

3.1.2 Perturbed climatologies

One possible approach to account for phenomena not represented in the empirical climatologies is to generate multiple fine-scale variations of the average descriptions (multiple statistical realizations) simulating fine-scale structure or gravity wave effects (e.g. Assink et al. 2014a; Smets et al. 2016). Here, however, we take a simpler parametric approach, incorporating the addition of a Gaussian wind jet (Jones 1986) centred at a specific altitude, which creates perturbations of the range-independent climatological (HWM14/MSIS2.0) descriptions. This simplified type of perturbation is intended to adequately compensate for missing fine-scale structure that is not included in either hybrid atmospheric specifications nor empirical climatologies, without having to resort to a full gravity-wave model parametrization and multiple ensemble realizations. This approach, which was further motivated by De Negri & Matoza (2023), aims to sufficiently “nudge” the profiles to form a duct when the observations imply that the duct existed, yet the unperturbed profiles do not predict arrivals.

3.1.2.1 Static perturbation

Originally, De Negri & Matoza (2023) defined the Gaussian perturbation parameters after manually analysing effective wind profiles for the cases studied. The parameters of altitude, width and amplitude that would increase the horizontal winds were set to avoid unrealistic arrivals. In this type of perturbation, the Gaussian enhancement is fixed at 40 km altitude by default, increasing the along-path winds with an amplitude, $c_0$⁠, that corresponds to a 30 per cent of the maximum along-winds (m/s) within the 30–50 km altitude limits (see eq. 1; De Negri & Matoza 2023).

where $\mathrm{△}h$ is the width (km), and $h_0$ is the central altitude (km) at which the perturbation is applied. That is, the amplitude (⁠$c_0$⁠) is centred at a fixed altitude (⁠$h_0$⁠), with a $\pm$10 km width (⁠$\mathrm{△}h$⁠). This perturbation (⁠$\tilde{c}$⁠) is applied only when positive along-path stratospheric winds are present, defaulting to raw (unperturbed) empirical climatologies when not. However, this choice of perturbation has likely underestimated the number of stratospheric arrivals (De Negri et al. 2022; De Negri & Matoza 2023), while the numerous observations from Yasur have been regarded as dominantly stratospheric (e.g. Le Pichon et al. 2005a; Antier et al. 2007; Assink et al. 2014a).

3.1.2.2 Adaptive perturbation

Attempting to approximate the higher amplitude and smaller oscillation lengths of typical gravity wave perturbations (e.g. Green et al. 2011; Drob et al. 2013; Smets et al. 2014), we tested several reasonable parameter modifications of eq. (1), to examine how the altitudes, widths and amplitudes of perturbations could increase stratospheric arrivals at IS22. After these preliminary experiments, we incorporated a second type of perturbation based on eq. (1). We decided to (1) increase the amplitude perturbation (⁠$c_0$⁠) to a 90 per cent of the maximum along-wind speeds from 30 to 60 km (i.e. up to $\sim$16 m s⁻¹ mean increase during along-wind days), (2) reduce the width (⁠$\mathrm{△}h$⁠) to 1 km and (3) centre the perturbation at a height (⁠$h_0$⁠) that follows the local maximum per day and time (see Fig. S14 and Table S1 in the Supporting Information). These perturbations, despite not being actual gravity-wave perturbations, better resemble the distortion length of one gravity wave oscillation and the maximum perturbation value it could have (e.g. Green et al. 2011; Drob et al. 2013; Smets et al. 2014). Additionally, we apply the adaptive perturbation proportional to the factor R (see eq. 2) to the zonal winds (u), and proportional to the factor $1-R$ to the meridional winds (v), thus guaranteeing a perturbation that follows the actual wind magnitude ratios (Fig. 3).

To summarize, in this study we will use “raw climatologies”, “hybrid models”, “static perturbed climatologies” (default perturbation) and “adaptive perturbed climatologies” (last introduced) to obtain the atmospheric descriptions for each modelled case (Fig. 3).

3.2 Propagation models

After inspecting the multiyear (2003–2022) set of detections (Figs 2a and b; also Figs S4 and S5 in the Supporting Information), we decided to calculate year-long infrasound propagation models corresponding to 2004 for Ambrym volcano (herein, “2004 Ambrym”) and 2019 for Yasur volcano (herein, “2019 Yasur”). These two years have numerous candidate detections for each case, presenting complete data sets to compare with the models.

We cover all the possible combinations with this methodology, calculating six different models per volcano and year: range-independent with raw climatologies, range-independent with static perturbed climatologies, range-independent with adaptive perturbed climatologies, range-independent with hybrid model, range-dependent with raw climatologies and range-dependent with hybrid model.

With each model, an iterative ray tracing search is automatically performed for every day of the year, every six hours [00:00, 06:00, 12:00 and 18:00 Coordinated Universal Time (UTC)]. The source–year–day–time runs try to estimate the back-azimuth deviation using arrivals near enough IS22 (⁠$<$5 km), with turning heights up to $\sim$120 km (thermospheric). We define the back-azimuth deviation as the mean difference between the true back-azimuth and the modelled back-azimuth of the considered arrivals, which mostly come from stratospheric and lower thermospheric ducts (De Negri & Matoza 2023).

Using the year-long sets of back-azimuth deviations for each model, we build year-long 1-D spline interpolations based on the daily means (Fig. 4, see green lines for each model). The interpolated curves will serve to compare models, and more importantly, to calculate the accuracy of each model when compared with the observations.

All models managed to produce a mostly complete set of arrivals to calculate the back-azimuth associated with each source–year–date–time combination. However, the range-dependent, raw climatology model for Ambrym volcano failed to produce ground intercepts (arrivals) near enough IS22 (Fig. 4k). We interpret this issue as the indication of an unstable propagation model, where the ray propagation becomes too sensitive to small changes in the launching parameters. However, further clarifications will be part of future research as we gathered enough results to complete the comparison with the data sets.

4 RESULTS

4.1 Regional winds and detectability

Yasur and Ambrym groups of detections have a marked yearly oscillating back-azimuth behaviour. This phenomenon correlates well with the yearly regional zonal wind patterns, which are clearly captured by climatology-based and hybrid atmospheric models (e.g. Fig. 5 for winds during 2019).

Infrasound propagating from Ambrym and Yasur to IS22 (Fig. 1) will encounter eastward stratospheric (⁠$\sim$30–60 km) zonal winds from April to October and westward winds from November to March (e.g. Figs 5a, c, e and g). Yasur infrasound will find a sustained propagation duct when westward zonal winds occur, whereas the eastward winds create upwind propagation conditions, making them almost undetectable (e.g. Figs 5b, d, f and h). At the same time, sources north from IS22 (i.e. Ambrym/Lopevi/Gaua/Ambae), will have no clear favourable stratospheric ducting conditions. Instead, they will suffer high crosswinds during prevalent eastward/westward zonal winds. However, at thermospheric altitudes (⁠$\sim$110–120 km), there is another theoretical duct that could enhance infrasound propagation to IS22 (Figs 5b, d, f and h). During wind regime transition periods (seasonal equinoxes), characterized by weak zonal winds (e.g. see wind profiles from March–April and September–November in Fig. 5), the minimal crosswinds could enable infrasound refracted at thermospheric heights to be more detectable, coincident with the times when most of the detections from Ambrym/Lopevi azimuths appear (Fig. 2f).

4.2 Observations: interpretations and shortcomings

4.2.1 Seasonal arcs and transient phenomena

The shared arc-like (sinusoidal) behaviour of the modelled detections (Figs 4 and 5) suggests that the candidate detections coming from Yasur and Ambrym back-azimuths (Fig. 2) are mainly due to the effect of the seasonal stratospheric zonal winds.

The near-continuous volcanic explosive signals from Yasur to IS22 shows as “bottom half sinusoids” or “capped arc” features due to low detectability during upwind conditions in winter (eastward winds; e.g. Fig. 5 for 2019). Transient stratospheric westwards ducts during dominant eastwards winds, linked with sudden stratospheric warming events in the region (Assink et al. 2014a) (e.g. Fig. 5b near July and September for 2019), allow detections to be seen for days to weeks in a row (e.g. Figs 4g and i for Ambrym and Figs 4h and j for Yasur near day 170).

Most of the detections of the “Ambrym group” occur during transition zonal wind regimes (March to April and August to October), which means that the bulk of the Ambrym group happens during a weakened or non-existent stratospheric duct. Indeed, the atmospheric effective sound speeds profiles for the region suggest that these candidate detections are enabled by the thermospheric duct (⁠$\sim$110–130 km), whereas the sinusoidal azimuthal deviations are caused by the crosswinds due to the seasonal zonal winds at stratospheric heights (e.g. see Fig. 5). As confirmed by the models, the crosswinds imposed on the propagating signals from Ambrym result in stronger azimuth deviations due to their North-to-South propagation path (Fig. 1). The peak-to-peak asymmetry in the modelled back-azimuths and the observed candidate detections indicates stronger westward (summer) than eastward zonal winds (winter) in the region.

It is likely that the “second arc” present in the folded observations (Figs 2d and f) is an overlap of candidate detections from Gaua and Ambae volcanoes. The models indicate that the daily variations (Fig. 4) are not higher than the back-azimuth differences from Ambrym/Lopevi and Gaua/Ambae (⁠$\sim$7–4${}^{\circ }$⁠). However, further modelling and comparison with other data sets are required to clarify this hypothesis.

4.2.2 Artifacts and frequency bands

At the start of 2018, we observe two groups of detections spanning from $\sim$10 to 45${}^{\circ }$ that appear to be artifacts (Fig. 2a; also Figs 2d–f at the end of January, February and the start of March). These detections are also present in the 26-band PMCC data set during the same time window (Fig. S4, Supporting Information), and were likely produced by temporary station issues, possibly due to a malfunction of one of the sensors. The data set also shows the formation of azimuthal bands of detections at $\sim {1.25}^{\circ }$ steps for the Yasur group of detections (Figs 2b and e) that appear to be artifacts. They are present in the 26-band PMCC data (Figs S4b, c, e and f, Supporting Information) but not the 15-band PMCC data (Figs S5b, c, e and f, Supporting Information) and are likely a processing artifact effect of the 26-band configuration.

4.3 Comparison with year-long and multiyear data sets

We aim to compare the back-azimuth interpolated model curves from year-long infrasound back-azimuth deviation estimations at IS22 (see Fig. 4), with the high-frequency (1–3 Hz) infrasound products (Hupe et al. 2021) for the years 2004, 2019, and the whole stacked (folded) multiyear time window (2003–2022). However, the ever-present non-volcanic overlapping sources (clutter), and propagation scatter effects make the comparison of a single interpolated curve difficult. Therefore, before calculating the error estimation of each model, we need to apply reduction techniques to the observations to represent the average behavior of the candidate (volcanic) detections.

4.3.1 Data reduction

To reduce clutter, we first select the high-frequency products (Hupe et al. 2021) that have a quality parameter $Q>0.4$ (Hupe et al. 2022). By default the products are composed of groups of detections with a size above 40 (Hupe et al. 2022), while for the other two data sets we set a minimum size threshold of 10 for each group of detections (families). Then, we consider only the detections with back-azimuths within $\pm {10}^{\circ }$ around the true direction of each volcano. We calculate the mean of the selected detections within the interval bounded by the 10th and 90th percentiles for the back-azimuth distribution in 7-d consecutive time windows, producing the values that we use to compare with the modelled interpolations. The 7-d means will be the reference data to compare with the model interpolations, allowing us to measure the performance of the models.

4.3.2 2004 Ambrym and 2019 Yasur comparison

The interpolated curves, centred at each back-azimuth, generally follow the most numerous candidate detection groups for both 2004 Ambrym and 2019 Yasur cases (Figs 6a and c). The 7-d 90 per cent confidence reduced observations are a much more clear indication of the fitting of each model to the observations (Figs 6b and d).

For 2004 Ambrym case, the reduced observations have a large negative deviation (⁠$\sim$$-{7.5}^{\circ }$⁠; Fig. 6b) at the start of the year, followed by a linear increase from days $\sim$40 to 120, and reaching values of $\sim {2.5}^{\circ }$ towards the middle of the year (day $\sim$150). The second half of the year (day $\sim$200–300) has fewer detections, which follow a negative slope in a more scattered manner than the first half. The match between the observations and the modelled results is good from days $\sim$40–120 for all models, both following a similar rate of back-azimuth change in time. However, the observed deviation values are underestimated by all models in the beginning of this period (days $\sim$40–120), while the hybrid model interpolations (blue lines) are the only curves that follow the higher slope of the transition from negative to positive deviations from days $\sim$100 to 120. As previously hypothesized, the back-azimuth variability of the detections is better portrayed by the hybrid predictions, which are able to follow detections with negative deviations during periods of positive deviations in the empirically constrained climatologies (e.g. day $\sim$180). The frequencies of the reduced observations remain near 1 Hz in most cases.

For 2019 Yasur case, the reduced observations show a good agreement with all the interpolated curves when the westward stratospheric duct is present (days 1 to $\sim$100 and $\sim$300 to 365; Fig. 6d). The hybrid models again seem to better follow groups of detections that are outliers from the seasonal predictions (e.g. day of year $\sim$180, with $\sim$1.8 Hz, and azimuth deviation of $\sim 0^{\circ }$⁠). Interestingly, from days $\sim$250–320, the back-azimuths are under-predicted by all the models. The reduced observations show distinct frequency changes through the year, periods of higher frequencies (⁠$\sim$1.8 Hz) happen at the beginning/end of the year and in the middle of the year, while lower frequencies (⁠$\sim$1.2 Hz) appear between these periods. Following the hybrid model results (Figs 4h and j), we interpret this effect as the relation between turning height and transmission loss. Favourable summer stratospheric ducts and transient stratospheric along-path winds during winter, which are likely driven by the QBO effects at these southern tropical latitudes (Anstey & Shepherd 2014; Assink et al. 2014b), allow higher frequencies to reach IS22 (shorter travel path). On the other hand, weak stratospheric ducting conditions force most ray paths to rely on the higher altitude thermospheric duct (⁠$\sim$110–120 km), lowering the mean frequencies of the arrivals at IS22.

4.3.3 Ambrym and Yasur multiyear (2003–2022) comparison

The models show general agreement with the multiyear observations again from days 1 to $\sim$100, and $\sim$250 to 365 (Fig. 7). However, the linear-like clusters of detections (i.e. frequency bands; see Section 4.2.2) clearly appear in steps of $\sim {1.25}^{\circ }$ around the back-azimuth of Yasur, skewing the overall statistics for this group (Fig. 7c). Fortunately, the model comparison with the 15-band 2003–2017 data set (Fig. S8, Supporting Information) shows that the most important characteristics of the back-azimuth changes are still well represented.

The hybrid models, despite their specific time dependency, better represent the slope of the transitions from negative to positive deviations (e.g. day 100–120).

Between day $\sim$120 and $\sim$250, the possible Gaua/Ambae group of detections (see Section 2.3) appears under Ambrym’s modelled deviations (Fig. 7a) from $\sim$−2.5${}^{\circ }$ to $-{7.5}^{\circ }$⁠. This group falls outside the interval bounded by the 10th and 90th percentile of the back-azimuth distribution in the 7-d rolling window, therefore it is erased from the reduced representation (Fig. 7b). Their frequency content (⁠$\sim$1.6 Hz) is slightly higher than most of the detections (⁠$\sim$1–1.4 Hz), with only one notable group that has clearly higher frequencies than the rest (⁠$\sim$1.8–2.2 Hz, near day 200). This suggests mixing of detections from transient stratospheric ducts with detections from thermospheric ducts (Fig. 7b). Additionally, if these detections are indeed from Gaua/Ambae, this could indicate a higher characteristic source frequency than Ambrym (⁠$\sim$1 Hz).

More outlier detections build up in the multiyear view for Yasur’s comparison (Fig. 7d; days $\sim$100 to $\sim$250). These detections show smaller deviations than predicted by the climatology-based models, only partially matched by the hybrid curves. During days $\sim$280 to $\sim$320 groups of less numerous, lower frequency detections (⁠$\sim$1.2–1.4 Hz) transition towards more numerous, higher-frequency groups (⁠$\sim$1.8–2.2 Hz). However, this frequency effect is clearly visible in the 15-band 2003–2017 stacked analogous comparison (see Fig. S8 in the Supporting Information), resembling the case of 2019 for Yasur (Fig. 6d).

4.3.4 Accuracy of the models

Up to this point, we have shown that the models are qualitatively comparable with the reduced data sets (Figs 6 and 7). Now, we formulate a simple monthly and yearly quantitative evaluation of the model-data discrepancy.

Following the 7-d time window of each reduced data point (see Section 4.3.1), we first calculate the mean of the modelled back-azimuth deviations in 7-d time intervals. Then, we calculate the monthly discrepancy using the Euclidean distances between the two sets of 7-d values (model versus observations), building means by month (Fig. 8). These values preserve the sign that corresponds to the difference between the observations and modelled back-azimuth deviations, reflecting the expected azimuthal direction of the discrepancy. Additionally, we calculate the year-long differences by averaging the magnitudes of their monthly discrepancies. These will be the final values representing the overall accuracy of each model with respect to the reduced observations (Table 1; also Figs S9 to S13 in the Supporting Information).

Ambrym monthly differences for 2004 (Fig. 8a) oscillate between a negative maximum of $\sim$−3${}^{\circ }$ to −4${}^{\circ }$ around December/January for all models, and a maximum positive peak of $\sim {1.5}^{\circ }$ around August for all but the hybrid models, which have a less-negative to maximum positive peak on April (⁠$\sim$−0.5${}^{\circ }$ and $\sim 1^{\circ }$ for the range-dependent and range-independent cases, respectively). The best agreement between models and observations occurs around November and March, with differences smaller than 1${}^{\circ }$⁠. This is coincident with the times when most of the detections appear, along the “main arc sides” pointed out in Fig. 2(f). The negative differences in the middle of the year are clearly caused by either clutter infrasound or other volcanic sources such as Gaua/Ambae. Despite being able to predict observations produced by transient ducts (e.g. QBO/SSW), the hybrid models accumulate the highest differences during this time.

The range-dependent, hybrid model shows largest difference, with a yearly mean of 1.59$\pm {1.07}^{\circ }$ (Table 1; also Fig. S9 in the Supporting Information). On the other hand, the range-independent, hybrid model has the smallest difference, with a yearly mean of 1.07$\pm {0.99}^{\circ }$⁠. Among the empirical climatology based models, the range-independent, static perturbed model follows with a 1.12 $\pm {0.82}^{\circ }$ mean difference.

The Yasur monthly differences for 2019 (Fig. 8c) are mostly below 1${}^{\circ }$ magnitude and negative. This is reflected in the yearly means, which have values from 0.43 to 0.58${}^{\circ }$⁠, and spreads from $\sim$0.3 to 0.4${}^{\circ }$⁠. This shows that the models have a better agreement with the observations in this case. However, in June and July there are disagreements between the climatology-based and hybrid models that differ in sign and have the highest magnitudes, marking a period within the “capped” part of the arc hypothesized in Fig. 2(e). These again are discrepancies caused by periods of lower numbers of detections likely driven by clutter infrasound, or transient ducts. In this case, the hybrid models seem to correctly predict a group of detections in July that could be associated with the effects of QBO westward equatorial winds in the second half of 2019 (see Figs 5a and b), but further research is required to clarify this assertion.

The model with the smallest yearly mean difference in this case is the range-independent, adaptive perturbed climatologies (0.43 $\pm {0.29}^{\circ }$⁠). The range-dependent, hybrid model has the largest relative difference (0.58 $\pm {0.39}^{\circ }$⁠). This mismatch increase can be mainly attributed to a peak of $\sim {1.5}^{\circ }$ difference in June, which is an overshoot effect of the interpolation.

Overall, the considerable spread of the mean euclidean differences for both cases shows that the models, rather than better or worse in absolute terms, are comparable. This indicates that the empirical climatologies are valid to understand the specific year-long, first-order behaviour of the infrasound observations for 2004 and 2019.

The comparison with the 2003–2022 stacked data set shows an increased yearly mean difference for all models (Table 1; Figs 8b and d; also Fig. S10 in the Supporting Information). As the multiyear comparison is suited for the empirical climatologies, the hybrid descriptions will be an upper limit reference of the monthly and year-long differences.

The monthly differences show that, when stacking the data, the overall discrepancies are negative and peak around August and May for Ambrym and Yasur cases, respectively (see Figs 8b and d). As previously shown, these negative peaks occur during periods of low numbers of candidate detections from either Ambrym or Yasur, caused by unfavourable ducting conditions, and non-seasonal detections that can be associated with transient favourable ducts (QBO or SSW events) and clutter infrasound. February to April, and October to November are overall the months with smallest differences for Ambrym, while September to February (cyclically) are the ones for Yasur. These months align with the arc-like features pointed out in Figs 2(e) and (f), during times of favourable ducting conditions.

Both multiyear mean differences show the same relative results (Table 1; also Fig. S10 in the Supporting Information). The range-independent, adaptive perturbed climatologies provide the minimum differences: 0.90 $\pm {0.66}^{\circ }$ for Ambrym while 0.86 $\pm {0.79}^{\circ }$ for Yasur. The range-dependent, hybrid descriptions go up to 1.77 $\pm {1.43}^{\circ }$ for Ambrym and 1.21 $\pm {0.94}^{\circ }$ for Yasur.

5 DISCUSSION

The accuracy measurements we present in Table 1 are an attempt to compare our rather simple models with a large amount of observations which have a much more complex behaviour. The mean differences represent a measure of agreement between the year-long back-azimuth deviation interpolated curves and the reduced data, for a specific year (2004 or 2019) or as a multiyear stack (2003–2022). We remark that these values serve to compare between the models we tested for this specific area and the times considered, but do not provide a general interpretation of the performance of each model scheme for the globe unless future experiments show similar results. Keeping this in mind, we see that all the models perform similarly well despite the considerable differences among the used atmospheric descriptions (i.e. climatologies versus hybrid; raw climatologies versus perturbed climatologies) and ray tracing model types (range-independent versus range-dependent). The mean differences remain in the range $\sim$0.6–1.8 deg, above IS22 nominal resolution at 1 Hz (⁠$\sim {0.5}^{\circ }$⁠; Le Pichon et al. 2005a). These differences range from 0.6 to 1.2 deg for Yasur (⁠$\sim$400 km), increasing to 0.9–1.8 deg for Ambrym (⁠$\sim$700 km). Such a difference increase from Yasur to Ambrym cases can be attributed to the intrinsic compounded error of ray propagation modeling with travel path length (⁠$\sim$400 to $\sim$700 km), and the increase in unmodelled atmospheric variability in the empirical climatologies at above mesospheric to lower thermospheric atmospheric altitudes (⁠$\sim$80–120 km; Emmert et al. 2020), which are particularly relevant for ray paths from Ambrym to IS22.

There is significant spread (⁠$\sigma$⁠) around the mean differences, representing 42–92 per cent of each value (see Table 1). It is important to note that this value does not provide an accurate interpretation of the actual expected error, as the error distribution tends to be one-sided. To better understand the expected error in the back-azimuth deviation, the monthly curves as Fig. 8 should be used, as they provide not only the magnitude but the sign of the back-azimuth deviation of the modelled infrasound in time. The variance of the candidate infrasound observations is highly influenced by overlapping non-volcanic sources (clutter) which remain present in the reduced data sets (Figs 6 and 7). Overall, these differences (⁠$\sim$1–3${}^{\circ }$⁠) seem to be in the same range of previous studies (Le Pichon et al. 2005a; Antier et al. 2007; Assink et al. 2014a; see Fig. 8).

The range-independent models provide comparable or better results than the range-dependent models (Table 1). This seems counterintuitive, but it tells us that the atmospheric descriptions in the area of study only change significantly in altitude at each modelled time. Indeed, the smallest horizontal resolution is 0.25 deg (ERA 5), or about 27 km, meaning that atmospheric changes in the order of $\sim$100 km or bigger can be mapped sufficiently well. That is, for a specific day and time, localized atmospheric changes (wind speed, temperature, etc.) at this scale are not large enough to produce arrivals that are much different from the along-path mean descriptions of a representative layer (range-independent). Considering the increased computation time that range-dependent ray tracing requires compared with range-independent ray tracing (10 times as mentioned in the infraGA manual, Section 3.2; Blom 2014), the range-independent models are the right choice with the level of resolution and the spatiotemporal scales we pursue for modelling the azimuth deviation effects for a robust and rapid calculation scenario (e.g. a near-real-time volcano monitoring tool).

The hybrid models are the only ones able to map shorter term wind reversals that happen in the middle of the year (see Figs 6 and 7). However, the monthly means (Figs 8a and b) indicate that the hybrid models predict back-azimuths that deviate significantly from those observed. A rule of thumb to calculate year-long predictions for a specific area would be to favour climatology-based models during periods with stable seasonal ducting, while switching to NWP-based atmospheric models when events as SSW or QBO-driven wind reversals have been observed in the past. Such analysis could also benefit from the characterization of seasonal trends of stratospheric waveguides and detectability on pre-defined areas or station locations (e.g. Blom et al. 2023).

We note that here we use the persistent Yasur and Ambrym volcanic sources to evaluate our method on continuous data. However, this method is generally intended for sudden, unexpected, large explosive eruptions without recent priors, such as the 2011 Puyehue-Cordón Caulle eruption, Chile; or the 2015 Calbuco eruption, Chile (De Negri & Matoza 2023). For well documented, persistent, infrasound sources such as Yasur or Ambrym volcanoes, other approaches based on data regression or machine learning could better perform in terms of their predictive accuracy. However, these approaches would not be suitable for volcanic observatories without long records of prior detections and new back-azimuths, real time atmospheric data access, and/or computational resources. That is, climatologies with ray tracing (predictions based on atmospheric structure only) are still needed to obtain critical information to characterize new transient infrasound during a volcanic emergency (De Negri & Matoza 2023).

5.1 Seasonal and daily changes of the models

Two interesting characteristics stand out at seasonal and diurnal time scales (Figs 9 and 10, for 2004 Ambrym and 2019 Yasur).

First, the asymmetry of the negative versus positive peak seasonal back-azimuth deviations is much more marked in the empirical climatologies than the hybrid models. This peak-to-peak asymmetry seems to be driven mainly by the 12:00 UTC arrivals, especially in the climatology-based results. For Yasur volcano (southwest propagation), the empirical climatologies show deviations that have maximum values of $\sim$$-5^{\circ }/+{2.5}^{\circ }$⁠, while the hybrid models have maximum values of $\sim$$-5^{\circ }/+5^{\circ }$ (Fig. 9). For Ambrym volcano (south propagation), a similar asymmetry is present for the empirical climatologies, but the maximum deviation values are $\sim$$-7^{\circ }/+4^{\circ }$⁠, while the hybrid models predict extremes of $\sim$$-5^{\circ }/+6^{\circ }$ (Fig. 10). This is related to the change rate of the back-azimuth deviation from the seasonal negative peak to the positive peak (and vice-versa) which is higher for the hybrid models. On average, the rates of back-azimuth change per day on the first 100 d for the Yasur models are about $7.5/100{[}^{\circ }/\mathrm{day}]$ for the climatologies (Figs 9a, b, c and e, right) while $10/100{[}^{\circ }/\mathrm{day}]$ for the hybrid models (Figs 9d and f, right). For Ambrym, these values increase to about $9/100{[}^{\circ }/\mathrm{day}]$ for the climatologies (Figs 10a–c, right) while $12/100{[}^{\circ }/\mathrm{day}]$ for the hybrid models (Figs 10d and e, right). As the empirical climatologies are based on multiyear interpolated observations back to $\sim$30 yr (Drob et al. 2015; Emmert et al. 2020), this may possibly indicate a shift over the last three decades to higher transition rate changes between seasons. However, this is only partially supported by qualitative comparison with two specific years in the area (2004 and 2019). To better understand this behaviour, further analysis should be performed by comparing a more complete set of year-long models (e.g. 2004–2020), and calculating the back-azimuth rates of the interpolated arrivals.

Second, all models present a similar semidiurnal behaviour, showing the biggest deviations at 00:00 and 12:00 UTC, and the smallest deviations at 06:00 and 18:00 UTC (e.g. days 1 to $\sim$100, and $\sim$300 to 365). The change in back-azimuth deviation at distinct times of the day suggests the thermospheric duct, whose temperature is affected by the Sun’s daily energy input, may have a visible influence in the total observed back-azimuth deviation. However, this semidiurnal effect is more marked for the climatology-based than the hybrid models, suggesting that stochastic atmospheric variations could reduce the importance of this contribution. In Section 5.2, we show that there is an observable difference in number of detections during the day at least for Ambrym case (Fig. 11).

The diurnal standard deviations can be above 1${}^{\circ }$ in periods of strong zonal winds (westwards; day $\sim$300 to 365 and 1 to $\sim$80; Figs 9 and 10 vertical bars around interpolations on right subfigures), declining to less than 0.5${}^{\circ }$ during weaker zonal wind conditions. These variations represent a secondary effect as the seasonal back-azimuth deviations can have magnitudes between $\sim$5 and 7${}^{\circ }$⁠, but could add significant differences in the back-azimuth source location for long-range propagation.

5.2 Feasibility of thermospheric arrivals

Our models use ground intercepts with turning heights at stratospheric (⁠$\sim$40 km) and thermospheric (⁠$\sim$110 km) altitudes (Fig. 4), maximizing the number of available (near enough IS22) arrivals. But infrasound that refracts at higher altitudes (thermospheric) suffers a higher transmission loss (Sutherland & Bass 2004), and could be undetectable at the recording sensors.

Interestingly, when arranging the multiyear observations from Ambrym direction by time of day (Figs 11b, d, f and h; also Fig. S15 in the Supporting Information), we see two 6-hr time windows with larger numbers of detections centred at $\sim$07:00 and $\sim$17:00 UTC (18:00 and 04:00 local time, respectively). These are followed by two periods of low detections, centred at $\sim$12:00 and $\sim$01:00 UTC (23:00 and 12:00 local time). On the other hand, the detections from Yasur direction (Figs 11a, c, e and g; also Fig. S15 in the Supporting Information) are mostly continuous and more numerous, with one low detection time window at $\sim$01:00 UTC. This low detection rates period has been attributed to daily increased local winds, which significantly reduce the signal-to-noise ratio at IS22 (Le Pichon et al. 2005b).

The diurnal variability is likely the effect of migrating solar tides in the mesosphere and thermosphere with dominant periods of 24, 12 and 8 hr (Le Pichon et al. 2005b). Throughout the year, the thermospheric duct turning height, mainly driven by changes in temperature, is consistently lower at 06:00 and 18:00 UTC (⁠$\sim$113–115 km), and higher at 00:00 and 12:00 UTC (⁠$\sim$117-121 km; see Fig. S16 in the Supporting Information). A shorter travel path (i.e. lower duct) could reduce the transmission loss by about 40 dB, allowing these waveforms to be detectable at IS22 (Sutherland & Bass 2004; Le Pichon et al. 2005b). Moreover, the infrasound propagation modelling results suggest that the bulk of detections at IS22 from the direction of Ambrym are thermospheric (Fig. 4). This means their detectability would be highly influenced by periods of lower transmission loss, matching the observations.

The thermospheric duct dependence of Ambrym detections can also explain why they appear during zonal wind transition periods (equinoxes) in the multiyear stacked data set (March–June; August–October; Fig. 2f). From the ray tracing modelling point of view, when the stratospheric duct weakens, the cross-winds reduce for North–South infrasound, making long-range infrasound propagation less sensitive to initial launching parameters. Then, the modelled rays rely on the thermospheric duct, which does not impose strong cross-winds, allowing the iterative search to find arrivals near IS22 at the cost of increased transmission losses (see Figs 2f and 5). On the other hand, the detections from Yasur are produced by East–West infrasound propagation, and will be mostly dependent on the westward stratospheric duct (January–March; November–December; Fig. 7). This means that the diurnal variations will be less relevant as most of the detections are not influenced by thermospheric ducting changes.

We further investigate the transmission loss versus time of day relation with our methodology. Using range-independent hybrid atmospheric descriptions, we calculate transmission loss (TL) estimations for 2004 Ambrym and 2019 Yasur every 30 d at 00:00, 06:00, 12:00 and 18:00 UTC from the medians of the arrivals that fall at less than 50 km from IS22. Our results show that, as hypothesized, smaller attenuations concentrate at 06:00 and 18:00 UTC for both volcanoes (Fig. 12; Table S2, Supporting Information). Based on the wealth of observations at IS22 and ray tracing simulations, there is good evidence to indicate that thermospheric arrivals from Ambrym have indeed been detected at IS22, at 06:00 and 18:00 UTC, aided by lower thermospheric duct altitudes.

6 CONCLUSIONS

We tested a methodology that aims to calculate robust, first-order estimates of long-range infrasound arrival characteristics, particularly the back-azimuth deviation, using empirical climatologies. We produced year-long models of infrasound propagation for arrivals at IMS station IS22 (Grande Terre, New Caledonia) from the Vanuatu archipelago region, including from Yasur (Tanna island, at $\sim$400 km and 43${}^{\circ }$N azimuth) and Ambrym (Ambrym island, at $\sim$670 km and 12${}^{\circ }$N azimuth) volcanoes.

Our goal was to investigate if the empirical climatologies were sufficient, and to what extent, to model the infrasound propagation path effects of the atmospheric changes in altitude of the horizontal winds, temperatures, densities and pressures during a year. We chose a case where the source locations can be assumed as known and nearly-continuous, hence leaving the source mislocation and detectability effects as induced by the atmospheric variability in time. We modelled such atmospheric parameters for infrasound propagation from the ground to 150 km with empirical climatologies (HWM14/MSIS2.0), and also considered two approaches for automatically perturbing the descriptions to enhance along-path stratospheric winds. Our simplified perturbation schemes are intended to mimic the effect of gravity-wave perturbations, without having to resort to a full statistical ensemble treatment and thus saving computation time and simplifying the implementation.

To evaluate the accuracy and limitations of the above climatology-based approaches, we also repeated the calculations using more accurate hybrid models that merge the ERA 5 descriptions (winds and temperatures) from the ground to $\sim$80 km with the climatologies for higher altitudes (⁠$>$80 km). To investigate the degree of model complexity required to match observations, we performed both range-independent and range-dependent 3D ray tracing modelling (infraGA/geoAC; Blom & Waxler 2012). In total, we compared six different atmospheric model parametrizations with the multidecadal (2003 up to 2022) observations. We used the 2003–2022 open-access high-frequency (1–3 Hz) PMCC products first released by Hupe et al. (2021).

We found that the range-dependent propagation models provide comparable results to the range-independent cases, supporting the use of a simplified, layered model of the atmosphere for regional and global calculations. From July to April, all of the models depict seasonal back-azimuth changes that are generally aligned with the year-long and stacked multiyear observations. From May to June, the reversal of the zonal stratospheric winds from prevalent westwards to eastwards winds shows a reduction in observations, larger differences with the models, and only transient groups of detections likely produced by regional scale wind perturbations and transient favourable ducts (e.g. Sudden Stratospheric Warming events). The mean difference per year between the models and the high-frequency PMCC products shows that the empirical climatologies can provide comparable results with the hybrid models. The match is generally better for Yasur volcano than for Ambrym, which is likely the combined effect of increased crosswinds and propagation path lengths, plus other possible eruptive sources from similar locations such as Gaua and Ambae volcanoes complicating the attribution of detections uniquely to Ambrym. For specific years, only the hybrid models were capable of reproducing detections from transient ducts during the worst matching period (April–June). The hybrid models also better characterize the back-azimuth changes during the equinox transition periods towards the stratospheric seasonal wind reversal (around April and October).

The results from our analysis of the IS22 persistent multiyear long-range volcanic-infrasound detections indicate that empirical climatologies are useful and adequate for a rapid infrasound propagation calculation scenario (such as for a volcano monitoring system). This represents a step forwards from assuming a single value of celerity and neglecting cross-winds, which has previously been assumed for rapid calculations (e.g. Matoza et al. 2017). In particular, we found that the range-independent, raw climatology models are well-suited to obtain robust first-order estimates of arrival back-azimuths for large-scale, year-long infrasound propagation in this region. However, the climatologies do not perform as well as the hybrid models when the stratospheric duct weakens (equinoxes) or when short-term (days) wind reversals allow favourable ducts to appear during upwind conditions (summer). In such cases, a more accurate model with hybrid descriptions should be performed at a later time, given more time to download the specifications and run the simulations. We also emphasize that propagation models based on climatologies can be pre-computed, since climatologies do not use any current atmospheric data. Pre-computation of infrasound propagation paths using climatologies could be useful for near-real-time infrasound monitoring applications such as in volcano-monitoring systems (De Negri & Matoza 2023).

Our proposed adaptive perturbation of the empirical climatologies, which is intended as a simplified alternative to a full gravity-wave realization approach (ensemble averaging), was successful at enhancing stratospheric predictions. This approach may also be convenient for a rapid scenario where only an approximation of expected backazimuth deviations is needed. As stratospheric gravity waves are known to considerably reduce infrasound transmission losses (i.e. enhance propagation and ducting conditions; Listowski et al. 2024), future enhancements to the adaptive approach could be explored.

Finally, we remark that none of the atmospheric models we considered, including the hybrid specifications, fully explain all of the observations (e.g. more negative deviations on the first $\sim$50 d of 2004 for Ambrym, and $\sim$250–350 of 2019 for Yasur on Fig. 6). Thus, our study supports the notion that continuous infrasound observations can be used to verify and potentially improve atmospheric specifications (e.g. Le Pichon et al. 2005b; Antier et al. 2007; Assink et al. 2014a; Vorobeva et al. 2023). To better resolve middle atmospheric dynamics in the mesospheric to lower thermospheric altitude range, other numerical weather prediction models could be used in place of the hybrid approach we implemented, such as the whole atmosphere community climate model with thermosphere and ionosphere extension (Liu et al. 2018) which is nudged with the modern-era retrospective analysis for research and applications (e.g. Vorobeva et al. 2023). However, this is beyond the aim of this study, as our goal was to assess a simplified low-computational cost methodology for first-order estimates.

This work focuses on volcanic infrasound detections at one station in the International Monitoring System. Generalizing these results will be important, both for other volcano-to-array paths globally and for other monitoring applications (e.g. global nuclear test monitoring), as mislocations due to unmodelled backazimuth deviations are a global issue. We suggest that additional long-range studies using other persistent volcanic source locations with large available detection records (e.g. Etna volcano) could be a step forward in this direction, further evaluating the applicability of this method beyond volcanic infrasound source location in near-real time. Additionally, considering other persistent anthropogenic or natural infrasound sources with large historical sets of detections (e.g. microbaroms, dams, etc.) could serve to understand to what extent this methodology could be applied to other infrasound generating phenomena.

ACKNOWLEDGMENTS

RSD thanks C. Scanlon for reviewing and improving the manuscript. All authors thank the CTBTO and the IMS station operators for guaranteeing the high quality of the infrasound data. The views expressed herein are those of the authors and do not necessarily reflect the views of the CTBTO Preparatory Commission. We also wish to thank Dr. David Green and Dr. Stephen Arrowsmith for their reviews, which helped us to improve the manuscript. This work was funded by National Science Foundation (NSF) Grant EAR-1847736.

DATA AVAILABILITY

All International Monitoring System data are available for scientific studies through the CTBTO vDEC platform (https://www.ctbto.org/specials/vdec/). The PMCC products used in this study can be accessed through the Geoportal of the Federal Institute for Geosciences and Natural Resources (BGR, Germany; https://geoportal.bgr.de) (see Hupe et al. 2022). For the hybrid models, we used the CDS API tools freely provided by European Centre for Medium-Range Weather Forecasts (ECMWF) to obtain the necessary ERA 5 reanalysis profiles, publicly available for academic research (https://confluence.ecmwf.int/display/CKB/How+to+download+ERA5). The climatology (HWM/MSISE) model codes are available from Drob et al. (2015) and Emmert et al. (2020), respectively. The 3-D ray tracing infrasound propagation algorithms are available for research purposes from the Los Alamos National Laboratory (LANL, USA) GitHub repository (https://github.com/LANL-Seismoacoustics/infraGA). The main code used in this project (ARCADE), is available for research purposes in the GitHub repository (https://github.com/rodrum/arcade) along with instructions and examples. The last stable code release can be accessed in Zenodo in 10.5281/zenodo.13383432 (De Negri 2024).

REFERENCES