A Three‐Dimensional Array for the Study of Infrasound Propagation Through the Atmospheric Boundary Layer

The Royal Netherlands Meteorological Institute (KNMI) operates a three‐dimensional microbarometer array at the Cabauw Experimental Site for Atmospheric Research observatory. The array consists of five microbarometers on a meteorological tower up to an altitude of 200 m. Ten ground‐based microbarometers surround the tower with an array aperture of 800 m. This unique setup allows for the study of infrasound propagation in three dimensions. The added value of the vertical dimension is the sensitivity to wind and temperature in the atmospheric boundary layer over multiple altitudes. In this study, we analyze infrasound generated by an accidental chemical explosion at the Moerdijk petrochemical plant on 3 June 2014. The recordings of the tower microbarometers show two sequential arrivals, whereas the recordings on the ground show one wavefront. This arrival structure is interpreted to be the upgoing and downgoing wavefronts. The observations are compared with propagation modeling results using global‐scale and mesoscale atmospheric models. Independent temperature and wind measurements, which are available at the Cabauw Experimental Site for Atmospheric Research, are used for comparison with model output. The modeling results explain the signal arrival times; however, the tower wavefront arrivals are not explained. This study is important for understanding the influence of the atmospheric boundary layer on infrasound detections and propagation.


Introduction
Infrasound, or acoustic waves with frequencies below 20 Hz, is used for monitoring natural and anthropogenic sources (Campus & Christie, 2010). These atmospheric sources displace large volumes of air, generating infrasound. Furthermore, infrasound is one of the techniques used for the verification of the Comprehensive Nuclear-Test-Ban Treaty (CTBT; Dahlman et al., 2009).
Microbarometers measure infrasound. Generally, they are set up in an array to determine source characteristic parameters: the direction and velocity of the wavefront. A global network of microbarometers is part of the International Monitoring System, which is used for the verification of the CTBT. The infrasound detections from the International Monitoring System are also used for research on, for example, volcanoes (Dabrowa et al., 2011), explosions (Christie et al., 2005), and meteors (ElGabry et al., 2017).
The atmosphere plays a central role in the study of infrasound. Local wind conditions determine the noise conditions near an infrasound microbarometer and therefore the detectability (Scott & Raspet, 1992). Moreover, the propagation conditions are determined by the state of the atmosphere. Infrasonic signals may propagate over long ranges due to low attenuation and the presence of atmospheric waveguides (Drob et al., 2003). The refraction of infrasonic waves, and therefore the formation of waveguides, is governed by vertical gradients in temperature and wind. Infrasonic waves refract upward in the presence of negative vertical gradients and vice versa. Throughout the atmosphere, one can distinguish tropospheric, stratospheric, and thermospheric waveguides (e.g., Evers et al., 2012;Smets et al., 2016;Waxler & Assink, 2019).
In this study, the focus is on local (tens of kilometers) propagation conditions in the atmospheric boundary layer (ABL). Near-surface (<1-km altitude) atmospheric variability has a large impact on local infrasound propagation. The ABL is the lowermost part of the troposphere whose structure is directly influenced by radiative heating and cooling of the Earth's surface (Stull, 1988;Wilson, 1996). Figure 1 presents the main features of the ABL during day (convective mixed layer) and night (stable nocturnal boundary layer) time. During daytime convection mixes the ABL. The ABL is then turbulent, which increases the wind noise levels. During nighttime turbulence is suppressed. Due to a laminar flow and lower wind noise levels, the signal-to-noise ratio (snr) of infrasound observations increases. The interplay between daytime upward refraction and nocturnal downward refraction has influence on signal level. Downward refraction increases the snr. Fee and Garcés (2007) show that ABL changes affect ambient noise levels, the snr, and infrasound detectability. Previous studies focused on low-frequency audible acoustic waves in relation to the ABL (e.g., Blom & Waxler, 2012;Talmadge et al., 2008;Waxler et al., 2006Waxler et al., , 2008. The local and regional variability of infrasound due to small-scale ABL variations is not fully examined. For practical reasons microbarometers are typically installed on the surface as two-dimensional arrays. This study focuses on improving the understanding of the influence of the ABL on infrasound detectability, using a three-dimensional infrasound array and ray propagation modeling. This array is located at the Cabauw Experimental Site for Atmospheric Research (CESAR) observatory (CESAR Observatory, 2019). The CESAR Infrasound Array (CIA) has a unique tower up to 200-m altitude. The tower microbarometers are depicted in Figure 1.
Infrasound signals are used from the Shell chemical plant explosions on 3 June 2014 near the city of Moerdijk in the Netherlands (Figure 2). Infrasonic arrivals at 40 km from the source indicate either tropospheric/ABL arrivals or the onset of a direct wave propagating to higher altitudes. The first reactor explosion occurred 20:48:26 UTC (local time at Moerdijk is UTC+2) and a second vessel explosion 20 s later at 20:48:46 UTC (Joustra et al., 2015) under evening ABL conditions. In the ABL a positive temperature gradient is expected  during the night and a negative gradient during daytime. The signals of the first explosion are studied to show the first three-dimensional analysis. Turbulence decreased since the Moerdijk explosions occurred 1 hr after sunset. This leads to favorable detection conditions during the Moerdijk explosions.
In this study, three-dimensional observations of regional infrasound signals are presented, showing the added value of recording infrasound in three dimensions. The raypaths and arrival times are modeled to explain the tower and ground recordings.
The paper is organized as follows. In section 2 the data acquisition, array setup and array processing will be described. In section 3 the two-dimensional array analysis of the Moerdijk explosion is discussed. The tower data will be evaluated and an analysis of the tower data will be given using a frequency-wavenumber (F-K) spectrum. Furthermore, wavefronts are modeled using ray theory as a first approximation, deriving travel times and raypaths. The weather models and ray propagation modeling results are described in section 4. Finally in section 5, a second event is discussed, and the paper research findings are discussed and summarized.

