Assessing middle atmosphere weather models using infrasound detections from microbaroms

The nonlinear interaction of ocean surface waves produces coherent infrasound noise—microbaroms—between 0.1 and 0.5 Hz. Microbaroms propagate through the atmosphere over thousands of kilometres due to low absorption and efficient ducting between the ground and the stratopause. These signals are globally and permanently detected by the International Monitoring System (IMS) infrasound network, which has been established to monitor compliance with the Comprehensive Nuclear-Test-Ban Treaty. At the International Data Centre (IDC) in Vienna, where IMS data are routinely processed, microbarom detections appear in overlapping frequency bands, and are treated as false alarms. Therefore, understanding the variability in microbarom detections is essential to support the IDC in the reduction of the false alarm rate. In this study, microbarom amplitudes and the direction of arrivals at the German infrasound station IS26 were modelled. For the simulations, the source was described by an operational ocean wave interaction model, and the signal amplitude was modelled using a semi-empirical attenuation relation. This relation strongly depends on middle atmosphere (MA; i.e. 15–90 km altitude) dynamics; however, vertical temperature and wind profiles, provided by numerical weather prediction (NWP) models, exhibit significant biases and differences when compared with high-resolution light detection and ranging instrument (lidar) soundings in altitudes where infrasound signals propagate. To estimate uncertainties in the modelled amplitude, a fully autonomous lidar for MA temperature measurements was installed at IS26. Temperature and wind perturbations, considering observed biases and deviations, were added to the operational high-resolution atmospheric model analysis produced by the European Centre for Medium-Range Weather Forecasts. Such uncertainties in horizontal winds and temperature strongly impact propagation conditions, explaining almost 97 per cent of the actual detections, compared to 77 per cent when using the direct output of the NWP model only. Incorporating realistic wind and temperature uncertainties in NWP models can thus significantly improve the understanding of microbarom detections as well as the detection capability of a single station throughout the year.


