Short-term spatial change in a volcanic tremor source during the 2011 Kirishima eruption

Abstract

Volcanic tremors are indicators of magmatic behavior, which is strongly related to volcanic eruptions and activity. Detection of spatial and temporal variations in the source location is important for understanding the mechanism of volcanic eruptions. However, short-term temporal variations within a tremor event have not always been detected by seismic array observations around volcanoes. Here, we show that volcanic tremor sources were activated at both the top (i.e., the crater) and the lower end of the conduit, by analyzing seismograms from a dense seismic array 3 km from the Shinmoedake crater, Kirishima volcano, Japan. We observed changes in the seismic ray direction during a volcanic tremor sequence, and inferred two major sources of the tremor from the slowness vectors of the approaching waves. One was located in a shallow region beneath the Shinmoedake crater. The other was found in a direction N30°W from the array, pointing to a location above a pressure source. The fine spatial and temporal characteristics of volcanic tremors suggest an interaction between deep and shallow conduits.

1. Introduction

Volcanic tremors are commonly observed around active volcanoes, and occur in the conduit system between the magma chamber and crater. The duration of most volcanic tremors is much longer than an ordinary tectonic earthquake, and the signal frequency comprises fundamental, and harmonic, oscillations (McNutt, 2005). Based on seismic observations, numerous studies have modeled the mechanical source of volcanic tremors (e.g., Julian, 1994; Chouet et al., 1997; Hellweg, 2000; Jellinek and Bercovici, 2011). Seismic array observations that use slowness analysis of seismic waves are a powerful tool for locating the origin of tremors. In particular, array analyses using the multiple signal classification (MUSIC) technique (Schmidt, 1986; Kikuma, 1999) have revealed the properties of mag-matic activity in volcanic regions (Goldstein and Chouet, 1994).

Shinmoedake crater, which is in the Kirishima volcano in Kyushu, Japan (Fig. 1), began a phreatic eruption on 19 January, 2011, with the eruption type changing to Vulcanian via sub-Plinian. According to geodetic analysis based on GPS data (Geographical Survey Institute of Japan, 2011), the estimated pressure source was located at a depth of 7.7 km, 6 km northeast of the crater (Fig. 1). The results also showed that the uplift of the volcano (ongoing since the previous year) had stopped, and subsidence of the volcano started after an eruption on 26 January, 2011. The largest explosion occurred on 1 February, with strong sonic waves affecting the residential area around the volcano. A volcanic plume over 3000 m high also erupted from the crater (Japan Meteorological Agency, 2011). Volcanic tremors occurred frequently during early February. The tremors had characteristics (i.e., long-duration harmonic oscillations) similar to those found at some other volcanoes. Here, we investigate the characteristics of volcanic tremors from an analysis of seismic array data.

Fig. 1.

Topographic map of seismic array location. Detailed configuration of the array is displayed lower left. Circles and the square indicate sensors and the broadband seismometer, respectively. The red circle and the pink dotted area show the pressure source from GPS data (Geographical Survey Institute of Japan, 2011), and the region with a velocity perturbation of −4% to −6% at a 3 km depth (Nishi and Kagiyama, 2002), respectively. The array response in the slowness domain at 2 Hz is shown at the lower right.

2. Slowness Analysis of Volcanic Tremors

After the start of the eruption, a seismic array was urgently deployed on 29 January, 2011, along a road 3 km from the Shinmoedake crater in the direction N115°W. Twenty-five 3-component seismometers with a natural frequency of 2 Hz were installed at a sensor interval of 20–40 m. The location and shape of the array, and array response functions at 2 Hz, are shown in Fig. 1. Data was sampled every 1 ms by a recording system controlled by a GPS-calibrated clock. In this study, we use seismograms recorded by 23 of the 25 sensors with good signal-to-noise (S/N) ratios. We also installed a broadband seismometer at the center of the seismic array (KU.KRSY in Fig. 1).

