Inertial effects due to eruption-induced atmospheric disturbances identified by superconducting gravimeter observations at Matsushiro, Japan

The violent eruption of Hunga Tonga-Hunga Ha’apai volcano on January 15, 2022 induced strong atmospheric disturbances, which traveled around the world as atmospheric Lamb waves. When this wave passed through the superconducting gravimeter station at Matsushiro, Japan, a large signal of gravity changes was recorded. Also, barometers installed around Matsushiro recorded wave trains of pressure changes. Analysis of the barometer data revealed that the atmospheric disturbances traveled as plane waves. Applying the theory of atmospheric loading for traveling plane waves, the observed gravity changes were well reproduced by a sum of three components of atmospheric loading, namely, Newtonian, free-air and inertial effects. In particular, the inertial effect of atmospheric loading, which is rarely observed, was clearly identified in the gravity data. From the theoretical modeling, an estimate of rigidity in the shallow region of the Earth was also obtained.


Introduction
Atmospheric loading has significant effects on surface geophysical observations using highly sensitive instruments including the superconducting gravimeter (Prothero and Goodkind 1968;Warburton and Goodkind 1977;Goodkind 1999;Hinderer et al. 2007). Usually, this effect is considered as a combination of the following two agents (Hinderer et al. 2014). For one thing, a positive change in atmospheric pressure causes larger upward Newtonian attraction by atmospheric mass, resulting in a decrease in gravity (Newtonian effect). For another, a positive change in atmospheric pressure causes downward displacement of the ground, resulting in an increase in the gravity measured on the surface of the Earth (freeair effect). In total, a decrease in gravity is observed for an increase in atmospheric pressure, because the former effect is larger in magnitude. These effects are commonly described using a single factor of proportionality representing the magnitude of a change in gravity acceleration for a unit change in surface atmospheric pressure. This factor, called atmospheric admittance, takes a value close to −3 × 10 −11 ms −2 Pa −1 (the minus sign means that a positive pressure change causes a negative change in gravity), and varies slightly depending on the location as well as on temporal and spatial scales (Hinderer et al. 2014). Considering the precision of the superconducting gravimeter, this is a really large effect, and many efforts have been made to improve the precision of atmospheric corrections. They include frequency-dependent admittances (Neumeyer 1995;Kroner and Jentzsch 1999;Crossley et al. 2002), 2D or 2.5D atmospheric models (Merriam 1992;Mukai et al. 1995;Boy et al. 2002) and 3D atmospheric models (Neumeyer et al. 2004;Boy and Chao 2005;Klügel and Wziontek 2009).
In addition to the Newtonian effect and the free-air effect, Zürn and Wielandt (2007) pointed out, in their theoretical considerations on the minimum of background seismic noise spectra in the vertical component, that there should be an inertial effect of atmospheric loading on surface acceleration measurements. Among the several models considered in that paper, they investigated the case of acoustic-gravity waves based on the work by Sorrells (1971), and provide necessary formulae for calculating gravity effects. Because the inertial effect originates from the acceleration of the ground to which the sensors are fixed, it does not show up at zero frequency, and can emerge only at higher frequencies.
Although Zürn and Wielandt (2007) did not demonstrate a clear observational evidence in real records, Zürn and Meurers (2009) showed an example of phase reversal in the atmospheric effects at millihertz band observed by the superconducting gravimeter C025 at Vienna, which implies the existence of the inertial effect of atmospheric loading. Forbriger et al. (2021) also investigated the inertial effect of atmospheric loading on observations with seismometers and gravimeters at Black Forest Observatory, based on the Bouguer plate model. Although daily variations in the atmosphere have spectral power in such frequency bands, they seldom appear as a clear signature of the inertial effect in the records of gravimeters or seismometers, because such variations are spatially random. Spatially coherent variations in the atmospheric pressure are needed so that the inertial effect of atmospheric loading can be identified by surface gravity or seismic observations.
On January 15, 2022, the Hunga Tonga-Hunga Ha'apai volcano in Tonga erupted violently. This event induced strong disturbances in the atmosphere which propagated Imanishi Earth, Planets and Space (2022) 74:54 several times around the world as atmospheric waves. The waves are considered to be atmospheric Lamb waves (Lamb 1910;Lindzen and Blake 1972;Gill 1982). The atmospheric Lamb waves are classified into a special branch of propagating waves in the stratified atmosphere, and propagate along the bottom boundary of the atmosphere (Watada 2009;Arai et al. 2011). The atmospheric Lamb waves are not dispersive, meaning that they do not change waveforms as they travel except for the effect of energy dissipation. After traveling about 8000 km from Tonga, the waves passed Japan about 8 h after the eruption. A local barometer network built around the superconducting gravimeter station at Matsushiro, Japan (Imanishi et al. 1997(Imanishi et al. , 2004 recorded wave trains of the pressure signals from this event. The recorded pressure signals were almost identical, indicating that the waves took the form of plane waves and therefore that highly coherent variations in the surface pressure had occurred in this region. Also, gravity changes associated with the passing of atmospheric disturbances were precisely recorded by a superconducting gravimeter at Matsushiro. In this paper, we show a clear evidence of the inertial effect of atmospheric loading, by analyzing both the gravimeter and barometer recordings. This analysis also led to rough estimation of rigidity of the Earth.

