# Classification of volcanic tremors and earthquakes based on seismic correlation: application at Sakurajima volcano, Japan

Theodorus Permana, Takeshi Nishimura, Hisashi Nakahara, Nikolai Shapiro

<figure><img src="image_1.png" /></figure>

## To cite this version:

Theodorus Permana, Takeshi Nishimura, Hisashi Nakahara, Nikolai Shapiro. Classification of volcanic tremors and earthquakes based on seismic correlation: application at Sakurajima volcano, Japan. Geophysical Journal International, 2022, 229, pp.1077-1097. 10.1093/gji/gab517. insu-03595126

HAL Id: insu-03595126

https://insu.hal.science/insu-03595126

Submitted on 3 Mar 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Classification of volcanic tremors and earthquakes based on seismic correlation: application at Sakurajima volcano, Japan.

Theodorus Permana<sup>1</sup>,<sup>2</sup> Takeshi Nishimura,<sup>1</sup> Hisashi Nakahara<sup>1</sup> and Nikolai Shapiro<sup>2</sup>

Accepted 2021 December 20. Received 2021 December 1; in original form 2021 June 23

## SUMMARY

Classification of volcanic tremors and earthquakes is an important part of volcano monitoring. Conventional classification technique relies on visual characterization of the amplitudes, frequency, or duration of seismic signals. In some cases, such classification is difficult to be maintained continuously, especially during high eruptive activities. In the present study, we develop a classification technique based on analyses of seismic wave correlation. The technique is based on a measure of spatial coherence using the distribution width of eigenvalues from the seismic covariance matrix that is referred to as 'spectral width'. We use the spectral width characteristics to extract the frequencies where the seismic signals are the most spatially coherent. We analyse 6-month continuous seismic data at six seismic stations deployed on Sakurajima volcano, Japan. The classification is performed every 10 min to recognize volcanic tremor and B-type earthquakes, both of which show unclear onsets and are dominant at around 1–4 Hz. Their frequencies of high spatial coherence are different from each other, providing a basis for automatic classification of both types of seismic events. Our classification results show that an increase in seismic activity, particularly volcanic tremors, occurred during high eruptive activity and is well matched with the events that are routinely determined by the Japan Meteorological Agency (JMA) based on visual examination of the data. We discuss the volcanic activities by combining our classification results with those from detections and location methods based on seismic correlation. Most of the source locations of volcanic tremor and B-type earthquakes are distributed at the same region with depths of less than 4 km beneath the active craters, which are almost consistent with the source locations of explosion earthquakes and B-type earthquakes determined from P-wave onsets in the previous studies. Also, the distribution of horizontal locations is spatially elongated at shallow depths above a dyke that is inferred to extend from Aira caldera at a depth of 7–9 km. We obtain deep and shallow sources for the tremor that occurred before and after an explosive eruption, respectively. We also obtain shallow sources beneath the active crater during a continuous tremor that accompanies increasing eruptive activity. The source locations of the classified volcanic tremor and B-type earthquakes may represent the movement of magma and gas from the shallow magma chamber at 5 km depth to the active craters.

**Key words:** Body waves; Volcano seismology; Wave propagation; Volcano monitoring.

## 1. INTRODUCTION

Monitoring of volcanic tremors and earthquakes are one of the basic and important monitoring tools for evaluating the volcanic activity. They provide insights into the underground physical processes that are related to the movements of volcanic fluids, such as magma, volcanic gas and hot waters (e.g. Chouet et al. 1987; Iguchi 1994; Nishimura & Iguchi 2011). A common type of volcanic earthquake is volcano-tectonic earthquakes (VTs) that are

generated by brittle failure due to magma pressure and/or regional stress. VTs generally show clear arrival times of *P* and *S* waves and rapid changes in seismic amplitudes that make them easily detected and located. On the other hand, volcanic tremors often show irregular variations of amplitude, emergent onset and long duration that ranges from minutes to days. Also, low-frequency earthquakes and B-type earthquakes often lack clear *P* and *S* waves. Since these events are considered to be directly or indirectly generated by volcanic fluids undergrounds, monitoring their sources over time

<footnote n="1"><sup>1</sup>Department of Geophysics, Graduate School of Science, Tohoku University, Sendai 980-8578, Japan. E-mail: theodorusp@sci.hokudai.ac.jp</footnote>

<footnote n="2"><sup>2</sup>Institut des Sciences de la Terre, UMR CNRS 5275, F-38041 Grenoble, France</footnote>
enables us to obtain information on volcanic fluids migration before and during an eruption (e.g. Ogiso et al. 2015; Kurokawa et al. 2016).

Conventional monitoring methods by volcano researchers and monitoring agencies involve the detection and classification of various kinds of seismic signals at active volcanoes by visual examination of the amplitudes, dominant frequency and/or duration of seismic signals, as well as related surface activities. Such procedures are difficult when the activity of volcanic tremors and earthquakes becomes very high, especially at the time of eruption. In addition, the classification criteria may not be unified between the researchers, staffs and agencies. Hence, automatic and/or semi-automatic methods using digital seismograms have been proposed. There are conventional methods that rely on clear and rapid changes in seismic amplitudes and spectra, such as the short-time average/long-time average (STA/LTA) method (Allen 1978) and the methods using the kurtosis and skewness of seismic signals (e.g. Saragiotis et al. 2002). For volcanic earthquakes with emergent onsets like volcanic tremor, signal coherence over several stations in a seismic network is used for detection. For example, the beamforming method (e.g. Rost & Thomas 2002) computes the slowness or apparent velocity of coherent signals, assuming the seismic source of interest is spatially localized at a distant region. Droznin et al. (2015) used stable and coherent cross-correlation functions (CCFs) as a 'fingerprint' for detecting volcanic tremors. Seydoux et al. (2016) measured spatial coherence by using the eigenvalues of the seismic covariance matrix that is obtained from seismograms at several stations, and detected volcanic tremors and earthquakes at Piton de la Fournaise volcano. Soubestre et al. (2018) used those eigenvalues to identify different episodes of volcanic tremors at the Klyuchevskoy group of volcanoes.

The detection methods using seismic wave coherence do not discriminate the types of seismic events. Classification is possible if the characteristics of each type of event such as waveform or spectral shape and duration are known. But, defining those characteristics requires a detailed examination of seismic activity at a target volcano. Unique waveform shapes of different types of events were extracted and used in template matching detection, such as the studies of non-volcanic low-frequency earthquakes in Japan (Shelly et al. 2007) and Taiwan (Tang et al. 2010). However, volcanic tremor amplitudes are often irregularly changed and do not have a specific shape. Also, several types of seismic events may have similar characteristics, such as the frequency contents of long-period events (LPs) and volcanic tremors (e.g. Fehler 1983). Several studies have proposed to use machine learning and artificial intelligence for automatic classification of volcanic earthquakes and tremors (e.g. Malfante et al. 2018; Titos et al. 2018).

In this study, we develop a classification technique based on seismic correlation analyses to discriminate volcanic tremors and earthquakes from continuous seismic records. Our technique is able to classify seismic events with unclear onsets with no *a priori* information. In Section 2, we describe the data recorded by a volcano monitoring network deployed at Sakurajima volcano with an interstation distance of about a few kilometres, which is a common scale at many monitored volcanoes. In Section 3, we extend the detection method of Seydoux et al. (2016) to classify volcanic tremors and earthquakes. For the discussion of the classified seismic events, we also apply the intercomponent correlation coefficient detection method (Journeau et al. 2020) and the CCF-based source scanning algorithm (SSA) location method (Permana et al. 2020). Both methods also involve the computation of CCFs. In Section 4, the

classification method along with the detection and location methods are applied to 6-month continuous seismic data at Sakurajima volcano, Japan, during which tremors and earthquakes are often observed during increasing eruption activity. We evaluate the detection, classification and location results and discuss two cases of volcanic activities, particularly when tremor is present. We compare the results with infrasound data and reported activities to understand the volcanic processes at Sakurajima. In Section 5, we evaluate the classification method by comparing with other methods, volcanic activity reports and studies on volcanic system at Sakurajima.

## 2. DATA AT SAKURAJIMA VOLCANO

Sakurajima volcano is a stratovolcano located at Kyushu Island in southern Japan. It is a former volcanic island located at the southern edge of the 20-km wide Aira Caldera in Kagoshima Bay. The island was connected with mainland Kyushu by the lava flow during the Taisho eruption in 1914. Today, it is a round-shaped land with a dimension of 12 km × 9 km and a maximum elevation of 1.1 km above sea level (Fig. 1). Sakurajima is one of the most active volcanoes in Japan. Since 1955, eruptive activity mainly occurred at Minamidake crater until 2006 when the activity shifted to Showa crater (Kobayashi et al. 2013). Recently, eruptive activity at the Minamidake crater has resumed (JMA 2018).

Sakurajima volcano has been closely monitored by the Japan Meteorological Agency (JMA), universities and the local government. Seismometers, infrasound sensors, tilt and strain meters, as well as ash samplers are installed on the volcano flanks. In the present study, we analyse seismic data at six JMA stations for 6 months from 2017 April to September (Fig. 1). Each station is equipped with a short-period (1 Hz) three-component seismometer. We preprocess the data by applying mean and trend removal and correcting for instrument response. We also use the reports of seismic and volcanic activities by JMA during 2017 April–September: the occurrence time of volcanic tremors and earthquakes, the occurrence time of eruptions, the location of eruption (crater), plume height and various information of volcanic activities obtained from instrumental and visual observations.