We applied the MUSIC analysis to the observed array data to determine the direction of volcanic tremor rays. The MUSIC spectrum is more sensitive to the directional variation of a beam approaching the array, compared with ordinary beam-forming analyses like frequency-wavenumber spectrum, semblance and slant stacking (Schmidt, 1986; Goldstein and Archuleta, 1991). In this study, MUSIC analysis is applied to a complex waveform comprising the seismic wave and its Hilbert transform in real and imaginary parts, respectively. A correlation matrix R is obtained from the waveforms in an array with M sensors:

((1))

where the parentheses indicate the ensemble average within a time window of the waveform, and the asterisk indicates a complex conjugate. An eigenvalue decomposition of the matrix determines both signal, and noise, subspaces composed of large and small eigenvalues. The eigenvectors in the noise subspace and the array response function provide the MUSIC spectrum as a function of the wavenumber vector k. The spectrum for a narrow frequency band of f Hz is defined as,

((2))

with

where e denotes eigenvectors in the noise subspace, and

a: is the array response function, and where the superscript H indicates a Hermitian matrix, L is number of signals comprising the waveform, and x _i is the positioning vector of the i-th sensor. We assume that the number of signals (L) is equal to 1 in this analysis, since the maximum eigenvalue of the correlation matrix dominates the others (see Fig. 2(c)). Many explosive eruptions and volcanic tremors occurred at the volcano until early February 2011. We explored tremor signals with durations longer than 30 min. The MUSIC analysis is carried out on the vertical component of the seismograms obtained by the seismic array. The MUSIC spectrum is calculated for a correlation matrix obtained from data in a time window of 4 s. The waveform data are passed through a band-pass filter of 1.5–2.5 Hz. The time window starts from the beginning of the tremor, and is shifted every 2 s until the tremor ceases. The slowness of the MUSIC spectra range from −1.2 to 1.2 s/km with a 0.01 s/km step in the north-south and east-west directions. Due to the array shape, spatial aliasing of MUSIC spectra in this analysis does not occur within the slowness range (see Fig. 1). In order to eliminate velocity heterogeneity just beneath the profile, we adopt static correction factors of the sensors to obtain the correlation matrix between waveforms, which is required for the MUSIC spectrum calculation. Here, the residual of fitting the plane wave to the direct P-wave in the seismograms of the explosion-quake serve as the static correction factors. Figure 2(a) shows sections of the MUSIC spectra in the time-azimuth domain for a fixed slowness value of 0.75 s/km. The normalized signal power is plotted here to identify the ray direction of the seismic waves. The color indicates the logarithmic strength of the spectrum, as shown in the scale bar. This tremor continued for about 30 min. The time-azimuth sections of the tremor, for the second time window in Fig. 2(a) for various slowness ranges, are displayed in Fig. 2(c). These figures clearly show that most of the waves propagated toward the − 115° azimuth (i.e., N115°W, which means that the waves came from the direction of the crater, N65°E). In addition, at several intervals, the waves arrived from N30°W, as indicated by the purple arrows in Fig. 2. This implies that the source of the seismic waves changed location during the tremor sequence.

Fig. 2.

MUSIC spectrum variations with lapse time. (a) Time-azimuth section of MUSIC spectrum at a slowness of 0.75 s/km. Each section is extracted for a time window indicated by a red bar with a numeral below the seismogram, at the top of the figure. The red dashed line indicates the crater direction. The thick arrow points to the spectrum from where the wave arrives, the N30°W direction. Warmer colors show the increasing strength of the spectrum; (b) Same as (a) except for the time-slowness section at the azimuth N65°E. The spectrum is masked by a dark shade at the lapse time when the wave does not come from the direction of the crater. The red dashed line indicates the reference slowness of the P-wave estimated from the explosion. (W), (X), (Y), and (Z) are wave types (see text); (c) Time-azimuth section of the spectra for various slowness values. The white dashed area indicates S-waves from N30°W. Most of the strong signals are composed of waves from the crater, with various slowness values. The eigenvalue of the correlation matrix for the lapse time-window variation in the MUSIC estimation is shown at the bottom of the figure. The line color indicates the order of the eigenvalue as shown in the legend. The seismograms displayed at the top of (a) and (c) are the waveforms in the vertical component at station KR25 (see Fig. 1).

