Microbarom radiation and propagation model assessment using infrasound recordings: a vespagram-based approach

This study investigates the use of a vespagram-based approach as a tool for multi-directional comparison between simulated microbarom soundscapes and infrasound data recorded at ground-based array stations. Data recorded at the IS37 station in northern Norway during 2014− 2019 have been processed to generate vespagrams (velocity spectral analysis) for five frequency bands between 0.1 and 0.6 Hz. The back-azimuth resolution between vespagrams and a microbarom model is harmonized by smoothing the modelled soundscapes along the back-azimuth axis with a kernel corresponding to the frequency5 dependent array resolution. An estimate of similarity between the output of a microbarom radiation and propagation model and infrasound observations is then generated based on the image processing approach of mean-square difference. The analysis revealed that vespagrams can monitor seasonal variations in the microbarom azimuth distribution, amplitude, and frequency, as well as changes during sudden stratospheric warming. The vespagram-based approach is computationally inexpensive, can uncover microbarom source variability, and has potential for near-real-time stratospheric diagnostics and atmospheric model 10 assessment.

well as for long-term assessment of model product uncertainties, particularly when applied to data from a global network of infrasound stations. A key aspect of this approach is that benchmarking between model and infrasound vespagrams considers all 55 back-azimuth directions rather than just the direction of the dominant microbarom source, as done in several previous studies (Garcés et al., 2004;Hupe et al., 2019;De Carlo et al., 2019;Smirnov et al., 2020;De Carlo et al., 2020). The microbarom soundscape at a station is typically a sum of components stemming from a wide spatial distribution of ocean regions, and recently den Ouden et al. (2020) demonstrated that an iterative decomposition of the array spatial covariance matrix using the CLEAN algorithm (Högbom, 1974) can be exploited to resolve the back-azimuth and trace velocity of the most coherent wave 60 front arrivals.
A long-term ambition is to exploit microbarom infrasound datasets to enhance the representation of stratospheric dynamics in atmospheric model products and hence increase the accuracy of both medium-range weather forecasting and sub-seasonal climate modeling (Büeler et al., 2020;Dorrington et al., 2020;Domeisen et al., 2020a, b). In addition to prospective numerical weather prediction improvements, the suggested vespagram-based approach may be applied in multi-technology studies of of a microbarom signal from the Atlantic Ocean, with a back-azimuth of 225 • and a 350 m/s apparent velocity typical of the stratospheric regime (Garcés et al., 1998;Whitaker and Mutschlecner, 2008;Nippress et al., 2014;Lonzaga, 2015). The base resolution of the array was taken to be the 1-sigma beam width of the Gaussian fitted to the array response at a constant velocity 4 https://doi.org/10.5194/angeo-2020-78 Preprint. Discussion started: 8 December 2020 c Author(s) 2020. CC BY 4.0 License. of 350 m/s (dashed line in Figure 1b) for each frequency band. The resulting resolution was found to be: 35 • , 23 • , 16 • , 13 • and 90 10 • for 0.1 − 0.2 Hz, 0.2 − 0.3 Hz, 0.3 − 0.4 Hz, 0.4 − 0.5 Hz and 0.5 − 0.6 Hz band, respectively. It should be noted that this estimate is based on the homogeneous medium plane-wave time-delays between the array elements only and does not take into account meteorological conditions at the station, noise, or other coherence loss mechanisms that may result in a wider beam width.
In array signal processing, separating coherent from incoherent parts of the recorded signal, as well as the separation between 95 different simultaneous arrivals are important concepts. When analyzing the wavefield in terms of a given horizontal slowness vector (e.g., described in terms of apparent velocity and back-azimuth), delay-and-sum beamforming (Ingate et al., 1985) is usually applied in combination with the underlying plane-wave model assumption. This method applies time-delays to the array sensor traces to focus on wave fronts arriving with a specific horizontal apparent velocity and a specific back-azimuth direction, hence amplifying wavefield components with the horizontal slowness of interest, while suppressing other components.

