Explaining global patterns of microbarom observations with wave action models

Abstract

We present a methodology to model the spatio-temporal variations of microbarom detections at a global scale. Our model combines the source term resulting from the non-linear ocean-wave interaction and a simplified description of the long-range infrasound propagation through the stratospheric waveguide. We compare model predictions with observations at infrasound stations of the International Monitoring System between 2008 and 2009. Our results show a first-order consistency between the observed and modelled trends of microbarom backazimuth detections for most stations. Taking into account stratospheric wind effect on the infrasound propagation systematically improves the fit between the observations the model predictions. However, correctly predicting patterns of weekly variation of detections turns out to be more challenging and would require further improving the source and the propagation models. Short-term and regional quantitative comparisons could then be carried out based on the metrics developed in this study.

Interface waves, Wave propagation, Acoustic properties

1 INTRODUCTION

Modern seismological and infrasound networks produce large amounts of continuous records that are dominated by background noise which has strong amplitudes near 0.15–0.2 Hz. The large amplitudes of background seismic and atmospheric waves, secondary microseisms and microbaroms, are generated by the interaction of ocean gravity waves with the seafloor and the atmosphere, respectively (e.g. Donn & Naini 1973; Rind 1980; Hedlin et al. 2002; Bowman et al. 2005) caused by the non-linear interference of oceanic waves with the same frequency propagating in opposite directions (Longuet-Higgins 1950). Microbaroms generation is directly proportional to oceanic waves interaction (Longuet-Higgins 1950; Waxler & Gilbert 2006; Ardhuin & Herbers 2013). In the case of microseisms, there is an additional term related to the coupling effect with bathymetry (e.g. Kedar et al. 2008; Ardhuin et al. 2011).

Applications based on the ambient seismic field have been developed during the last decade (e.g. Campillo et al. 2011a,b), leading to new methods of passive seismic imaging that were applied at different scales (e.g. Shapiro et al. 2005; Ritzwoller et al. 2011; Mordret et al. 2013) and for monitoring temporal changes in the crust (Brenguier et al. 2008a,b). For applying these methods, distribution of noise sources is a critical issue (e.g. Froment et al. 2010). Studies targeting the origin of microseisms have a long history (e.g. Hasselmann 1963; Webb & Cox 1986; Webb 1992). Their global source distribution is reasonably well modelled by coupling the Longuet Higgins theory to numerical wave action models (e.g. Kedar et al. 2008; Ardhuin et al. 2011). The resulting patterns are generally compared with observations of excited surface and body waves (e.g. Stehly et al. 2006; Stutzmann et al. 2009; Landès et al. 2010; Hillers et al. 2012).

Recent developments of infrasound networks at global and continental scales facilitate the analysis of acoustic waves for probing atmospheric structure and providing new insights into the physics of large atmospheric phenomena (e.g. Hedlin et al. 2012). Both, localized events (e.g. Walker et al. 2011) and background noise (e.g. Le Pichon et al. 2006) provide useful data set for inversion studies. Similar to seismology, atmospheric noise has a significant potential advantage because they are continuously produced and recorded allowing interferometric approaches (Haney 2009) to be applied for sounding the middle and upper atmosphere.