I N T RO D U C T I O N
The Atmospheric dynamics Research InfraStructure in Europe (ARISE) project has established an infrastructure of complementary observation technologies, such as light detection and ranging instruments (lidars) and infrasound recordings, for enhancing knowledge of middle atmosphere (MA) dynamics (Blanc et al. 2018). The in-frasound technology therefore provides an appropriate monitoring technique because low-frequency sound can propagate through the atmosphere over large distances due to minor attenuation (Evers & Haak 2010) at a range of altitudes where routine observations are lacking (e.g. Le Pichon et al. 2015).
The infrasound monitoring technology is one of the verification techniques used for compliance with the Comprehensive Nuclear-Test-Ban Treaty (CTBT). The International Monitoring System (IMS) has been established to provide worldwide coverage for the detection of nuclear explosions with a minimum yield of 1 kt TNTequivalent. When fully implemented, 60 stations make up the IMS infrasound network, which is dedicated to globally monitoring the atmosphere. At the International Data Centre (IDC) in Vienna, the IMS data are routinely processed to detect coherent, low-frequency pressure waves (e.g. Marty 2019). Among these are microbaroms, a quasi-continuous natural infrasound source generated by standing ocean surface waves (Donn & Naini 1973). They are regularly and globally detected, covering a frequency range between 0.1 and 0.5 Hz. In the context of the CTBT, microbaroms produce a high false alarm rate in automatic detection lists of the IDC (Arrowsmith 2018) and thus extend the time required to provide the Reviewed Event Bulletins. Understanding the seasonal variations in microbarom sources and atmospheric propagation parameters allows a quicker identification and rejection of these detections.
The focus of this study was to better understand seasonal variations in the characteristics of the microbarom detections, such as the amplitude and direction of arrival, as recorded at IMS station IS26 in southern Germany (48.85 • N, 13.71 • E). The modelling was carried out with regard to the spatio-temporal variability of the source term and propagation conditions in the MA. The spatio-temporal evolution of the sources was represented by the ocean wave interaction model that was developed by the French Research Institute for Exploitation of the Sea (IFREMER) for estimating seismic noise (Ardhuin et al. 2011). Prediction of the azimuthal distribution was adapted from Landès et al. (2014). Using MA specifications (vertical profiles of wind speed and temperature) derived from numerical weather prediction (NWP) model data of the European Centre for Medium-Range Weather Forecasts (ECMWF), Landès et al. (2014) explained the seasonal variation of microbarom signals.
In this study, the signal amplitude at the station was modelled by applying the attenuation relation proposed by Le Pichon et al. (2012). This accounts for the effects of the source frequency, propagation range and along-path effective sound speed as a measure for atmospheric propagation conditions. Microbaroms can potentially propagate through the atmosphere over large distances due to low absorption at 0.1-0.5 Hz (Sutherland & Bass 2004) and efficient atmospheric ducting between the ground and the stratopause (Drob et al. 2003;Landés et al. 2014). Ducting is dependent on the 3-D wind and temperature field and the propagation direction. The main characteristics of sudden stratospheric warming (SSW) events have been successfully derived from directional microbarom amplitude variations resulting from changes in stratospheric propagation conditions (e.g. Smets & Evers 2014). As infrasonic waves propagate into the MA, significant features of the vertical structure of temperature and wind are reflected in the signal detected on the ground (e.g. Kulichkov 2010). Unresolved fluctuations in the temperature or wind profiles that form the waveguide may significantly affect the received signals.
Uncertainties in the ECMWF products used were quantified using recent high-latitude observations. A temperature bias in the ECMWF high-resolution (HRES) operational analysis (L137) part of the Integrated Forecast System (IFS) has been highlighted by comparing it to measured temperature profiles from MA lidar systems (Ehard et al. 2017;Hildebrand et al. 2017). For quantifying biases and deviations at mid-latitudes, and their potential effects on amplitude predictions at IS26, a portable lidar for measuring temperature profiles between 30 and 90 km altitude was colocated to the station for 7 months in 2016. The distribution of the effective sound speed profile was here adjusted using measured differences between ECMWF analyses and lidar profiles. This was then used as an input for calculating the attenuation of the source pressure. A range of uncertainty for the predicted amplitude at IS26 was inferred, showing that uncertainties in NWP models can strongly affect the predicted amplitudes. The benefit of measuring wind profiles at IMS stations to better understand the variability of the observed signal amplitude and station detection capability is essential in the context of the CTBT.

Infrasound propagation and attenuation in the atmosphere
For the eight-element infrasound array IS26, the Progressive Multi-Channel Correlation (PMCC) algorithm (Cansi 1995) was used for detecting coherent waves with frequencies up to 5 Hz. Focusing on the frequency range of microbaroms (0.1-0.5 Hz), the detections follow a seasonal variation in azimuth (see Supporting Information Fig. S1). During winter, strong predominant signals are detected from northwesterly directions; whereas during summer, signals are more scattered and less pronounced, originating from northerly or southeasterly directions. This is related to prevailing winds in the stratosphere forming the essential waveguide between the ground and stratopause altitudes (e.g. Pilger et al. 2017).
In order to estimate the waveguide, the effective sound velocity (v eff ) was introduced: where c is the adiabatic sound speed equal to 20.05 √ T with temperature T, v h is the horizontal wind vector projected in the direction of wave propagation (in the following, denoted by azimuth α) described by vector n. The effective sound speed ratio (v eff −ratio ) is the quotient of the temperature-and wind-dependent sound speed at about 50 km altitude relative to the effective sound speed at the ground. Ratios greater than or equal to 1 indicate favourable conditions for acoustic ducting in the ground-to-stratosphere waveguide. Le Pichon et al. (2012) parametrized factors controlling attenuation for different frequencies (f) and values of v eff −ratio representative of realistic atmospheric conditions. The following formula was proposed for calculating the attenuation coefficient: at a distance R (in km) from a reference distance to the source of 1 km. The first term of eq. (2) describes the near-field attenuation in the shadow zone, where a parametrizes the dissipation of direct waves. The second term accounts for far-field attenuation in the geometrical acoustic duct region, and comprises three parameters: b defines the geometrical spreading and attenuation of waves, and s is a scaling distance representing attenuation in the shadow zone, the width of which is given by d.
These parameters were taken from tabulated values (Le Pichon et al. 2012), whereas b had to be additionally interpolated with regard to v eff −ratio . In order to account for geometric acoustic ducting in the stratopause region, the maximum of v eff −ratio was calculated between 40 and 60 km altitude for all backazimuths (β = α -180 • ), with a resolution of 1 • . As input for eq. (1), the temperature and wind vertical profiles were extracted from the ECMWF HRES operational analyses, provided at 137 vertical levels up to 0.01 hPa.