The standard deviation of slowness in the MUSIC spectrum estimation, as defined by Goldstein and Archuleta (1991) gives

((3))

where δt is the deviation of the measured phase arrival time from the expected plane wave arrival time, M is the number of sensors, and Δx is the station separation of the sensors. In the second term of Eq. (3), S _n is the signal-to-noise ratio of the seismogram, n is the number of samples for the estimation, Γ is the aperture of the array, and f is the frequency. The first term of Eq. (3) relates to the heterogeneous structure just under the surface below the array. The second term corresponds to the estimated MUSIC spectrum. Substituting the parameters for the deviation in this study, which are M = 23, n = 4000 (4 s × 1 kHz sampling), S _n > 2, Γ = 0.4 km, f = 2 Hz, we obtain a second term value of less than 10⁻³ s/km; therefore, the accuracy of the MUSIC spectrum estimation relating to the array specification is very high. On the other hand, we estimated the first terms for P- and S-waves from the plane wave fitting to the direct P- and S-waves in the seismograms of six natural earthquakes, giving P- and S-wave terms of 0.038 and 0.116 s/km, respectively. The value of the first term is much larger than the second term, which implies that surface heterogeneity strongly affects the MUSIC estimation. However, this also implies that the slowness change in the lapse time windows is reliable, because the amount of change is larger than the standard deviation.

3. Estimated P-, and Surface, Wave Slowness

The P-wave velocity structure in this area has been estimated using tomographic analysis (Nishi and Kagiyama, 2002). However, the ray direction should be estimated from the velocity profile below the array. A reference event whose location is known improves the estimation. As such a reference event, we took an explosive eruption that occurred on 1 February, 2011, at the Shinmoedake crater. The seismograms recorded by the array and the MUSIC spectra derived from them are shown in Fig. 3. These are P-, and surface, waves which, according to the slowness analysis, arrive at the array with a slowness of 0.42 and 1.0 s/km from the crater, respectively. The wave types can be identified from the particle motion of the ground velocities. The P-, and surface, waves are from the direction of the crater. The radial motion and retrograde orbit imply P-, and surface, waves (Rayleigh waves), respectively. In addition, we applied polarization analysis to the seismogram every 0.5 s, with a 1-s window. Azimuths of the eigenvectors with maximum and moderate eigenvalues at the surface wave arrival are similar to those for P-wave arrival (i.e., comparison value in the green window in Fig. 3(c) with an orange window). In Fig. 3(b), the normalized ratio of e₁ by (i.e. e₁/e_tot) and its product with e₂ (i.e. e₁e₂/e_tot), where e₁, e₂, e₃ are maximum, moderate, and minimum eigenvalues, are plotted for the lapse time variation. The e₁/e_tot at the surface wave arrival (values in the orange window) becomes smaller than those at the P-wave arrival (green window). On the other hand, the e₁e₂/e_tot in the orange window is large compare with that in the green window. These results show that the particle velocity at the surface wave part traces in the radial-vertical plane. This also suggests that the later-arriving wave is a Rayleigh-type wave. The slowness values for P-, and surface, waves are used as a reference for estimating the direction of the tremor source. Namely, the incident angle of the seismic waves to the array can be estimated relative to that of the explosion, based on the discrepancy between their slowness values. Waves with slowness values that are greater than the P-wave slowness (0.42 s/km) are likely to be composed of S-, or surface, waves. Tremors with a slowness of 1.0 s/km can be interpreted as surface waves. The direction of the tremor’s source is determined relative to the reference event.

Fig. 3.

MUSIC spectrum analysis applied to the explosion on 1 February, 2011. Top: seismograms of the vertical component observed by the array. Center: time-slowness section of the spectrum from N65°W used to obtain the reference slowness of P-, and surface, waves (0.42 and 1.0 s/km, respectively). Bottom: particle motion analysis applied to the broadband seismogram (KU.KRSY): (a) three-component seismogram; (b) ratios of e₁ to e_tot, and e₁ × e₂ to e_tot, where e₁ and e₂ are maximum and moderate eigenvalues. e_tot is defined in the text; (c) azimuths of the eigenvectors; (d) right and left particle motion are extracted from P-, and surface, wave arrivals in the green and orange time windows, respectively. Both orbits indicate the ground velocity in a radial direction from the crater.