At Sakurajima, volcanic tremors and earthquakes are observed before, during and after an eruption. Examples of one-minute waveforms of tremors and earthquakes are shown in Fig. 2. Volcanic earthquakes are commonly classified into A-type, B-type and explosion earthquakes (e.g. Iguchi 1994), following the classification by Minakami (1974). A-type earthquakes show clear *P*- and *S*-wave onsets and are typically dominant at high frequencies of about 10 Hz or larger (Figs 2a and b). Brittle shear faulting in the rocks has been thought of as the mechanism of A-type earthquakes. B-type earthquakes often show unclear onsets of *P* and *S* waves and dominant frequencies of 1–5 Hz (Figs 2c and d). This type of earthquake is attributed to the fluctuations in the magma and gas in the conduit. Although the onsets are usually unclear, some B-type earthquakes have been located beneath the craters at depths of about 3 km using the onsets of *P* wave (Iguchi 1989; Iguchi 1994). Explosion earthquakes, which sometimes show clear onsets and dominant frequencies of 1–5 Hz (Figs 2e and f), are associated with an explosion that rapidly ejects volcanic material (Iguchi 1994). JMA catalogue during the period of 2017 April–July reported 55 A-type earthquakes, 1792 B-type earthquakes and 14 explosion earthquakes. The JMA catalogue also lists volcanic tremors: we classify them into 312 tremor
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Map of Sakurajima volcano and the seismic network. Colour denotes elevation from sea level in km. Blue lines represent the outlines of craters. Black triangles represent the JMA seismic stations used in this study. Map in the upper left shows the location of Sakurajima volcano in Japan (red square).</figcaption></figure>

episodes with durations ranging from 1 to about 10 min (referred to as 'short-duration tremor'), and 22 tremor episodes with longer durations (referred to as 'continuous tremor'). The tremors show dominant frequencies at around 1–4 Hz and occur continuously or intermittently, usually during eruptive activity. Those tremors often show unclear onsets and their amplitudes vary irregularly (Figs 2g and h). Volcanic tremors are also attributed to the activity of magma and gas in the conduit.

Reported eruptions in JMA catalogues are commonly classified into non-explosive and explosive eruptions that are referred to as 'eruption' and 'explosion', respectively. An 'explosion' is an eruption that generates significant explosion earthquakes and/or infrasound signals while an 'eruption' does not accompany those phenomena (JMA 2017).

# 3. METHODS

## 3.1 Analysis of continuous seismic data

We analyse the seismic records at Sakurajima volcano for consecutive, non-overlapping 10-min windows. This time window length is set to be sufficient to obtain a good performance of the detection, classification and location for volcanic tremors occurring at Sakurajima volcano. Each window is divided further into smaller overlapping subwindows. Fig. 3 shows an illustration of seismic velocity records, $u(t)$, of $N$ stations that are divided into 10-min time windows and their subwindows. We change the length and overlap of the subwindows according to the purposes of the analyses in the following sections. To avoid unwanted artifacts in the seismic spectra, the data at each subwindow are tapered by using a 5 per cent tapered cosine window. Therefore, the overlap of the subwindows is necessary to avoid losing information from the tapered part of the data.

## 3.2 Detection characteristics of volcanic tremors and earthquakes using network covariance matrix analysis

We distinguish the seismic signals of volcanic tremors and earthquakes from ambient noise by evaluating the spatial coherence of recorded signals over a seismic network. Volcanic tremors and earthquakes are assumed to be generated at a localized region, while seismic signals of ambient noises consist of the waves generated by many sources distributed around the stations. Following Seydoux et al. (2016), we assume that a seismogram consists of seismic signals from a number of independent sources which is represented by the rank of the covariance matrix. The largest eigenvalue of this matrix represents a coherent signal while the smaller eigenvalues correspond to noises. The width of eigenvalues distribution at a specified frequency is called the 'spectral width', where a smaller spectral width indicates higher spatial coherence. Averaged spectral widths over a frequency band are used for 'detection' by identifying those showing high spatial coherence. A brief description about the method and the preliminary results using our Sakurajima data are provided in Appendix A. We use seismic data from the vertical component of the seismometers, and divide the 10-min windows into 25 s subwindows with 80 per cent overlap.

The network covariance matrix analysis does not discriminate between different types of seismic events. Therefore, we further develop a method to classify volcanic tremors and earthquakes. First, we investigate the spectral width characteristics of different types of events by computing the spectral width every 0.02 Hz in a frequency band of 0.5–10 Hz. We select only 10-min windows that contain only one type of events according to their occurrence times in the JMA catalogue. We only use the time windows with a minimum spectral width of 1 or smaller that correspond to the seismic signals characterized by high spatial coherence. We average the spectral width as a function of frequency over the selected 10-min windows for each type of events and refer to them as the
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Examples of one-minute waveforms (left-hand column) and their amplitude spectra (right-hand column) of volcanic earthquakes and tremors that occurred in 2017 at Sakurajima volcano: (a) and (b) A-type earthquake, (c) and (d) B-type earthquake, (e) and (f) explosion earthquake and (g) and (h) volcanic tremor.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Schematic illustration of the analyses. Seismic records $u(t)$ from $N$ stations are divided into 10-min time windows and further divided into subwindows. The subwindows are shown with no overlap for simplicity.</figcaption></figure>

characteristic spectral widths. A total of 1090 windows are used to obtain the characteristic spectral width of continuous tremor, 366 windows for short-duration tremor and 1292 windows for B-type earthquake. However, there are only 26 selected windows for A-type and 5 for explosion earthquakes. These small numbers of windows are not sufficiently representing their spectral widths characteristics. In the following analysis, we give a focus to the classification of volcanic tremors and B-type earthquakes.

These characteristic spectral widths are shown in Figs 4(a), (d) and (g) for continuous tremor, short-duration tremor and B-type earthquake, respectively. For each 10-min window, we extract the

frequency of the 10 smallest spectral width minima. The distribution of all extracted minima is then presented as a histogram by using 0.02 Hz bins. The spectral width minima are shown in Figs 4(b), (e) and (h) and their histograms in Figs 4(c), (f) and (i).

The shape of the characteristic spectral width for the two categories of tremors is very similar, indicating that both tremors are generated by the same volcanic process. High coherence (i.e. low spectral width) is generally observed at the frequencies of 0.7–7 Hz, which is wider than the expected dominant frequency of volcanic tremors. Their histograms show clear peaks at similar frequencies between 0.9 and 2 Hz. The spectral width minima of continuous tremor are more concentrated at this frequency band than those of short-duration tremor (Figs 4b and c). On the other hand, the characteristic spectral width and the histogram for B-type earthquakes show a dominant minimum at about 1.25 Hz. The histogram shows low and high peaks in the frequency band of 0.9–3 Hz. The most dominant peaks are found at 1.25 and 1.8 Hz.

These frequencies having high spatial coherence are approximately equal to the dominant frequencies of each type of the events (Fig. 2). Moreover, the histogram peaks represent a unique set of frequencies where the seismic wave is the most coherent for each
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Spectral width characteristics of (a)–(c) continuous tremor, (d)–(f) short-duration tremor and (g)–(i) B-type earthquakes. (a), (d) and (g) the spectral width versus frequency for each 10-min window (grey lines) and their average (black line). (b), (e) and (h) 10 smallest spectral width minima from each time window. (c), (f) and (i) the histogram of the spectral width minima.</figcaption></figure>

type of event. We use these frequencies to improve the detection method by classifying the seismic events.

## 3.3 Classification method of volcanic tremor and B-type earthquake

To examine whether a time window contains the signal from volcanic tremor or B-type earthquake, we may compare the spectral width as a function of frequency with the characteristic spectral width of both types of seismic events. The spectral widths characteristics of continuous and short-duration tremors are similar to each other so that we do not analyse them separately, although we may classify them by examining the duration time of signal amplitudes. Here, we propose to compare the minima of the spectral width instead of using all spectral width information. For a 10-min time window, we compute a detection measure specific to each

type of events:

$$ \varepsilon_h = \frac{\sum_f w_h(f) [\sigma(f)^{\text{env}} - \sigma(f)^{\text{smooth}}]}{\sum_f w_h(f)}, \quad (1) $$

where the subscript h represents different types of seismic events: T for volcanic tremor and B for B-type earthquakes. σ(f) is the spectral width as a function of frequency that is obtained using the procedure described in Appendix A (Fig. 5a). σ(f) smooth is the smoothed spectral width that is obtained from the convolution of σ(f) with a 0.2 Hz-wide Hann window. σ(f) env is the upper envelope of the smoothed spectral width which is obtained by applying spline interpolation on the local maxima of σ(f) smooth (Fig. 5b). The σ(f) is transformed to σ(f) env − σ(f) smooth to reduce local variations in the frequency domain and to focus only on the spectral width minima instead of the shape of the spectral width curves (Fig. 5c).
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Example of $\sigma(f)^{\text{env}} - \sigma(f)^{\text{smooth}}$ calculation for a 10-min window and the weight functions for calculating $\varepsilon_h$. (a) Spectral width $\sigma(f)$. (b) Smoothed spectral width $\sigma(f)^{\text{smooth}}$ (blue line) and its envelope $\sigma(f)^{\text{env}}$ (magenta line). (c) $\sigma(f)^{\text{env}} - \sigma(f)^{\text{smooth}}$. (d) Weight function for tremor $w_T(f)$. (e) Weight function for B-type earthquakes $w_B(f)$. Note that the x-axis values in (b) is different from (a) because $\sigma(f)^{\text{smooth}}$ is a product of a convolution of $\sigma(f)$ with a Hann window, which involves multiplication and summation of the two data.</figcaption></figure>

