Identification of the infrasound signals emitted by explosive eruption of Mt. Shinmoedake by three-dimensional ray tracing Open Access

Infrasound is inaudible sound of frequencies lower than 20 Hz and can propagate for a long distance with little attenuation. This characteristic is suitable for the detection of distant huge events such as tsunamis, typhoons, tornados, volcanic eruptions, and nuclear tests [e.g., Bedard and Georges (2000), Campus (2004), and Headlin et al. (2002)]. The Comprehensive Nuclear-Test-Ban Treaty Organization (CTBTO) has been establishing the International Monitoring System (IMS) consisting of 60 sites of infrasound sensors around the world (Christie and Campus, 2010) to estimate the source position and the scale of geo-hazards as well as nuclear tests.

It is predicted that Japan is threatened by a Nankai Trough earthquake of M8–M9 within 30 years (Headquarters for Earthquake Research Promotion, 2013) with a probability of ∼70%. In the course of this earthquake, it is also predicted that a tsunami with a height of more than 10 m will hit a large area along the Pacific Coast longer than 1000 km in the southwest region of Japan. For the mitigation of the tsunami disaster, we have developed an infrasound sensor, ADXII-INF01, in cooperation with SAYA Inc., and we have constructed a dense network of infrasound sensors in Kochi, Japan (Yamamoto and Yokota, 2015). The ADXII-INF01 covers the frequencies of $10−3$ to 6.25 Hz with a data polling rate of 4 Hz. The sensor network located more than 200 km away from Mt. Shinmoedake detected infrasound signals emitted by its explosive eruption on 10 March 2018. To clarify the nature of the observed wave signals, we have done an analysis of the observed signals by ray tracing calculations, and we have identified waves refracted from the troposphere, the stratosphere, the mesosphere, and the thermosphere.

In Sec. II, we show the observational data obtained from the explosive eruption of Mt. Shinmoedake. Section III describes a ray tracing formulation based on the Hamilton equations; we present equations for calculating the ray paths and formulas for travel time and travel distance of infrasound to the observation sites. In Sec. IV, we compare the results of the ray tracing calculations with the observed travel time and distance of the infrasound and identify the nature of the observed signals. Summary and discussion are given in Sec. V.

Mt. Shinmoedake is a part of the Mt. Kirishima cluster of volcanos on Kyushu Island, Japan. The volcano is located at 31.91° N, 130.88° E, and its elevation is 1412 m. It is recorded to have erupted in the years 1716, 1717, 1822, 1959, 1991, 2008, 2009, 2011, 2017, and 2018. A recent explosive eruption associated with atmospheric pressure pulses and a seismic tremor was observed on 6 March 2018. The explosive eruptions occurred that year more than 40 times in March and once in April and June, respectively.

During the explosive eruptions in March, atmospheric pressure pulses stronger than 100 Pa were observed several times by the microbarometer MB2005 installed by the Earthquake Research Institute (ERI), University of Tokyo. At the time of eruption, the microbarometer was located 3 km away from the crater of Mt. Shinmoedake and was located at 31.89° N, 130.86° E. The microbarometer covered a frequency range of 0.01–4 Hz, and its sampling rate was 100 Hz. The eruption of Mt. Shinmoedake observed at 1:54 (Japan Standard Time, JST) on 10 March 2018 produced the most significant atmospheric pressure pulses in the eruptions that occurred in 2018. Figure 1 shows atmospheric pressure pulses received by the sensor of ERI, indicating that the observed highest pressure was 610 Pa at 01:54:43 (JST) (Ichihara, 2018).

Our infrasound sensor network located more than 200 km away from Mt. Shinmoedake also detected and recorded the signals. Figure 2 shows the locations of our infrasound sensors on Shikoku Island, Japan. They are located in the region of $204<r<343$ km from Mt. Shinmoedake and $25 °<φ<37 °$⁠, where $φ$ is the azimuth angle of the sensors measured counterclockwise from eastward, taking the origin at Mt. Shinmoedake. Table I summarizes the distance r of the sensor sites measured from Mt. Shinmoedake and their azimuth angle $φ$⁠.

Figure 3 shows bandpass-filtered waveforms received by the sensors at the polling rate of 4 Hz at sites 1–6 in the frequency range of 0.05–0.3 Hz, which implies a spatial resolution of 6 km for each of the sites. The highest atmospheric pressure pulse detected by each of the sensors was between 1.0 and 1.5 Pa.

