Takanori Matsuzawa ¹ , Takumi Matsumoto ¹ and Toshikazu Tanada ¹ ¹ National Research Institute for Earth Science and Disaster Resilience Address: 3-1 Tennodai, Tsukuba, Ibaraki 305-0006, JAPAN (Received: Sep.1, 2016 Accepted: Oct.21, 2016) Short Running Title: Seismic inversion of Kuchinoerabujima eruption
Abstract
We observed long-period seismic waves excited by eruptions in 2014 and 2015 at Kuchinoerabujima Volcano, Japan. The propagations of seismic waves were clearly recognized in F-net seismograms in southwest Japan. The excitation force of three eruptions was analyzed using F-net broadband seismometers. The analysis of source time functions of three-component single force shows that vertical force is dominant in the eruptions at 12:24 on Aug. 3, 2014, and at 12:09 on May 29, 2015 (local time). In addition to the vertical force, lateral force in the north-south direction was also dominant in the eruption at 9:59 on May 29, 2015. This suggests that the mechanism of this eruption is different from the other two eruptions. The lateral force may be related to the momentum exchange by the emission of deposits inclined to the north- south direction. Keywords: Kuchinoerabujima Volcano, Eruption, Single force, Seismic Inversion
INTRODUCTION
1.
Long-period seismic wave excitations are sometimes observed in volcanic eruptions. To reveal the mechanism of eruptions, the excitation force is estimated by a seismic inversion technique. For example, lateral force is estimated in the 1980 St. Helens eruption, and attributed to a landslide followed by a lateral blast (Kanamori and Given, 1982). In the 1989 Izu-Oshima volcanic eruptions, vertical single force and vector dipoles are estimated, and interpreted as the force caused by the sliding down of the lava, and the volumetric change in cracks or extended magma reservoirs, respectively (Takeo et al., 1990). In the 2010 eruption with caldera formation on Miyake Island, the very-long-period seismic signal is explained by a piston model (Kumagai et al., 2001). In the eruptions of Asama Volcano in 2004, downward and subsequent upward vertical force are interpreted as a result of removal of capping and viscous drag force by ascending magma, respectively (Ohminato et al., 2006). As in these studies, analysis of the seismic excitation source is useful for discussing the eruption process. Kuchinoerabujima Volcano is located in the south of Kyushu Island, Japan (Fig. 1). This volcano consists of several edifices. Shintake and Furutake Volcanos in Kuchinoerabujima have repeatedly erupted in the last 10,000 years. Lithic tephra indicates that these are Vulcanian-type and phreatomagmatic eruptions (Geshi and Kobayashi, 2006). Historical records since 1841 show that eruptions occurred at or around Shintake crater, and repeated at intervals of tens of years (Geshi and Kobayashi, 2006; Japan Meteorological Agency (JMA), 2013). In very recent years, an explosive eruption from the Shintake crater was observed at 12:24 on Aug. 3, 2014 (local time; i.e., UTC+9) (JMA, 2014). In the following year, an explosive eruption at Shintake crater occurred at 9:59 on May 29, 2015 (JMA, 2015a). Due to this eruption, the residents of Kuchinoerabu Island were evacuated from the island. The plume height of this eruption reached more than 9000 m above the rim of the crater (JMA, 2015a). At about 12:09 on the same day, another eruption was observed (Japan Coast Guard, 2015). Eruptions were also observed at 12:17 and 16:31 on Jun. 18, 2015, and at 9:43 on Jun. 19, 2015 (JMA, 2015b) at Kuchinoerabujima Volcano. Hereinafter, we call the eruptions at 12:24 on Aug. 3, 2014, at 9:59 on May 29, 2015, and at 12:09 on the same day as Eruption A, B, and C, respectively. We analyze these three eruptions in this paper.
2.
SEISMIC DATA
The National Research Institute for Earth Science and Disaster Resilience (NIED) operates a nation-wide broadband seismic network called F-net (Okada et al., 2004). At the F-net stations, Streckeisen STS-1, 2, or 2.5 broadband sensors are installed in tunnels. Figure 2 shows the seismic data of F-net at the time of eruptions. In Eruption A, B, and C, propagations of long-period (more than 10 s) seismic signals are recognized even at a distance of hundreds of kilometers from Kuchinoerabujima Volcano (Fig. 2(a), (b), and (c)). We note that all seismic data in the analysis and figures are prepared after correction of the instrumental response. KYKF station is located on Yakushima Island next to Kuchinoerabu Island (Fig. 1). The distance from the crater to KYKF is about 20 km. Seismic signals of eruptions were clearly observed in the broadband seismograms. Figure 3 shows the seismograms and the particle motions in a horizontal plane at KYKF. Particle motions of Eruption A and C are almost in the radial (east-to-west) direction. The particle motion of Eruption B is dominant in the transverse (north-to-south) direction. This suggests that the direction of the applied force in Eruption B is different from the other eruptions. Seismograms of Eruption A are slightly noisy in the whole of Japan even before the eruption (Fig. 2(a)). In Fig. 3, the linearity of the particle motion of Eruption A is less clear than those of the other eruptions. An Mw 6.9 earthquake occurred at 00:22 (UTC) in the Micronesia region before Eruption A. The seismograms of Eruption A, as in Fig. 2 and 3, seem to be affected by the coda part of surface waves excited Figure 2 UD-component seismograms of F-net at the time of eruptions. Seismograms are band-pass filtered between 10 s and 20 s. The vertical position of each trace is based on the distance from Kuchinoerabujima Volcano. The scale of the amplitude of waveforms is the same in all panels. Each panel shows the eruption at (a) 12:24 on Aug. 3, 2014 (Eruption A), (b) 9:59 on May 29, 2015 (Eruption B), (c) 12:09 on May 29, 2015 (Eruption C), (d) 12:17 on Jun. 18, 2015, (e) 16:31, Jun. 18, 2015, and (f) 9:43 on Jun. 19, 2015. The propagations of seismic waves were recognized in (a), (b), and (c), as indicated by black tilted arrows. Aso long-period events, which were detected by a matched filter technique (Matsuzawa et al., 2015), are indicated by gray tilted arrows. In the left panels, the locations of F-net stations are shown by squares. The locations of Kuchinoerabujima Volcano and Aso Volcano are indicated by a black and a gray triangle in the same map, respectively. Propagations of small seismic waves are also found in Fig. 2 (b), (c), and (e) at the time without eruptions. These did not originate from the Kuchinoerabu region, but from the Aso volcanic region. They are interpreted as long-period events in Aso Volcano (Kaneshima et al., 1996). As indicated by the gray arrows in Fig. 2, long-period events are detected at 10:07, 12:08, and 12:13 on May 29, 2015, and 16:30 on Jun. 18, 2015 by a matched filter technique using F-net data (Matsuzawa et al., 2015). These detected times correspond to the start of propagations of seismic waves from the Aso region. After Eruption C, three eruptions have been reported as described in the previous section. Figure 2 (d), (e), and (f) show the waveforms around the time of these reported eruptions. However, clear propagations from Kuchinoerabujima Volcano are not recognized in these three eruptions at F-net stations. It is difficult to analyze such very small seismic sources of these eruptions using only F-net data. Therefore, we analyzed the eruptions with clear seismic signals (i.e., Eruption A, B, and C) in this paper.
3.  INVERSION METHOD
We estimate the source time functions of excitation force using a waveform inversion method to reveal the mechanisms in eruptions and adopt a discrete wavenumber integral method (Takeo, 1985; Takeo et al., 1990) to calculate synthetic waveforms. The same laterally layered seismic structure is assumed as that in F-net (Kubo et al., 2002). In this analysis, source time functions of three-component single force are estimated simultaneously by least squares inversion with a smoothing constraint. Each component of source time functions is expressed by a first order B-spline function that is spanned by half-overlapped basis functions (Fig. 4). duration is 100 s. The number of basis functions,  N , is 100 in our result except for the case in Fig. 5(c). In a discrete form, the synthetic waveform of the  j -th component at the  i -th station by the single force applied in the  k -th direction is expressed as
g iz ( t ) is the synthetic waveform, when the source time function is a unit basis function. s x ( t ) is a source time function and indicates the amount of applied force. K is the number of force directions (e.g., K = 3 for three-component force). N is the number of points of the source time function. This is easily extended to the case of the moment tensor source, replacing   s t ( t ) t ) with the moment rate function of the k -th basis moment tensor and the synthetic waveform excited by a unit basis source time function of the corresponding moment tensor, respectively. In addition, waveforms by a mixed source of single force and moment tensor source are expressed by a simple summation of two equations in these two cases. In this inversion analysis, the sampling interval of source time function is the same as the decimated interval of the observed waveforms ( Δt =1 s). Therefore, t m  and  τ 1  are given as
(2). M is the number of data points in each decimated observed waveform, and 180 in this analysis except for the case in Fig. 5(c). We assume that g ijk (t) is zero for t < 0 or t > L t （L = 80 in this study). Waveform data and source time functions are written in a vector form as and
respectively. I  and J  are the number of used stations and number of components of each sensor, respectively. Using these vectors, the observation equation is written as
Here, G is an IJM × KN size matrix. The matrix coefficients of G are given by (1). In this analysis, a smoothing constraint is introduced, as the second order difference is zero. This is expressed in an equation as
For convenience,   s k ( τ 0 )  and   s k ( τ N +1 )  are assumed to be zero in the above equation. The set of Equation (6) for all   k   and l   is written in matrix form as Fs = 0
Combining this constraint, the observation equation is modified to
In this equation, &  is a weighting parameter of the smoothing constraint (i.e., hyper parameter). The observation equation with the constraint (8) is solved by a least squares method. As an optimized hyper parameter, we choose  α , which minimizes Akaike's Bayesian Information Criterion (ABIC) (e.g., Yabuki and Matsu’ura, 1992). The source of the synthetic waveform is assumed as a point source at a depth of 0.5 km in the Shintake crater. Observed and synthetic waveforms are band-pass filtered between 10 and 100 s.
4.
RESULTS
Figures 5, 6, and 7 show the results of estimated source time functions of Eruption A, B, and C, respectively. Both the north-south and vertical components are dominant in Eruption B, while only the vertical component is dominant in Eruption A and C. This means that seismic waves in Eruption A and C are mainly excited by the vertical force. Seismic waves in Eruption B are also excited by the lateral force in the north-south direction. In Eruption B, the peak amplitudes of the lateral single force are similar to that of the vertical force. In addition, the rise of the vertical force precedes the lateral force. The results of the force directions and the delayed lateral force of Eruption B suggest that the excitation force of Eruption B is different from the other two eruptions and is complex. component single force. (b) Comparison of synthetic (thick dashed lines) and observed waveforms (thin lines). In Eruption A, force in the north-south direction is dominant in the longer period (more than 30 s), especially after 30 s, while vertical force is dominant in the shorter period (10–20 s) within the first 30 s. The propagating surface waves by Eruption A are dominant at the period from 10 to 20 s. Surface waves by the far-field earthquake are dominant above 20 s at the time of Eruption A. Therefore, the longer-period lateral force after 30 s apparently arises from the effect of the far-field earthquake. Figures 5(c) and (d) show the result when the source duration is limited to the first 30 s. The characteristic waveform at KYKF is well explained without the large lateral single force. In addition, the peak value of the estimated single force in Fig. 5(c) is larger than that in Fig. 5(a), and also larger than those of Eruption B and C. The peak value in Fig. 5(a) seems to be underestimated by the disturbance by the far-field earthquake. This is supported by the peak amplitude of the seismic wave at KYKF that is larger than those of Eruption B and C (Fig. 3).
5. DISCUSSION
Our results show that both the lateral and vertical single forces are dominant in Eruption B. However, it is also possible that other processes (e.g., earthquake, or crack opening) may occur in the eruption. As such mechanisms are represented by moment tensor, therefore, we examined a seismic source with vertical single force and six-component moment tensor. Figure 8 shows the result in this case. As the variance reduction (VR) is 80.2%, waveforms are relatively well explained. However, this is smaller than the VR of 84.5% in the case of the three-component single force, although the number of model parameters in Fig. 8 is more than double. This suggests that the contribution from the lateral single force is suitable for explaining the seismic waveforms. In terms of the geological observation, it is suggested that the eruptions in 2014 and 2015 at Kuchinoerabujima Volcano are directed blasts (Kobayashi, 2015; Kobayashi and Yamaguchi, 2015). Kobayashi and Yamaguchi (2015) suggested that Eruption B was a lateral blast. Force in the north-south direction in our result supports the occurrence of lateral blast in Eruption B, as the lateral mass movement is suggested from the lateral single force. In addition, a blast in the lateral direction is reported in the eruptions at the other volcanos (e.g., Kieffer, 1981; Kanamori and Given, 1982). In Eruption B, the rise of the source time function in the vertical direction precedes that in the north-south direction by several seconds. The delayed lateral force is also supported by the camera observed by JMA, which shows that Eruption B is in the vertical direction at first, and then in the inclined direction after several seconds (e.g., Jiji Press, 2015). The excitation force of Eruption C was dominant in the vertical direction. This means that the excitation force changed in each eruption even at the interval of two hours. Kobayashi and Yamaguchi (2015) suggest that blocking material on the crater caused a lateral blast in Eruption B. Perhaps, the material has been removed from the crater by Eruption B. Then, the direction of Eruption C may be vertical again as in Eruption A. In addition, the source time function of Eruption B is more complex and has a slightly higher dominant frequency than those in the other eruptions. Such features may be related to the coexistence of vertical and lateral blasts that are caused by the removal of blocking material at the crater. In Eruption B, small long-period force is found both in the north-south and vertical direction from 60 s to 80 s (Fig. 6(a)). This may also be attributed to the eruption. Kanamori and Given (1982) reported that a large-scale landslide excited long-period seismic waves after a lateral blast in the 1980 St. Helens eruption. However, this mechanism does not seem to be applicable to the excitation force of Eruption B, as no large landslides were reported in this eruption. Deceleration or landing of laterally emitted deposits may cause such momentum exchange from 60 s to 80 s, as found in the camera observation by JMA (e.g., Jiji Press, 2015). Further study, including the location of sources, might reveal whether the source at each stage is around the surface or in the deeper part. Eruption B was larger than Eruption C, as confirmed by the satellite images (JMA, 2015c), while the amplitudes of seismic excitation are similar. As single force excitation is interpreted as the momentum exchange by mass movement, this suggests that a large fraction of mass movement in Eruption C occurred below the surface. Eruption B or dike opening, which preceded Eruption C, might cause a low-pressure region in the vent or magma chamber. Then, fluid movement or degassing from magma might flow into that region from the deeper part. Further study on the depth and the detailed excitation mechanism would be helpful to discuss this point. Our obtained source time functions show oscillating behaviors. In Eruption C, for example, maximum downward force is found after small downward force and moderate upward force. This may be caused by the uncertainty of the assumed Green's function and microseismic noise, as the dominant period of source time function is relatively short (typically 10–20 s). Data observed at close stations are likely to be less affected by such effects due to the short path length. In a future study, an inversion analysis with close stations might be suitable for discussing the details of the source time function (e.g., Ohminato et al., 2006). In that study, the topography of the island should be considered, while a laterally layered structure is assumed in this study that uses distant stations. Another cause of such oscillation is the effect of the applied band-pass filter. Nakano et al. (2008) pointed out that the obtained source time function is band-pass-filtered, when the filter is applied only to the observed seismic waveforms. In this analysis, this effect can be avoided, as the same band-pass filter is applied to both the observed and synthetic waveforms (Yagi, 2014). However, Yagi (2014) showed that the applied band-pass filter indirectly affects the solution in the case of damped least squares inversion. As our smoothing constraint is different from that in Yagi (2014), we briefly discuss the effect of band-pass filter, replacing the term of constraint in (12), (13), and (14) in Yagi (2014). For simplicity, we assume continuous functions instead of a discrete form as in Yagi (2014). In the time domain, the inversion problem is attributed to minimizing
Here, * and B(t) mean convolution in the time domain and a response function of the band-pass filter, respectively. The smoothing constraint of the second order difference (Equation (6)) is expressed by the second order derivatives of the source time function. From Parseval's theorem, this is equal to minimizing
in the frequency domain. Here,  ·(ω)  denotes the Fourier transformed function. Therefore, the source time function in the frequency domain is given by
means a complex conjugate. As B*(ω)B(ω)<<1 out of the pass band, In the above equation, becomes very large in the high frequency band. Therefore, the high frequency component in the solution was reduced indirectly by the effect of the band-pass filter. The strength of this effect is controlled by the hyper parameter   α . Such effect may also cause the small apparent oscillation in our result due to the lack of the high frequency component. Observations close to the volcano are sometimes unavailable due to natural hazards or other factors (e.g., human error). This risk increases when an eruption occurs. Some observations at stations close to Kuchinoerabujima Volcano were disrupted by the eruptions (e.g., JMA, 2016). However, continuous monitoring and analysis with F-net data are available even in such cases, as F-net sensors are widely distributed in Japan and installed in places with a low-noise level that enables us to analyze small signals. Analysis with distant seismic stations is important as an alternate monitoring method of volcanic activity, when immediate analysis or monitoring is required, especially, in the case of lack or failure of data from close stations.
CONCLUSIONS
6.
When an eruption on Aug. 3, 2014 (Eruption A) and two eruptions on May 29, 2015 (Eruption B and C) occurred in Kuchinoerabujima Volcano, the propagations of seismic waves were clearly recognized in southwestern Japan in the F-net broadband data (Fig. 2). However, the propagations of seismic waves were not clear at the time of the three eruptions on Jun. 18 and 19, 2015. This means that the excitation force of these three eruptions was smaller than those of Eruption A, B, and C.We analyzed the seismic waves in Eruption A, B, and C, and estimated the excitation force by the least squares inversion method with a smoothing constraint using F-net broadband data. The excitation force of Eruption A and C was mainly in the vertical direction. In Eruption B, however, lateral force was also applied just after the vertical force. This is consistent with the emission of deposits observed by a JMA camera (Jiji Press, 2015), and the suggestion from a geological study (Kobayashi and Yamaguchi, 2015)
Acknowledgements
We used the Generic Mapping Tools (Wessel and Smith, 1998) to prepare figures. This work is partially supported by MEXT/JSPS KAKENHI, Grant-in-Aid for Special Purposes. We appreciate two anonymous reviewers and a sub editor for their helpful comments.
References:
Geshi, N. and T. Kobayashi, 2006. Volcanic activities of Kuchinoerabujima volcano within the last 30,000 years, Bulletin of the Volcanological Society of Japan, 51, 1–20. Japan Coast Guard, 2015. Press release of Coordinating Committee for Prediction of Volcanic Eruption No. 132 on Jun. 15, 2015, Japan Meteorological Agency. Japan Meteorological Agency, 2013. National catalogue of the active volcanoes in Japan, The fourth edition, Japan Meteorological Business Support Center, 1500pp. Japan Meteorological Agency, 2014. Press release of Coordinating Committee for Prediction of Volcanic Eruption No. 130 on Oct. 23, 2014, Japan Meteorological Agency. Japan Meteorological Agency, 2015a. Press release of Coordinating Committee for Prediction of Volcanic Eruption No. 132 on Jun. 15, 2015, Japan Meteorological Agency. Japan Meteorological Agency, 2015b. Press release of Coordinating Committee for Prediction of Volcanic Eruption No. 133 on Oct. 21, 2015, Japan Meteorological Agency. Japan Meteorological Agency, 2015c. http://www.jma-net.go.jp/sat/data/web89/himawari8_sample_data.html, and http://www.jma- net.go.jp/sat/data/web89/parts89/himawari8_sample_data/original/20150529_TRC_kuchinoerabu.wm v (accessed on Nov. 16, 2016). Japan Meteorological Agency, 2016. Press release of Coordinating Committee for Prediction of Volcanic Eruption No. 134 on Feb. 17, 2016, Japan Meteorological Agency. Jiji Press, 2015. https://www.youtube.com/watch?v=WtMC2GYKkY0 (accessed on Aug. 23, 2016). Kanamori, H. and J. W. Given, 1982. Analysis of long-period seismic waves excited by the May 18, 1980, eruption of Mount St. Helens –A terrestrial monopole?, J. Geophys. Res., 87, 5422–5432. Kaneshima, S., H. Kawakatsu, H. Matsubayashi, Y. Sudo, T. Tsutsui, T. Ohminato, H. Ito, K. Uhira, H. Yamasato, J. Oikawa, M. Takeo, and T. Iidaka, 1996. Mechanism of phreatic eruptions at Aso volcano inferred from near-field broadband seismic observations, Science, 273, 642–645. Kieffer, S. W., 1981. Blast dynamics at Mount St Helens on 18 May 1980, Nature, 291, 568–570. Kobayashi, T., 2015. Reconsideration of volcanic eruptions associated with blast, Chikyu Monthly, 37, No. 6, 220–227. Kobayashi, T. and K. Yamaguchi, 2015. "Pyroclastic flows" generated in the 2014 and 2015 eruptions of Shindake, Kuchinoerabujima Island, Abstract of the Volcanological Society of Japan, 2015 Fall Meeting, p.13. Kubo, A., E. Fukuyama, H. Kawai, and K. Nonomura, 2002. NIED seismic moment tensor catalogue for regional earthquakes around Japan: quality test and application, Tectonophys., 356, 23–48. Kumagai, H., T. Ohminato, M. Nakano, M. Ooi, A. Kubo, H. Inoue, and J. Oikawa, 2001. Very-long-period seismic signals and caldera formation at Miyake Island, Japan, Science, 293, 687–690. Matsuzawa, T., T. Matsumoto, and T. Tanada, 2015. Detection of long period events of the Aso volcano applying a matched filter technique to F-net broadband seismic data, Japan Geoscience Union Meeting 2015. Nakano, M., H. Kumagai, and H. Inoue, 2008. Waveform inversion in the frequency domain for the simultaneous determination of earthquake source mechanism and moment function, Geophys. J. Int., 173, 1000–1011, doi:10.1111/j.1365-246X.2008.03783.x. Ohminato, T., M. Takeo, H. Kumagai, T. Yamashita, J. Oikawa, E. Koyama, H. Tsuji, and T. Urabe, 2006. Vulcanian eruptions with dominant single force components observed during the 2004 Asama volcanic activity in Japan, Earth Planet Space, 58, 583–593. Okada, Y., K. Kasahara, S. Hori, K. Obara, S. Sekiguchi, H. Fujiwara, and A. Yamamoto, 2004. Recent progress of seismic observation networks in Japan –Hi-net, F-net, K-NET, and KiK-net–, Earth Planet Space, 56, 15–28. Takeo, M., 1985. Near-field synthetic seismograms taking into account the effect of anelasticity –The effect of anelastic attenuation on seismograms caused by a sedimentary layer–, Papers in Meteorology and Geophysics, 36, No.4, 245-257. Takeo, M., H. Yamasato, I. Furuya, and M. Seino, 1990. Analysis of long-period seismic waves excited by the November 1987 eruption of Izu-Oshima volcano, J. Geophys. Res., 95, 19377–19393. Wessel, P. and W. H. F. Smith, 1998. New, improved version of Generic Mapping Tools released, EOS Trans. AGU, 79, No. 47, 579. Yabuki, T. and M. Matsu'ura, 1992. Geodetic data inversion using a Bayesian information criterion for spatial distribution of fault slip, Geophys. J. Int., 109, No. 2, 363–375. Yagi, Y, 2014. Filtering effect in waveform data inversion for seismic source process, Zisin (Journal of Seismological Society of Japan), 2nd.ser., 66, 147–149, doi:10.4294/zisin.66.147.