100
However, the slowness vector models are not always accurate (Gibbons et al., 2020). In particular, the actual shape of the wave front arriving at infrasound arrays may differ from a theoretical plane-wave due to meteorological conditions and turbulence at the station, which make the underlying assumption of a locally homogeneous effective sound speed invalid. In this case, the beamforming is less efficient and the reduced array gain results in lower stack amplitude and signal distortion (Rost and Thomas, 2002).

105
To determine an unknown slowness vector component and to study the spatial structure of the wavefield over time, one can use the vespa (velocity spectral analysis) processing. This not only enhances the signal as the beamforming does, but also allows one to determine either the direction or apparent velocity of incoming signal. The vespa method estimates the power of the signal either for a fixed apparent velocity with varying back-azimuth or for a fixed back-azimuth with varying apparent velocity.
The result of the vespa processing is usually presented as an image displaying the power of incoming signal as a function of 110 time and back-azimuth (or apparent velocity) called vespagram. Despite that vespa is a widely applied in seismological array data studies (e.g., Davies et al., 1971;Kanasewich et al., 1973;Muirhead and Datt, 1976;McFadden et al., 1986), it has not previously been exploited in peer-reviewed microbarom infrasound studies.
The vespa processing procedure described below is applied to each analyzed time window and frequency band: 1) For each sensor n of an array, we extract signal recording x n (t) that corresponds to the time window of interest. The 115 analysis here is done for an 1h moving time window, evaluated every 30 min. In general, the time series recorded at sensor n at the location r n can be written as where y(t) represents a plane wave front signal, and s hor is the horizontal component of the slowness vector.
3) Apply a Butterworth bandpass filter to recordings. Calculations are performed for five equally spaced frequency bands 120 that are within the microbarom frequency range (see Figure 1b). x n (t + r n · s hor ). (2) In this study, classical linear vespa processing (Davies et al., 1971) is applied where the noise suppression is proportional to square root of N (Rost and Thomas, 2002). A beam is generated at each 1 • in back-azimuth, for the fixed apparent velocity of 350 m/s, which is within stratospheric arrival regime (Garcés et al., 1998;Whitaker and Mutschlecner, 2008;125 Nippress et al., 2014;Lonzaga, 2015). That allows to estimate signals coming from all directions but from approximately the same height corresponding to stratospheric altitudes. 5) Calculate mean squared pressure (power) of each beam to get an estimate of incoming signal strength as a function of back-azimuth and time.
Steps (1) -(5) are applied to all analyzed years of data.