Three-Dimensional Setup
In this study, a three-dimensional infrasound array up to 200-m altitude is used. Table 1 shows the locations of the microbarometers. The location of the explosions near Moerdijk is at 51.68579 • N and 4.56430 • E ( Figure 2). Ten microbarometers are placed on the Earth's surface and five on the tower. The tower microbarometers are located at an altitude of 60-, 100-, 140-, 180-, and 200-m altitude (Table 1). A schematic array setup is shown in Figure 3. During the time of the Moerdijk explosions, all 10 ground microbarometers were operational. To optimize the detections for higher frequencies, a smaller aperture array was chosen (CIA01-CIA06), located closest to the tower (Campus & Christie, 2010). The microbarometers are developed at the KNMI (Royal Netherlands Meteorological Institute) and able to resolve pressure fluctuations over a wide dynamic range. After correcting for the instrument response, larger meteorological pressure disturbances with amplitudes in orders of hectopascals can be resolved. The differential pressure microbarometers are capable of measuring frequencies between 0.001 and 20 Hz (Mentink & Evers, 2011).
The wind leads to pressure disturbances that are incoherent over a small area when compared to acoustics. This property is used to filter out wind-induced noise by sampling over a larger area. Porous hoses are used for the ground microbarometers as a noise reduction technique (Hedlin et al., 2003). The tower The angle between the apparent velocity (c app ) and the group velocity (c g ) is indicated by , the incidence angle.
microbarometers have been mounted without porous hoses. The microbarometers sample the pressure field at 40 Hz.