The weight function $w_h(f)$ is uniquely defined for each type of event based on the peaks of the histogram of spectral width minima. We first average the histograms of short-duration and continuous tremors. Then, we normalize the histograms of tremor and B-type earthquake by the maximum peak and select the frequencies between 0.9–2 Hz with normalized peak value larger than 0.5. In the next step, we iteratively examine each selected frequency to determine its weight: for a selected frequency $f_s$, an initial weight function is defined with an initial weight of $w_h(f_s) = 1$ and 0 at other frequencies. The $\varepsilon_T$ and $\varepsilon_B$ at $f_s$ are calculated for each time window of tremors (continuous and short-duration) and B-type earthquakes, respectively. Then, we average $\varepsilon_T$ from a total of 1456 windows (1090 continuous tremors and 366 short-duration tremors) and $\varepsilon_B$ from 1292 windows of B-type earthquakes. The calculation is repeated to obtain the average $\varepsilon_T$ and $\varepsilon_B$ for all of the selected frequencies. The average $\varepsilon_T$ values are then normalized by the largest value in all of the selected frequencies and are then used as the weight function $w_T(f)$. Similarly, we normalize the average $\varepsilon_B$ values to obtain $w_B(f)$. We further set the weights of <0.5 to

zero. The final weight functions are shown in Figs 5(d) and (e) for volcanic tremor and B-type earthquake, respectively.

In Figs 6(a)–(c), we show an example of $\varepsilon_T$ calculation using the spectral widths for a volcanic tremor on May 12. The value of $\sigma(f)^{\text{env}} - \sigma(f)^{\text{smooth}}$ (Fig. 6b) is high at several frequencies that correspond well with those of non-zero weights (Fig. 5d). When $\sigma(f)^{\text{env}} - \sigma(f)^{\text{smooth}}$ at those frequency is large, the obtained $\varepsilon_T$ comes to be large and suggests an occurrence of tremor (Fig. 6c). Figs 6(e)–(g) show an example of $\varepsilon_B$ calculation for a 7-hr period on June 17 when many B-type earthquakes occurred. The detection measure $\varepsilon_B$ is high especially when the spectral width minima at around 1.25 Hz is present (Fig. 6g).

To judge which type of events occurs, we introduce a classification value $L$ that is computed by subtracting the detection measure of B-type earthquake $\varepsilon_B$ from that of volcanic tremor $\varepsilon_T$:

$$L = \varepsilon_T - \varepsilon_B. \qquad (2)$$

A positive $L$ value indicates the characteristics of volcanic tremors (Fig. 6d) while negative value indicates that of B-type earthquakes (Fig. 6h). For the time windows when both volcanic tremor and B-type earthquakes are present, the $L$-value is closer to zero than the case when only volcanic tremor or B-type earthquake is present. $L$ is also close to zero for A-type earthquakes, explosion earthquakes and noise.

## 3.4 Detection using intercomponent correlation coefficient

We use the intercomponent correlation coefficient method as a complementary detection method to discuss the classified volcanic activity. For the description of the method, please refer to Journeau et al. (2020). This method assumes that if the source of volcanic tremor is stable during a given time window, the intercomponent cross-correlation of the seismic signal at a single station is also stable in time and serves as a ‘fingerprint’ of that tremor. We analyse seismic data at 1–4 Hz using 100 s subwindows with 50 per cent overlap, where we compute the CCFs from three available components: north–south, east–west and vertical components. We further average the correlation coefficients over the three component pairs and then over all stations. The preliminary analysis using Sakurajima data show that the correlation coefficients are obtained within the range of 0.2–0.4. We assume that a higher correlation coefficient represents a more stable seismic source in a 10-min window and may indicate the detection of volcanic tremor or earthquake.

## 3.5 Location using CCF-based source scanning algorithm

The source location of seismic signals in a 10-min window is determined using the CCF-based SSA method described in Permana et al. (2020). The method uses the information on the traveltime difference between seismic waves recorded at two stations. This traveltime difference is represented by the lag time of the peak amplitude of a CCF. The seismic source is determined using a grid search scheme to find a location where the theoretical traveltime differences match the observed ones from the CCF peaks. The CCF-based SSA extends a previous location method (Droznin et al. 2015) to 3-D space by using ray tracing in heterogeneous velocity model to estimate the theoretical traveltime differences. A minimum of three station pairs is required, where more pairs will improve the accuracy. To properly locate the source, direct seismic waves are expected to be dominant. Seismic scattering may reduce the location accuracy through the diffraction of these direct waves. We suppress
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Example of $ε_h$ and $L$ computation for 7-hr data of (a)–(d) volcanic tremor on May 12 and (e)–(h) B-type earthquakes on June 17. (a) and (e) Spectral width $σ(f)$, (b) and (f) $σ(f)^{env} - σ(f)^{smooth}$, (c) and (g) $ε_h$ values and (d) and (h) $L$ values. $ε_B$ for May 12 and $ε_T$ for June 17 are not shown.</figcaption></figure>

the incoherent part in the CCF by computing and stacking CCFs over several subwindows (e.g. Bensen *et al.* 2007) instead of following other location methods that use CCF, such as using double correlation (Li *et al.* 2017b), higher-order correlation functions (Li *et al.* 2017a), or transforming CCF envelope into probability density (Li & Gudmundsson 2020).

We analyse vertical component seismic data at 1–4 Hz using 20 s subwindows with 75 per cent overlap. The CCFs are computed at each subwindow and stacked over a 10-min window for all available station pairs. Permana *et al.* (2020) have shown the advantage of using stacked CCF of simulated tremor and obtained location errors of about 1 km in 3-D directions. Smooth envelopes of the CCFs are then computed with a smoothing parameter equivalent to 1 s to reduce location error caused by inaccuracies in traveltime prediction. We prepare a 3-D grid with a dimension of 11 km × 14 km × 8 km, covering the area shown in Fig. 1. The gridpoint intervals are set to be 90 m × 80 m × 10 m in latitude, longitude and vertical directions, respectively. For the computation of theoretical traveltimes, we assume that direct S wave is dominant, and we prepare a 3-D S-wave velocity model where the velocity values are based on 1-D P-wave structures from Miyamachi *et al.* (2013) and Ono *et al.* (1978). The uppermost value of the 1-D velocity is set at the topography. More details on the velocity model and the preliminary results using our Sakurajima data are provided in Appendix B.

# 4. RESULTS

## 4.1 Overview of detection, classification and location results

The spectral width, classification value $L$, correlation coefficient and the source location are computed for each 10-min time window

from 2017 April to September. We use the averaged spectral width at 1–4 Hz for the detection using network covariance matrix analysis. We show the results together with the occurrence time of volcanic tremors and B-type earthquakes from the JMA catalogue in Figs 7 and 8 on 2017 May and June, respectively, and in [Figs S1.1–S1.4 of] Supporting Information S1 for the other months.

In general, continuous tremors are well detected, classified and located. The periods when continuous tremors are reported well matched with the time windows of low averaged spectral width values of about <1.35 and correlation coefficients of >0.3, as shown at time periods with pink colour in Figs 7(f) and (g). Positive $L$ values for tremors are shown as time periods with purple colour in Fig. 7(e). On the other hand, when short-duration tremors and B-type earthquakes intermittently occurred on May 13–23, the large averaged spectral width, small correlation coefficients and $L$ values around zero indicate that short-duration tremor and B-type earthquake are poorly recognized in this period. However, $L$-value can clearly distinguish B-type earthquakes in some time windows, such as when B-type earthquakes occurred and no other type of events are reported, for example, on June 8 and 17 (Fig. 8). When the activity of B-type earthquakes is high, several time windows may contain more than one earthquake, increasing the length of coherent signals. Also, there are several cases when B-type earthquakes are well classified with low $L$ but the averaged spectral width is large, such as on May 24 (Fig. 7). This is because the information on the frequencies where high coherence is expected is removed in the averaged spectral width but is used by the $L$-value. The correlation coefficient is generally less sensitive to detect short-duration tremors and B-type earthquakes compared to the averaged spectral width.

During the occurrence time of continuous tremors, the source locations are more spatially constrained. The horizontal locations during this period are mostly distributed at the flank of Sakurajima
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Detection, classification and location during 2017 May. (a) JMA catalogue of volcanic tremor and B-type earthquake occurrences. (b) The latitude (in degrees), (c) longitude (in degrees) and (d) depth (in km) of the source location. The grey scale colour scheme denotes the averaged spectral width value. Dashed line and dashed-dotted line denote the location of Minamidake and Showa crater, respectively. (e) Classification value L. (f) Averaged spectral width. (g) Correlation coefficient.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Detection, classification and location during 2017 June. Figure description is the same as Fig. 7.</figcaption></figure>
volcano, and the depths are distributed at less than 4 km. A map showing the source location of continuous tremors from 2017 April to July is provided in [Fig. S2.1 (Supporting Information S2)]. The variations in the vertical location are larger than those in the horizontal one. On May 4 and 7, tremor source locations are spatially concentrated over a longer time period than the tremor duration reported by JMA (Fig. 7). This indicates that coherent signals are captured by the CCFs over a longer time period than the reported tremors. The source locations of short-duration tremor and B-type earthquake tend to be more scattered. However, some B-type earthquakes with low averaged spectral width in June 17 (Fig. 8) show a similar source location distribution to that of continuous tremors. These B-type earthquakes show high signal-to-noise ratio (SNR) and there is no occurrences of the other type of events. On June 2–5, when tremors are reported, the source locations are distributed in a wider region compared to those in May. The averaged spectral widths are consistently low, but the correlation coefficients during these tremors are not as high as those obtained during the other continuous tremors. Inspection of the seismic waveforms reveals that this tremor episode actually consists of short-duration tremors with a duration of 1–15 min that occurred intermittently with an interval between 30 s and 20 min.

