A Comprehensive Study of Infrasound Signals Detected from the Ingolstadt, Germany, Explosion of 1 September 2018

The explosion at the Ingolstadt oil refinery was widely recorded at seismic and infrasound stations deployed throughout Central Europe, to distances of several hundred to a thousand kilometres. This study focuses on the wealth of data recorded at infrasound stations in Central and Eastern Europe, while from the many detecting seismic stations within 400 km range, only seismic and seismo-acoustic arrivals at the close-in Gräfenberg array are considered here. Most of the infrasound stations are acoustic arrays enabling us to apply array processing techniques to determine relevant wave field parameters, such as backazimuth and slowness (resp. trace velocity). These parameters not only confirm the source direction, but also put constraints on the observed arrivals’ propagation modes. Wave field parameters suggest that we observe tropospheric arrivals to about 150 km and stratospheric and/or thermospheric returns for longer distances. 1D, 2D and 3D ray tracing predict tropospheric arrivals to westerly directions up to distances of 100 km, beyond which only thermospheric returns are obtained azimuth-independent beyond 250–300 km. Stratospheric returns do not follow from any of the increasingly complex ray tracing models. Parabolic equation propagation modeling however suggests that in a number of cases stratospheric ducting may be possible. However, neither the tropospheric seismo-acoustic arrivals at the Gräfenberg array nor the various arrivals at IMS station IS26 could be modeled. Therefore, the Ingolstadt explosion along with the observed infrasonic phases provide an excellent test bed to investigate our ability in realistically forecasting atmospheric wave propagation with existing algorithms and available atmospheric models.


Introduction
Probing the atmosphere by the propagation of infrasound signals has now become an important scientific field with the advent of the Comprehensive Nuclear-Test-Ban Treaty. There, a network of infrasound sensors within the International Monitoring System (IMS) is surveilling the globe for atmospheric pressure disturbances that could originate from nuclear testing, especially in the atmosphere, and thus would constitute a violation of this treaty. As this infrasound sensor network is sparse and one primary objective within the verification tasks for the Treaty is the accurate location of such a source, an accurate model of atmospheric propagation conditions is required. Only this can ensure the precise interpretation of the observed signals from the associated atmospheric waves, which can travel closely along the Earth's surface within the tropospheric duct or at greater altitudes in the stratospheric or thermospheric ducts (Drob et al. 2003). If the stratospheric duct, at altitudes of 40-60 km, is well developed depending strongly on the wind and temperature field, it is usually very efficient in the transport of acoustic energy to greater distances, as the attenuation within this duct is not nearly as strong as at higher altitudes. There, at altitudes of more than 100-120 km, the thermosphere allows for infrasound returns due to the large effective sound speeds based on the prevailing high temperatures. However, such thermospheric returns are not as widely observed as stratospheric returns, and usually only for rather low signal frequencies, since the attenuation in the thermosphere is much higher than in the stratosphere (Sutherland and Bass 2004).
In this study we investigate the infrasound signals observed from the 1 September 2018 explosion on the outskirts of the city of Ingolstadt, when at the oil refinery in Vohburg a leakage in a tank generated a large cloud which ignited at 3:11:45 UTC and generated seismic and acoustic waves that were recorded at local and regional seismic sensor networks within 400 km distance of the source (Fuchs et al. 2019). This ground truth information is available from the media (e.g. Donaukurier 2018), where also a later second explosion is mentioned. The closest sensors at distances of some 15-20 km from the explosion site belong to the Gräfenberg seismic array (GRF), where not only strong seismic signals were observed but also later seismo-acoustic phases from the atmospheric waves. Fuchs et al. (2019) looked at the seismic data from the temporary AlpArray network (Hetenyi et al. 2018) and found, besides the corresponding seismic phases, different seismo-acoustic signal groups. They attributed these to signals that propagated in the tropospheric duct with celerities of around 332 m/s, on average, to distances of up to 200 km, in the stratospheric duct with a mean celerity of 292 m/s and observed between 150 and 325 km distance, or most likely in the thermosphere with a mean celerity of about 250 m/s identifiable at distances of 250-350 km from the source. Except for the strong signals in the GRF data, we focus here exclusively on observations from infrasound sensor stations throughout Central Europe to distances of 1000 km and more, by incorporating additional information from a number of infrasound stations in Eastern Europe that were presented recently (Czanik et al. 2018).