Characterization of atmospheric disturbances
We start with characterizing the atmospheric disturbances observed at and around the Matsushiro station (see Figure 1). At Matsushiro, the superconducting gravimeter iGrav #028 has been in operation since February 2019. In an effort to improve the precision of atmospheric correction for the superconducting gravimeter, we built a local network of barometer observations around Matsushiro. Four of such barometer stations are still in operation, as shown in Figure 1 and Table 1. The sampling interval of the barometers and the gravimeter is 1 s. Figure 2 shows the barometer data from the five stations including Matsushiro. Since the barometer stations are located at different heights, absolute values of barometric pressure are different from station to station (Table 1). Aside from the difference in the absolute values and arrival times, it is evident that the barometers at the five stations have recorded very similar wave trains of atmospheric signals from the Tonga volcano. To extract the temporal changes in atmospheric pressure from the volcano, a linear trend was estimated and removed from the data for each station. Extracted signals were subject to calculation of cross-correlations, from which differences in arrival times with respect to Matsushiro were estimated. The result was − 115 s (advance) for Komoro and + 83 s (delay) for Ohmachi. Here we excluded Matsumoto and Nagano stations, because the clocks of the recording systems at these stations were incorrect at the time of the Tonga event, due to lost connection to internet time servers. We verified that the wave trains of the atmospheric disturbances from the three stations were almost identical if the data from the Komoro and Ohmachi stations were shifted by the estimated time differences.  Table 1. The square indicates the area in which Newtonian attraction by the atmospheric mass is calculated. Atmospheric waves from the Tonga volcano travel from the south-east direction through the north-west direction Given this result, we can conclude that the observed waves propagated as plane waves. The assumption of plane waves will be justified considering the distance (about 8000 km) from the Tonga volcano and the size (about 100 km × 100 km) of our study area. The velocity and azimuth of the waves can be estimated using the difference in arrival times between stations. Figure 3 shows the time difference at each station with respect to Matsushiro, plotted against the distance from the volcano. Fitting a linear function to these data while taking into account the uncertainties (1 s) in time difference estimation gives (3.49 ± 0.02) s km −1 as the slowness of the waves, in other words, (286 ± 2) ms −1 as the propagation velocity. This value is close to, but significantly lower than, 310 ms −1 , the value typically quoted as the propagation velocity of the atmospheric Lamb waves (e.g., Nishida et al. 2014). One possible explanation of this discrepancy is the effect of background eastward winds in the layer of the atmosphere where Lamb waves occupy, which is particularly strong around Japan in winter seasons. For the present, we identify the observed waves as atmospheric Lamb waves without going into meteorological details. This does not have serious influence on the following analysis of the atmospheric loading.
The azimuth was estimated to be N45W. This means that the waves came from the south-east direction and traveled to the north-west direction.