Short-duration tremors and B-type earthquakes often show larger averaged spectral width, smaller correlation coefficients and *L* values close to zero. Also, their source locations are more scattered. These are mainly because their duration times are much shorter than 10 min. In such case, the accuracy of the detection, classification and location is decreasing due to increasing noise contamination. Although it is possible to reduce the time window length for analysis, we show in Appendix C that a shorter time window produces smaller differences between the averaged spectral widths of seismic events and that of noise. These smaller differences may cause some difficulties during detection and subsequently affects the classification method. The similar effect is also observed during correlation coefficient computation. In addition, the location accuracy deteriorates when using shorter time windows.

In the following sections, we discuss two volcanic activities where volcanic tremors and B-type earthquakes are observed on May 6–8 and August 22–25. The two time periods of activities show different styles of eruption. We use the averaged spectral width, *L*-value and the correlation coefficient to select ‘reliable’ source locations. We discuss the temporal changes in the source locations by comparing with seismic and infrasound measurements and JMA reports. Then, we interpret the volcanic processes associated with the obtained source locations. We calculate the spectrograms of seismic signals at V.SFT2 station and infrasound signals at V.SKD2 station every 10-min windows with a 2-min shift. Spectrum at each time window is normalized by its maximum value. All time information is presented in local time.

## 4.2 Results on volcanic activity on 2017 May 6–8

Fig. 9 shows the analysis results, waveform data and their spectrograms for a volcanic activity on May 6–8. A tremor started on May 6 with unclear onset, then its amplitude increased at around 03:00 on May 7. The tremor maintained its amplitude before abruptly decreasing at 09:10. A low seismicity period then lasted for about 2 hr, during which several B-type earthquakes are observed. An explosion occurred at 11:08 on May 7 at Showa crater. Then, no significant seismic events are observed for about 30 min and another small eruption occurred at 11:42, followed by a tremor that shows

rapidly varying amplitudes. JMA catalogue reported the tremor after the explosion and eruption, but not the tremor before the explosion.

During the tremors before and after the explosion, infrasound signals are recorded. (Fig. 9a). This strongly suggests that the tremors are associated with eruptions effusing volcanic ash and/or gas. The infrasound spectra are dominant at 0.8–3 Hz both before and after the explosion (Fig. 9b). Seismic spectra show two peaks at about 1.25 and 2.6 Hz (Fig. 9d). The 1.25 Hz peak is wider before the explosion, while the 2.6 Hz peak does not show significant changes.

Averaged spectral width of around 1.35 or less, correlation coefficients of around 0.3 or more, and *L* of around 0.35 or more are observed during the tremors. However, higher averaged spectral width and lower correlation coefficient during the tremor before the explosion indicate lower spatial coherence (Figs 9e and f), during which the source locations are less concentrated in space at 2–7 km depths and south of the craters (Figs 9h–j). During the low seismicity period, B-type earthquakes are poorly detected by the averaged spectral width and correlation coefficient but successfully classified with *L* (Fig. 9g). The explosion is detected by the averaged spectral width, but the correlation coefficient is low and the *L*-value is around zero. The source locations during the occurrence of B-type and explosion earthquakes are spatially scattered. For the tremors after the explosions, from around 14:00 to 22:00 on May 7, the sources are mostly located at shallower depths of <1 km. Tremor amplitude decreases shortly after 22:00, and their source locations are scattered.

We select ‘reliable’ tremor source locations by setting some reliability criteria: averaged spectral width of <1.35, *L* of >0.35, and correlation coefficient of >0.3 (Fig. 10a). Selected sources are mostly distributed in the south of the craters. Before the explosion, most of the tremor sources do not satisfy the reliability criteria, except for a few of them that are located at 4–5 km depths. On the other hand, the reliable tremor sources after the explosions are determined at 2 km depth or less. We use less restrictive reliability criteria to show the source locations of B-type earthquakes (Fig. 10b): averaged spectral width of <1.7, correlation coefficient of >0.2 and *L*-value of < −0.35. Because the location accuracy is low, their source locations are distributed deeper but approximately at a similar area with that of volcanic tremors.

We validate the difference of depths between the tremor before and after the explosion by averaging the normalized CCF envelopes in four 2-hr tremor windows: two before the explosion and two after the explosion (Fig. 11). We indicate the traveltime differences from two pseudo-sources beneath the Showa crater. One is a deep source at 4.5 km depth which represents the ‘reliable’ source locations obtained before the explosion, and another is a shallower source at 0.5 km elevation. We compare the CCF envelopes for station pairs V.SKA2-V.SKRD (Figs 11a–d) and V.SKD2-V.SKRD (Figs 11e–h), which have an interstation distance of 8.06 and 1.28 km, respectively. Overall, the envelope peaks before the explosion are well matched with the traveltime differences of the deep source. However, by closely looking at the envelopes before the explosion, we see that the traveltime differences of the shallow source seem to be matched with another peak of the envelopes (Fig. 11e) or at high values close to the peaks (Figs 11a, b and f). On the contrary, the envelopes after the explosion match with the predicted traveltime differences in the opposite way (Figs 11c, d, g and h). These results suggest that the seismic signals from shallow and deep parts may be present before and after the explosion. The spatial coherence and source location accuracy for the tremor before the explosion are
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Results on volcanic activity on May 6–8. (a) Infrasound signal at VSKD2 station. Red stars denote the explosion and eruption from JMA reports. (b) Infrasound spectrogram. (c) Seismic signal at V.SFT2 station. (d) Seismic spectrogram. (e) Averaged spectral width. (f) Correlation coefficient. (g) Classification value L. (h) Latitude, (i) longitude and (j) depth of the source locations, respectively. Source location colours represent the time shown in the colour bar at the bottom. Dashed line and dashed-dotted line denote the location of Minamidake and Showa crater, respectively. Grey shaded areas denote the time of 'reliable' source locations shown in Fig. 10.</figcaption></figure>

noticeably lower. This may be because the strength of the shallow and deep sources are not so different. Since the infrasound signals were observed during the tremor episodes before and after the explosions, the shallow signals are possibly related to ash and/or gas emission from the crater.

Differences in the source depths before and after the explosion may be related to the changes in the seismic spectra, particularly the

1.25 Hz peak. These differences can be explained by the changes in the source time function of volcanic tremor and/or the seismic generation region inside the conduit. Various studies have discussed several volcanic fluid-related processes that can generate seismic waves, such as fluid-driven crack mechanism (Aki et al. 1977; Chouet 1988), gas flow (Seidl et al. 1981), gas bubbles in volcanic fluids (Leet 1988) and the interaction of fluid flow with
<figure><img src="image_13.png" /><figcaption>**Figure 10.** 'Reliable' source locations of (a) volcanic tremor and (b) B-type earthquakes on May 6–8. Source location colour represents the time according to the colour bar at the bottom of Fig. 9.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Average normalized CCF envelopes (black lines) and their standard deviation (grey-shaded areas). Each column shows a time period and each row shows a station pair. Red and green vertical lines denote the computed traveltime difference from a theoretical deep and shallow source, respectively.</figcaption></figure>

barriers in the conduit (Honda & Yomogida 1993). We found that the 2.6 Hz seismic spectral peak is not always observed during tremor occurrences. The 2.6 Hz spectral peak may be also related to the generation of seismic waves deeper in the conduit, possibly at 4–5 km depth where the tremor sources are estimated. Multiple peaks in the seismic spectra can also be explained by seismic wave excitation in the conduit by several resonators (Hurst & Sherburn 1993). On the other hand, stable infrasound spectra indicate that there are no significant changes in the processes that generate the infrasound. Although the conduit and/or vent size can affect the tremor generation and amplitude (e.g. McNutt 1994; McNutt & Nishimura 2008), the vent size and the pressure variation at the vent during the ejection of volcanic gas and ash in this episode may be stable regardless of the amplitude variations. We refer to the tremors preceding and following the explosion as 'eruption tremor' by considering that stable infrasound signal is observed during the

tremors. The period of low seismicity may represent a period of gas trapping and pressure buildup in the sealed conduit that leads to the explosion (Nishimura et al. 2013).

## 4.3 Results on volcanic activity on 2017 August 22–25

We present another example of volcanic activity on August 22–25 which shows changing style of tremor occurrences associated with changing eruption activity. On August 22 at around 22:30, a continuous tremor occurred and lasted for about 12 hr (Fig. 12c). Tremor amplitudes gradually started to increase at 00:00 on August 23. The tremor shows stable amplitudes from about 05:00 and ends at 10:50, followed by a spike of infrasound signals indicating an explosion (Fig. 12a) and a short-duration tremor with sudden onset and larger seismic amplitudes. During the following two days,
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Results on volcanic activity on August 22–25. The red stars in (a) denote the explosion time from JMA reports. The rest of the figure description is the same as in Fig. 9.</figcaption></figure>