Array Processing
The array processing of infrasonic recordings uses the signal coherency, in order to extract events of interest. Two processes are part of data processing: signal detection and plane wave parameter estimation using beamforming. The coherency between the recordings in the time domain (Melton & Bailey, 1957) and frequency domain (Smart & Flinn, 1971) can be quantified using the Fisher ratio (F). The probability of a detection can be estimated through the framework of Fisher statistics (Evers, 2008). It follows that the snr relates to the Fisher ratio (F) through Here, N indicates the number of microbarometers. The waveforms are preprocessed before beamforming. The data are detrended and band-pass filtered between 0.4 and 10 Hz. This filter is chosen to avoid interference of low-frequency coherent noise from, for example, microbaroms (∼0.2-Hz) and for improving signal detections of interest. In the time domain, resampling is done to 500 Hz to permit a more precise fitting to the slowness grid. The source parameters are estimated using a plane wave model to characterize the wavefronts, that is, back azimuth and apparent velocity. A plane wave is described by an incidence angle ( ) and azimuth clockwise with respect to the north ( ) (Figure 3). Wavefront propagation is described by a vector normal to the wavefront and is called the slowness vector ⃗ p = (p x , p , p z ).
The back azimuth is indicated by b = + 180 • and directs to the event origin. The apparent velocity (c app ) is the horizontal component of the group velocity (c). The group velocity is the velocity with which the envelope of the wave propagates. |⃗ p| H is the horizontal component of the slowness vector and is the reciprocal of the apparent velocity. The group velocity and the apparent velocity are related through the incidence angle with respect to the horizontal (equation (3)).
The recordings are split in time segments of 5 s with 90% overlap. A horizontal slowness grid is used. For every grid point the recordings are delayed and summed. For each slowness beam (p x , p y ) and time segment, the Fisher ratio is evaluated. The slowness beam with the highest Fisher ratio is stored for event analysis.
The slowness grid is arranged for apparent velocities and back azimuth values of interest: apparent velocities range between 300 and 500 m/s, with steps of 1 m/s, and back azimuth values between 100 • and 300 • (resolution of 1 • ), for the time domain analysis. The frequency domain analysis is focused on the event; therefore, the parameters are chosen closely to the events parameters retrieved from the time domain analysis: a back azimuth range (217-219 • with steps of 1 • ), frequencies between 0.4 and 10 Hz (band stepping of 0.01 Hz and band averaging of 0.1 Hz), and an apparent velocity range (330-400 m/s, with 1-m/s steps) is used.
As mentioned above, the ground array is used for the beamforming technique. To analyze the tower observations a F-K analysis is applied, using the vertical components of the wavefield. This analysis gives a solution of the complete slowness vector. The method is fully explained in, for example, Capon (1969) and Denholm-Price and Rees (1999). A planar wavefield is assumed.
A F-K spectrum is calculated by applying a two-dimensional Fourier transform and is presented in equation (4).
with frequency f , travel time t, and the vertical wavenumber k z at position z. The wavenumber vector ⃗ k gives the spatial vector perpendicular to the wavefront. The F-K analysis is applied on band-pass-filtered tower data (between 0.4 and 2 Hz) assuming interelement distance of 40 m, in the frequency domain. Negative k z values in the F-K spectrum indicate downgoing wavefronts and upgoing wavefronts have positive k z values.
The slowness vector relates to the wavenumber vector by the angular frequency ( ): ⃗ k = ⃗ p · = ⃗ p · 2 . Substituting these equations in equations (2) and (3) gives Fitting equation (6) to the F-K spectrum, the incidence angle ( ) on the tower sensors can be estimated.
Here, the group velocity equals the adiabatic sound speed (340 m/s), since the wind speed is close to 0.

Results: Array Processing and Waveform Analysis
This section discusses two-dimensional and three-dimensional array analyses of the Moerdijk explosion signals. The differences between ground observations and tower observations are described in section 3.1, and time and frequency domain beamforming results are discussed in section 3.2.