Identifying sources of microbaroms detected in southern Germany
Microbarom source terms were modelled using noise generation theory due to nonlinear ocean wave interaction, as developed by Ardhuin & Herbers (2013). Input to their model is data from NASA's wave action model, WAVEWATCH-III and ECMWF surface winds (Ardhuin et al. 2011). The output covers a global 0.5 • × 0.5 • grid up to ±78 • latitude, with a temporal resolution of 3 hr.
The source region considered herein was the area between 20 • N and 78 • N, within a longitudinal range between 60 • W and 80 • E. As shown in Fig. 1, the source amplitudes are strongest in the North Atlantic Ocean, where the dominant frequency is between 0.2 and 0.3 Hz. This corresponds to the mean backazimuths of microbarom detections at IS26, ranging from 240 • to 360 • . Along the coastlines, in shallower water, and in spatially limited seas, the dominant frequencies are higher due to reduced wavelengths compared to the open ocean (Fig. 1a). Further analyses in this study referred to the dominant frequency range between 0.15 and 0.35 Hz. In addition, filtering parameters were applied to PMCC processing results in order to focus on detections with high significance (see Supporting Information Table S1). These filters resulted in 6893 detections remaining for the considered 3 yr period, from 2015 to 2017, of which more than 95 per cent originated from northwesterly directions. For these dominant directions, the model's maximum acoustic source pressure, determined for each δβ ± 1 • , is shown in Fig. 2(a). Since the model of Ardhuin et al. (2011) was developed for estimating seismic noise, a scaling factor, F, needed to be applied to transfer the source pressure into infrasound amplitudes. The exact value of F was not fully elaborated, however. In this study, we simply compared the observed amplitude to the modelled one on a relative scale, the scaling factor was F = 1:2000. Further developments should account for the source directivity (Brekhovskikh et al. 1973) and the source's geometric characteristics-the fact that several gridpoints contribute to the detected signals since standing ocean waves are a surface-like, rather than a point-like, independent source of infrasound (Smets & Evers 2014).
To model the amplitude, each gridpoint of the source model was assigned the azimuth, α, and distance, R, to IS26. For each α (equal to β -180 • ), the daily mean value for v eff −ratio was obtained by averaging the ECMWF HRES analyses of temperature and horizontal winds at the receiver at 00 and 12 UTC (Fig. 2b). Using R and v eff −ratio , the attenuation coefficient, A P (eq. 2 with f set to 0.2 Hz), was calculated for each potential source gridpoint; the attenuation of the maximum sources is shown in Fig. 2(c). At all gridpoints, the source pressure was then multiplied by the attributed A P and the scaling factor F to obtain the residual pressure. The expected amplitude at the receiver was estimated at the maximum of the residual pressure for each δβ ± 1 • (Fig. 2d). PMCC detections are superimposed in Fig. 2(d). Obviously, the modelled amplitudes reflect both source strength and prevailing winds (Fig. 2b) in their seasonal variations throughout the year.