tremor activity has changed from continuous to intermittent occurrences. Each tremor occurrence has a duration of 1–30 min with an interval ranging from 10 min to more than one hour. JMA reported nine explosions from 18:00 on August 23 to 09:00 on August 24 (JMA 2017). The activity then significantly decreased during which no explosions are reported, but several tremors with small amplitudes are observed. Six explosions are reported between 18:00 on August 24 and 02:00 on August 25. Infrasound data indicate several smaller explosions may have occurred. Several volcanic tremors are

observed shortly after the explosion. From August 21 to 28, JMA reported intermittent small eruptions that ejected volcanic blocks and lights from Showa crater during the night that caused by hot gas and lava. Several explosion sounds were reported on August 23–25, and daily amount of volcanic gas ($SO_2$) in August was the highest at the end of the month (JMA 2017).

We are able to detect both continuous and intermittent occurrence of tremors using the averaged spectral width and correlation coefficient, which show similar values with those for May 6–8 activity
<figure><img src="image_16.png" /><figcaption>**Figure 13.** 'Reliable' source locations of (a) volcanic tremor and (b) B-type earthquakes on August 22–25. Source location colour represents the time according to the colour bar at the bottom of Fig. 12.</figcaption></figure>

(Figs 12e and f). A seismic event with large amplitudes at 14:40 on August 24 shows a low averaged spectral width and low correlation coefficient. The event is classified as a B-type earthquake.

The continuous tremor initially shows source locations at 1–4 km depths. From around 22:30 to 00:30 on August 23, the tremor shows high averaged spectral width and low correlation coefficient, suggesting low spatial coherence probably due to lower SNR (Figs 12h–j). The tremor sources are then consistently located at shallow depths above sea level just beneath the Showa crater. During the intermittent occurrences, the sources are scattered but mostly located at 0–2 km depths. Between 06:00 and 23:00 on August 24, the volcanic activity decreases and the sources are mostly unreliable because the seismic records mainly consist of noises. Then, the source locations become concentrated again at less than 2 km depth between 00:00 and 03:00 on August 25, following several explosions. Fig. 13(a) shows 'reliable' source locations of volcanic tremor selected using the same criteria with those for the May 6–8 activity, which are mainly those from the continuous tremor. Since the continuous tremor sources are estimated at shallow depths beneath Showa crater, it is inferred to be an eruption tremor that is associated with the continuous emission of volcanic ash and gas. There are many source locations of intermittent tremors that do not satisfy the reliability criteria despite showing high coherence. The distribution of reliable source location is found to be vertically elongated at depths of 0–4 km. Fig. 13(b) shows 'reliable' source locations for B-type earthquakes with the same criteria as those for tremor but with $L < -0.35$. The sources are located at depth of 0–2 km beneath Minamidake and Showa craters.

Seismic spectra show that a single peak at 1.25 Hz is recognized, particularly during the continuous tremor. Infrasound spectra reveal that the frequencies of 1–1.5 Hz become more dominant compared to the tremors on May 6–8. This suggests that vent condition and/or properties of volcanic fluids have changed accompanying more frequent explosions and ejection of volcanic material. This more intense eruptive activity is reflected by the shallow source locations.

During some intermittent occurrences of tremors, both seismic and infrasound spectra show multiple peaks that change gradually with time within the 1–4 Hz frequency band (Figs 12b and d).

It is difficult to classify these tremors, because the frequencies of high coherence change with time. As a result, some tremors are not clearly classified or misclassified as B-type earthquakes (Fig. 12g). Gradual changes of seismic spectral peaks with time have been previously observed during volcanic tremors that follow Vulcanian eruptions at Sakurajima (e.g. Maryanto et al. 2008). The possible process that may cause such spectral changes is the resonance of the gas pocket in the conduit at shallow depths.

# 5. DISCUSSION

## 5.1 Comparison with other methods

We compare the classification results with the detection results using intercomponent correlation coefficients, averaged spectral width, the STA/LTA method, and the daily number of volcanic tremors and earthquakes as well as eruption plume height from JMA reports (JMA 2017). Note that the number of volcanic earthquakes includes A-type, B-type and explosion earthquakes, and that the number of volcanic tremors and earthquakes reported by JMA have been verified by human observers.

JMA reports show that volcanic tremor activity is concentrated during two time periods: from the end of April to early June with a maximum number of 25, and in August–September with a maximum number of 78 (Fig. 14a). The daily number of volcanic earthquakes is less than 100 in early April (Fig. 14b), 0–40 from the end of April–June, and larger than 100 in August–September when tremor activity also increases. The eruption plume height data in Fig. 14(h) indicates that the eruptive activity is high during the two time periods of high volcanic tremor activity.

The STA/LTA method (Allen 1978) is applied to the vertical component seismic data at 1–4 Hz. The STA/LTA method raises a detection when a ratio of average amplitude in a short time window over that in a long time window exceeds a threshold called trigger. We use 1 s short time windows, 30 s long time window and a trigger ratio of 3. The daily detections of STA/LTA method are similar to those of volcanic earthquakes reported by JMA,
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Comparison of daily number of events during 2017 April–September: (a) volcanic tremor from JMA, (b) volcanic earthquakes (including A-type, B-type and explosion earthquakes) from JMA, (c) STA/LTA method, (d) the network covariance matrix analysis, (e) intercomponent correlation coefficient, (f) classification of B-type earthquakes and (g) classification of volcanic tremor. (h) Eruption plume height (JMA). Red and black circles denote the explosion (explosive eruption) and eruption (non-explosive eruption), respectively.</figcaption></figure>

except on June 2–5 where sudden amplitude changes during the intermittent occurrence of tremors produce high amplitude ratios that trigger the detection (Fig. 14c).

For the detection based on seismic correlation, we count the daily number of 10-min windows which satisfy the following criterion for each parameter: an averaged spectral width of <1.3 and a correlation coefficient of >0.28. For the classification results, we first compute the average of *L* values at time windows where volcanic tremor is reported by JMA. Similarly, we also compute the averaged *L* values at time windows of reported B-type earthquakes. Those average values are used as the criteria to obtain the number of daily events. That is, the criteria of *L* > 0.55 and *L* < -0.70 are used to classify volcanic tremor and B-type earthquake, respectively.

The averaged spectral width detects both volcanic tremor and earthquakes, and it better detects continuous tremors (Fig. 14d). The number of detections in May are larger than those reported by the JMA, but the difference is probably because the JMA counts the number of tremor occurrences while the averaged spectral width counts the number of time windows that contain the tremors. The intercomponent correlation coefficient produces similar results with the averaged spectral width (Fig. 14e) but is less sensitive for short-duration tremors and earthquakes whose duration times are shorter than 10 min.

The classification using *L*-value also shows the highest activity of volcanic tremor during the above-mentioned two time periods (Fig. 14g). The number of daily detections in April–June when continuous tremor is dominant is larger than that in August–September,

in contrast with JMA reports. On the contrary, when short-duration tremor is dominant in August–September, some tremors are not well detected. The classification of B-type earthquakes (Fig. 14f) shows some similarity with JMA reports during April–June, where many earthquakes occurred in early April followed by low activity until the end of June. Similar with the result from volcanic tremors, the daily number of classifications in August–September is much lower than that from JMA reports.

We examine the similarity of classification results with the JMA catalogue. We count the number of time windows when volcanic tremor is reported in JMA catalogue and the *L*-value is positive and larger than a specified *L* threshold in a 10-min time window. The total number of the time windows is then presented as a percentage with respect to the total number of time windows when tremors are reported in the JMA catalogue. Similarly, we count the time windows for the case of B-type earthquake that associated with negative *L*-value that is smaller than a threshold. Fig. 15(a) shows the percentage versus the absolute value of the *L* threshold (note that threshold for B-type earthquake is negative). When the threshold is close to zero, which means that the sign of *L*-value classifies the event, we obtain 80.9 per cent successful classification for volcanic tremor and 74.7 per cent for B-type earthquake. The similarity with the JMA catalogue is decreasing with increasing absolute threshold. We also compute the percentage of the time windows in which a specified *L* threshold is satisfied but no event in the JMA catalogue (Fig. 15b). The percentage is rapidly decreasing with increasing absolute *L* threshold. For the previously determined *L* > 0.55 and
<figure><img src="image_18.png" /><figcaption>**Figure 15.** (a) Similarity of classification with the JMA catalogue and (b) classification of time windows that do not contain events in JMA catalogue. Values are in percentage with respect to the number of time windows containing events in JMA catalogue and presented as a function of the absolute *L* threshold. For an absolute threshold *x*, the true criteria are *> x* for volcanic tremor and *< -x* for B-type earthquake.</figcaption></figure>

$L < -0.70$, we obtain 67.5 per cent more volcanic tremor and 92.5 per cent more B-type earthquake than the JMA catalogue. Future investigations using seismic stations deployed very close to the active craters may clarify whether or not these newly detected events are really volcanic tremors and B-type earthquakes.

## 5.2 Comparison with volcanic activity