Infrasound Observations
The signals measured on the ground (lower six traces in Figure 4) show similar waveforms with peak overpressure of ±12 Pa in the raw data and around 3-Pa peak amplitude for the band-passed-filtered (0.4-2 Hz) data. The data are filtered for the frequency range with the highest-energy content (section 3.2). A frequency difference can be noted between the tower and ground microbarometers, which is likely related to the high-frequency attenuation due to the tower (see also Figure S1 in the supporting information for spectra of the data). The effect of the tower on the pressure observations will be studied in the future. The tower wavefronts arrive in arrival pairs. This is interpreted as wave interference of two wavefronts and plays an important part in the tower observations. The superposed wavefronts will be discussed in sections 3.2 and 5.  In total, about 120 s of coherent infrasound is detected from the direction of Moerdijk. This is longer than the time of the two explosions. Smaller impulsive events are observed with peak amplitudes of about 0.13 Pa (Figure 5b). A possibility, interpreted from the apparent velocities, is that the smaller high signal-to-noise periods are shorter explosions and the long-lasting coherent energy is due to the fire. Furthermore, multipathing can be part of the signal with higher apparent velocities.

Array Analysis
The tower observations are different from the ground observations. A hypothesis to explain the tower observations by wave interference is that two wavefronts are propagating along the tower microbarometers. To test this hypothesis, a model of two superposed waveforms recorded by the ground array is fitted to the real data ( Figure 6a).
The best fit between the real data and the simulated waveforms is when two ground waveforms are added with delay times of 0.975, 0.65, 0.47, 0.179, and 0.1 s for stations CIA11-CIA15, respectively. This shows a decrease in delay time with increasing altitude.
Furthermore, the frequencies and vertical wavenumbers are presented in a F-K spectrum (Figure 6b). Upgoing and downgoing wavefronts are observed with wavelengths between ∼100 and 1,000 m. The downward directed wavefront has slightly higher amplitudes, which implies a convergence of waves moving downward (e.g., Waxler et al., 2008). A line is fitted through the spectrum, using equation (6). The best fit is for a ±9 • incidence angle ( ). Using the apparent velocity from the ground array analysis ( Figure 5) and equation (3), the incidence angle is arccosarccos(340∕344) = 8.7 • , which is similar to the incidence angle determined from the tower observations. In conclusion, using the F-K analysis, the incidence angle can be estimated for an incoming, refracted, or reflected waveform. The third dimension in the array analysis provides insight in the wavefront orientation when arriving at the array.

Ray Propagation Modeling and Numerical Weather Prediction Models
In this section, ray propagation modeling based on ray theory is briefly introduced. Geometrical acoustics is used to understand the sequential arrivals recorded in the tower. Section 4.1 introduces the atmospheric weather models and independent wind and temperature measurements. In sections 4.2 and 4.3, ray theory is explained and the approximations are discussed. Furthermore, the raypaths are shown by making use of various weather prediction models.