4. Origin of the Volcanic Tremor

As seen in Fig. 2, seismic waves approaching the array propagate from the direction of the crater, except for the waves around the intervals indicated by the arrows. Firstly, we focus on the waves that do not originate from the crater. Figure 4 shows the time-slowness sections at azimuths 65° and 150° around the interval indicated by the arrow in section (2) of Fig. 2, and the section in the two-dimensional slowness domain. The wave direction originating from the crater suddenly changed to N30° W; it remained so for about 30 s and then returned to the crater. The slowness value was 0.63 s/km, which is larger than that of the P-wave; as discussed above, this implies that these are S-waves. Interestingly, this direction almost coincides with that of the location of the pressure source from the geodetic data (Geographical Survey Institute of Japan, 2011). The slowness analysis can only reveal the direction of the source. The possible location of the tremor source may be determined from the ray path approaching the array. The reliability of the source location could be expressed by spatial probabilistic parameters (e.g., Scherbaum et al., 1997; Saccorotti et al, 1998; Almendros et al., 2001). In this study, we estimate the spatial likelihood distribution of a tremor source as described in Scherbaum et al. (1997). The spatial likelihood at a point (x) in the medium can be calculated from the discrepancy between an optimum slowness vector (s _mx ) from MUSIC analysis and a slowness vector (s(x)) of the seismic ray from the point to the array. For simplicity, we assume a Gaussian error distribution of slowness with a standard deviation of σ as given by Eq. (3). The likelihood is defined as

((4))

where s(x) is determined using a velocity structure, x and the array location. Unfortunately, we do not have information about the S-wave velocity structure; therefore, we adopt a simple depth-dependent S-wave velocity structure (v _s (z) = 1.7 + 0.7z) with the assumptions: (1) 90% of the S-wave velocity at the surface is the surface wave velocity obtained above; and (2) an S -velocity structure calculated from P-waves by Tsutsui et al. (1996) with V _p /V _s = 1.73. We then determined the assumed structure after smoothly connecting the two velocity structures. The spatial distribution of likelihood at the wave from N65°E (T1 in Fig. 4) and N30°W (T2) are shown in Fig. 5. Two vertical cross-sections indicate the distributions of the likelihood on the planes shown in Fig. 5. One is the vertical plane which passes through the positions of the crater and the pressure source. Another plane includes the locations of the array and the crater. As can be seen from Fig. 5, the likelihood becomes large along the ray path of waves that approach the array in the optimum direction obtained by the MUSIC estimation. The wave from the direction of the crater at the earlier lapse time reached the array from a source located within the high-likelihood area, which includes the subsurface of the crater. During this eruption, Shinmoedake was the only crater emitting volcanic material. Therefore, we conclude that the tremor source generating waves originating from the direction of the crater is located just beneath the crater. On the other hand, the pressure source inferred from the geodetic data is located within the high-likelihood area, although the area extends over several square kilometers. It is considered that at least the pressure source is one of the possible sources of the tremor. The ray direction of N30°W also corresponds to a low-velocity region, as estimated by a previous tomographic study (Nishi and Kagiyama, 2002). The first important characteristic of the observations is that the estimated source locations changed during the tremor, which persisted for only a few minutes. The findings indicate that there might be at least two major tremor sources—one located at/below the crater, and another around the pressure source. This suggests that the tremor occurs at both ends of the conduit connecting the magma source to the crater, and was activated in a tremor sequence.

Fig. 4.