In Fig. 16(a), we show all source locations from 2017 April to September that satisfy the reliability criteria defined in Section 5.1. Most of the sources are horizontally distributed in SW-NE direction from the south of Minamidake to Showa crater, and vertically at depths from 4 km to the ground surface. B-type earthquake sources are mainly distributed at shallow depths of <2 km beneath the craters, in a similar region with that of volcanic tremors. A similar source distribution for volcanic tremor and B-type earthquakes confirms that both types of events are generated in the conduit beneath the craters. In our analysis, source locations at 7 km depth or more may not be reliable. This depth limit of 7 km is obtained from a test described in Supporting Information S3 using the CCFs at different depths [(Figs S3.1–S3.3)].

As shown in Fig. 16(b), previous studies suggested that the magma is supplied from a reservoir at 11 km depth beneath Aira caldera to the shallower one at 5 km depth beneath Sakurajima (e.g. Eto & Nakamura 1986; Yokoyama 1986; Kamo 1989). The depths of 0–2 km beneath the craters have been suggested as a volcanic conduit from the hypocentres of earthquakes associated with eruptions such as BL-type and explosion earthquakes (Ishihara 1990;

Iguchi 1994), which are determined from P-wave onsets. Note that the source locations of explosion earthquakes are not well determined using the CCF-based SSA method, even when their SNRs are high. Our estimated sources are not consistent with the hypocentres beneath the craters which were determined from arrival times of P waves (Ishihara 1990; Tameguri et al. 2001). This is possibly because their waveforms are complex. Their generation mechanism has been explained by expansion and contraction in the conduit followed by the generation of air-shock just beneath the crater (Tameguri et al. 2002). In addition, their typical duration is about 1 min, shorter than most short-duration tremors. Also, there is no case in our data where several explosions are observed in a single time window, in contrast to B-type earthquakes.

A small magma pocket at 2 km depth beneath the summit is inferred from the study of the seismicity gap of VTs (Ishihara 1988). Such seismicity gap is not observed in our results. This is probably because our estimated tremor sources are associated with magma and/or gas movement from the shallow reservoir to the craters through this magma pocket, and/or because the observation periods are different.

The regional stress at Aira Caldera and Sakurajima has been described as dominated by extension in WNW-ESE direction (Kimura 1985). Hidayati et al. (2007) determined a tensile fault that extends in SW-NE direction at depths of 7–9 km from ground deformation data. This tensile fault elongates horizontally beneath the Showa crater and ends at the SW of Sakurajima where VTs are observed at 7–8 km depths. The tremor source distribution is well consistent with the tensile fault, although the fault is located much deeper.

## 5.3 Application and computational cost

The analyses in this study can be divided into two major steps: analysis of known seismic events to obtain the weight function for each type of events; and the detection, classification and location of continuous data. A considerably large number of known events may be needed to extract their spectral width characteristics. Different types of events need to show spectral width characteristics that are distinct to each other, otherwise the classification will be difficult. At Sakurajima, we use the narrow frequency bands showing high spatial coherence to define the weight function. We noted that such narrow frequency bands may not be observed in some volcanoes. In the cases where different event characteristics are observed, some modifications on the details of computation shown in Fig. 5 may be necessary.

In this study, we perform a ‘thorough’ analysis that requires extensive computations by applying detection, classification and location method to each 10-min window. We measure the computation time for a 10-min window using a processor with four cores at 3.6 GHz speed and multithreading capability. It takes less than 5 min to compute the spectral width, *L*-value, correlation coefficient and the CCFs. On the other hand, the grid search of the location method requires more computational time, especially for higher spatial resolution. We may only focus on the ‘reliable’ results by first defining some reliability criteria based on the threshold of the averaged spectral width, correlation coefficient and *L* values for volcanic tremor and B-type earthquake. Only the time windows that satisfy those criteria may proceed with source location determination. A two-step grid search approach can be performed faster without reducing the spatial resolution: (1) we obtain an initial source location in a 3-D grid with a larger grid interval; (2) we refine the location by using a smaller subgrid around the initial source
<figure><img src="image_19.png" /><figcaption>**Figure 16.** (a) Source locations of volcanic tremors and B-type earthquakes that satisfy the criteria described in Section 5.1. (b) Schematic illustration of the volcanic system from previous studies overlaid to the source locations in (a). The generation depth for explosion (Ex), B-type (B) and A-type earthquakes (A) is based on Iguchi (1994). The location of magma pocket is inferred from Ishihara (1988). A magma chamber at 5 km depth is proposed by Eto & Nakamura (1986). The tensile fault model is proposed by Hidayati et al. (2007).</figcaption></figure>

location and a smaller grid interval. We can also combine the grid search with another method such as the gradient method (Tarantola 1987).

# 6. CONCLUSIONS

We developed a method for the classification of volcanic tremors and earthquakes based on correlation analyses of seismic records. We extended a detection method that measured the spatial coherence and calculated the distribution width of the eigenvalues extracted from the seismic covariance matrix called the spectral width. We analysed continuous data during 2017 April–September from six seismic stations at Sakurajima volcano, Japan. We extracted unique frequencies of high spatial coherence that define the weight functions to classify volcanic tremors and B-type earthquakes. The results show that about 80 per cent volcanic tremors and 75 per cent B-type earthquakes reported in JMA catalogues are successfully classified by our method. The sources of the classified seismic events were located and evaluated using the cross-correlation-based methods. Most sources are estimated between 4 km depth and the ground surface, from the south of the active craters to beneath the craters. B-type earthquakes are mostly located at less than 2 km depth at the same region as the tremor sources. The classification result is consistent with the reported volcanic activities, and the distribution of estimated source locations is consistent with generation depth of volcanic earthquakes and magma supply system proposed by various studies. Our technique can be used as an automatic tool for the detection and location of volcanic tremors and earthquakes at active volcanoes.

# ACKNOWLEDGEMENTS

Permana was supported by the International Joint Graduate Program in Earth and Environmental Sciences (GP-EES) of Tohoku University. This study was supported by the Ministry of Education, Culture, Sports, and Technology (MEXT) of Japan, under its The Second Earthquake and Volcano Hazards Observation and Research Program (Earthquake and Volcano Hazard Reduction Research). We thank JMA for providing various data at Sakurajima. Ryohei Kawaguchi kindly supported the use of seismic activity at Sakurajima. The work of Shapiro was supported by the European Research Council (ERC) under the European Union Horizon 2020 Research and Innovation Programme (grant agreements 787399-SEISMAZE). We are grateful for the reviews and comments by two anonymous reviewers that helped us improve this paper. Extensive computations for this study were performed using a high-performance computing infrastructure at the Institut des Sciences de la Terre (ISTerre), Grenoble. There is no conflicting interest regarding the content and publication of this paper. Permana, Nishimura and Shapiro designed the study and analysis method. Permana performed the computations with supervision from Nishimura, Nakahara and Shapiro. Permana prepared the manuscript with inputs from all co-authors.

# DATA AVAILABILITY

Seismic and infrasound data from the JMA are available through the Data Management Center of the National Research Institute for Earth Science and Disaster Resilience (NIED) at http://www.hinet.bosai.go.jp. Data of eruption activity at Sakurajima, monitored by
the JMA, can be accessed at http://www.jma-net.go.jp/kagoshim a/vol/kazan_top.html. Summarized reports on volcanic activity by the JMA are publicly available at https://www.data.jma.go.jp/svd/v ois/data/tokyo/STOCK/monthly_v-act_doc/monthly_vact.php. Cat- alogue of seismic activity at Sakurajima is not publicly available.

REFERENCES

Aki, K., Fehler, M. & Das, S., 1977. Source mechanism of volcanic tremor: fluid-driven crack models and their application to the 1963 Kilauea erup- tion, *J. Volc. Geotherm. Res.*, **2**, 259–287.

Allen, R.V., 1978. Automatic earthquake recognition and timing from single traces, *Bull. seism. Soc. Am.*, **68**, 1521–1532.

Bensen, G.D., Ritzwoller, M.H., Barmin, M.P., Levshin, A.L., Lin, F., Moschetti, M.P., Shapiro, N.M. & Yang, Y., 2007. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements, *Geophys. J. Int.*, **169**(3), 1239–1260.

Chouet, B., 1988. Resonance of a fluid-driven crack: radiation properties and implications for the source of long-period events and harmonic tremor, *J. geophys. Res.*, **93**, 4375–4400.

Chouet, B., Koyanagi, R. & Aki, K., 1987. The origin of volcanic tremors in Hawaii: part II, *U.S. Geol. Surv. Prof. Pap.*, **2**(1350), 1259–1280.

Droznin, D.V., Shapiro, N.M., Droznina, S.Ya., Senyukov, S.L., Chebrov, V.N. & Gordeev, E.I., 2015. Detecting and locating volcanic tremors on the Klyuchevskoy group of volcanoes (Kamchatka) based on correlations of continuous seismic records, *Geophys. J. Int.*, **203**(2), 1001–1010.

Eto, T. & Nakamura, S., 1986. Ground deformation around Sakurajima volcano—1974–1982—, in *5th Joint Observation of Sakurajima Volcano*, pp. 11–21, ed. Kamo, K., Sakurajima Volcanological Observatory, in Japanese.

Fehler, M.C., 1983. Observation of volcanic tremor at Mount St. Helens volcano, *J. geophys. Res.*, **88**(B4), 3476–3484.

Hidayati, S., Ishihara, K. & Iguchi, M., 2007. Volcano-tectonic earthquakes during the stage of magma accumulation at the Aira Caldera, Southern Kyushu, Japan, *Bull. Volcanol. Soc. Japan*, **52**(6), 289–309.

