A simple method to evaluate the air-to-ground coupling efficiency: a tool helping the assessment of seismic/infrasonic energy partitioning during an eruption

A volcanic eruption transmits both seismic and infrasound signals. The seismo-acoustic power ratio is widely used to investigate the eruption behaviors and the source dynamics. It is often the case that seismic data during an eruption are significantly contaminated or even dominated by ground shaking due to infrasound (air-to-ground signals). To evaluate the contribution of infrasound-originated power in the seismic data, we need a response function of the seismic station to infrasound. It is rare to obtain a seismo-acoustic data set containing only infrasound signals, though it is ideal for calculating the response function. This study proposes a simple way to calculate the response function using seismo-acoustic data containing infrasound and independent seismic waves. The method requires data recorded at a single station and mainly uses the cross-correlation function between the infrasound data and the Hilbert transform of the seismic data. It is tested with data recorded by a station at Kirishima volcano, Japan, of which response function has been constrained. It is shown that the method calculates a proper response function even when the seismic data contain more significant seismic power (or noise) than the air-to-ground signals. The proposed method will be useful in monitoring and understanding eruption behaviors using seismo-acoustic observations.


Introduction
A propagating acoustic wave in the atmosphere induces local ground oscillation (Ben-Menahem and Singh 1981;Sabatier et al. 1986). Therefore, it is recorded not only by infrasound sensors but also by seismic sensors. Such signals in seismometer data are called 'ground-coupled air waves' in volcanology (e.g., Johnson and Malone 2007;De Angelis et al. 2012;Fee et al. 2016;Smith et al. 2016;McKee et al. 2018;Matoza et al. 2019). However, the term 'ground-coupled air waves' was originally used to indicate the pressure oscillation generated by propagating seismic waves (Donn and Posmentier 1964). The bidirectional usage of the term has confused the discussion. Matoza and Fee 2014 distinguished air-ground coupling an ground-air coupling. In this manuscript, we call the former in seismometer records as an air-to-ground signal and the latter in acoustic records as a ground-to-air signal to clarify the coupling direction.
The generation efficiency of the air-to-ground signal is ∼ 0.1 − 10 µm/s/Pa (Ichihara 2016; Novoselov et al. 2020), which is much larger than that of the groundto-air signal ( ∼ 0.0003 Pa/(µm/s)) (e.g., Kim et al. 2004;Watada et al. 2006;Ichihara et al. 2012;Kurokawa and Ichihara 2020). The former is usually more significant in the seismo-acoustic observation during volcanic eruptions. The air-to-ground signals recorded by seismometers are used to investigate infrasound when few Open Access *Correspondence: ichihara@eri.u-tokyo.ac.jp 1 Earthquake Research Institute, University of Tokyo, Tokyo 113-0032, Japan Full list of author information is available at the end of the article infrasound sensors are available (Johnson and Malone 2007;Ichihara et al. 2012;De Angelis et al. 2012;McKee et al. 2018). On the other hand, they can produce significant power in the seismometer records and disturb the analyses of seismic signals associated with eruptions (Nakamichi et al. 2013;Ichihara (2016) and other surface phenomena, such as snow avalanches (Heck et al. 2019;Marchetti et al. 2020).
The Volcanic Acoustic-Seismic Ratio (VASR) that is the infrasonic power over the seismic power (Johnson and Aster 2005), has been used to investigate the energy partitioning between seismic waves and infrasonic waves. The energy partitioning allows us to infer the change of eruptive behaviors (Sciotto et al. 2006;Johnson and Aster 2005;Ichihara 2016;Palacios et al. 2016;Fee et al. 2020). When the same eruption mechanism is kept, the seismic and infrasonic eruption tremors tend to vary their amplitudes with the eruption parameter, such as the plume height, keeping their amplitude ratio constant (Ichihara 2016;Fee et al. 2016;Haney et al. 2018;Sciotto et al. 2019). In such situations, it is particularly important to distinguish whether the VASR or the amplitude ratio calculated from data represents the seismic-infrasonic energy partitioning or the ground response to the incident infrasound (Ichihara 2016). More essentially, we need to distinguish if the recorded waveforms mainly consist of seismic or infrasonic signals. It is achieved by the propagation velocity when seismic arrays are available (Nakamichi et al. 2013;Heck et al. 2019;Marchetti et al. 2020). Otherwise, we should refer to the response function of the seismic station to infrasound. Ichihara (2016) obtained the ground response to infrasound as a function of frequency at a station in Kirishima volcano, Japan ( Fig. 1). This case is unique in that a good infrasound source is available: Sakurajima, about 42 km away ( Fig. 1c), frequently transmits explosion infrasound (Fig. 1a). Because of the distance, seismic waves from the source only rarely reach the stations at Kirishima, or if they do, they are well separated in time from the infrasound signals. The response function, consisting of the amplitude ratio (Fig. 1d) and the phase shift of the vertical ground velocity to the pressure data ( Fig. 1e), was obtained using the initial 10 s containing a clear infrasound signal, and those from 15 events were stacked (Ichihara 2016). In general, such an ideal infrasound source is rare in the field. Kurokawa and Ichihara (2020) used the seismic and infrasonic spectral ratios of airplane signals as the reference of the air-to-ground coupling efficiency to investigate the origin of signals associated with volcanic activity at a remote island, Ioto, Japan. The airplane signals were the only known acoustic source at their station, though their frequency components and the incident angles were much higher than volcanic signals.
At seismo-acoustic stations near an active volcano, clear infrasound signals may be recorded during eruptions. However, the seismic data should always contain seismic waves as well as infrasound waves. This paper presents a convenient method to evaluate the response function using data containing both seismic and infrasonic signals. We propose to use the response function to examine the contribution of the air-to-ground signals in the VASR of data to assess the seismic-infrasonic energy partitioning more appropriately. The method requires only an infrasound sensor and a collocated seismometer without significant wind noise. If there is an additional infrasound station in a short distance (<1 km), we can also correct the wind noise effect.

Basic theory
We denote vertical ground velocity recorded by a seismometer as d v and pressure change recorded by an infrasound sensor as d p . For simplicity, we assume the incident waves are dominated by a single seismic wave ( v s ) and a single infrasonic wave ( p a ). We assume that d v consists of v s , an air-to-ground signal ( v a ) generated by p a , and wind-induced ground oscillation ( v w ), while d p consists of p a and wind noise ( p w ). The contribution of the ground-to-air signals is assumed negligible, as mentioned above. These assumptions are represented by A cross-correlation function between two time series, d 1 and d 2 , in a given time window [t, t + T ] is represented by where τ is the time delay of d 1 to d 2 . The corresponding cross correlation coefficient Hereafter, we omit t that specifies the time window. We call E(d) as the power of d, but its unit is the square of the unit of d instead of J/s.
Among v s , v a , v w , p a , and p w in Eq.
(1), we assume no pair except v a and p a has a correlation. Although v s and p a may have a correlation if their source is common, their correlation should be found with a time shift significantly different from that of v a to p a , considering the velocity difference between the seismic waves and infrasound (Ichihara et al. 2012). The correlation between v w and p w is small if the distance between the infrasound sensor and the seismometer is larger than the correlation length of wind noise. For example, the distance of 5 m is enough to suppress the wind-noise correlation at frequency > 1 Hz in ordinary wind conditions (Shields 2005;Ichihara et al. 2012). Therefore where W p represents the power ratio of the wind noise to the infrasound. Under the same wind condition, the infrasound data are usually more significantly affected by wind noise. Therefore, we neglect , and E(d p ) are calculated from the observed data. Our aim here is to estimate E(v a )/E(p a ) from the observed data, which is useful to obtain the seismo-acoustic power ratio Infrasound from an explosion, at 22:38 (hh:mm JST) on 13 January 2011, of Sakurajima volcano recorded by an infrasound sensor (a) and a broad-band seismometer (b) at the SMN station near Kirishima volcano, 42 km to the north of Sakurajima in southern Kyushu, Japan (c). The vertical dashed lines in a and b indicate 0, 10, 130 s after the arrival of the explosion signal. The amplitude (d) and phase (e) response of the ground velocity to incident infrasound at the SMN station. The red dots represent the response function at 1-7 Hz used by Ichihara (2016). All the plots are reproduced from Ichihara (2016), but the phase delay in e has been converted from the range of [−π, π ] to [−2π , 0] It is assumed that p a is propagating along the ground surface in the direction of the x-axis. The component of an angular frequency, ω , of p a is where α is the sound speed in the atmosphere. When the ground is a homogeneous elastic half space, and its seismic speeds are much larger than α , the vertical ground velocity ( v aω ) induced by p aω , is where ( , µ) are the Lamé's constants of the ground (Ben-Menahem and Singh 1981).
Equation (9) indicates that v aω has a phase shift of −π/2 with respect to p aω , and this phase shift does not depend on the frequency. The actual ground is not homogeneous or flat, and the incident infrasound wave is not necessarily along the ground surface. Nevertheless, the phase shift of −π/2 is observed in many places (Ichihara et al. 2012;Ichihara 2016). Although the sensor distance can cause a time shift of correlation, it can be corrected (McKee et al. 2018). On the other hand, the phase delay significantly deviates from the theoretical value and fluctuates at high frequency (Fig. 1e), possibly due to the ground structure. Besides, H ps depends on the frequency, especially in the low frequency (Fig. 1d). Therefore, it is difficult to estimate the response function theoretically. Here we derive an equation relating H ps and the observed data, assuming the phase shift is −π/2 . The relation is tested with data in the following sections.
First, the phase of d v is shifted by π/2 , using the Hilbert transform, h i (d v ) , as where Imag[·] takes the imaginary part. The phaseshifted time series are denoted with a subscript h, namely Second, we calculate R(d h v , d p ; τ ) , and search its maximum around τ = 0 . We denote τ that gives the maximum as τ max . It should be consistent with the propagation time of infrasound between the infrasound sensor and the seismometer. When any signals in d v , other than v a , correlate with p a , CC(d h v , d p ; τ ) tends to have a maximum significantly deviated from the expected time delay (Ichihara et al. 2012). Such data sets are not adequate to use in the current method.
Third, we approximate The last approximation holds if the wind noise is negligible. We use Eqs. (12) and (5)  It is noted that this equation is applicable on the condition that an infrasound signal p a exists. Equation (14) indicates that we should refer to H 2 ps for interpreting the observed power ratio, E(d v )/E(d p ) , the reciprocal of which is associated to VASR. The seismic power E(v s ) increases the power ratio above H 2 ps . Equation (14) is particularly useful when the air-to-ground signal is comparable with the seismic power. We can also identify the condition of E(v s ) ≪ H 2 ps , in which the observed E(d v )/E(d p ) cannot provide the meaningful seismicacoustic energy partitioning.
When the wind noise is significant, we need another infrasound station close by to evaluate the relation between E(d p ) and E(p a ) . The calculation with wind noise is given in Additional File 1. With wind noise, Eq. (14) becomes is significantly contaminated by the wind noise power, and thus E(d v )/E(d p ) (or VASR) is not useful to investigate the seismic-acoustic energy partitioning, either.

Method to estimate H ps from data
We use the data presented in Fig. 1a and b as d p and d v , respectively. They are infrasound from Sakurajima recorded by an infrasound sensor (Hakusan, SI102, 0.1-1000 Hz) and a seismometer (Nanometrics, Trillium 120P) at Kirishima, respectively. The sampling rate is 100 Hz for both sets of data. Ichihara (2016) has shown that the part of d v containing the air-to-ground signal is well reproduced by d p with the response function at 1-7 Hz shown by the red dashed line in Fig. 1d and 1e. As the infrasound signal is larger than the background wind noise (Fig. 1a), we assume W p = 0 in the following analysis. On the other hand, there is noticeable oscillation before the arrival of infrasound in the seismograph (Fig. 1b), indicating that the background seismic signal is comparable with the air-to-ground signal. The background signal might be the volcanic tremor associated with the Shinmoe-dake activity at Kirishima (Ichihara and Matsumoto 2017), but we do not discuss its origin in detail in this study. Here we treat all signals other than v a in d v as v s .
We use the four frequency bands: 1-3.5, 3.5-7, 7-12, and 12-18 Hz. For the selection of the frequency bands, we refer to the feature of the known response function shown in Fig. 1d and e. The phase shift of d v to d p is constant around −π/2 at 1-7 Hz, while the amplitude ratio H ps is relatively stable in 7-18 Hz. In general cases without known response functions, we may arbitrarily select the frequency bands.
To check that the data include infrasound signals, we examine the cross-correlation coefficient, R(d v , d p ; τ ) . We apply a zero-phase-shift bandpass filter to d p and d v and calculate R(d v , d p ; τ ) , in the individual bands. Figure 2 shows the results, in which the time origin is taken at the onset of the infrasound signal. In each of Fig. 2a-d, the lowest frame shows R(d v , d p ; τ ) , calculated in a 2.56-s time window sliding every 0.5 s. We can clearly see the typical correlation pattern of air-to-ground signals, that is a positive peak in τ > 0 and a negative peak in τ < 0 with a node at τ = 0 (Ichihara et al. 2012), continuously after the arrival of the infrasound. It confirms that the data include p a and v a in the period in all the frequency bands.
To obtain H ps applying Eq. (13), we use the bandpass-filtered data from 10 s to 130 s after the arrival of the infrasound as d p and d v . The initial 10-s window is excluded, because it includes the strong pulse used Fig. 2 Cross-correlation analyses of the data in Fig. 1 at a 1-3.5 Hz, b 3.5-7 Hz, c 7-12 Hz, and d 12-18 Hz. In each panel, the top and middle frames show the band-passed infrasonic data ( d p ) and seismic data ( d v ), respectively. The bottom frame shows the cross-correlation coefficient between the seismic data, d v , and infrasonic data, d p , with the vertical axis representing the time delay, τ , of d v to d p . The horizontal dashed line indicates τ = 0 to obtain the response function in Fig. 1. First, we shift the phase of d v forward by π/2 to obtain d h v . Second, we calculate R(d h v , d p ; τ ) . Third, we take the maximum of R(d h v , d p ; τ ) , denoted as R max , in τ = ±0.01 . This span, corresponding to one data point, is allowed, considering the time delay between the two sensors separated by ∼ 5 m and possible phase delay due to noise. Fourth, we calculate E(d p ) and E(d h v ) by the mean square of the data. Finally, we calculate H ps , using Eq. (13). Figure 3a compares the value of H ps obtained by the current method and by Ichihara (2016). They are consistent. It is noted that the previous value was estimated from the beginning 10 s of 15 events, while the current value is from a single event in the subsequent 120 s. Figure 3b shows the corresponding R(d h v , d p ; τ ) for each band. They have the maximum values within τ ± 0.01 s, confirming that the phase shift of −π/2 has been recovered successfully.

Results and discussion
The data used in this analysis are ideal to evaluate H ps , because d p and d v mainly include p a and v a , respectively. We do not always have such data sets. Here we examine the effect of v s with the assumption that v s is independent Fig. 3 a Amplitude response of the ground velocity to infrasound, H ps , obtained by this study is shown for 1-3.5 Hz (red line), 3.5-7 Hz (green line), 7-12 Hz (blue dashes), and 12-18 Hz (purple dashes). The results are compared with the response function obtained by Ichihara (2016) (black line; Fig. 1d). b Corresponding cross-correlation coefficient between the phase-shifted seismic data, d h v , and the infrasonic data, d p . The colors indicate the frequency bands as in a Fig. 4 Effect of background seismic signals, v s , in the seismic data, d v , on the estimation of H ps . The colors indicate the frequency bands: 1-3.5 Hz (red), 3.5-7 Hz (green), 7-12 Hz (blue), and 12-18 Hz (purple). One-hundred random functions are bandpass filtered and used as v s . The horizontal axes in a, c, and d are Ŵ , the square-root of power ratio of v s to the original d v that is approximated as v a . It is noted that the actual power ratio E(v s )/E(v a ) is larger than Ŵ 2 = E(v s )/E(d v ) as d v itself contains a significant power of background signals. The lines and error bars show the means and the standard deviations of H ps (a), the maximum correlation coefficient, R max (d), and the time delay, τ max , that gives R max (c) for the hundred random functions. The closed circles on the vertical axies are the values without v s . The relation between the individual R max and H ps for all the cases are shown with circles in b of p a . We artificially add v s made in the following way. We generate a random function that is centered at zero and has the same length as d v . We apply the same band-pass filter to this function and call it as d bg . We normalize d bg to have the same power as d v , by multiplying it by . Then, we amplify it by an arbitrary factor Ŵ to make a hypothetical v s . Namely should be larger than Ŵ , because d v contains some background signal other than v a . We may see it by comparing the amplitudes before and after the infrasound arrival in d v and d p (Fig. 2). It is also supported by Fig. 3b that R max < 1.
Including v s into d v , we follow the same procedure as the previous paragraph to calculate H ps . For each frequency band and each Ŵ , we tested with 100 random functions. The mean and the standard deviation of H ps , τ max , and R max are presented as functions of Ŵ in Fig. 4a, c, and d, respectively. We also show the relation of the individual H ps and R max in Fig. 4b. As the background signal power, E(v s ) , increases, the cross-correlation coefficient decreases (Fig. 4d). On the other hand, H ps stays around the expected value ( Fig. 4a and b). The average values are almost the same as those without artificial noise, H 0 ps (the closed circles on the vertical axis in Fig. 4a). A similar result is obtained from a test using actual volcanic tremor data instead of the random functions (Additional File 2). These results support that we can evaluate H ps properly using a single event signal, even when the background signals in the seismometer data (seismic waves or other noise) are comparable with the air-to-ground signals. In cases with more significant contaminations ( R max < 0.3 ), we may estimate H ps reasonably when we have repetitive Fig. 5 Result of an analysis using the seismo-acoustic eruption tremor data in a 1-3.5 Hz, b 3.5-7 Hz, c 7-12 Hz, and d 12-18 Hz. In each panel, the top two frames show the acoustic and seismic data in the frequency band ( d p and d v , respectively), and the third shows the maximum correlation between d p and d h v . In the bottom frame, the open circles show the calculated H ps using this data in a 120-s time window moving every 60 s, the colored dashed line indicates the reference value from Fig. 3, and the crosses represent E(d v )/E(d p ) in the 120-s windows or persistent infrasound signals and can average in a long time window or many windows.
Finally, we analyze data of seismo-acoustic eruption tremors recorded by the same station during the 2011 eruption of Shinmoe-dake at Kirishima volcano. Three subplinian events occurred, one in the afternoon of 26 January 2011 and two on the next day. Here we use the data from 0:00 to 4:50 on 27 January 2011, associated with the second event (Ichihara 2016). We apply the bandpass filter and calculate H ps using Eq. (13) in a 120-s time window moving every 60 s. Figure 5 shows the results in the frequency bands (a) 1-3.5 Hz, (b) 3.5-7 Hz, (c) 7-12 Hz, and (d) 12-18 Hz. In each of Fig. 5a-d, the top two frames show the bandpass-filtered acoustic and seismic waveforms, and the third shows R max . The calculated H ps is presented with open circles in the bottom frame. They distribute close to the colored dashed line, which indicates the expected value, H 0 ps (the circle of the same color in Fig. 4a). It means that the current method reasonably estimates H ps with the eruption tremor data. The deviation and scattering become larger when R max is smaller, which is expected from the result in Fig. 4. The crosses in the same bottom frame show the root-mean-square amplitude ratio, E(d v )/E(d p ) . They are always larger than H ps , particularly after 3:00 in the climatic phase of the second subplinian event, above 3.5 Hz (Fig. 5b-d). Namely, the seismic eruption tremors ( v s ) dominate the air-to-ground signals ( v a ) generated by the acoustic eruption tremors.

Conclusions
We developed Eq. (13) that relates the ground response to infrasound ( H ps ) and the observed seismo-acoustic data ( d v and d p ). We have confirmed that it provides a proper value of H ps using a single infrasonic signal when the seismic data contain background seismic signals (noise) as large as the air-to-ground signal. Even if the background noise is more significant, the proposed method allows us to estimate H ps reasonably by averaging the values obtained from a large number of signals or time windows. Constraining H ps is essential for interpreting the observed seismic-acoustic power ratio to the energy partitioning between seismic waves and infrasound at the source.