U N C E RTA I N T I E S I N I N F R A S O U N D AT T E N UAT I O N M O D E L L I N G U S I N G L I DA R DATA
Overall, Fig. 2 shows reasonable agreement between the PMCC detections and the predicted backazimuths. The observed seasonal variations in amplitude are also consistent with the source and propagation models. In particular, large amplitudes and the increased number of detections in winter agree with the modelling, as well as the significantly reduced number of microbarom detections during summer. This indicates that infrasound detections originating from ocean swell in the North Atlantic are well understood with respect to infrasound propagation in the atmosphere; however, it has recently been shown that temperature and wind profiles provided by ECMWF analyses differ from measurements in the MA (e.g. Hildebrand et al. 2017). Since ECMWF analyses are commonly used for determining the effective sound speed, it is important to assess how these differences would affect microbarom modelling. This can be essential during the equinoxes and summer, when v eff −ratio is close to, or even below, 1. Investigation of such uncertainties may lead to a better understanding of outliers in the detections and, thus, false alarms at the IDC.
The Compact Rayleigh Autonomous Lidar (CORAL; for a short description, see Kaifler et al. 2017) was colocated to IS26 between May and November 2016. More than 400 hr of observations, over 83 nights, were used to quantify bias in ECMWF temperatures between 30 and 70 km altitude. All the night-time mean differences between CORAL and the ECMWF analyses are shown in Fig. 3. The overall mean shows a cold bias in ECMWF analyses from 45 km upwards, amounting to 5 K at 50 km altitude (deviations within 2σ ∼ 4.5 K) and 12 K at 58 km altitude (2σ ∼ 6 K). This is in accordance with observations made by Ehard et al. (2017) and Hildebrand et al. (2017). Note that Rapp et al. (2018) already considered monthly mean observations from the same CORAL observations and found significant biases between ECMWF analyses and both radio occultation and CORAL measurements (Rapp et al. 2018, figs 7 and 8).
The following approach estimated the effect of biases in NWP models on amplitude modelling. Considering the model output times, t mod , given at intervals of 3 hr: for each time period, t mod ± 3 h, the maximum observed amplitude in the direction β obs was  For nightly CORAL measurements between May and November 2016, the nightly mean ECMWF temperature was calculated between 18 and 06 UTC, with the lidar mean temperature corresponding to the respective observation duration. 83 nightly mean differences between ECMWF and CORAL (grey) resulted in a cold ECMWF mean bias (red) above 45 km. This increased up to 12 K at 58 km altitude, where deviations of 6 K are within the 2σ interval. compared with the maximum predicted pressure at IS26 in the direction β obs ± 5 • . Therefore, two predicted amplitude ranges were computed (note that, if no detection was found within a 6 hr time window, interpolation considered the direction of the detection that is nearest in time). (i) One deterministic amplitude value was obtained for a v eff −ratio based on the ECMWF temperature and wind profiles at the receiver (cf. Section 2). Using this value, an amplitude range was calculated in order to reflect short-time variations in source strength. Upper and lower limits were determined as 2 d sliding maxima and minima of the deterministic values, respectively.
(ii) v eff −ratio was randomly perturbed, given the range of the observed temperature deviations measured during the CORAL campaign in 2016. ECMWF wind profiles were combined with the perturbed temperature profiles incorporating the measured mean bias (up to 12 K) and deviations from ECMWF analyses (up to 6 K within 2σ ). The range of uncertainty of the modelled amplitude was obtained by first excluding amplitude deviations outside the 2σ confidence interval, and then applying the 2 d sliding maxima and minima, as per (i).
According to eq. (1), a temperature bias alters the speed of sound, c, but wind is another important term contributing to the effective sound speed. Although no obvious mean bias was found by Hildebrand et al. (2017) and Rüfenacht et al. (2018) for the horizontal wind components in the ECMWF analyses, differences within the 2σ confidence interval amounted to ±40 ms −1 . The temperature bias and deviations found using the CORAL data qualitatively agree with the ones found using the RMR lidar at the ALOMAR observatory in northern Norway in March and August 2016. Therefore, ECMWF wind differences derived from the same measurement campaign were applied to the simulations. Analogous to estimating amplitude uncertainties by temperature, as described in (ii), 100 calculations for perturbations of meridional and zonal ECMWF winds, incorporating observed 2σ deviations (±40 ms −1 ), were performed. The temperature effect and combined effect of temperature and wind uncertainties are shown in Fig. 4.