Honda, S. & Yomogida, K., 1993. Periodic magma movement in the conduit with a barrier: a model for the volcanic tremor, *Geophys. Res. Lett.*, **20**(2), 229–232.

Hurst, A.W. & Sherburn, S., 1993. Volcanic tremor at Ruapehu: charac- teristics and implications for the resonant source, *New Zealand J. Geol. Geophys.*, **36**(4), 475–485.

Iguchi, M., 1989. Distribution of the initial motions of volcanic mi- croearthquakes (B-type) at Sakurajima volcano, *Annu. Disas. Prev. Res. Inst.*, *Kyoto Univ.*, **32B-1**, 13–22 (in Japanese with English abstract).

Iguchi, M., 1994. *A Vertical Expansion Source Model for the Mechanisms of Earthquakes Originated in the Magma Conduit of an Andesitic Volcano: Sakurajima, Japan*, PhD thesis, Kyoto University, Kyoto.

Ishihara, K., 1988. Geophysical evidence on the existence of magma reser- voir and conduit at Sakurajima volcano, Japan, *Ann. Disast. Prev. Res. Inst.*, *Kyoto Univ.*, **31B1**, 1–15 (in Japanese with English abstract).

Ishihara, K., 1990. Pressure sources and induced ground deformation asso- ciated with explosive eruptions at andesitic volcano: sakurajima volcano, Japan, in *Magma Transport and Storage*, pp. 335–356, ed Ryan, M.P., John Wiley & Sons.

JMA, 2017. *Volcanic activity of Sakurajima volcano (2017)*, annual report, Japan Meteorological Agency, Tokyo, in Japanese.

JMA, 2018. *Volcanic activity of Sakurajima volcano (2018)*, annual report, Japan Meteorological Agency, Tokyo, in Japanese.

Journeau, C., Shapiro, N.M., Seydoux, L., Soubestre, J., Ferrazzini, V. & Peltier, A., 2020. Detection, classification, and location of seismovolcanic signals with multicomponent seismic data: example from the Piton de la Fournaise volcano (La Réunion, France), *J. geophys. Res. Solid Earth*, **125**(8), e2019JB019333, doi: 10.1029/2019JB019333.

Kamo, K., 1989. A dialogue with Sakurajima volcano, in *Proc. Kagoshima International Conference on Volcanoes*, pp. 3–13, Nippon Institute for Research Advancement.

Kimura, M., 1985. Back-arc rifting in the Okinawa Trough, *Mar. Pet. Geol.*, **2**, 222–240.

Kobayashi, T., Miki, D., Sasaki, H., Iguchi, M., Yamamoto, T. & Uto, K., 2013. *Geological map of Sakurajima volcano*, 2nd edn, Geological Survey of Japan.

Kurokawa, A., Takeo, M. & Kurita, K., 2016. Two types of volcanic tremor changed with eruption style during 1986 Izu-Oshima eruption, *J. geophys. Res. Solid Earth*, **121**(4), 2727–2736.

Leet, R.C., 1988. Saturated and subcooled hydrothermal boiling in ground- water flow channels as a source of volcanic tremor, *J. geophys. Res.*, **93**, 4835–4849.

Li, K.L. & Gudmundsson, O., 2020. A probabilistic tremor lo- cation method, *Geophys. Res. Lett.*, **47**(4), e2019GL085538, doi: 10.1029/2019GL085538.

Li, K.L., Sadeghisorkhani, H., Sgattoni, G., Gudmundsson, O. & Roberts, R.G., 2017a. Locating tremor using stacked products of correlations, *Geophys. Res. Lett.*, **44**, 3156–3164.

Li, K.L., Sgattoni, G., Sadeghisorkhani, H., Roberts, R. & Gudmundsson, O., 2017b. A double-correlation tremor-location method, *Geophys. J. Int.*, **208**(2), 1231–1236.

Malfante, M., Mura, M.D., Mars, J.I., Metaxian, J.-P., Macedo, O. & Inza, A., 2018. Automatic classification of volcano seismic signatures, *J. geophys. Res. Solid Earth*, **123**, 10,645–10,658.

Maryanto, S., Iguchi, M. & Tameguri, T., 2008. Constraints on the source mechanism of harmonic tremors based on seismological, ground defor- mation, and visual observations at Sakurajima volcano, Japan, *J. Volc. Geotherm. Res.*, **170**, 198–217.

McNutt, S.R., 1994. Volcanic tremor from around the world: 1992 update, *Acta Volcanol.*, **5**, 197–200.

McNutt, S.R. & Nishimura, T., 2008. Volcanic tremor during eruptions: temporal characteristics, scaling and constraints on conduit size and pro- cesses, *J. Volc. Geotherm. Res.*, **178**, 10–18.

Minakami, T., 1974. Seismology and volcanoes in Japan, in *Physical Vol- canology*, pp. 1–27, eds. Civetta, L., Gasparini, P., Luongo, G. & Rapolla, A., Elsevier.

Miyamachi, H. et al., 2013. Shallow velocity structure beneath the Aira Caldera and Sakurajima volcano as inferred from refraction analysis of the seismic experiment in 2008, *Bull. Volcanol. Soc. Japan*, **58**(1), 227– 237.

Nishimura, T. & Iguchi, M., 2011. *Volcanic Earthquakes and Tremor in Japan*, Kyoto University Press.

Nishimura, T. et al., 2013. Mechanism of small vulcanian eruptions at Suwanosejima volcano, Japan, as inferred from precursor inflations and tremor signals, *Bull. Volcanol.*, **75**(12), doi: 10.1007/s00445-013-0779-1.

Ogiso, M., Matsubayashi, H. & Yamamoto, T., 2015. Descent of tremor source locations before the 2014 phreatic eruption of Ontake volcano, Japan, *Earth Planets Space*, **67**(1), doi: 10.1186/s40623-015-0376-y.

Ono, K., Ito, K., Hasegawa, I., Ichikawa, K., Iizuka, S., Kakuta, T. & Suzuki, H., 1978. Explosion seismic studies in south Kyushu especially around the Sakurajima volcano, *J. Phys. Earth*, **26**(Suppl.), S309–S319.

Permana, T., Nishimura, T., Nakahara, H., Fujita, E. & Ueda, H., 2020. Reliability evaluation of volcanic tremor source location determination using cross-correlation functions, *Geophys. J. Int.*, **220**(2), 1300–1315.

Rost, S. & Thomas, C., 2002. Array seismology: methods and applications, *Rev. Geophys.*, **40**(3), 2–1-2-27.

Saragiotis, C.D., Hadjileontiadis, L.J. & Panas, S.M., 2002. PAI-S/K: a robust automatic seismic P phase arrival identification scheme, *IEEE Trans. Geosci. Remote Sens.*, **40**(6), 1395–1404.

Seidl, D., Schick, R. & Riuscetti, M., 1981. Volcanic tremors at Etna: a model for hydraulic origin, *Bull. Volcanol.*, **44**(1), 43–56.

Seydoux, L., Shapiro, N.M., de Rosny, J., Brenguier, F. & Landès, M., 2016. Detecting seismic activity with a covariance matrix analysis of data recorded on seismic arrays, *Geophys. J. Int.*, **204**(3), 1430–1442.

Shelly, D.R., Beroza, G.C. & Ide, S., 2007. Non-volcanic tremor and low- frequency earthquake swarms, *Nature*, **446**, 305–307.

Soubestre, J., Shapiro, N.M., Seydoux, L., de Rosny, J., Droznin, D.V., Droznina, S.Y., Senyukov, S.L. & Gordeev, E.I., 2018. Network-based detection and classification of seismovolcanic tremors: example from the Klyuchevskoy volcanic group in Kamchatka, *J. geophys. Res. Solid Earth*, **123**(1), 564–582.
Tameguri, T., Iguchi, M. & Ishihara, K., 2001. Reexamination of moment tensors for initial motion of explosion earthquakes using borehole seismograms at Sakurajima volcano, Japan, *Earth Planets Space*, **53**, 63–68.

Tameguri, T., Iguchi, M. & Ishihara, K., 2002. Mechanism of explosive eruptions from moment tensor analyses of explosion earthquakes at Sakurajima Volcano, Japan, *Bull. Volc. Soc. Japan*, **47**(4), 197–215.

Tang, C.-C., Peng, Z., Chao, K., Chen, C.-H. & Lin, C.-H., 2010. Detecting low-frequency earthquakes within non-volcanic tremor in southern Taiwan triggered by the 2005 Mw8.6 Nias earthquake, *Geophys. Res. Lett.*, **37**(16), L16307, doi: 10.1029/2010GL043918.

Tarantola, A., 1987. *Inverse Problem Theory*, 1st edn, Elsevier.

Titos, M., Bueno, A., Garcia, L. & Benitez, C., 2018. A Deep Neural Networks Approach to Automatic Recognition Systems for Volcano-Seismic Events, *IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens.*, **11**(5), 1533–1544.

Yokoyama, I., 1986. Crustal deformation caused by the 1914 eruption of Sakurajima volcano, Japan and its secular changes, *J. Volc. Geotherm. Res.*, **30**, 283–304.

## SUPPORTING INFORMATION

Supplementary data are available at GJI online.

**Figure S1.1.** Detection, classification and location during 2017 April.

**Figure S1.2.** Detection, classification and location during 2017 July.

**Figure S1.3.** Detection, classification and location during 2017 August.

**Figure S1.4.** Detection, classification and location during 2017 September.