Table I summarizes the standard deviation σ for all sites. The standard deviation of the ambient sound pressure at each site was calculated on the basis of the observed data for 24 h on 10 March 2018. We regard the observed sound pressure higher than the 3.0 σ level to be the signals. To clarify the wavepackets in Fig. 3, we assume the signals detected intermittently whose time interval is shorter than 10 s to be one group; the time interval corresponds to a frequency of 0.1 Hz. At least two wavepackets are identified in Fig. 3. Furthermore, we identify the signals of the $3σ$-confidence level at sites 1–6 in time interval 3 as indicated by red arrows in Fig. 3. The durations of time intervals 1–3 at sites 1–6 are summarized in Table II. We shall examine these features in Sec. IV by the ray tracing of the infrasound emitted from Mt. Shinmoedake.

In the present observation, the propagation distance of the wave to the sites (200–340 km) is much longer than the wavelength (⁠$<6$ km). In such a situation, one can introduce the concept of a ray in discussing its propagation and can apply geometrical acoustics [e.g., Landau and Lifshitz (1959)], in which the wave is propagated along a ray whose tangent is parallel to the wavevector at any point.

The formulation discussed here follows that presented by Landau and Lifshitz (1959) and is analogous to calculating the motion of a particle using Hamilton equations [see also Blom (2019), Jones and Bedard (2018), and Lighthill (1978)]. Analogous Hamilton equations for a wave of angular frequency ω and wavevector k are given by

where r denotes the position vector. Note that ω corresponds to Hamiltonian and k to momentum in particle mechanics. The first equation gives the group velocity of propagation of the wave, and the second one determines the time variation of the wavevector. If the atmospheric conditions can be regarded to be in the steady state (⁠$∂ω/∂t=0$⁠) during wave propagation, we see from Eq. (1) that ω is kept constant along a ray path, since

Considering acoustic wave propagation for no wind in the atmosphere, one has the dispersion relation given by

where c is sound speed and $k=|k|$⁠. A wave is expressed by

with the use of the dispersion relation given by Eq. (3), where i is the imaginary unit. This dispersion relation is applicable under the condition that the effect of gravity and the Earth's rotation can be ignored [e.g., Scott et al. (2017)]. Namely, it holds for frequencies of $ω/2π≫g/2πc∼5×10−3 Hz≫Ω/π∼2.3×10−5 Hz$⁠, where $g=9.8 m/s2$ is the standard acceleration of gravity and $Ω=7.3×10−5 rad/s$ is the angular speed of Earth's rotation [see Jones and Bedard (2018)].

Equation (3) is a dispersion relation for no wind in the atmosphere. The dispersion relation in the presence of wind is derived in the following manner. Let K be a wind-blowing inertial system and K′ be another inertial system moving with wind velocity u relative to K. The position vector $r′$ and the wavevector $k′$ in K′ are given (Landau and Lifshitz, 1959) by

There is no wind in K′; thus, the wave observed in K′ is given by $ψ∝ exp [i(k′·r′−ck′t)]$⁠, where c is sound speed in the system with no wind. In consequence, the wave expressed in terms of the quantities in K is expressed by

indicating that the dispersion relation in the wind-blowing system K is given by

Let us examine the situation in which the wind velocity u is horizontal and varies with altitude z in the region of concern. We assume that u is uniform in the horizontal plane at each altitude z. Namely, the wind velocity field is written as $u=(ux(z),uy(z),0)$⁠, where we assume that u is steady during the propagation of the wave. The dispersion relation (7) is then expressed by

and the results of the Hamilton equations (1) are given by

Equation (12) indicates that the wavevector projected onto the xy-plane is constant in time.

We express constant k_x and k_y in terms of the spherical coordinates as

where k₀ is the wavenumber at t = 0, and θ₀ and $ϕ0$ are the initial zenith and azimuth angles, respectively. The dispersion relation (8) is then expressed by

where

with

being a horizontal unit vector parallel to $k0$ projected onto the xy-plane. Note that U is independent of θ₀.