Atmospheric Specifications and Independent Measurements
In this study, three weather models are used for ray propagation modeling. The European Centre for Medium-Range Weather Forecasts (ECMWF) high-resolution model, the High Resolution Limited Area Model (HIRLAM), and the HIRLAM ALADIN Research on Mesoscale Operational Numerical Weather Prediction in Euromed (HARMONIE) model are used. These models are a smooth and deterministic representation of the true atmosphere, assimilating various ground-and satellite-based measurements (Bauer et al., 2015). Numerical weather prediction models have a horizontal resolution on the order of a few kilometers. This means that small-scale atmospheric variations are not resolved (Jakob, 2010). The small-scale structure in the atmosphere can be relevant, as it contributes to phenomena such as shallow waveguide formations and scattering of infrasound (Norris et al., 2009). This has been studied previously on local and regional scales (Kim et al., 2018).
Sound propagation depends on the temperature and wind conditions (Drob et al., 2003). These parameters can be obtained from these weather models. ECMWF's operational forecasts (Version Cycle 38r2) are part of the Integrated Forecast System (European Centre for Medium-Range Weather Forecasts, 2018). Throughout this paper, these models are further referred to as ECMWF, HIRLAM, and HARMONIE. Differences in parameterization and assimilation may result in different profiles for the same location and time. Hydrostatic models assume that vertical accelerations are negligible. ECMWF and HIRLAM are hydrostatic. HARMONIE is nonhydrostatic; therefore, it can resolve nonhydrostatic features, such as thunderstorm updrafts, downdrafts, and deep convection. A summary of important model characteristics is given in Table 2.
An analysis product is not produced by data assimilation alone. The observations are also used to correct for errors in the short forecast of the previous analysis outcome. HIRLAM is a regional model and uses the boundary conditions of ECMWF (e.g., Undén et al., 2002). For HARMONIE, a special HIRLAM model run transforms hourly boundary conditions from ECMWF to hourly boundary conditions for HARMONIE. The horizontal resolution is 2.5 km (Bengtsson et al., 2017). HARMONIE is therefore permitting smaller convection scales.
In addition to the weather models, independent temperature and wind speed measurements are analyzed at the same location as CIA. Kim et al. (2018) show that in situ measurements (radiosonde sounding data) could provide a more accurate prediction, for example, a waveguide was pronounced in the radiosonde data that was not pronounced in the weather prediction model.
The tower measurements, up to 200-m altitude, are closest to the real atmosphere. Every 10 min the in situ average wind speed and temperatures are measured on the tower. The wind speed values are measured by a cup anemometer that samples at 4 Hz. The temperatures are measured by a Pt 500 element that samples every 12 s. For higher altitudes, a LAP3000 wind profiler is used for the wind speed measurements. Temperature profiles are retrieved from the microwave radiometer (Rose et al., 2005). The LAP3000 wind profiler measures the wind speed and direction from 137-m altitude up to an altitude of about 7 km. The wind profiler data are postprocessed, and they represent a 30-min average (20:30-21:00 UTC; Gaffard et al., 2006). The absolute temperatures are measured every half hour up to an altitude of 10 km and represent 5-min averages (measured 20:50 UTC). These independent measurements are compared to the numerical weather models.

Propagation Modeling
In a horizontally layered atmosphere, ray propagation can, as a first approach, be described by the effective sound speed (Godin, 2002). For local propagation and horizontally propagating waves, the effective sound speed (c eff ) in direction of propagation can be used instead of the group velocity. The effective sound speed is defined as a combination of the sound speed (c T ) with the along track (w a ) wind speed (c eff = c T + w a ),' where c T = √ RT(z) and w a = ⃗ u·n. is the ratio of specific heat and R is the gas constant. ⃗ u·n represents the horizontal wind speed (⃗ u) in the propagation directionn. Here, the propagation direction is approximately 38.4 • with respect to the north. Upward refraction occurs when the vertical effective sound speed gradient is negative and refraction back to the surface (downward refraction) may occur from regions where the effective sound speed becomes larger than its surface value due to a positive gradient. In the latter case, the waves are trapped in between the ground and a layer aloft.
Given an atmospheric state, sound propagation can be simulated by the infrasound wave equation. As a first approach for the description of a sound ray, the focus is on the ray theoretical solution of the infrasound wave equation, also called the geometrical acoustic approximation. This is the high-frequency asymptotic solution to the wave equation. The derived Eikonal equation is used for computing raypaths and arrival times of the first arrival. Details are available in literature (Blom & Waxler, 2017;Brekhovskikh & Godin, 1999;Pierce & Smith, 1981). Ray theory does not describe interference and diffraction from small-scale structures but can be used to understand propagation conditions to the first order. Raypaths and travel times between Moerdijk and CIA are simulated by a three-dimensional ray tracer using a one-dimensional atmosphere. A review of numerical methods for prediction of infrasonic signal characteristics is given by Waxler and Assink (2019). Figure 7 presents the wind in direction of propagation, the absolute temperature, and the effective sound speed at location and time of the Moerdijk explosions according to ECMWF (9-hr forecast), HIRLAM (3-hr forecast), and HARMONIE (nowcast). They vary in time and space resolution ( Table 2). Additionally, independent wind speed and temperature measurements during the explosion time, located at CIA, are plotted.