**Figure S2.1.** Estimated source locations of continuous volcanic tremors during 2017 April–July. All source locations are plotted regardless of the averaged spectral width, correlation coefficient or L values.

**Figure S3.1.** Average CCF envelopes as a function of depth for station pair V.SKA2-V.SKR.C. (a) Location of areas A and B. (b) Average of CCF envelopes in 0.1 km depth bins beneath area A and (c) area B. Red points denote the peak of the envelopes.

**Figure S3.2.** Average CCF envelopes as a function of depth for station pair V.SKD2-V.SKR.C beneath (a) area A and (b) area B.

**Figure S3.3.** Average CCF envelopes as a function of depth for station pair V.SKD2-V.SKR.D beneath (a) area A and (b) area B.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

## APPENDIX A: NETWORK COVARIANCE MATRIX ANALYSIS

Following Seydoux *et al.* (2016), for a given subwindow *m*, we define the matrix of seismic spectra at frequency *f* as:

$$ \mathbf{u}_m(f) = \begin{bmatrix} u_{m1}(f) \\ u_{m2}(f) \\ \vdots \\ u_{mN}(f) \end{bmatrix}, \qquad (A1) $$

where $u_{mi}(f)$ is the complex spectrum of the seismogram $u_{mi}(t)$ at the *i*th station. The network covariance matrix is expressed as:

$$ \mathbf{V}(f) = \frac{1}{M} \sum_{m=1}^{M} \mathbf{u}_{m}(f) \mathbf{u}_{m}^{\dagger}(f), \qquad (A2) $$

<figure><img src="image_20.png" /><figcaption>**Figure A1.** Example of eigenvalue distribution and spectral width from two 10-min windows of (a) volcanic tremor and (b) ambient noise.</figcaption></figure>

where *M* is the number of subwindows in the 10-min time window. The matrix **V**(*f*) contains the average of the direct product between **u**<sub>*m*</sub>(*f*) and its Hermitian transpose **u**<sub>*m*</sub><sup>†</sup>(*f*). Because **V**(*f*) is an *N* × *N* Hermitian matrix, we can decompose it into its complex eigenvectors **v**<sub>*i*</sub> and the corresponding real positive eigenvalues $\lambda_i$:

$$ \mathbf{V}(f) = \sum_{i=1}^{N} \lambda_i(f) \mathbf{v}_i(f) \mathbf{v}_i^{\dagger}(f). \qquad (A3) $$

The eigenvalues $\lambda_i$ are sorted in descending order with the first eigenvalues being the largest. The spectral width is represented as:

$$ \sigma(f) = \frac{\sum_{i=1}^{N} (i-1) \lambda_i(f)}{\sum_{i=1}^{N} \lambda_i(f)}. \qquad (A4) $$

A large number of subwindows ensure a robust estimation of the spectral width. The seismic data are decimated from 100 to 25 Hz to shorten the computation time. Then, spectral whitening is applied. The spectral widths are estimated every 0.02 Hz.

For spatially coherent signals, which may be generated by volcanic tremor source, the first (largest) eigenvalue can be significantly dominant compared to the others. As a result, we observe a small $\sigma(f)$ (Fig. A1a). Contrary, larger $\sigma(f)$ is obtained for the data consisting of incoherent noises (Fig. A1b).

Since volcanic tremors and B-type earthquakes are abundant during the period of 2017 April–July, the spectral width is averaged for the frequency band of 1–4 Hz, where both types of events are known to be dominant. We remove low averaged spectral widths of <0.3 that are caused by missing data. Fig. A2 shows the averaged spectral widths during the period of 2017 April–July. We indicate the occurrence time of volcanic tremors and A-type, B-type, explosion earthquakes, short-duration tremors and the longer continuous tremors based on JMA catalogue.

During 2017 April–July, the estimated averaged spectral widths range from 0.5 to 1.9. The averaged spectral width is fluctuating when B-type earthquakes are observed, particularly in early April and throughout July. Such fluctuation is also observed when short-duration volcanic tremors occurred, for example, on May 13–23.
<figure><img src="image_21.png" /><figcaption>**Figure A2.** Averaged spectral width at 1–4 Hz computed every 10-min for 2017 (a) April, (b) May, (c) June and (d) July. The upper panels show the occurrence time of the A-type, B-type and explosion earthquakes, as well as short-duration (shown as ‘short’) and continuous tremors, according to JMA report. The averaged spectral widths are shown by black lines in the lower panel and also shown by colour scale to compare them with the seismic activity.</figcaption></figure>

Low averaged spectral widths of 1.1–1.3 are obtained over a period of several hours on May 7–12 and 28–30, and also for a period of more than two days on June 2–5. These time periods of low averaged spectral width in May and June are consistent with the occurrence time of continuous tremors. The averaged spectral widths of B-type earthquakes and short-duration tremors are generally larger than that of continuous tremor, as shown in the period of May 13–23.

A-type earthquakes show large averaged spectral width due to their higher frequency content. Explosion earthquakes also show large averaged spectral width which makes it difficult to detect them using spectral width analysis. The averaged spectral width tends to be lower during night time and higher during the daytime. This daily variation may be related to human activities, and its amplitudes are smaller than the changes caused by volcanic tremor or B-type earthquakes.

## APPENDIX B: APPLICATION OF THE CCF-BASED SSA

We prepare a 3-D velocity model by first combining 1-D *P*-wave velocity structures from two previous studies: *P*-wave velocity values at depths of <4 km are taken from a refraction seismic study at Sakurajima volcano and Aira Caldera (Miyamachi et al. 2013) and those at depths of >4 km from a seismic exploration using artificial sources around the Sakurajima volcano (Ono et al. 1978). Assuming that tremor and B-type earthquake signals are dominated by *S* wave, we calculate *S*-wave velocities using $V_p/V_s$ ratio of $\sqrt{3}$. A smooth 1-D structure is then generated by using a polynomial fitting of order 4. For each gridpoint, we assign the smoothed *S*-wave velocities by setting the topography as the uppermost value of the 1-D velocity structure. Fig. B1(a) shows the 1-D velocity structure combined from the two previous studies and the smoothed one. Fig. B1(b) shows a vertical cross-section of the 3-D velocity model.
<figure><img src="image_22.png" /><figcaption>**Figure B1.** (a) Grey line represents the 1-D velocity structure with combined velocity values from Miyamachi *et al.* (2013) and Ono *et al.* (1978). The black line is the smooth model used in the 3-D velocity model. (b) A vertical cross-section of the 3-D S-wave velocity model at the latitude of the maximum elevation. The velocity profile in (a) is extended in depth from the ground surface.</figcaption></figure>

We then compute the theoretical traveltimes from each station to all gridpoints through the 3-D velocity model.

Fig. B2 shows the estimated source locations in 2017 May, during which many continuous and short-period tremors are reported by JMA. The estimated sources are more concentrated in space over the time periods when continuous tremors occurred (Fig. B2a). Contrary, for the time periods of the occurrence of short-duration tremors and A-type, B-type and explosion earthquakes, the estimated sources are not well determined and more scattered in space. Similarly, the estimated sources are scattered for the time windows containing only incoherent seismic noise. As an attempt to distinguish the seismic events from noise, we examine the values of the brightness function (Permana *et al.* 2020). We calculate the brightness range at each time window as $B_{max} - B_{min}$, where $B_{max}$ is the maximum brightness and $B_{min}$ is the minimum one (Fig. B2d). The

brightness range is high not only when the source locations are more spatially constrained but also for the other time windows. Therefore, we may not be able to distinguish the noise using the brightness function. Higher brightness range during nighttime and lower values during daytime indicate the effects of human activity and artificial noise in our data.

## APPENDIX C: SELECTION OF TIME WINDOW LENGTH

We examine the effect of time window by calculating the averaged spectral width, the classification value *L* and source location of tremors as well as B-type earthquakes using different time window lengths: 10, 5, 2 and 1 min. The subwindow lengths and overlap amounts for the 5- to 1-min window are proportionally scaled-down from those used in our analysis using 10-min window, to ensure an equal number of subwindows. Fig. C1 shows the results on 6-hr data. During continuous tremor on August 23 (Fig. C1a), the averaged spectral width is consistently low for all time window lengths (Fig. C1b). When the tremors occur intermittently, the difference between averaged spectral width of tremor and that of noise becomes smaller for shorter time windows. The intercomponent correlation coefficients also show a similar effect, although it is not shown here. Fig. C1(c) shows that the vertical source locations are deteriorating with shorter time windows. Similar observation is also obtained during the occurrences of B-type earthquakes on June 16 (Figs C1d–f). When identifying short-duration or intermittent events, we may not observe stable averaged spectral widths or source locations over several time windows, in contrast with continuous tremor. Although it is possible to use shorter time windows, the computation time becomes longer and the detection of seismic events from noise needs to be carefully examined.
<figure><img src="image_23.png" /><figcaption>**Figure B2.** Source location every 10-min window during 2017 May. (a) Seismic activity by JMA, (b) latitude (in degrees), (c) longitude (in degrees) and (d) depth (in km) of the source locations. (e) Brightness range.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure C1.** Comparison of results using different time window lengths for (a)–(c) tremor on August 23 and (d)–(f) B-type earthquakes on June 16. (a) and (d) Filtered waveform at station V.SFT2, (b) and (e) average spectral width and (c) and (f) vertical location. All analyses are performed at 1–4 Hz.</figcaption></figure>