D I S C U S S I O N
Of the observed amplitudes, 76.8 per cent are in the range of the modelled amplitudes based on ECMWF analyses of wind and temperature. This rate, however, varies with season, as does the number of detections at IS26. During summer (May-August), in only 12 per cent of all time windows at least one signal is detected; 29.4 per cent of the maximum amplitudes can be predicted by solely using ECMWF analyses. During winter (November-February), these rates are significantly higher (47 and 81.4 per cent, respectively). . The modelled amplitude range, based on ECMWF data, was determined from sliding 2 d minimum/maximum deterministic values. This accounts for the fluctuation in source strength, which is also seen in the detections (black dots). Amplitude modelling uncertainties have been estimated by ECMWF temperature and wind perturbations based on CORAL measurements at IS26 and on RMR lidar data from the ALOMAR observatory in northern Norway (Hildebrand et al. 2017;Rüfenacht et al. 2018), respectively.
As seen in Fig. 4, warming the stratopause by the determined mean bias of up to 12 K, and incorporating 6 K uncertainty, weakly impacts downwind propagation during winter for arrivals from prevailing northwesterly directions. Nevertheless, amplitude predictions improve by 10.3 per cent. When combining temperature perturbations with wind deviations of up to 40 ms −1 , even 96.5 per cent of detections are predicted. Remaining discrepancies can be seen in the aftermath of SSWs, which occurred in January 2015 (Manney et al. 2015) and 2017 (Xiong et al. 2018); however, these are well resolved by the ECMWF analyses (Fig. 2b) and the predicted amplitudes, which match observations half of the seasonal mean value (16 mPa). In December 2016, large amplitudes were modelled, despite a markedly reduced number of detections (60 and 56 per cent less compared to 2015 and 2017, respectively). The uncertainty is significantly increased, so that amplitudes can range from 1 to 90 mPa, whereas ECMWF analyses only cover amplitudes between 9 and 40 mPa. The uncertainty indicates that propagation conditions might have changed, possibly because of a minor SSW.
Obviously, temperature and wind perturbations can markedly affect the propagation conditions. As a result, uncertainties in the modelled amplitude increase during the equinox periods (Fig. 4), when the mean zonal wind circulation of the stratosphere weakens and reverses. The number of time windows that show at least one detection are comparable to those during winter (41 per cent), and the amplitude prediction rates are slightly higher (+4 per cent). In summer, the low prediction rate using ECMWF analyses (29.4 per cent) is caused by amplitudes modelled below the station's noise threshold (about 2.5 mPa; see Matoza et al. 2013). These were either not detected by the PMCC algorithm or filtered out (Supporting Information Table S1); however, 12 per cent of the time windows still provide observations, exhibiting amplitudes near the threshold. These observations can only be predicted when accounting for wind and temperature deviations within the 2σ confidence interval. The prediction rate is significantly improved, to 87.2 per cent (temperature only: 40.9 per cent). This implies that atmospheric conditions had temporarily changed, weakening attenuation at the receiver, especially in 2015 (Fig. 4).
Overall, modelling the amplitudes while considering uncertainties, resulting from MA wind and temperature models, significantly improves the ability to explain microbarom amplitudes. The improvement amounts to 8.7 and 19.9 per cent when accounting for uncertainty effects of temperature and combined wind and temperature, respectively.
Considering large propagation ranges for North Atlantic sources (more than 3000 km; see Supporting Information Fig. S2), thermospheric returns are unlikely (Sutherland & Bass 2004), as opposed to arrivals from North Sea sources involving shorter propagation paths. Thus, part of the unresolved variability observed in the amplitude of the microbaroms from the North Atlantic can provide useful integrated information about dynamical properties in the stratospheric waveguide. Occasionally, at shorter ranges, microbaroms likely originate from the Mediterranean (β = 165 • )-for example, in summer 2017-despite unfavourable propagation conditions in the stratospheric waveguide. Remaining discrepancies might be caused by errors in the tabulated parameters a, b, d and s, used in eq. (2). These are mean values, determined from synthetic simulations (Le Pichon et al. 2012), and possibly do not fully describe all atmospheric states.
In this study, ECMWF wind and temperature profiles were taken at the receiver, instead of the entire propagation path (e.g. Landès et al. 2014). Since, for the first time, a lidar was located at an IMS infrasound array, this allowed direct comparison of NWP products (temperature data) with lidar profiles. Although propagation modelling might be improved by considering atmospheric conditions over the entire path, the approach used here has already provided results that are consistent with the observations (Fig. 2). Further studies should explore how far this assumption is valid within the studied region, dependent on the time of the year. Moreover, colocating a lidar to IS26 that is capable of measuring wind profiles would be beneficial in estimating the impacts of local uncertainties in ECMWF wind profiles on amplitude modelling; however, given the agreement of temperature biases found in the lidar data, using deviations from comparable wind observations at the ALOMAR observatory is the best available approximation to date.
Further investigations are also needed to more rigorously determine the scaling factor, F, of the microbarom source term in order to pursue quantitative comparisons between observations and modelling results. Additional improvements would include integrating the microbarom source term along each great circle path, assuming independent radiating monopole sources (Waxler & Gilbert 2006), rather than considering the maximum source contribution.
Despite comparable predicted amplitudes in different directions (Fig. 2d), PMCC detections do not completely cover the expected range of backazimuths. The detections are rather concentrated on specific directions. Besides the fact that atmospheric specifications are not considered over the entire propagation path, the reduced spreading of the backazimuths could be explained by intrinsic limitations of the PMCC algorithm. PMCC solely detects the most prominent coherent sources in a given time and frequency window (Cansi 1995). In future applications, alternative array processing techniques, such as the multiple signal classification algorithm (Schmidt 1986), could overcome this limitation by providing detections from different directions at the same time.