Result of Propagation Modeling
Differences in effective sound speed profiles in the ABL are observed (Figure 7c). Regional and global models have difficulties resolving the stable boundary layer, in particular when the wind speed values are low (Holtslag et al., 2013). Other differences can be explained by the different model parameters (e.g., dynamical core) and resolutions as mentioned in Table 2.
The Moerdijk explosions occurred 1 hr after sunset. A temperature inversion is observed for all weather models, with its maximum around 60-m altitude. Compared to the in situ measurements, the weather prediction models show lower maximum temperatures and temperature gradients (Figure 7b). The positive temperature gradient contributes to downward refraction. A negative wind speed gradient in the HIRLAM and HARMONIE profiles in direction of propagation is shown near the surface (Figure 7a). The in situ wind speed measurements show low wind speed values near the surface (<200 m) with a positive gradient. Of all considered models, ECMWF is in closest agreement with the in situ data.
All weather models predict the existence of an acoustic waveguide within the ABL. However, the characteristics of the waveguide vary among the models (Figure 8). The waveguide is formed between the ground and 1,500-, 700-, and 600-m altitude for HARMONIE, ECMWF, and HIRLAM, respectively. This difference has implications for the arrival times and shape of the propagating wavefronts in the ABL (Figure 9). The strength of the waveguide can be attributed to the difference between the effective sound speed at the surface and at the turning altitude. For ECMWF this difference is 1.9 m/s, for HIRLAM 1.4 m/s, and for HARMONIE 1.8 m/s. The ECMWF waveguide is therefore strongest. Furthermore, ECMWF is the only model that predicts downward refraction of the raypaths in the first couple of 100 m. For HIRLAM and HARMONIE first upward refraction and subsequent downward refraction are modeled. Multiple ground bounces on the Earth's surface produce upgoing and downgoing wavefronts arriving at the tower. A positive vertical slowness (blue colored) indicates an upgoing wavefront and a negative vertical slowness (red colored) a downgoing wavefront. The modeled wavefronts using HARMONIE and HIRLAM do not explain the tower observations. Instead of the observed decrease in delay time between the wavefronts with increasing altitude (∧-shape), an increase in delay time (∨-shape) is modeled. This is because the modeled wavefront reflects at the surface when arriving at the CIA tower. The modeled travel times are ∼118.0, ∼119.3, and ∼118.9 s using ECMWF, HIRLAM, and HARMONIE, respectively. The modeled travel time using ECMWF is closest to the observed travel time of ∼118.0 s.
In conclusion, none of the selected models describes the atmospheric state in such a way that ray theory can simulate the observed wavefront structure to explain the tower observations. The atmospheric weather models are an average representation of the atmosphere at the event time and location. However, the ABL is dynamically changing in smaller temporal and spatial resolution and models are not suited to resolve these small-scale variations. The results show the importance of various descriptions of the atmosphere and its influence on the wavefront arrivals in the ABL. The tower observations are an added value to help understand the ABL influence on ray propagation.