Time-slowness section of the spectrum at azimuths N65°E and N150°E. The time window is taken as shown in the green box in the KR25 seismogram. This part corresponds to the lapse time indicated by the arrow in section 2 of Fig. 2(a). T1 and T2 indicate the lapse times when the spatial likelihood distributions are plotted in Fig. 5. The two-dimensional slowness section for the time window in the pink square in the left sections is shown in the right figure. The green dashed line indicates the direction of the pressure source in Fig. 1. The slowness of the wave indicated by the green arrow in Fig. 2(a) is estimated to be 0.63 s/km from N30°W, which implies an S-wave.

Fig. 5.

Vertical cross-sections of the likelihood distribution at the lapse times T1 and T2, shown in Fig. 4. A hotter color corresponds to a higher likelihood value, as shown by the scale bar. Left and right sections are displayed along the red and blue lines in the inserted map at the upper right, respectively. Red and yellow triangles show the location of the crater and the seismic array. The pink circle is the location of the pressure source obtained by the Geographical Survey Institute of Japan (2011).

Next, we investigate the variation in slowness within the part of the seismic ray originating from the direction of the crater. Figure 2(b) displays the cross-section in the slowness-time domain for a fixed azimuth of 65°. Generally, the slowness varies from 0.15 to 1.0 s/km. We divide the seismogram into five types, according to the detailed slowness and ray directions of the tremor. These are wavelets from the crater with slowness values (s) of approximately (W) s < 0.42 (P-wave slowness); (X) 0.42 ≤ s < 1.0; (Y) s = 1.0 s/km (surface wave slowness); and (Z) from the direction N30°W. The waves (X) are S-waves because the slowness exceeds that for P-waves. The slowness variation corresponds to the change in the source depth and/or the wave type. The waves categorized as (W) and (X) probably travelled from the source at depths at-, and below, the location of the explosion source. The source for the waves (Y) is probably located close to the surface of the crater because the waves propagated at surface wave speeds. As seen in Fig. 2(b), the sequence of the tremor occurrence has an interesting feature, in that wave type (Y) seems to follow the occurrence of the tremor with types (W) and (X).

5. Discussion and Summary

The characteristics found in the Kirishima volcano suggest a model for the occurrence of volcanic tremors. Once a tremor occurs at either edge of the conduit, an instability, such as shaking or wobbling, can be generated in the material inside the conduit due to flow and pressure gradients. Then tremor generation at an edge could induce tremor at another point because of the instability. A change in wave type can be caused by a phenomenon in which magma bubbling up at/below the crater generates seismic waves (e.g., tremor type waves (W) and (X)). Then, as the magma rises towards the crater, the oscillations within the magma create tremor type waves (Y). Such a magma flow at the crater could reduce the pressure inside the conduit, probably generating a magma flow wave at/around the pressure source (e.g., activating the source of waves (Z)). This qualitative model can explain the characteristics of the tremor observed in the present study. The models proposed in the previous studies can explain the phenomena of the short-term change in source location. If the conduits that can generate tremors by the interaction between magma instability and the conduit wall (e.g. Julian, 1994) are formed just beneath the crater and above the pressure source, a tremor-source-change behavior could occur. Figure 6 shows the temporal frequency change of the tremor discussed above. The frequency composition of the dominant signal does not change much with variation of the source location. Only a few ten percent change in the predominant frequency are found. In this case, a very special condition in Julian’s model is needed to explain the observational fact that the spectral feature is not sensitive to the depth variation of the source. Such a feature might correspond to the model by Jellinek and Bercovici (2011).

Fig. 6.

Temporal variation of the spectra of the tremor analyzed in Fig. 4. The spectrum is taken from a time window of 10 sec. A hotter color shows the increasing strength of the power spectrum. The green box indicates the same as that in Fig. 4.

Volcanic tremors are phenomena commonly observed worldwide. Magma behavior is of great importance for monitoring and predicting volcanic activity. Our results have revealed part of this behavior, implying that dense seismic array observations with a high sampling rate can contribute to our understanding of volcanic activity and contribute to the prediction of eruptions. However, a multi-array analysis, which we were unable to conduct in this study, is required to determine the exact location of the tremor source and the conduit geometry. In addition, geodetic observation is also needed to monitor magma behavior. Such comprehensive observations will enable further understanding of the behavior of volcanic eruptions.