Evaluation of atmospheric loading effects on gravity
Since we have collected necessary information on the Lamb waves from the Tonga volcano, now we can calculate each of the three (Newtonian, free-air and inertial) terms of atmospheric loading that would have affected the gravity observed at Matsushiro. For the calculation, we use the data of atmospheric pressure recorded at Matsushiro as representative of the whole study area. Zürn and Wielandt (2007) gives formulae for calculating the effects of atmospheric loading for several theoretical models. Here we adopt the "acoustic-gravity wave model" (Sect. 4.3) of that paper, because we deal with traveling plane waves. For a single angular frequency ω, the three equations (13), (15) and (16) of Zürn and Wielandt (2007) give the relevant Fourier components of Newtonian, free-air and inertial effects, respectively. Therefore, the necessary steps to be taken here will be (i) Fourier transform the time domain data of surface atmospheric pressure, (ii) multiply the Fourier components by the factors in the right-hand sides of the three equations and (iii) inverse Fourier transform the frequency domain data back to the time domain. These procedures are quite simple to implement. However, we found that applying this method to our data gives an unrealistically large value (larger than 9 × 10 −9 ms −2 ) for the Newtonian effect. The cause of this phenomenon is likely the assumption of the half-space adopted in the theory. Because the Newtonian attraction slowly decays with distance, a flat Earth model with infinitely wide spatial extension may lead to a serious overestimation of the total force. On the other hand, the free-air and the inertial effects are related with local displacement fields of the ground, which will be affected to a lesser degree by the half-space treatment compared with the Newtonian effect. Therefore, here we choose to adopt a less sophisticated grid-method for calculating the Newtonian effect, whereas the method of Zürn and Wielandt (2007) is applied to calculate the free-air and the inertial effects.
For calculating the Newtonian effect, we consider a square area of 100 km × 100 km size around Matsushiro (Figure 1). The whole area is divided by a 360 × 360 grid. The size of each cell is about 286 m, corresponding to the propagation length of the Lamb waves in 1 s (the sampling interval of the barometer data). Within one cell, atmospheric pressure is taken to be uniform. Each cell is further divided into 29 × 29 sub-cells, each having the size of about 10 m. Vertically, 10,000 layers with 5-m spacing are used. The total height of the atmosphere considered is 50 km. It is noted that the atmospheric Lamb waves are evanescent (see Fig. 5 of Arai et al. (2011)). Following equation (10) of Zürn and Wielandt (2007), we assume that density perturbations for the traveling waves depend on the height z as exp(−z/Hγ ) , where H is the scale height and γ is the specific heat ratio of the atmosphere. We assume that H = 8 × 10 3 m and γ = 1.4 . We ignore topography and curvature of the ground. For each epoch, Newtonian attraction due to density perturbations from all the sub-cells is summed up to give the total effect at station Matsushiro.
For calculating the free-air and the inertial effects, equations (15) and (16) of Zürn and Wielandt (2007) are applied. Converting the horizontal wavenumber into angular frequency ω in the two equations, the Fourier components for the free-air and the inertial effects are given by and respectively. and µ are the Lamé constants and c is the wave velocity. p ω is the Fourier component of temporal pressure changes for angular frequency ω . We assume that = µ , so that the factor ( + 2µ)/( + µ) is equal to 1.5. The vertical gravity gradient δg/δz is taken to be 3.08 × 10 −6 s −2 . Note that magnitudes of these two effects are inversely proportional to the value of rigidity µ . Because the Fourier component for the inertial effect is proportional to the angular frequency, it is necessary to suppress higher frequencies in the calculation. After some trials, we found that limiting the frequency band up to 0.0035 Hz gave a good result when compared with observed gravity residuals. The same cutoff is applied also to the Fourier component of the free-air effect.
For the observational data of gravity to be compared with theoretical calculations, we apply tidal correction using the parameters determined from more than ten years of gravity observations at Matsushiro. We do not apply the usual atmospheric correction with an admittance ( −3.3 × 10 −11 ms −2 Pa −1 in the case of Matsushiro). The residual gravity time series are then lowpass-filtered. Again, after some trials, we found that choosing a cutoff frequency of 0.0025 Hz gave a reasonable result. The difference in the optimal cutoff frequencies for pressure and gravity is likely due to the different frequency ( f ) dependence of the power spectral densities of these data at the millihertz band ( ∼ f −1 for pressure and ∼ f 0 for gravity). Figure 4b shows the calculation results of the three terms of atmospheric loading, using the atmospheric pressure data at Matsushiro shown in Figure 4a. The attraction term (blue) looks like a smoothed version of the atmospheric pressure. Because we assume that the atmospheric disturbances take the form of plane waves traveling in one direction, spatially integrating attractions from the cells over the study area is similar in effect to applying a lowpass filter to the time domain data. When the change in atmospheric pressure takes its maximum value of 174 Pa at 11:37 (UTC), the gravity decrease from the attraction term reaches 7.3 × 10 −9 ms −2 . If we simply calculate the ratio between these numbers, we obtain −4.2 × 10 −11 ms −2 Pa −1 as an apparent admittance. The reason why this is larger than usually seen comes from the fact that in the case of Lamb waves the atmospheric admittance due to Newtonian attraction is theoretically about 1.4 times as large (in the absolute sense) as that for the atmosphere in static equilibrium (Zürn and Wielandt 2007). During the initial positive pressure pulse of the Lamb waves (from 11:18 to 11:46), the negative changes in gravity due to the Newtonian effect are partly cancelled by positive changes due to the free-air (green) and the inertial (red) effects. The rigidity µ is taken to be 40 GPa. The free-air effect is temporally smoother than the inertial effect, because the latter is derived by differentiating the former twice with respect to time for each Fourier component. The sum of the three terms (magenta) is shown in Figure 4c. It reproduces well the observed gravity changes (cyan). The calculation reproduces also the latter parts of the wave train (from 11:50 to 12:30), indicating that our model also fits those later phases of atmospheric waves. The good agreement between observations and theoretical calculations proves that the theory of Zürn and Wielandt (2007) is correct. In particular, incorporating the inertial effect of atmospheric loading which is often neglected is the key feature to understand the observed gravity signals. In other words, this result may be regarded as a clear identification of the inertial effect of atmospheric loading by surface gravity observations. If we subtract the Newtonian and the free-air terms from the observed gravity, we can isolate the inertial term. Figure 4d shows the inertial effect thus extracted, to be compared with the calculated one shown in Figure 4b. This is the "signal" of the inertial effect of atmospheric loading, identified by precise gravity observations with a superconducting gravimeter. This identification was made possible because the atmospheric disturbances were (1) energetic enough and (2) spatially coherent. Although this is not the first report of the signal of this kind (Zürn and Meurers 2009), the Tonga event has provided a rare opportunity for studying such a phenomenon with sufficient signal-to-noise ratio. Fig. 4 a Atmospheric pressure changes recorded at Matsushiro. b Newtonian (blue), free-air (green) and inertial (red) terms of atmospheric loading, calculated using the barometer data shown in a. c Sum of the three terms of atmospheric loading (magenta), compared with the observed gravity changes (cyan). d The observed gravity changes with the calculated values of Newtonian and free-air effects subtracted The rigidity µ is an undetermined parameter in the theoretical calculations. After some trials, we found that the observation is well explained when we choose µ = 40 GPa, not 30 GPa or 50 GPa. In other words, a rough estimate of rigidity has also been obtained from our analysis (Wang and Tanimoto 2020). The theory we rely upon in this study was developed on the basis of an elastic half-space with homogeneous elastic properties. Therefore, it is not easy to specify the effective depth to which the present estimate of rigidity is sensitive. Considering that the wavelength of the Lamb waves used in our analysis ranges from a few tens to a few hundreds of kilometers, the estimated rigidity may represent the property in the upper part of the Earth with similar vertical scales. Indeed, rigidity of 40 GPa corresponds to the depth of 15-20 km according to PREM (Dziewonski and Anderson 1981).

Conclusion
The atmospheric disturbances induced by the eruption of the Hunga Tonga-Hunga Ha'apai volcano on January 15, 2022 were recorded as clear signals in gravity by the superconducting gravimeter at Matsushiro, Japan. We have successfully modeled the gravity changes by atmospheric loading based on the theory of Zürn and Wielandt (2007). In particular, the inertial effect due to atmospheric loading was clearly identified, almost exactly as predicted by the theory of Zürn and Wielandt (2007).
In addition, analysis of the data has enabled rough estimation of rigidity of the Earth, maybe corresponding to the physical property in the shallow region up to several tens of kilometers depth. From a gravimetric point of view, the Tonga event served as a natural experiment of measuring "how rigid the Earth is" by pressing a wide area of the ground simultaneously while making precise gravity observations on the Earth's surface.