If $ż$ turns to be $ż<0$ at a certain altitude, the wave emitted upward from z = 0 in the direction $(θ0,ϕ0)$ turns its path downward at this altitude. The turning altitude $zt$ is characterized by $ż=0$⁠, namely, $kz(zt)=0$ from Eq. (11), in consequence $k(zt)=k0 sin θ0$⁠. Equating ω given by Eq. (15) at $z=zt$ and the source altitude z₀ in view of the invariance of ω along a ray path [see Eq. (2)], one obtains the equation that determines the turning altitude $zt$⁠, namely,

The wavenumber k at altitude z is calculated by equating ω given by Eq. (15) at the altitudes of z and $zt$⁠. The result is given by

Since $k=kz2+kx2+ky2=kz2+k02 sin2θ0$⁠, Eq. (19) leads to

The time required to reach the turning altitude $zt$ is calculated with the use of Eq. (11), which yields

Substituting k given by Eq. (19) and k_z given by Eq. (20) expressed as functions of z into Eq. (21), we obtain

where the plus and minus signs should be adopted for the upward (dz > 0) and downward (dz < 0) propagation of the wave, respectively. The following equations [Eqs. (23) and (25)–(28)] are those for upward propagation. From the time reversal symmetry, the travel time from $z=z0$ to $z=zt$ is equal to that from $z=zt$ to $z=z0$⁠. Expressing the time for the wave to reach the turning altitude $zt$ from z = z₀ by $τ/2$⁠, we have

The x-coordinate of the turning point is calculated from Eqs. (9) and (11), which yield

with the use of Eq. (14). Substituting Eqs. (19) and (20) into Eq. (24), we obtain

The x-coordinate $X/2$ of the turning point at $z=zt$ of the wave starting from a source at the origin is given by

In a similar manner, the y-coordinate $Y/2$ is obtained to be

The end point $R=(X,Y,z0)$ of a path of the ray emitted from the source $(0,0,z0)$ is expressed by

with the use of Eq. (17).

In the following, we set $z0=0$ unless otherwise stated.

We calculate ray paths and travel time for the infrasound emitted from the source at Mt. Shinmoedake. We employ the Mass Spectrometer and Incoherent Scatter Radar Extended model (NRLMSISE-00) (Picone et al., 2002) for the vertical profile of the temperature. Figure 4(a) shows the vertical profile of the temperature T, and Fig. 4(b) shows sound speed c calculated by using the temperature profile. For the wind profiles, we use Horizontal Wind Model 2014 (HWM14) (Drob et al., 2015) and Modern-Era Retrospective analysis for Research and Applications, version 2 (MERRA-2) (Gelaro et al., 2017). HWM14 is an empirical model based the observations from the ground and satellites for 50 years and reproduces long-term averaged climatology of the atmosphere (Drob et al., 2015). MERRA-2, developed by the Global Modeling Assimilation Office of NASA (Gelaro et al., 2017), is an atmospheric reanalysis dataset based on satellite and ground-based observations. MERRA-2 provides reanalysis data with time resolution of 3 h and a spatial resolution of $0.5°×0.625°$ in longitude and latitude (Global Modeling and Assimilation Office (GMAO), 2015). Figures 4(c) and 4(d) show the vertical profiles of zonal (positive to the east) wind u_x and meridional (positive to the north) wind u_y. We interpolated u_x and u_y given by MERRA-2 at 0:00 and 3:00 JST to deduce those at 2:00 on 10 March 2018 in the area of Mt. Shinmoedake (31.9° N, 130.9° E). The HWM14 profile is shown in blue and the MERRA-2 profile in red for comparison.

In the following, we adopt HWM14 for the altitude of $z≤200$ km, which is set to include all rays falling onto all of the sites. For the rays confined in $z≤76$ km, we use MERRA-2.

Figure 5 shows examples of paths of the rays emitted toward the site 5 (⁠$ϕ0=37°$⁠) for the initial zenith angles $θ0=30°, 65°$⁠, and 85°. Here, we used HWM14 for the wind profile and set $z0=0$⁠, since the turning altitudes are much higher than the source altitude as will be seen below. Figure 5(a) indicates the ray paths projected onto the xy-plane together with those for the windless case by a dashed line for reference. The deviation of the end points of the sound rays from the line connecting the source and site 5 less than 3°. The deviation angles for the other sites are around 3° as well. Figure 5(b) shows paths of those rays projected onto the xz-plane, indicating that the rays return to the ground for $θ0≳30°$⁠, for which the turning altitudes are $zt∼110$ km.