The detection capability of the International Monitoring System (IMS) network is a key concern for the Comprehensive Nuclear-Test-Ban Treaty (CTBT). The microbaroms interfere with the detection of explosive events (e.g. Evers & Haak 2001; Stevens et al. 2002), so their understanding is essential for the treaty verification. However, observations of microbaroms strongly depend on the stratospheric wind and temperature structure along the propagation path (e.g. Garcés et al. 2004; Le Pichon et al. 2006; Brachet et al. 2010). The propagating medium is more dynamic and varies on much shorter time scales compared to seismology. Therefore, the spatial and temporal variability of the atmosphere complicate microbarom source studies using remote observations (Drob et al. 2003). Using 5 yr of continuous records from 40 infrasound stations of the global IMS network, Landès et al. (2012) provided averaged seasonal patterns of global microbarom source locations. Although the source mechanism of the microbarom emission as described by Waxler & Gilbert (2006) is well accepted, modelling global microbarom observations remains challenging. At local scales and for specific events like Hurricanes, numerical models have been proposed to explain observations (e.g. Stopa et al. 2011). In our study, we model global microbarom observations with specific interest to the seasonal variations of backazimuths and number of detections. Following the source mechanism of Waxler & Gilbert (2006), we consider the wave–wave interaction model from Hillers et al. (2012) together with a simplified propagation model that uses wind forecasts from the ECMWF atmospheric circulation model (http://www.ecmwf.int). Simulations are then systematically compared with observations at all considered IMS stations.

2 OBSERVATIONS OF MICROBAROMS

The IMS includes a global network of infrasonic stations designed to detect atmospheric explosions to verify the CTBT We focus on microbarom detections at 40 IMS stations between 2008 January and 2009 December (Fig. 1). The spatial distribution of these sites covers a wide range of latitudes, atmospheric and oceanic conditions. All infrasound stations are composed of four or more microbarometers and include wind-filtering systems and communication facilities (Christie & Campus 2010). To detect coherent plane waves crossing the array, we use the Progressive Multi-Channel Correlation (PMCC) algorithm (Cansi 1995) with a configuration to 15 log-scaled frequency bands between 0.01 and 5 Hz (Le Pichon et al. 2010; Matoza et al. 2013). A grid search is performed over successive overlapping time windows and frequency bands. A coherent arrival in a particular time window and frequency band is registered as a ‘pixel’. Pixels with similar wavefront properties (frequency, time, azimuth and trace velocity) are then grouped into detections, so-called ‘families’. We select detections with an average frequency ranging in the 0.1–0.4 Hz microbarom band. Nine stations (IS11, IS14, IS43, IS45, IS50, IS51, IS56 and IS57) with poor microbarom detections (noisy sites, instrumentation failures…) are not considered for this study.

Figure 1.

Location map of the 40 IMS stations used in this study.

Fig. 2 shows the azimuthal distributions of detections at two stations (IS05 and IS26) averaged over 1 week periods. Observation during 2 yr reveals clear seasonal patterns. For example, at station IS05, in austral summer, the dominating azimuth is ∼125°, while in winter it switches to ∼225°. Similar seasonal variations in dominating azimuths are observed at most stations (Fig. 3). Moreover, opposite dominating azimuths are observed in the northern and southern hemispheres. This observed seasonal pattern reveals seasonal changes in the microbarom source intensity, and is caused by the seasonal variations of the prevailing stratospheric winds (e.g. Le Pichon et al. 2006). Microbarom detections also contains signals propagating through the thermosphere (e.g. Garcés et al. 2004). In this study, since we essentially focus on first-order seasonal variations and consider long propagation range, thermospheric propagation is not taken into account in our simplified modelling.

Figure 2.

Distribution of observed microbarom detections at stations IS05 (a) and IS26 (b). The grey areas indicate observation periods without data. The colour scale corresponds to the logarithm of the number of detections. Left-hand small panels plot the normalized azimuthal distribution of detections.

Figure 3.

Weekly distribution of the number of detections at 33 IMS stations between 2008 and 2009. Numbers near the station name correspond to their latitudes. Stations are sorted by latitude. The colour corresponds to the weekly averaged azimuth. The height of bars indicate detection number (log scale). For each station, the upper limit of the y-axis corresponds to 15 000 detections.

We follow the approach proposed by Landès et al. (2012) and average backazimuths of microbarom detections to estimate dominant source regions (Figs 4a and b). The reconstructed source regions are compared with maps highlighting regions where oceanic waves interact according to the Longuet-Higgins's theory [Figs 4c and d, model from Hillers et al. (2012) described in the next section]. We note that the observed and predicted source regions are consistent and follow the global atmospheric circulation model.

Figure 4.

Microbarom source location based on cross-bearing azimuth of PMCC detections for January and July (Landès et al. 2012) (a, b). Wave–wave interaction terms averaged between 0.1 and 0.3 Hz for 2008 January and 2008 July (c, d).

3 PREDICTING THE NUMBER OF MICROBAROM DETECTIONS

We propose a statistical approach to model the temporal and azimuthal dependence of the number of microbarom detections. In order to capture variations from seasonal down to weekly scales, the microbarom source model is defined with an 1-d temporal resolution. Our model uses the daily averaged ECMWF wind model and the wave interaction model developed by Kedar et al. (2008). Following the global scale approach of Hillers et al. (2012), this model estimates the directional spectral density function F(f, θ) of the wave interaction intensities using Wave Action Model (NOAA Wavewatch III, http://polar.ncep.noaa.gov/waves). The Wavewatch III results cover latitudes between −78° and 78° on a 1.25° × 1° longitude–latitude grid with a directional resolution of 15°. There are 25 logarithmically spaced frequency bands between 0.362 and 0.037 Hz (3.3–27.2 s). Sea ice data were acquired from the NOAA assimilation system archives. We use a 3-hr temporal resolution between 2008 January and 2009 December. For comparison with observations, the interaction terms are averaged every day. This model was used to evaluate the consistency between the observed and predicted excitation patterns (Hillers et al. 2012). As opposed to other more recent source models (e.g. Ardhuin et al. 2011), wave reflections along coastlines are not here taken into account.

The microbarom source mechanism is described by the nonlinear interaction of oceanic waves with the same frequency travelling in opposite directions (Longuet-Higgins 1950). Waxler & Gilbert (2006) applied this theory and described the interaction term at the ocean–atmosphere interface. In addition, Ardhuin & Herbers (2013) provide a complete theory of noise generation in the solid Earth, oceans and atmosphere and find similar equations with no bathymetry dependence for the microbarom generation. The wave–wave interaction ϕ at a source point r_s for a frequency f is given by:

\begin{eqnarray} \phi (r_{\rm s},f) \ = \ \int _0^{\pi } F\left(r_{\rm s},\frac{f}{2},\theta \right) F\left(r_{\rm s},\frac{f}{2},\theta +\pi \right)\mathrm{d}\theta , \end{eqnarray}

(1)

where F is the directional spectral density function of oceanic waves, f denotes the frequency and θ the azimuth. We compute daily averages of ϕ between 0.1 and 0.3 Hz to obtain the microbarom source density P(r_s) (Figs 4c and d). With this description, a surface element dS(r_s) emits isotropically infrasound with intensity proportionals to P(r_s) dS (r_s).

We adopt a simplified approach to model infrasound propagation through the atmosphere assuming that microbarom signals mainly propagate along great circle paths in the stratospheric waveguide. We neglect thermospheric paths for two reasons: (1) attenuation of thermospheric phases becomes large for propagation distances greater than 1000–2000 km and frequencies near 0.2 Hz (Le Pichon et al. 2012) and (2) we focus our modelling on the first-order seasonality of microbarom detections caused by the stratospheric wind reversal.

We introduce the propagation term Γ_w(r_s, r_r) indicating the presence of a stratospheric waveguide from the ground to 50 km altitude along the great circle connecting the source (r_s) and the receiver (r_r). The efficiency of the waveguide is described by the dimensionless parameter C_eff-ratio defined by the ratio between the effective sound speed at 50 km altitude (defined from the daily averaged ECMWF wind model) and the sound speed at the ground level. Near 0.2 Hz, long range propagation is favoured when C_eff-ratio > 0.9 (e.g. Le Pichon et al. 2012). At a distance of 1000 km, the attenuation is about 60 dB considering a reference distance of 1 km from the source. Here, propagation condition between r_s and r_r is approximated by

\begin{eqnarray} \Gamma _{\rm w} (M,S) = \left\lbrace \begin{array}{ll} 1 &\quad \textrm {if}\ C_{{\rm eff-ratio}} > 0.9\ \textrm {for each point along the}\\ &\qquad \textrm{great circle path between}\ S\ \textrm {and}\ M.\\ 0 &\quad \textrm {otherwise} \end{array} \right. \nonumber\\ \end{eqnarray}

(2)

We then incorporate attenuation coefficients derived from multiple numerical simulations using parabolic formulation (Le Pichon et al. 2012). At large distances, far beyond the first stratospheric bounce, the following empirical expression is used for the energy attenuation:

\begin{eqnarray} \Gamma _{\rm a} (r_{\rm s},r_{\rm r}) \ = \ \left( \frac{r_0}{\Delta _{{\rm sr}}} \right) ^ {2\beta }, \end{eqnarray}

(3)

where r₀ is a reference distance equal to 1 km, Δ_sr is the distance between the source and the receiver, and the transmission loss parameter β is set to 0.8 in the 0.1–0.3 Hz frequency band.

Combining the source and propagation terms, the microbarom energy P(r_s, r_r) emitted by a surface element dS(r_s) and received at r_r is given by:

\begin{eqnarray} P(r_{\rm s},r_{\rm r}) \ = \ \Gamma _{\rm a} (r_{\rm s},r_{\rm r}) \Gamma _{\rm w} (r_{\rm s},r_{\rm r}) P(r_{\rm s}) {\rm d}S(r_{\rm s}). \end{eqnarray}

(4)

To compare model predictions and observations, we consider that the number of microbarom detections at stations r_r from the source region r_s is proportional to the radiated energy P(r_s, r_r). Another assumption is that no signal can be detected from sources whose intensity is below the threshold P_t. This parameter is fixed to P_t = 10⁻⁴ Pa² Hz⁻¹ at all stations. It represents the minimum power spectral density level from which recorded microbarom signals are considered as detections. Further studies could improve this approximation by taking into account the minimal noise level observed at each station.

For a given station , the sum of contributions from all sources within the azimuth range [θ − 4°, θ + 4°] and with P(r_s, r_r) > P_t is a proxy for the number of detections associated with θ.

Repeating this computation for the range of azimuths and days yields an azimuth-time distribution of number of detections N(θ, t) (where t is time) as shown in Fig. 5. We compute these functions for the 33 selected IMS stations between of 2008 and 2009. The method used for a quantitative comparison between model predictions and observations is described in the following section. Note that the value of P_t is chosen after trials taking into account the trade off between detections of dominant source regions and multiple detections (possibly noisy) of weaker source regions.

Figure 5.

Distribution of predicted microbarom detections at stations IS05 (a) and IS26 (b). The colour scale corresponds to the logarithm of the number of detections. Left-hand small panels plot the normalized azimuthal distribution of detections.

4 COMPARING OBSERVATIONS AND PREDICTIONS

The qualitative analyses presented in Figs 3 and 6 indicate a good consistency between the modelled and observed annual trend of microbarom detections, more especially in the southern hemisphere. For stations IS05 and IS26, Figs 2 and 5 present consistent incoming arrivals of ∼200–250° in the boreal summer season, and ∼50–150° in winter. However, for some stations, the agreement is less successful. For example at station IS30 (Fig. 7), the model does not predict correctly signals with backazimuths between 170° and 250°, corresponding to sources in the Indian Ocean.

Figure 6.

Weekly azimuth of the 33 selected IMS stations between 2008 and 2009 obtained from the modelling. Numbers near the station name correspond to their latitudes. Stations are sorted by latitude. The colour corresponds to the weekly averaged azimuth. The figure has to be compared with the Fig. 3 associated to observations.

Figure 7.

Comparisons between observed (a) and predicted (b) microbarom detections at station IS30. The grey areas indicate observation periods without data. The colour scale corresponds to the logarithm of the number of detections. Left-hand small panels plot the normalized azimuthal distribution of detections.

To quantify the accuracy of microbarom predictions at different stations, we compute correlations between the observed and the predicted seasonal patterns (such as shown in Figs 2 and 5). C_az is defined as the correlation coefficient of the marginals of the number of detections along the backazimuths axis.

\begin{eqnarray} M^{\rm obs}_{\rm az}(\theta ) &=& \int N^{\rm obs}(\theta ,t) \,\mathrm{d}t \nonumber\\ M^{\rm pred}_{\rm az}(\theta ) &=& \int N^{\rm pred}(\theta ,t) \,\mathrm{d}t \nonumber\\ C_{\rm az} &=& C_{\rm corr} \big[M^{\rm obs}_{\rm az}(\theta ), M^{\rm pred}_{\rm az}(\theta ) \big], \end{eqnarray}

(5)

where C_corr denotes the correlation coefficient, indexes ^obs and ^pred refer to the observed and predicted values, respectively.

We further introduce S_corr to compare the temporal variations of the azimuthal distribution between the observed and predicted microbarom signals:

\begin{eqnarray} S_{\rm corr} = C_{\rm corr} [ N^{\rm obs}(\theta ,t), N^{\rm pred}(\theta ,t) ]. \end{eqnarray}

(6)

We compute C_az and S_corr for all stations (Fig. 8). The overview of C_az (Fig. 8a) shows a good correlation between the main observed and predicted backazimuths, except for stations IS27, IS34 and IS39. The distribution of S_corr is characterized by lower similarities (Fig. 8b). For stations where S_corr is lower that 0.3, the modelling fails to explain the temporal variation of the observed detections.

Figure 8.

Distribution of C_az (a) and S_corr (b) for 33 IMS stations.

To evaluate the wind effect in the model, we compute the azimuth-time distribution of detections N(θ, t) without the Γ_w term (eq. 2) and we compare the new distributions C_az and S_corr with those previously computed (Fig. 9). We can see that introducing stratospheric wind effects significantly improves predictions at almost all stations. Fig. 10 shows another comparison displaying the global azimuthal distributions of microbaroms and emphasizing the importance of the wind term in the modelling. It can be seen, that without accounting for wind effects, many modelled sources are not observed. Introducing the wind term corrects this strong discrepancies in most cases.

Figure 9.

Comparison of C_az (a) and S_corr (b) for 33 IMS stations with and without propagation effects.

Figure 10.

Normalized azimuthal distributions M_az(θ) of microbaroms for 33 IMS stations. Observations are in grey. Predictions with and without including propagation effects are superimposed (red and green lines, respectively).

5 DISCUSSIONS AND CONCLUDING REMARKS

We have developed a methodology to model the global temporal variations of microbarom detections. It considers both microbarom source energy resulting from ocean wave interaction models (e.g. Hillers et al. 2012) and long-range infrasound propagation through the stratospheric waveguide. Despite its simplicity, this approach explains to first order the dominant trends of microbarom observations. The correlation between the observed and the modelled azimuthal distributions is generally good, more specifically for middle latitude stations where stratospheric wind prevails. We have also shown that taking into account the effects of the stratospheric winds on the propagation systematically improves the agreement with observations.

Discrepancies observed for some stations in Fig. 8 may have multiple origins. First, the accuracy of the wave action model is limited. In particular, the model is not defined in polar latitudes. Also, it does not account for coastal wave reflections that can cause strong oceanic wave interference and microbarom generation in some areas. Another source of discrepancies can be explained by our simplified propagation modelling which ignores thermospheric phases. This effect becomes non negligible for stations located close to source regions where thermospheric attenuation is relatively low. Finally, the ECMWF wind model used in this study has known uncertainties due to the lack of assimilated measurements in the stratosphere (e.g. Charlton-Perez et al. 2013) and is still a research topic for ongoing international projects (e.g. SPARC http://www.sparc-climate.org/ and ARISE http://arise-project.eu).

Improvements of the microbarom predictions at smaller temporal scales (from weekly to daily fluctuations) would require improved source and propagation models. The description of the infrasound propagation through the atmosphere used in our study is approximate and all terms from eq. (4) can be formulated more accurately. Also, we expect that tropospheric and mesospheric winds would also play a non negligible role, especially at high latitude and equatorial regions. One of the main shortcoming of the used source model is non-accounting for coastal reflections. This limitation has already been mentioned by Hillers et al. (2012). As shown by Ardhuin et al. (2012) and Obrebski et al. (2012), wave interaction induced by coastal reflections may play an important role in the generation of microseisms and we expect a similar effect for microbaroms. Considering the source models that include coastal reflections (e.g. Ardhuin et al. 2011) would improve microbarom predictions in the future. The accuracy of these models could then be assessed by comparing observations with modelling results using metrics similar to those developed in our study.

We would like to thank Fabrice Ardhuin for his useful comments and discussions. The research leading to these results has partly been performed within the ARISE project (http://www.arise-project.eu) and was supported by funding from the FP7 ERC Advanced grant 227507 (WHISPER), by the European Community's Seventh Framework Programme (FP7/2007-2013) under grant agreement 284387, and from the CEA, DAM, DIF, France.