Infrasound Observations at Regional Infrasound Sensors
The oil refinery at Vohburg near Ingolstadt, Bavaria, is located at the geographical coordinates of latitude 48.766°N and longitude 11.6°E and depicted in Fig. 1. North of the explosion site the 13 stations of the GRF seismic array (Buttkus 1986) are situated at distances from 15 to 110 km, showing rather strong seismo-acoustic signals at short distances, but they can be tracked to distances of 100 km. Only the farthest stations (GRA1, GRA3) lack any noticeable acoustic arrival. We only include these data in this study, firstly, for completeness, since we used the seismic phases at GRF for event location, especially to confirm and to specify an origin time of the event, and secondly, to infer the pattern for tropospheric infrasound propagation towards the North. For seismo-acoustic arrivals at other seismic stations, namely the AlpArray network with some 400 stations within 400 km of the explosion site, we refer to Fuchs et al. (2019) who investigated this data in greater detail.
The data that were analysed here are from acoustic sensors at different locations within about 600-700 km from Ingolstadt. Most of them consist of sensor arrays enabling us to apply array analysis techniques (Cansi 1995;Stammler 1993) to extract waveform parameters (as summarized in Table 1), except for single-sensor stations BFO and HLG/ HLGM1 and for infrasound array DIA (Fig. 1). For the latter we were only able to obtain data from 3 elements, with one of these elements providing data not suitable for analysis. Therefore, these stations yielded only celerity estimates from arrival time readings, while no backazimuth or slowness estimates could be made. The waveforms for all these stations are shown in Fig. 2, also including the GRF array data.
The closest infrasound array IS26 in the Bavarian forest is 155 km to the east from the explosion site. The waveforms show three signal groups, the first arriving with a celerity of 327 m/s, merging directly into the second group consisting seemingly of two distinct arrivals with celerity of 275-285 m/s. For these two signal groups we used an analysis window of about 1 min. A third group follows significantly later, i.e. more than 5 min and with considerably lower amplitude, but with a similar frequency content, with a celerity of a little less than 200 m/s within a window of 2.5 min. The frequency-wavenumber (F-K) analyses (Stammler 1993) of these three waveform groups are presented in Fig. 3a-c. The waveform parameters for these signal groups from F-K analysis in terms of backazimuth and slowness, as well as translated trace velocity, are given in Table 1. While the F-K analyses in Fig. 3a, b display isolated peaks in the wave-number plane, this does not apply for the third waveform group, where the peak seems blurred (Fig. 3c). Analyzing this longer waveform group in successive time windows of shorter length shows a steady increase of the slowness from 303 to 326 s/°(or decrease of the trace velocity from 367 to 341 m/s).
However, a note is in order with respect to this last signal. According to information from the media (Donaukurier 2018) and confirmed by officials of the oil processing facility (M. Raue, personal communication, 2019), there has been a secondary event that took place about 5-6 min after the initial explosion. It consisted of the deflagration of a gas cloud, which might not have produced significant seismic energy to be seen at the adjacent GRF stations, but enough acoustic energy to be detected at IS26. From the timing of this secondary source the last signal group at IS26 could correspond to the equivalent acoustic Map of Central Europe centered around the Ingolstadt explosion site (red star) displaying the locations of the infrasound stations (larger cyan triangles) within 700 km distance and the Gräfenberg seismic array (smaller black triangles) used in this study. Concentric circles depict the distance ranges to 700 km in intervals of 100 km Table 1 Infrasound array stations and the associated signal parameters of the signal groups analyzed for the Ingolstadt explosion (see F-K-analyses in Figs waves as observed from the initial explosion, by the earlier two signal groups. An additional F-K analysis was carried out on the seismic data of the co-located GERES array, with 10 elements operational on 1 September 2018, for the second waveform group, the largest in amplitude at the IS26 station, but still not visible in the seismic waveform data. The results of this analysis are shown in Fig. 3d, with estimated wave field parameters given in Table 1 and matching the results obtained for IS26 well. The consistency of our results strongly argues that seismic arrays may well be used in an attempt to identify acoustic waves even when the signal-to-noise ratio is definitively less than 1. As Fuchs et al. (2019) point out seismic stations (in their case the AlpArray network) are indeed useful to supplement infrasound stations in infrasound studies by using seismo-acoustic arrivals, as has been done e.g. by Hedlin et al. (2010), which studied acoustic propagation from a large bolide across the dense US Array seismic network.
The IS26 data were further processed with the PMCC method (Cansi 1995) to confirm the results obtained by interactive F-K analysis, with the outcome of this processing shown in Fig. 4. We obtain with PMCC again two signal groups over a time span of some 150 s corresponding to those shown in Fig. 3a, b with consistent parameters for back azimuth and trace velocity, with the latter parameter showing a steady increase from about 330-360 m/s. With PMCC we were, however, not able to retrieve any detection for the third signal group (see Fig. 3c) that follows with some 5 min delay. This fact is somewhat surprising since the signal has a signal-tonoise ratio of about 2-3 at frequencies of about 1-2 Hz.
Five other infrasound arrays at distances between 500 and more than 600 km are located to the northwest (IGADE, CIA, DBNI, EXLI) or to the east of Ingolstadt (PSZI) for which we were able to obtain waveform data. The F-K analysis of IGADE data for a 3 min window (190 s) is given in Fig. 5a showing consistency between theoretical and estimated backazimuth (see Table 1). The slowness (trace velocity) is found as 326 s/°(341 m/s). The Dutch stations CIA, DBNI and EXLI at a similar distance range but about 15°-30°in backazimuth further west from the source provide consistent azimuths for slightly shorter analysis windows (* 150-160 s), but the slowness values differ quite strongly ( Fig. 5b-d, Table 1). For the larger aperture CIA a slightly higher slowness value was found than for IGADE, while the small aperture arrays DBNI and EXLI provided low slowness estimates. Station PSZI in Hungary exhibits a fairly impulsive arrival and therefore F-K analysis was carried out for only a 80 s time window. It yields a consistent backazimuth, and typical slowness and trace velocity values (Fig. 5e, Table 1). For this station at 620 km range, infrasound station PCVI in the Czech Republic at shorter distance (* 300 km) as well as three infrasound stations in Romania (I67RO, BURARI, IPLOR) at longer distance (900-1200 km), Czanik et al. (2018) reported infrasound detections that they associated to the Ingolstadt explosion. Some of their estimated parameters will be used below for comparison with the results obtained in this work.
PMCC analysis was carried out as well for the IGADE data highlighting a rather long signal (* 300 s) with very stable estimates for arrival direction and trace velocity (Fig. 6). The signal identified by PMCC starts even earlier than the large waveform signal used in the analysis shown in Fig. 5a. Azimuth and trace velocity estimates from both methods are again in good agreement indicative for the Ingolstadt source.

Infrasound Propagation Modeling
Propagation modeling is carried out with a variety of approaches for the interpretation of the signals observed at the network of European infrasound stations. Most of these stations show very consistent b Figure 2 Record section of the available waveform data with the (seismo-)acoustic signals observed at the GRF as well as the infrasound arrays. The data are bandpass-filtered between 1 and 5 Hz. Traces are ordered according to increasing distance from bottom to top. The vertical line marks represent the arrivals times for phases with celerity of 330 m/s (GRF, GEC2, IS26), and 300 m/s otherwise. For IS26 the data windows for the three different analyzed phases (I1, I2, I3) are sketched and labeled  Figure 3 F-K-analyses for different signals groups at IS26 for a the initial emergent onset; b the second high amplitude onset; and c the late emergent arrival trailing the second arrival by about 5 min (see Fig. 2). For the time window of the second arrival at IS26 we analyzed the seismic data recorded at GERES (d) showing consistent results with respect to the infrasound data. Lines show backazimuth segments of 30°clockwise from north (top), circles quantify slowness in 100 s/°steps. The signal parameters extracted for the signal peak (yellow spots) are summarized in Table 1 4234 K. Koch and C. Pilger Pure Appl. Geophys.
pattern. These latter codes may provide additional insight when the methods for simpler modeling approaches may not have succeeded in predicting some or any acoustic arrivals at a particular station. Only a subset of results from the applied modeling is presented within this study, highlighting the ducting conditions into eastern and northern directions (using 2D ray tracing and 1D PE for the case of PSZI and IGADE, respectively). In the exemplary presentations for different stations we aimed at emphasizing basic wave propagation features but also the potential for particular stratospheric wave guides (using 2D PE towards IGADE and DBNI). And for the PE calculations we selected a frequency of 1 Hz which is relevant for the dominant frequencies of 0.5 Hz (IGADE) to 2 Hz (GRF) found for the data shown in Fig. 2, depending on the distance from the source.
Propagation modeling was carried out on the basis of the European Centre for Medium Range Weather Forecasts (ECMWF) atmospheric background model of the 6-h analysis data. It provides physical parameters (pressure, temperature, wind speeds) for altitudes up to 60-70 km, beyond which the atmospheric model is extended by smoothly merging it with the HWM07 and MSISE00 climatologies (Drob et al. 2008;Picone et al. 2002). Within the 2D ray tracing modeling a range-dependent sound speed profile is established from interpolating the atmospheric specification with a resolution of 200 m both in range and altitude. For the 3D modeling with GeoAc we defined a spatial grid of profiles with 1°9 2°spacing in latitude and longitude from 44°to 56°N and 2°to 30°E, respectively. The specific ECMWF model taken was for 06:00 UTC on 1 September 2018, where for the 1D PE modeling we selected the mid-range profile between source and receiver. In the case of 2D PE calculations the grid points were chosen with a resolution on the order of tens of meters horizontally (35 m) and vertically (15 m).
As an example for the atmospheric models used within the propagation modeling, the effective sound speed profiles from Ingolstadt towards different directions on 1 September 2018 are discussed (Fig. 7): to the north (meridional direction) towards the infrasound station HLG (and IGADE as well), to the east (zonal direction) towards the Hungarian infrasound station PSZI, and towards the opposite direction (west) to station BFO. These profiles, given at various ranges along the path, show sound speeds in the lower troposphere of about 325-335 m/s which is higher than the maximum speed in the stratosphere reaching up to 320 m/s at altitudes around 40 km or in the mesosphere/lower thermosphere between 90 and 105 km altitude. This renders effective sound speed ratios (c_eff_ratio), i.e. the ratio of the sound speed at altitude to the sound speed near the ground, of less than 1. Thus, no acoustic wave returns are favoured in any direction from these height regions. The ratio reaches the highest values, close to 1, to the west (BFO), intermediate values of 0.95 and more towards the north (HLG), and the lowest values (0.95 or less) to the east (PSZI). Only at altitudes larger than 120 km, the effective sound speed exceeds the surface sound speed, yielding c_eff_ratios well above 1. Hence, at larger distances expected infrasound signals should only result from thermospheric ducting.
The results of the propagation modeling to station IGADE (at a range of about 550 km), towards a northerly directions (i.e. azimuth 338°) from Ingolstadt is displayed in Fig. 8 using the 2D ray tracing and the 1D PE methods. The model exhibits larger sound speeds near the ground (\ 10 km altitude), a rather thick layer of increased velocity between 40 and 80 km altitude and the highest velocities in the thermosphere above heights of 110 km. From this sound speed layering a quite simple propagation pattern emerges in which only rays corresponding to thermospheric arrivals are possible. Interesting in this case is also the fact that two arrivals can be expected,  Fig. 3, but for the infrasound signals identified at a IGADE; b CIA; c DBNI; d EXLI; e PSZI. For the extracted signal peak parameters representative for the center of the yellow spot see Table 1 Vol. 177, (2020) A Comprehensive Study of Infrasound Signals Detected from the Ingolstadt, Germany 4237 one from a single bounce ray turning at an altitude of about 115 km and a second, steeper double bounced ray with a turning altitude of about 125 km. In the 1D PE code attenuation is taken into account (Sutherland and Bass 2004)-and from this modeling we see a clear reduction of the infrasound energy with distance and altitude. While the energy for thermospheric propagating has high amplitudes out to nearly 300 km, it decreases rapidly with distance, so that at distances more than 400 km, we obtain only amplitudes that have attenuated by about 100 dB. On the other hand, we obtain amplitudes beyond these distances having attenuated less, only to about 80 dB (i.e. 20 dB less), that extend only to stratospheric heights. Thus PE modeling suggests that some of the acoustic energy is indeed refracted at 40-50 km altitude, where an increase of the effective sound speed is observed. So in contrast to the ray tracing suggesting only thermospheric ducting, PE modeling argues for some odds of observing stratospheric rather than thermospheric arrivals.
For the travel path from Ingolstadt towards the east, i.e. station PSZI in Hungary, at a distance of some 600 km, the atmospheric model (see also Fig. 7) shows again higher velocity in the tropospheric layer, a rather thin band for the stratospheric layer between 40 and 50 km with moderate velocity and the thermospheric duct above 125 km altitude with the highest effective sound speeds. Therefore the situation is not significantly different than for the case Figure 7 Along-path effective sound speed profiles on 1 September 2018 06:00 UTC for different propagation directions from Ingolstadt (first profile) to the respective stations (last profile) used in the 2D modeling: a towards PSZI in an easterly direction; b towards HLG (and IGADE) in the northerly direction; c towards BFO in a westerly direction. Gray vertical lines indicate effective sound speed ratios of 0.9, 0.95 and 1.0 (with respect to the effective sound speed at the ground level). Only towards the west a ratio of 1 is nearly reached for stratospheric levels at height around 50 km, for which stratospheric ducting is favored most 4238 K. Koch and C. Pilger Pure Appl. Geophys. of IGADE (Fig. 8) in that ray tracing (Fig. 9) yields a thermospheric arrival, but no stratospheric returns. When modeling infrasound propagation with the PE method a similar pattern to IGADE is seen here as well, i.e. the wave field energy indicating thermospheric propagation attenuates even more rapidly, so that beyond about 200-250 km the major remaining energy is carried through a stratospheric duct.
Since we performed ray tracing in the 2D model we also applied a 2D PE algorithm for the case of IGADE and the Dutch station DBNI. The corresponding results are displayed in Fig. 10 showing again that the acoustic energy from the source is efficiently propagated between the ground and the stratospheric height of about 50-60 km. The energy reaching the thermosphere beyond 120 km is in contrast heavily attenuated as has been observed for the 1D case before. These additional modeling exercises thus suggest that the prediction of thermospheric arrivals at intermediate distances may need to be interpreted with due diligence and that in these cases stratospheric returns may be observed in reality. This is in agreement with findings by Le Pichon et al. (2012) and Landes et al. (2014), as will be discussed in the next section.
Since the modeling with 1D and 2D methods resulted in some tropospheric ducting at shorter distances of up to about 100 km distance and otherwise in thermospheric ducting where the rays reach turnings height of more than 110-120 km, we also applied the 3D ray tracing package GeoAc (Blom and Waxler 2012) to investigate whether this more sophisticated approach would allow additional rays within a stratospheric duct, as tentatively suggested by the PE calculations. The results from this modeling (Fig. 11) does not produce such additional propagation paths, but merely confirms that wave propagation though the ECMWF model for 1 September 2018 and Central Europe with ray tracing methods does not allow other than tropospheric and Propagation modeling for a northerly path towards station IGADE using (a top) ray tracing (no attenuation model included) showing two thermospheric eigenrays at or near the station and (b bottom) 1D PE modeling with considerable energy being leaked into a duct at stratospheric heights; thermospheric ducting is not apparent due to the attenuation model. The station location is marked by a triangle, the effective sound speed in (a top) is color-coded in the background Vol. 177, (2020) A Comprehensive Study of Infrasound Signals Detected from the Ingolstadt, Germany 4239 thermospheric paths, as indicated by the pattern of bounce points for the rays from Ingolstadt to any place in Central Europe. In particular this methods emphasizes that at some places about 500-600 km to the West and North multipathing occurs with one ray turning below 120 km and the other one turning above 120 km altitude.

Discussion
In this study we have analyzed the available infrasound signals from a heavy explosion at an oil processing facility near Ingolstadt on 1 September 2018 recorded at a considerable number of infrasound sensors that are deployed throughout Central Europe, namely Germany, the Netherlands and Hungary. These results will be augmented by the data that were published in work by Czanik et al. (2018) for additional infrasound stations in the Czech Republic and in Romania. From array processing methods of the signals recorded at the infrasound array stations (see Table 1), we find that they arrive from the direction of the explosion site, with the estimated backazimuth agreeing within 5.5°of the theoretical azimuth. For the estimated values for slowness or trace velocity, respectively, the signals show a larger range covering values of 339 s/°to 240 s/°in slowness translating into trace velocities between 327 and 463 m/s. These trace velocity values normally represent the sound speeds of the atmospheric waves at the altitudes where they turn down again. While most values are typical for sound speeds at stratospheric heights (Table 1), where often the values are alike 330-340 m/s, the higher values observed for DBNI and EXLI and the late arrival at IS26 could indicate thermospheric arrivals. While for IS26 it may be an option, it does not apply for DBNI and EXLI being at a similar range (and azimuth) as other stations like CIA and IGADE. We note here the very small Propagation modeling for an easterly direction towards station PSZI for a ray tracing and b 1D PE modeling. As in Fig. 8, with ray tracing supporting the interpretation of the observed infrasound signal as a thermospheric arrival, while PE modeling suggesting stratospheric ducting 4240 K. Koch and C. Pilger Pure Appl. Geophys. aperture of these arrays on the order of tens of meters, at least as the elements are concerned which were operational and for which we could retrieve waveform data. For these cases local wind anomalies could seriously distort the arriving wavefront and therefore impact on the estimated parameters, as was argued by Negraru et al. (2010) for high phase velocity measurement obtained in their study. A good explanation for the very high trace velocities cannot be given, especially, as the backazimuth estimates are in general rather close to the expected values. The celerity, i.e. the mean velocity referred to a path along the ground surface, calculated from the travel time divided by the station distance, has been used to characterize the potential propagation mode of acoustic waves. As Negraru et al. (2010) pointed out propagation within different atmospheric ducts can be associated with specific celerity ranges (see also Kulichkov 2000): (seismo-) acoustic arrivals with celerity values between 310 and 330 m/s are most likely related to the tropospheric duct, those with celerities between 280 and 320 m/s are indicative for stratospheric returns, while thermospheric returns show celerities between 180 and 300 m/s, depending on the range from the source. For all stations, listed in Table 1 or considered by Czanik et al. (2018), the estimated celerity values are shown in Fig. 12, where for most of them several data points are given, when the observations exhibited distinct pulses. For example, as shown in Fig. 2, the IMS infrasound station IS26 showed three such onsets, an emergent initial signal, a later arrival with two high amplitude onsets and a late arrival with low signal-tonoise ratio. In this case the corresponding celerity values are 330 m/s, 285-288 m/s and 270 m/s, resp., as displayed in Fig. 12.

Figure 10
Results of 2D PE calculations for two stations towards the north with azimuth difference of about 25°(IGADE vs. DBNI). Here too an atmospheric attenuation model is included by which thermospheric energy is highly reduced, while tropospheric and stratospheric energy is distinctly enhanced. The respective station location is marked by a triangle above the abscissa Vol. 177, (2020) A Comprehensive Study of Infrasound Signals Detected from the Ingolstadt, Germany 4241 The closest stations to the explosion site are eleven elements of the Gräfenberg seismic array (GRF) showing clear seismo-acoustic arrivals, and IMS infrasound station IS26 with acoustic signals, where the earliest onsets, and hence those with the highest celerities above 320 m/s as shown in Fig. 12, must then be associated with tropospheric ducting. For more distant stations, up to 600-700 km, the earliest arrivals are associated with celerities between 300 and 320 m/s. This fact argues strongly for stratospheric ducting, while the celerities associated with thermospheric arrivals obtained from 3D ray tracing suggest much lower values. In these cases, mostly related to the infrasound stations in Northern Germany and the Netherlands, the earliest arrival seems therefore to be of stratospheric origin, supporting the findings obtained from PE modeling. Only for the infrasound stations beyond 900 km, located to the east in Romania, observed celerities compare favorably with those obtained from ray tracing. The associated arrivals are therefore quite likely associated with thermospheric ducting, but multiply bounced stratospheric arrivals, as is also suggested for PSZI by Fig. 9, cannot be ruled out. As the ray tracing allows only for thermospherically refracted rays and is close to or within a shadow zone between the first and second thermospheric bounce (Fig. 11) this direction to the east could allow for favorable propagation conditions for stratospheric waves.
In general, ducting in the stratosphere (or thermosphere) is possible, if the effective sound speed ratio is larger or equal than 1, or is impossible otherwise. In terms of numerical modeling this is strictly the case for the ray tracing method. For wave theory methods, such as PE modeling, Le Pichon et al. (2012) and Landes et al. (2014) however studied various ranges of this ratio and found that stratospheric ducting may already occur at a ratio [ 0.9, depending on frequency, hence the stratospheric ducting mechanisms appearing in Figs. 8, 9 and 10 may indeed be realistic. Bounce point pattern from 3D ray tracing through a grid of ECMWF models covering Central and parts of Eastern Europe-green symbols indicate bounce points from tropospheric propagation; red and blue symbols are for thermospheric propagation with ray turning heights between 110 and 120 km or above 120 km, respectively. The locations of the source and the infrasound stations as given in Fig. 1  For some stations, in particular IS26, the ray tracing did not produce results indicating any potential propagation path, and this applies as well to PE modeling. The ECMWF atmospheric model for 1 September 2018 seems not to describe the atmospheric conditions accurately enough in order to produce relevant arrivals in the modeling. Of course this could be related to the rather coarse 6 h sampling of the weather prediction model of ECMWF used here. However, similar findings have been obtained by Fuchs et al. (2019), who were also unable to model any stratospheric arrivals even though they found strong evidence of a stratospheric propagation branch for distances between 150 and 325 km from the Ingolstadt source. In their case they applied the same 3D ray tracing scheme as used here and an even higher-resolution ECMWF model.
We argue here as well that the Ingolstadt explosion has occurred at a time near the autumnal equinox, where the winds in the atmosphere are variable and change rapidly from one dominant period to the other and which may not be captured adequately be numerical weather prediction models (Le Pichon et al. 2015). A similar example of potential stratospheric returns being difficult to model, if at all, from an event at equinox times, is available with the North Korean nuclear explosion of 3 September 2017, as reported by Assink et al. (2018) and Gaebler et al. (2019).
As laid out above, a secondary deflagration was reported to have occurred some 5-6 min after the explosion. From the arrival time of the last (third) signal group at IS26 and assuming it would not be a thermospheric return (observed celerity 194 m/s) from the initial explosion, but rather match the same propagation mode as the earlier arrivals, with a celerity of 325-330 m/s for this phase it can be estimated that the secondary source followed the explosion by approximately 5:30 min.

Concluding Remarks
The explosion on 1 September 2018 near Ingolstadt has been studied here and also by Fuchs et al.

Figure 12
Observed versus modeled [from 3D ray tracing (RT)] celerity values for all stations analyzed, complemented with results obtained by Czanik et al. (2018), plotted against distance and annotated with the station code. The parameters from modeling are shown as orange squares (two symbols per station represent multipathing (see also Fig. 11 and text). Multiple symbols for the observed values (blue diamonds) represent different onset times within an infrasound signal (e.g. IS26-see text and discussion on Fig. 3) Vol. 177, (2020) A Comprehensive Study of Infrasound Signals Detected from the Ingolstadt, Germany (2019). The data analyses in both studies have provided a clear set of acoustic arrivals at seismic and infrasound stations throughout Central Europe. While the seismic data of the AlpArray network analysed by Fuchs et al. (2019) could only be interpreted in terms of celerity based on apparent arrival time of the seismoacoustic phase, we have focused here on the infrasound array stations for which we were able to apply array processing techniques in order to obtain additional waveform parameters such as backazimuth and apparent (trace) velocity (or equivalently slowness) of the waveform signal. For the backazimuths we see a strong agreement with the theoretical values, with deviations of only a few degrees. For the trace velocity we find typical values representative for tropospheric or stratospheric ducting, i.e. values between 325 and 345 m/s, with the exception of the late arrival at IS26 and at the small aperture arrays DBNI and EXLI in the Netherlands. The deviation at IS26 seems to be related to a potential overlapping of wave groups, while the large differences for the two other stations remain unknown. Any array configuration errors seem not likely, as we obtain reasonable backazimuth estimates. Due to the rather small aperture of these latter arrays, on the order of tens of meters, the F-K analysis of this kind of data may lead to intrinsically inaccurate slowness values, i.e. a rather extended peak in the wavenumber plane (Fig. 5c, d). If local wind condition, as also referred to by Negraru et al. (2010) for their slowness observations, near the array may have an effect in the case of such apertures, is speculative and thus not further considered here. Considering the celerity estimates for the observed signals, we find that a majority of signals may indeed be stratospheric returns, which is in agreement with the stratospheric branch postulated from observations at AlpArray stations by Fuchs et al. (2019). As we are both not overly successful in modeling these returns with ray tracing methods through the available atmospheric model, the Ingolstadt explosion appears an excellent opportunity for assessing our ability of realistically predicting infrasound arrivals through established or newly developed algorithms and/or sound speed models, especially at times close to the equinoxes.