The observed signals include the rays that turn at some altitudes. Turning altitude $zt$ for the sound emitted toward the direction (⁠$θ0,ϕ0$⁠) is given by the solution of Eq. (18), namely, the solution of the equation

where

Figure 6 illustrates $f(z;ϕ0)$ toward site 5 (⁠$ϕ0=37°$⁠). For a turn of a ray to occur, θ₀ must be larger than a certain value because $f(z;ϕ0)$ as a function of z is bounded at the altitudes z of interest, i.e., $0<z<200$ km. The turning altitude $zt$ of the sound emitted toward the initial azimuth angle $ϕ0$ is given by the value of z at the intersection of the curve $f(z;ϕ0)$ and the vertical line at $1/ sin θ0$ in Fig. 6. If there are multiple solutions that satisfy Eq. (29), the lowest $zt$ gives the actual turning altitude within the framework of geometric acoustics. Figures 5, 6, and 8 indicate that there are four possible turning altitudes. They are in the lower thermosphere [$zt>110$ km; Figs. 5(b) and 6(a)], the mesosphere [$60<zt<66$ km; Fig. 6(b)], the stratosphere [$42<zt<45$ km; Fig. 6(b)], and the troposphere [$zt∼10$ km; Figs. 6(b) and 8].

Figure 7 shows the relation of the travel time τ and horizontal range r of the rays that fall on the line $ϕ=φi$ (⁠$i=1,..,6$⁠). The red and blue shaded areas indicate time intervals 1 and 2, respectively. Here, we define the “end point” of the ray by a circle of radius 3 km in view of the limitation on the spatial resolution of the wave of frequency $∼0.05$ Hz. The intervals of the initial azimuth and zenith angles are taken to be $25°<ϕ0<40°$ and $0°<θ0<90°$ with their increments of 0.1° in drawing the figure of τ versus r.

We can expect that the signals detected in time interval 1 are a direct wave propagated horizontally near the ground surface (⁠$z0=0$⁠). Apparent sound speed estimated from the distance and observed travel time between the source and each of the sites is 317 m/s, which corresponds to a surface temperature of $−25°C$⁠. On the other hand, the sound speed is calculated to be 334 m/s if we adopt the surface temperature of $3.5°C$⁠, which is the temperature averaged over five places near Mt. Shinmoedake and the observation sites (observation data by the Japan Meteorological Agency). This result suggests that the signals detected in time interval 1 are not simply sound waves propagated horizontally.

To reveal the reason why the observed speed of the wave is slower than the calculated one, we investigate the wave propagation in the troposphere in detail. In the propagation in the troposphere, the source altitude cannot be set to be $z0=0$ km because the elevation of Mt. Shinmoedake cannot be ignored compared with the altitude of the tropopause. We carried out ray tracing calculations by varying the source altitude in $0<z0<6$ km. We found that the rays arrive at all sites and are propagated near the ground surface for the source at $z0=2.5$ km. Figure 8(b) shows some of the paths of the rays emitted from the source at $z0=2.5$ km toward the direction of the site 5 $(ϕ0=37°)$ for $θ0=75°$⁠, $80°, 86°$⁠, and $88°$ together with $c(z)+U(z;ϕ0)$ in Fig. 8(a). The rays emitted at $86°<θ0<90°$ are confined in $2<z<10$ km, those at $85°<θ0<86°$ are confined in $0<z<10$ km and near the ground surface, those at $78°<θ0<85°$ turn once at $z∼10$ km and fall to the ground surface, and those at $θ0<78°$ do not turn at the troposphere. The figure indicates that the ray is not propagated straight horizontally but undergoes multiple refraction, which occurs because the z-component of the wavevector tends to be propagated downward at altitudes of $dc/dz>0$ and $dux/dz>0$ and vice versa [see Eq. (13)]. In the present case, the turns at the high altitude are mainly the result of the tropospheric jet stream around $z∼10$ km (see Fig. 8), and those at the low altitude are due to the negative temperature gradient (⁠$dT/dz<0$⁠). Net horizontal propagation speeds deduced from the travel times are 316–320 m/s for sites 1–6. These net speeds agree with the observed apparent sound speed. This result suggests that the signals detected in time interval 1 could be the waves refracted multiple times in the troposphere.