Discussion
In this study, recordings of a three-dimensional microbarometer array are used. The differences between two-and three-dimensional observations are analyzed using the infrasonic signals from the accidental Shell Moerdijk explosion on 3 June 2014. The signals are interpreted to be from the Moerdijk explosions, since the ground array beamforming showed that the signals originate from the Moerdijk direction with infrasonic frequencies and velocities.
The Moerdijk explosions occurred after sunset (local origin time is 22:50:24), with a positive temperature gradient below 200 m. This means that the ABL was in transition to become stable. Figure S2 shows the transition of the ABL on 3 June, derived from the meteorological sensors at CESAR. Figure S3 shows the vertical profiles around the time of the explosion. Here, downward refraction of infrasonic waves is observed due to the combined temperature inversion and low wind speed values. Downward refracted low-frequency audible sound in a stable ABL has been studied previously by Chunchuzov et al. (2017) and Talmadge et al. (2008).
It appears difficult to simulate these observed wavefront arrivals using global and regional numerical weather prediction models. Weather prediction models are an average representation of the reality and getting the small scale right is a field of active research (Heus et al., 2010;Siebesma et al., 2003). Differences are observed between the models near the surface. ECMWF vertical wind speed gradients are most similar to the in situ wind speed measurements. The ECMWF profile explains best the travel times. Both ∧and ∨shape wavefront arrivals are modeled; hence, ECMWF explains part of the wavefront structure (∧-shape) that is observed in the data. HARMONIE and HIRLAM solely model ∨-shape arrivals. Therefore, none of the models are suited to explain the observed upgoing and downgoing wavefront structures.
Ray theory is exact when atmospheric scales are large compared to infrasonic wave lengths ( ≪ L; Brekhovskikh & Godin, 1999). This means that for the signals with the lowest frequencies, generated by the Moerdijk explosions, it is disputable if the raypaths as shown here are affected by the fine-scale temporal and spatial atmospheric structure. For future studies, frequency-dependent full wave modeling of the wavefronts can be used (Assink et al., 2017;Waxler & Assink, 2019). Full wave modeling would be useful to understand the dispersion of infrasound in the ABL to higher order. The full wave modeling will still depend on the atmospheric weather prediction models that lack these small-scale resolution needed. Instead of weather prediction models that parameterize the ABL, a model can be used that solves the ABL (e.g., large eddy simulation).
It is shown that the tower waveforms can be represented as two sequential ground waveforms with decreasing delay time with increasing altitude (∧-shape). This means that approaching the surface, the delay time between the wavefronts increases. However, observations at ground level do not show continuation of this trend. Only one wavefront arrival is detected.
Besides the direct source-receiver path, refracted raypaths arrive at the receiver. Due to the time difference between the different paths, a single explosion can result in multiple arrivals. One hypothesis to explain the Moerdijk observations is that an upward propagating direct wave is recorded at the ground and in the tower. The downward refracted wave in the ABL is subsequently recorded, solely at the tower microbarometers. Since it can be trapped in an elevated waveguide with the lower boundary between 0-and 60-m altitude, the downward propagating wave is not reaching the surface microbarometers. Similar observations are presented in Lees (2015, 2017). A high-altitude balloon was used, where elevated acoustic waveguides are observed.
This hypothesis explains why there are two ground wavefronts detected in the tower arriving near horizontally (∼9 • ). Furthermore, it explains why there is an upgoing wavefront followed by a downgoing wavefront, which is not recorded at the surface microbarometers. To verify this hypothesis at tropospheric levels, more research is needed. Figure 10 shows results from three-dimensional ray tracing using the range-dependent HARMONIE model and the GeoAc ray tracer from Blom (2014). In contrast to the atmospheres used in Figure 8, here it is assumed that the atmosphere is not horizontally homogeneous. While the figure suggests that individual waveguides can be distinguished, it should be noted that the validity of ray theory for such small waveguides Figure 11. Band-pass-filtered tower observations of (a) the first Moerdijk explosion and (b) a sonic boom event over the North Sea on 24 January 2017. The red lines show the waveforms produced by addition of two ground wavefronts with increasing time between them with increasing altitude. The ∧-shape (a) and ∨-shape (b) are visible in the tower observations. Note the difference in y axis scale (pressure). is limited. However, this result illustrates the complexity of the acoustic waveguides that may be present in the ABL at any given moment in time. Nevertheless, multiple waveguides present in the ABL explain part of the tower and ground observations, which could explain why the refracted waveforms do not fully reach the surface.
Next to ∧-shape wavefronts (Figure 11a), an increase in delay time with increasing altitude has been observed as well, for different infrasonic events. This occurs when the wavefront reflects at the surface (∨-shape).
An example is a sonic boom from the North Sea (approximated distance of 240 km from CIA) on 24 January 2017. Here porous hoses were installed. The tower recordings are presented in Figure 11b. The same hypothesis as the Moerdijk explosions is tested for the sonic boom event: a model of two superposed ground waves is fitted to the real data. The best fit between the real data and simulated waveforms is by adding two infrasonic ground recordings with delay times of 0.17, 0.26, 0.35, 0.45, and 0.53 s for stations CIA11-CIA15, respectively ( Figure 11b). First, the highest station (CIA15) is reached and lastly the ground stations, from where the wavefront reflects upward. The tower observations provide this additional information that would not have been observed using ground observations.

