3.1 Determination of volcanic eruption time
Numerous international media on the Internet reported that the Tonga volcano erupted at about 17:10, but few reported more precise eruption time or how many noticeable eruptions occurred in the eruption sequence, resulting in the lack of many eruption details. Considering that the travel time error of seismic phase is generally less than 0.5min in a not very large distance, the travel time of seismic wave is used to determine the time of each eruption event. Some stations whose data can be downloaded from the IRIS seismic database are retrieved, in which Iris-AFI and Iris-MSVF are the nearest stations we obtained. Iris-MSVF is located 757km northwest of the Tonga volcano, and Iris-AFI is located 832km northeast of the Tonga volcano. Additionally, according to the pictures of volcanic eruptions reported on the Internet, the coordinate points (20.558 ° S, 175.382 ° W) are obtained by comparing the terrain of Google Earth. The vertical components of the two seismic stations are decomposed into 11 layers by Meyer wavelet, and then the lower parts of high-frequency components of d9, d10 and d11 layers (a11 and d11 corresponding to the lowest low-frequency and lowest high-frequency components) are extracted for reconstruction. Thus, we are able to obtain the obvious Rayleigh surface wave and some body waves with small amplitude caused by volcanic eruption (Fig. 5).
In simulation of the seismic wave, the full wave-field synthetic seismograms of QSEIS program are obtained from Professor Wang Rongjiang of GFZ in Potsdam, when he taught at Chinese Academy of Sciences (CAS) in 2017. The QSEIS program is based on the semi-analytical wave field theory of orthogonal decomposition. It can be applied effectively on the multi-layer medium model of local and teleseismic earthquakes. Given that the empirical theory is extremely complex with excessive formulas, we do not introduce it in detail; instead, we list the key parameters involved in the program setting in the appendix for readers, and the main theories can be referred to Wang (1999) and Wang et al. (2017).
According to the synthetic seismograms and the actual record, it can be seen there are two groups of Rayleigh waves at 04:18:39UTC and 04:21:23 UTC at the Iris-AFI station (Fig. 5). Then, by subtracting the travel time of Rayleigh waves, the corresponding eruption time of the first group of surface wave is 04:13:55UTC. The corresponding time speculated by Iris-MSVF station is about 2 seconds earlier, revealing very close results. However, the arrival of the body wave is about 20 seconds later than the actual record. Considering it is more ordinary that the travel time error of surface wave is more than ten seconds caused by the velocity model at the epicenter distance of about 800km, it is generally speculated that the first strong eruption is around 04:13:35UTC, that is, 20 seconds forward, which is more appropriate for the travel time of body wave. The second eruption time is about 240 seconds after the first (i.e. 04:17:35 UTC).
For the largest eruption, according to the similarity of waveform, it is speculated that the surface wave arrival is about at 04:31:30UTC on Iris-AFI and the deduced first arrival of body wave can be seen at about 04:28:20UTC at the two stations. Therefore, we conjecture that the starting time of the maximum eruption is about 04:26:26UTC. By the way, these three speculated eruption times are very close to the results that surmised by Yuen et al. (2022). In addition, from the abundance and complexity of phases in simulation or the actual Rayleigh waves, we further speculate there are more eruptions in the duration of the latter two eruptions, but with smaller recorded amplitudes. Meanwhile, the phenomenon possibly indicates that the top depth of magmatic activity at the second eruption is shallower than that at the first eruption. Namely the confining pressure is smaller for the second eruption without larger Rayleigh amplitude but within a longer duration. Similarly, the top depth of the magmatic activity of the third (largest eruption) may be much shallower. Maybe in the process of the largest eruption, some magma has lost the volcanic dome cover, or the crater diameter has been expanded several times, corresponding to the much smaller amplitude of surface wave and much longer duration. According to the latest research (O’callaghan, 2022), in May 2022, Professor Cronin led a group to sail over the volcano’s caldera, and indeed found that the four-kilometer-wide caldera had dropped about 650m. This may be a direct reflection as the deficit of the calderas (Mcbirney, 1990). Figure 5 also shows many intermittent eruptions in the subsequent time after the largest eruption.
3.2 Infrasound simulation and determination of tremor phases
In online reports around the world a few days after the Tonga volcano erupted, most of them employed the keyword "shockwave" to describe the significant short-term atmospheric pressure disturbance. However, from the analysis of the pressure disturbances recorded in Chinese Mainland, it can be found that their shape of initial onset is not sharp. Unlike the air wave generated by the explosion in the near field (Li, et al. 2018), it is more likely an infrasound wave, or a more professional terminology - Lamb wave (Francis, 1973). In order to reliably identify their attributes, the program QSSP of full wave-field synthetic seismograms provided by Professor Wang Rongjiang in 2017 is applied, which comprehensively considers the structures of atmosphere, ocean, mantle, core and self-gravity. The source depth is set at 3km above the Earth's surface (In the actual simulation, several different heights are used for calculation, in which the result of placing the source in the upper part of the stratosphere is reversed from the observed waveform) to analyze the infrasound wave generated by a model of spherical explosion. Considering the lengthy theory of the QSSP, we do not introduce them in detail in this paper, but the key parameters involved in the program setting are listed in the appendix for readers, and the relevant theories and detailed derivation of formulas can be further referenced (Alterman et al., 1959; Wu, 1997; Wang et al., 2017).
After a series of trials and errors of nearly 300 combinations in the past three months, the infrasound source time function composed of 14 eruption events has been accomplished. Why to select the source time function? Two factors are needed to be considered. First, the calculated synthetic waveform must be highly similar to the most of clearly recorded infrasonic waveforms in China, which includes amplitude, period and coda shape. Second, it should be also quite consistent in the actual travel time. Admittedly, there may be more small eruptions in the actual process, as reflected in the nearby monitoring of the eruption process of the Redoubt Volcano by McNutt et al. (2013). Vergoz et al. (2022) also mentioned that during the eruption of the Tonga volcano, some seismic stations had recorded a small number of earthquakes with mb4 ~ 5, although their SNRs are much lower. Meanwhile, witnesses 90km away from the volcano said they heard cracking and noise-like artist fire. Professor Cronin further speculated that there are some large eruptions and hundreds of small-scale eruptions (O'callaghan, 2022). Thus, the source time function composed of 14 eruption events is only a rough approximation. The atmospheric model used in the simulation is based on the American standard atmospheric static model (NOAA, 1976; Lewis, 2022) that is embedded in the instructional examples of QSSP, with a top height of 86km and a division in vertical with one layer per 1km (Fig. 6).
From the input parameters of QSEIS and QSSP programs, it can be seen that the focal depths setting for calculating seismic waves and infrasound waves are different. The focal depth of seismic waves in the crust is 3km underground, according to the depth estimation of the small eruption of the Changbaishan volcano in Northeast of China by Wei (2010) and Guo et al.(2015) and the speculated burial top depth of magma in the Tonga volcano by Cronin et al.(2017). However, the focal depth of infrasound wave propagation is 3km above the ground. They are cautiously set based on the practical observation. Although the sources of the two kinds of waves are both generated by the Tonga volcanic eruption, the direct sources should be different. The source of the crustal surface wave is still the mechanical vibration caused by the interaction between the magma chamber and the solid rock (which can also be caused by the counterforce of the big air explosion, but the direct source is still in rock), and the source of the air wave should be in the air near the ground and outside the volcanic vent. In addition, although the time of the three major eruptions is inferred from the seismic waves, there are still some different stages in the complex evolution, which includes these magma gushes from the underground to the surface, the high-temperature magma interacts with the low-temperature atmosphere, and finally the interaction evolves into infrasound wave propagation sources of different scales in the air. Thus, the start time of outward propagation of the two propagation sources (in solid and in air) may contain a more or less underlying time difference.
It should be noted that there are maybe more or less subjective assumptions in the above analysis. Fortunately, the difference between the focal depth and the possible actual depth is not too large, and the influence on the propagation travel time can be ignored. Based on the above considerations, the source time function settings in QSEIS of seismic wave simulation and QSSP of infrasonic wave simulation are also different, but the difference is in the order of tens of seconds. That is, in QSSP, the respective eruption time does not completely follow the time interval inferred from the seismic wave, but there is only a deviation of tens of seconds between each other (Fortunately, the sampling interval of the barometer is 1min that is difficult to distinguish the difference). Finally, the ultimate judgment is essentially based on the better fitting of the far-field pressure records. Moreover, it can be seen from Fig. 5 that there are still some differences between the original seismic recorded waveform (gray) and the part (black) after wavelet extraction, and the waveforms extracted by the two stations are also not highly consistent. These two features reveal that the underground eruption process of the Tonga volcano is indeed extremely complex. The initial eruption on the upper part may trigger the re-distributing of pore pressure and confining pressure in lower magma, and then trigger more anisotropic eruption events in the longitudinal and transverse directions (Heap, et.al, 2015). Incidentally, Cronin et al. (2017) previously carried out a bathymetric survey showing dozens of large or small craters that distribute around the Tonga volcano within 50km along the northeast direction. This survey further illustrates that there is a connection between the deeper magma of these volcanoes, which can result in the complexity of the source of observed seismic waves. If the seismic stations can be jointed surrounding the Tonga volcano within 100km, more inspired analyses will be obtained in future work.
The specific results of infrasound simulation are shown in Fig. 7 ~ 8 and table 1. In Fig. 7, the propagation snapshots of 1、3、5、7、9 and 11 hours are drawn on the Earth's surface in an axisymmetric manner. Figure 8 shows the time records at several different epicentral distances, among which the actual epicentral distance of Songpan, Jixian and Shantou stations (see Fig. 10) are used in calculation for subsequent comparison with the recorded waveforms. From the early propagation stage of the simulation, in the near-field region (i.e. within 3000km), the source effects of 14 sub-events can be distinguished, leading to much more high-frequency compositions and complexity in the waveform. Then, the different wave groups gradually show attenuation phenomenon in high-frequency and continuous merger. About 5 hours later, two predominant phases (Phase1, Phase2) are separated, and with the increase of epicenter distance, the travel time difference of these two phases also gradually increases. Indeed, there is the second phase in the actual record (Phase2' in Fig. 8 and Phase2 in Fig. 3), but the second phase is later than the theoretical travel time, especially at Jixian station, 70km east of Beijing, which may be delayed half an hour. How does this delay generate? We speculate that the Phase2 is actually a wave of Lamb (Francis, 1973) that alternately transmits forward between the top of the mesosphere and the Earth's surface with a high velocity (Fig. 6). The typical velocity of Lamb wave is 254m/s; however, it obviously has been affected by the winter westerlies in the mid latitudes of the Northern Hemisphere. Affected by the backward wind of the westerly belt, their delay can be explained. At the same time, the attenuation effect at high altitude is more obvious (Reed, 1972; Francis, 1973), thus the amplitude is smaller. In contrast, Phase1 propagates along the troposphere. As there are more obvious difference of tropospheric transversal temperature in each region, it can be seen that the deviation between the actual travel time and the theoretical travel time is significantly negatively correlated with the near-surface air temperature (Fig. 11&12). As for the greater delay of Jixian station, we speculate that it is further related to the sharp cooling in the adjacent areas of Beijing, especially within 200km northwest of Beijing since the noon of January 15 (A cooling of 12 ℃ appeared in a half day6; Haidian Meteorological Bureau even released7 a blue warning of strong wind of force 7, Fig. 13b), and the strong northwest wind. Phase3 in Fig. 3 should be the infrasound waves transmitted around the other side of the Earth, which can be deduced from the sum of travel time (~ 36.1h) of Phase 1 and Phase3 that travel just one circle of the Earth. Therefore, the measurement of the parameters of each pressure disturbance in Fig. 4a is actually the measurement of the group velocity of infrasound waves propagating in the troposphere, which reflect a result of the comprehensive superposition of multiple eruption events. Among them, the component of the largest eruption accounts for the largest proportion, and thus the time starting point of measurement starts from the time of the largest eruption.
Indeed, it should be noted that the seismic moments used in seismic wave simulation and atmospheric infrasonic simulation are different. The 0.15e + 18N•m is equivalent to an MS5.4 earthquake, and 0.63e + 18N•m is equivalent to an MS5.8 earthquake. They are the magnitude estimation of the two foreshocks before the maximum eruption. Since the downloaded IRIS data are recorded as Count Number and lack of sensitivity parameters, it is prior to compare with the waveform similarity but not rely on the actual physical amplitude. Fortunately, our aim of using seismic wave simulation is merely to constrain the eruption time. However, the 14 seismic moments such as 0.0311e + 19N•m (equivalent to an MS5.6 earthquake) and 0.1693e + 19N•m (equivalent to an MS6.1 earthquake) used for infrasound simulation have clear physical significance, because the calculation results are very close to the 2hPa amplitude of the actual atmospheric pressure record. By comparison, then, the magnitudes of these seismic moments are less than the magnitude of 500 ~ 1000 Hiroshima atomic bombs estimated on the Internet (the magnitude reflected by the cumulative intensity of all infrasonic sources is about MS6.6, that is, 1.2e + 19 joules and less than MS7.0; The equivalent magnitude of Hiroshima atomic bombs is referred to about MS5.3 according to Chinese traditional report), but larger than fourteen times of the estimation of 200MT TNT given by Vergoz et al. (2022). Interestingly, the online reports are greatly different in the magnitude estimation of the Tonga volcano released energy. Our results are between the two aforementioned values.
Table 1 The maximum and minimum pressure of simulation at different epicentral distances
(The value with * is the actual record)
Δ/°
|
4
|
12
|
20
|
28
|
36
|
44
|
52
|
60
|
68
|
79.492
|
86.968
|
93.803
|
79.492*
|
86.968*
|
93.803*
|
Max/Pa
|
5629.7
|
1228.0
|
539.4
|
377.8
|
297.3
|
232.6
|
195.7
|
178.9
|
158.5
|
144.0
|
140.3
|
133.6
|
151.7
|
142.7
|
141.7
|
Min/Pa
|
-3270.4
|
-693.1
|
-229.3
|
-244.5
|
-206.3
|
-153.3
|
-128.8
|
-114.7
|
-102.3
|
-80.5
|
-72.4
|
-67.4
|
-83.3
|
-92.1
|
-57.8
|
3.3 Determination of infrasound wave parameters
Based on the previous analyses, the estimated time of the maximum volcano eruption at 04:26:26 UTC is obtained, and the attributes of tremor phases are also analyzed. As a result, the travel time of the strongest pressure disturbance to each barometric station can be easily calculated. Two methods are used to measure the wave speed in the following. One is to use the slope of distance and travel time (Fig. 9a). The other is scanning by 0.5m/s steps within a certain speed range to obtain the speed by the smallest residuals between the theoretical travel time and the actual travel time (Fig. 9b). In Fig. 9a, the linear regression approach (Chatterjee&Hadi, 1986) is applied to all stations or part of stations (overlapped by the red line), and the wave speeds obtained are 296.1m/s and 320.0m/s respectively. Their average values are close to 308.5m/s, which corresponding to the smallest absolute error of travel time difference in Fig. 9b. The convergence speeds of the square root and square root of the absolute value (i.e. the residual is weighted differently) are also around 308m/s; hence, the propagation speed of the primary infrasound wave is recommended as 308m/s. Figure 9c shows the comparison between the distance curve obtained using a 308.5m/s speed and the distance curve calculated with coordinates, and Fig. 9d shows the corresponding specific time difference of the 462 stations.
Figure 10 shows the spatial distribution of the arrival time of infrasound waves at each station, and Fig. 11 shows the spatial distribution of time deviations between the theoretical and actual times. Figure 12 shows the average near-surface air temperature of ten minutes before and after the arrival of the peak of infrasound waves, Fig. 13 shows the period distribution, and Fig. 14 shows the amplitude distribution. According to the time deviation distribution, bounded on the line along the cities of Kashgar, Lanzhou, Xi'an, Hefei and Shanghai, clear late arrivals in the north and early arrivals in the south can be observed. Among them, arrival times of about 9 stations are still not consistent with the surrounding stations, with a difference of approximately 10 minutes. Two factors may account for the outliers of the deviations. On the one hand, the deviations may be attributed to the obvious wind interference on the propagation. On the other hand, the system built-in timing of a few stations may have a deviation of 1 to 2 minutes. In addition, the stations in the northwest of Beijing present obviously late arrivals.
Period statistics indicate that the average period of the 462 stations is 70min. Except for the adjacent regions of Beijing, most of the large and small periods are generally staggered with each other. In the estimation of period by Welch periodogram, it can be inferred that the period stability is obviously affected by the SNR. Furthermore, the high-frequency disturbance near the peak of the infrasound wave will reduce the period, whereas the high-frequency disturbance on both sides of the peak slope will increase the period estimation. In addition, from the view of the mathematical stability of spectrum estimation, the period estimation accuracy of a single bulge is prominently not as good as the repeated signals. Therefore, the results in Fig. 13 are only for references, and the average value has more practical significance.
The amplitude in Fig. 14 is measured from reconstructed data with two modified steps. The first step is Meyer wavelet decomposition and the second step is 3-point smoothing. In the specific measurement process, we select the arrival time as the time reference point, expand it for 40min before and after to extract the maximum and minimum values of the smoothed curve, and then subtract the maximum and minimum values to obtain the amplitude. The distribution results show that the average value is about 2.0hPa, and the amplitude of infrasound wave is relatively smaller in the area with lower background pressure (Fig. 15). The collective and larger amplitude at stations within 200km around Beijing is obvious, but there is no clear raise of the background pressure in this area (these mentioned stations correspond to the red points in Fig. 15). The schematic diagram of Fig. 16 is used to illustrate the phenomenon that when the phase_a' encounters the atmospheric cold front near Beijing propagating along the troposphere, the period and propagation speed decrease, but the amplitude increases. The travel time deviation of the Phase_b' and Phase_b propagating along the stratosphere gradually increases with the increase of propagation distance, due to the westerly influence of the westerly belt with the velocity of 20 ~ 40km/s (Liu, et al., 2010). As a result, Phase2 and Phase2' in Fig. 8 do not appear at the same time. In contrast, there is a small temperature raise on the line from Lanzhou to Xi'an; thereafter, the period increases and the infrasound amplitude decreases.
[6] https://tv.cctv.com/2022/01/15/VIDE4zRukfRt84aEfIQUQ328220115.shtml
[7] https://www.tianqi.com/alarmnews/2201141611010841.html