In time interval 2, Fig. 3 indicates that the signals are detected at all sites. We choose the first and second maxima of the sound pressure detected at sites 1–6 in time interval 2. The travel time calculated for sites 5 and 6 reproduces the observed ones within an error of 10 s as shown in Table III. These signals detected at sites 5 and 6 are identified to be the rays turning from the stratosphere and the mesosphere as seen in Fig. 7. For site 4, the observed first signal is identified to be the ray turning from the stratosphere, whereas its travel time of the ray turning from the mesosphere is just outside time interval 2, and its end point is just outside site 4. The signals observed at sites 1–3 in time interval 2 have not been reproduced in the present calculations. One of the reasons why the calculations do not reproduce some of the observed signals is an underestimate in the stratospheric and mesospheric wind speeds in the wind models of HWM14 and MERRA-2 because of the space and time variation of wind. See also Garcés et al. (2004) and Le Pichon et al. (2005).

In time interval 3, the observation shows that signals are detected at all sites as indicated by the red arrows in Fig. 3. The observed travel times of the signals at sites 1–3 are less than 40 s apart from the end of time interval 2, while those at sites 4–6 are more than 150 s apart from the end of time interval 2. According to the ray tracing, the signals detected at sites 4–6 are identified to be the rays turning from the thermosphere, and no ray turning from the thermosphere arrives at sites 1–3 as shown in Fig. 7. This result is consistent with the observation. However, the observations indicate that two signals arrived at the sites 5 and 6 and one signal was detected at site 4. The reason why only one signal was detected at site 4 is unclear at present. As for the observed signals close to time interval 2 at sites 1–3, the calculations predict that no ray arrives at sites 1–3.

On the basis of the dataset obtained by the densely distributed infrasound sensors on Shikoku Island in Japan, we have identified characteristic impulsive signals coming from the volcanic eruption at Mt. Shinmoedake on Kyushu Island on 10 March 2018 at all of the six observation sites, which are located between 200 km (site 1) and 340 km (site 6) away from Mt. Shinmoedake. We applied the ray tracing calculations to the analysis of the infrasonic waves in the frequency range between 0.05 and 0.3 Hz (wavelength range between 6 and 1 km) using the model atmosphere including the temperature profile provided by NRLMSISE-00 and the wind profiles provided by HWM14/MERRA-2.

We have shown that the observed signals can mostly be identified with one-to-one correspondence. The first waveform group having an observed propagation speed around 317 m/s is likely to be the waves propagated undergoing multiple refraction around an altitude of 10 km in the troposphere. On the other hand, the second waveform group detected at sites 4–6 is identified to be the waves turning from the stratosphere around 42 km and the mesosphere around 60 km. We have identified the third waveform group detected at the sites 4–6 to be the waves refracted from the lower thermosphere around 110 km. However, there are several signals that have not been reproduced by the ray tracing calculations such as the waveforms in the second and third groups detected at sites 1–3. One of the possibilities is that the stratospheric and mesospheric wind speeds were higher at the time of the event than those given by the wind model HWM14/MERRA-2 used in the present study. Another possibility is the limitation of the ray tracing method, which ignores wave behavior such as diffraction and scattering [see, e.g., Millet et al. (2007), Sabatini et al. (2019a), and Sabatini et al. (2019b)]. The present result demonstrates the importance of a dense and widely distributed infrasound observation network for the unambiguous identification of various wave signals of the infrasound emitted by not only the volcano. Such observations are possibly useful also for the retrieval of subtle structures of wind and temperature, combined with ray tracing calculations and, if necessary, more detailed ones [e.g., Sabatini et al. (2019a) and Sabatini et al. (2019b)].

We thank Dr. M. Ichihara for providing the data obtained by the microbarometer set by ERI near Mt. Shinmoedake and Dr. S. Watanabe and Dr. Y. Kakinami for fruitful discussion. The data provided by MERRA-2, HWM14, and NRLMSISE-00 are acknowledged. This work was partly supported by Ministry of Education, Culture, Sports, Science and Technology/Japan Society for the Promotion of Science (MEXT/JSPS) KAKENHI Grant No. 17H02062 and SECOM Science and Technology Foundation and ERI Joint Usage/Research Program (JURP) Grant Nos. 2014-B-15 and 2018-B-04.