Concluding Remarks and Future Studies
In this study, infrasound observations of a three-dimensional microbarometer array are shown in nocturnal boundary layer conditions. Here, the focus is on the influence of evening ABL conditions with low wind speed and shallow waveguides on local infrasound propagation. Clear infrasound arrivals of the Moerdijk explosions are observed in the ground recordings and two sequential arrivals in the tower recordings. The two sequential arrivals are interpreted to be interfering upgoing and downgoing wavefronts. The numerical weather prediction models are not suited to solve the small temporal and spatial ABL scales used here. Differences appear between the regional models, global models, and in situ wind and temperature measurements near the surface.
Future studies are needed on the effect of the ABL processes (e.g., turbulence) on infrasound waves and to understand the effect of the tower structure on the microbarometer observations. The focus will be on detections during the day, when convection and turbulence processes contribute significantly to the wind noise at infrasonic frequencies. The detectability decreases during high turbulent periods. With more information available on the influence of the ABL, future microbarometer array locations and configurations can be improved and hence the infrasound detections. Such knowledge is of importance for a successful verification of the CTBT, for which infrasound is one of the verification techniques. Moreover, future studies will focus on the use of infrasound as a remote sensing technique of the ABL (e.g., Johnson et al., 2012).
In summary, this work indicates that the added value of the third dimension is the sensitivity to the atmosphere. The state of the ABL, around the explosion time (20:50 UTC), induced relatively low noise levels and downward propagation due to the small wind speed and positive temperature gradient. The ray propagation modeling results, using all three weather prediction models, show downward refracted waves trapped in a near-surface (<1.5 km) waveguide. The simulated travel times correspond to the observed travel times. However, the upgoing and downgoing wavefronts were not explained by the modeled wavefront arrivals.
Concluding from the tower observations, a hypothesis is found for the Moerdijk case implying that there is a possible second waveguide at higher altitudes in which the downward propagating wave is trapped and not able to reach the surface. This would not have been observed from solely the ground observations and one-dimensional ray modeling using ray theory.
This study shows the complexity of modeling local infrasound raypaths due to smaller unresolved ABL scales, which was used for explaining the first vertical infrasound observations. contributions are funded through a VIDI project from the Netherlands Organisation for Scientific Research (NWO), Project 864.14.005. The data will be available in the KNMI database; see RDSA.knmi.nl/network/ NL for details. Data are available thanks to Gert-Jan van den Hazel and Jordi Domingo-Ballesta. FDSN open data webservices can be used to obtain data. Special thanks to Henk Klein Baltink (KNMI), Jordi Vila (Wageningen University), Sander Tijm (KNMI), Thaisa van der Woude (University of Technology Delft), and Mathijs Koymans (KNMI). The figures are made with Matplotlib v.2.1.0 (Hunter, 2007) and Generic Mapping Tools (Wessel et al., 2013). The GeoAc range-dependent ray tracer is used from Blom (2014). Independent measurements (Wind profiler measurements, tower data, and MWR data) are from the Ruisdael observatory.