Microbarom source and propagation modeling
In this section we summarize the approach applied to get directional spectrum of microbarom soundscape as a function of time.
The procedure is as follows. This microbarom model allows prediction of the location and intensity of the microbarom sources when applied to the 145 Hasselmann integral. The Hasselmann integral is derived from the output of the wave model and establishes a relationship between the source spectrum and the spectral densities of counter propagating waves for a given frequency (Hasselmann, 1963). The output of this step is an acoustic spectrum for each cell of the wave model. After applying these steps, and integrating over the frequency bands, we get an estimate of microbarom amplitude as function  (−95 ± 16), −21 ± 14 (−15 ± 8) and 26 ± 6 (34 ± 7).
A similarity index (SI), taken from an imaging processing approach, is introduced as where MSE is a mean squared error between normalized smoothed model output, P model (t, θ), and normalized vespagram,  195 In winter, SI for lower frequencies is stable and has values ∼ 1, with exceptions corresponding to increased noise level in vespagrams or to SSW events that are discussed below. Relatively low SI for higher frequencies can be explained either by spurious apparent sources corresponding to array response function side-lobes ( . The similarity estimation proposed allows detection of inconsistencies between the microbarom model and the vespa processing which might be used for identifying biases in atmospheric models. This is especially promisingly for low frequencies where side-lobes of array response do not appreciably affect analysis.

Examination of major sudden stratospheric warmings
Although this is not the main objective of the current study, in this section we examine the ability of the vespagrams to detect https://doi.org/10.5194/angeo-2020-78 Preprint. Discussion started: 8 December 2020 c Author(s) 2020. CC BY 4.0 License. SSWs usually occur in wintertime and are, in general, associated with a sudden and short increase in stratospheric temperature and mesospheric cooling at high / middle latitudes (Shepherd et al., 2014;Butler et al., 2015;Limpasuvan et al., 2016;Zülicke et al., 2018). SSWs are often classified into minor and major warmings, depending on whether there was a weakening or reversal this event is classified as minor, but in a number of studies it is referred to as major (Xiong et al., 2018;Conte et al., 2019).
Vertical dashed lines in Figures 5 -6 correspond to the onset days listed above when SSW criteria were met.
The infrasound signature reported by Donn and Rind (1971) and Evers and Siegmund (2009), which showed a significant change in direction of the infrasound arrival due to a change in favorable stratospheric waveguide, can be seen in Figure 6 for all that the signature appears late in the model data and its duration is much shorter than in vespagram, analogous to study by Smets et al. (2016). For higher frequencies ( Figure 3) the duration of a change from eastward to westward pattern is longer and continues until late March -early April that corresponds to reanalysis data (Manney and Lawrence, 2016).

235
Another feature revealed is a significant decrease in similarity index between model and observations during SSWs ( Figure   5) which is characteristic for all events under consideration. The smallest discrepancies in the direction of the dominant wave front between the model and infrasound data during SSWs reach about 5 • − 7 • , but the largest reach as much as 90 • − 100 • ( Figure 6). This may be caused by the following factors. The back-azimuth change during SSW usually appears earlier in the vespagrams than in the model with the difference of 3 to 24 hours. Similar results were previously obtained by Smets and

Discussion and Conclusions
In this study, we compare observed and predicted microbaroms soundscapes using a vespagram-based approach. Analysis is 255 performed based on calculation of microbaroms power as a function of time and back-azimuth at constant apparent velocity of 350 m/s. Note, however, that the vespagram-family of time-dependent microbarom data visualizations can be constructed also using other array processing techniques that estimate power as function of the slowness of the wave front, e.g., using robust estimators as explored by Bishop et al. (2020), or adaptive high-resolution approaches like Capon's method (Capon, 1969).
recorded infrasound data for all directions simultaneously, as opposed to methods where only the back-azimuth direction of maximum power is considered (e.g., Hupe et al., 2019;Smirnov et al., 2020). Since the vespa processing is computationally low-cost and able to track variations in microbarom parameters over extended periods spanning one or several years, it can be utilized for near-real time assessment of atmospheric model products and for developing infrasound-based stratospheric diagnostics. It can also be used when assessing changes in infrasound signatures over shorter time windows, e.g., during extreme which can account for both range-dependent atmospheric models and cross-wind effects (e.g., Smets and Evers, 2014;Smets et al., 2016). However, the inherent high-frequency approximation of the ray-theory can limit the modelling of diffraction and scattering effects (Chunchuzov et al., 2015) that can be important for the low-frequency microbaroms.
A more elaborate microbarom propagation model could also allow for an estimate of the full microbarom wavefield impinging an infrasound station, hence providing an estimate of its power within the full horizontal slowness space of plane wave front 280 directions (or a selected relevant region). This way, we could benchmark the modelled and recorded microbarom field at an infrasound array for each sliding time window without restricting the analysis to the region around a fixed apparent velocity as carried out in the current study. Notably, such "f-k plots" of modelled and recorded microbaroms are also (time-varying) images which can be assessed and compared using the versatile ecosystem of image processing and image comparison algorithms.
Future developments can include compilation of long-term time-dependent statistics of similarity between model and 285 infrasound recordings for multiple stations on global and regional scales, in order to define anomaly flag criteria which would indicate that there is unexpected inconsistency between model and observations due to, for example, biases in atmospheric model products. Moreover, we suggest to apply the presented approach in global assessment and comparisons of ocean wave-action model products, as well as in validation and further refinement of microbarom radiation estimation algorithms.