C O N C L U S I O N S A N D O U T L O O K
The combined use of the operational ECMWF HRES analyses and the IFREMER ocean wave interaction model, together with a semiempirical attenuation relation, allowed modelling of directional microbarom amplitude variations resulting from changes in stratospheric propagation conditions at infrasound station IS26. With this approach, the seasonal variation in microbarom amplitudes is well explained. In addition, accounting for known uncertainties in such NWP models, regarding MA winds and temperature, improves predictions especially during summer. In this context, using the atmospheric conditions at the receiver led to good agreement between modelling results and observations. Installing lidars along the paths between infrasound stations and potential sources would not be practicable, in terms of infrastructure requirements and cost effectiveness, due to the multitude of conceivable propagation pathways. Therefore, it can be anticipated that measuring wind and temperature profiles at IMS infrasound stations will improve the understanding of microbarom observations. Such efforts are essential in the context of the CTBT because microbaroms complicate verifications by producing high false alarm rates in routine data processing at the IDC. Better understanding the accuracy of NWP models, and how this applies to infrasound propagation in the MA, will help to reduce the high false alarm rates. Our results imply that adding a systematic MA temperature offset to the commonly used ECMWF temperature profiles through lidar measurements could improve predictions of the IMS infrasound network detection capability on shorter timescales by up to 10 per cent; however, further studies are needed, including more than one station. For more continuous operation, autonomous lidars, such as CORAL, or other technologies, such as wind radiometers (Rüfenacht et al. 2018), are highly appropriate, along with developments of lidar prototypes that are also capable of wind measurements. Capitalizing on such scientific and technical advances should reinforce the potential benefit of a routine and global use of natural infrasound for civil and scientific applications.

A C K N O W L E D G E M E N T S S U P P O RT I N G I N F O R M AT I O N
Supplementary data are available at GJ I online. Figure S1. Infrasound detections in the 0.1-0.3 Hz range at IS26. The effective sound speed ratio is colour-coded, indicating favourable (red) and unfavourable (blue) propagation conditions. Figure S2. Similar to Fig. 2, but for the entire backazimuth range, and also including the distance of the maximum sources (e). Between 30 • and 90 • , no potential microbarom sources exist in the geographic area of interest. Table S1. Filters applied in this study to consider PMCC detections with high significance only.
Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.