Abstract
We present the combination of permutation entropy (PE) and power spectral density (PSD) analysis on continuous seismic data recorded by short-period seismic stations during the 2010 Merapi volcano eruption. The calculation of PE aims at characterizing the randomness level in seismic noise, while the PSD parameters use to detect the background noise level in various frequency bands. It was previously observed that a significant reduction of randomness before the volcano eruption could be indicated as one of the short-term precursors due to the lack of high frequencies (>1 Hz) in the noise wave-field caused by high absorption losses as the hot magma uprises to the upper crust. The results show no significant reduction in signal randomness before the eruption series. The characteristic of events during the preeruptive period and the crisis tends to be chaotic (PE in the range 0.9 to 1). Further calculations show that the standard deviation in PE decreased in four days before the first eruption onset on 26 October. PE was stable at the highest values (very close to 1) and gradually returned to the previous fluctuation after the eruption onset. The level of background noise in the low- and high-frequency bands appeared to have the same tendency. The two main eruptions correspond to the two highest peaks of noise levels.
1. Introduction
Merapi volcano is located at the center of the Island of Java, Indonesia, about 30 km north of Yogyakarta. This volcano is categorized as one of the most dangerous volcanoes in the world because of the high level of eruption activities and the density of the surrounding population [1]. The last major eruption that occurred from October-November 2010 significantly impacted the social, environmental, health, and material aspects of the society around the volcano. This eruption is the first large explosive eruption (VEI ~4) observed by the modern multiparameter stations of the CVHGM (Center of Volcanology and Geological Hazard Mitigation), including six broadband seismic stations and four short-period seismic stations installed in the flank of Merapi (Figure 1). Merapi eruptions that occurred in the previous period with VEI 1 to 2 was happening in 1984, 1986, 1992, 1994, 1997, 1998, 2001, and 2006 were preceded by the growth of the lava dome (Merapi-type), while the appearance of the lava dome did not precede the October-November 2010 eruptions. In other words, the 2010 eruption is not the typical eruption type of Merapi, but a transition between effusive and explosive style [2].
One of the significant challenges for volcanologists is to predict the occurrence of major events during the volcanic unrest period. An essential part of the prediction effort is the determination of short-term precursors as an indicator. In Merapi, some precursors usually observed before the eruption, i.e., the increase of volcanic earthquakes (volcano-tectonic) events, a more concentrated smoke coming out from the crater, the change in the composition of gas and water, the increasing of temperature in the crater, and the deformation of the volcano edifice [3, 4].
In recent years, ambient seismic noise has been widely used along with the continuous observation of volcanic signals. For example, it can be used to monitor the physical property change of material below the volcano, and hence it can be used to predict the activity of a volcano [5–7]. Ambient seismic noise is characterized by the randomness of strength and position of the source. This noise also has characteristics related to the heterogeneity of the nature of the medium where the seismic signal is passed through [8, 9]. The variation in ambient seismic noise usually occurs earlier than other precursor signals, such as surface deformation or volcano-tectonic events due to subsurface rock fracture processes.
One new approach in the analysis of ambient seismic noise related to periods of volcanic unrest is the application of permutation entropy (PE) calculations. This method can be used to analyze the randomness of time-series signals based on the number of matches between the seismic signal pattern and permutation pattern on a particular embedding dimension [10]. A significant reduction in the temporal variation of PE approximately was observed ten days before the onset of the 1996 Gjálp eruption (Iceland) and indicated as a short-term precursor of volcanic eruption [11]. Further analysis with the calculation of the temporal variation in the dominant and centroid frequencies [12] showed that the reduction was due to the lack of high-frequency events (>1 Hz) associated with a high absorption losses due to the increase in magma moving up the earth’s surface.
PE calculations can be categorized into the “one-component single station” analysis method [13]. One of the advantages of this method is the minimum use of equipment. The calculation only requires at least one component (usually vertical) of seismic data from a station close to the summit of a volcano. Another method in the same category is power spectral density (PSD) calculations. PSD is the result of the Fourier transform of a signal autocorrelation [14]. This parameter can be used to show the intensity of the seismic signal in a specific frequency band. This method was previously used in the near-real-time analysis of several active volcanoes, including Villarrica volcano (Chile), Tungurahua volcano (Ecuador), and Teide volcano (Spain) [15]. The Teide volcano showed increased activity in early 2003 after 30 years of quiet periods [15, 16]. PSD calculation results show the intensity of the seismic signal in several frequency bands related to increasing of the volcano activity. The advantage of this method compared to the most popular way in volcanic monitoring, i.e., Real Seismic Amplitude Measurement or RSAM [17], is the ability to describe the temporal changes in the body of a volcano [18].
This paper presents the results of a combined analysis of the temporal variations of PE and PSD on seismic data recorded before and during the Merapi eruption crisis in October-November 2010. A comparison between these two methods with the RSAM method, along with the daily event calculations, is analyzed for the further possibility of the technique toward volcano events precursory signals in Merapi.
2. Materials and Methods
2.1. The 2010 Merapi Eruption
Merapi volcano overlies the Java subduction zone and is composed mainly of basaltic-andesite tephra, pyroclastic flow, lava, and lahar deposits. Eruptions during the twentieth century typically recurred every 4 to 6 years and produced viscous lava domes that collapsed to form pyroclastic flows and subsequent lahars. These eruptions were relatively small, with typical eruptive volumes of and magnitudes or volcanic explosivity index (VEI) of 1–2 [19–22], where magnitude is given by [].
2.1.1. Seismic Network
There were six broadband seismic stations (GRW, GMR, PAS, LBH, WOR, and CRM) and four short-period seismic stations (DEL, KLA, PLA, and PUS) that recorded activities before and during the October 2010 eruption. Broadband stations were equipped with Güralp CMG-40T/30s seismometer, while the short-period stations were equipped with L-22/2 Hz and the L-4C/1 Hz (Mark Products) seismometers. The availability of seismic data at each station is shown in Figure 2. During the active phase of the volcano, data transmission from most broadband stations on the network was disrupted, and there were a lot of data gaps as the consequences. The data also does not cover the eruption series. Fortunately, the short-period stations were relatively stable and recorded the eruption series events with a narrower frequency band. This paper shows the temporal variations of PE and PSD on seismic data records from PUS, DEL, and PLA short-period stations. The station selection based on the data completeness before and after the eruption is also the location of the station. The PUS station is a short distance from the crater, while the DEL and PLA stations are in the middle range (<6 km). The results are expected to represent the ambient noise characteristics in the other locations.
2.1.2. Pre- and Coeruptive Seismicity
The seismicity at Merapi volcano during the 2010 crisis shows that all types of events previously identified at Merapi [23] were represented in the 2010 activity: Volcano-Tectonic (VT) events, Low-Frequency events (LF), tremor, “Multiphase” events (MP), “guguran” = rock falls (RF), and Very-Long Period events (VLP). Early indications of increased seismic activity included four swarms of volcano-tectonic (VT) events on 31 October 2009, 6 December 2009, and 10 June 2010. The period from 20 September 2010 until the initial explosive eruption on 26 October 2010 was marked by a dramatic increase in all monitored parameters, including the seismicity rate [24]. The number of both volcano-tectonic (VT) events corresponding to shear fracturing in the edifice and multiphase events (MP, also called “hybrid” events) corresponding to magma movement increased exponentially in October 2010. Besides the sharp increase of VT and MP events, the number and magnitude of rockfalls (RF) also intensified before the eruption. From 1 to 18 October, more than 200 very long-period (VLP) signals were recorded at summit stations, with some large VLP events recorded at all broadband stations.
An explosive eruption began at 10 : 02 UTC on 26 October and ended at ~ midnight UTC. A period of relative calm ensued on 26–28 October and was followed by smaller explosive eruptions on 29 October (~17 : 10–19 : 00 UTC), 31 October (~7 : 30 and~8 : 15 UTC) and 1 November (~3 : 00 UTC). More than 150 large low-frequency (LF) events (with dominant frequencies ~2 Hz) occurred between 29 October and 3 November. On 3–4 November, the eruptive intensity increased again with a series of explosions. Early on 3 November, data from close-range seismic stations became saturated (>0.025 mm/s) due to increased intensity of tremor. On 4 November, the tremor increased again and was felt as Mercalli intensity 2–3 shaking at 10–20 km from the volcano. All four proximal real-time seismic stations were completely saturated (>0.025 mm/s). Seismic amplitudes from the distal Imogiri station at the time of the climactic explosion (4 November at ~17 : 01 UTC) were 5 times larger than signals associated with the 3 November explosion. The intermittent and sometimes sustained explosive eruptions during the night (local time) of 4–5 November, included the climactic eruption.
On 6 November, tremor amplitude decreased slowly in parallel with decreased explosive activity. Later on, 7–8 November, RSAM increased again and remained at relatively high levels for another 2 days, which prompted CVGHM to quickly rebuild parts of the seismic monitoring system that were heavily damaged during the climactic eruption on 5 November. Due to danger, new stations were temporarily set-up more than 10 km from the summit. The destroyed stations (PUS, DEL, and KLA) were rebuilt after the eruption ended. After 8 November, seismic activity (mainly tremor and some volcano-tectonic events probably associated with stress readjustment after the massive eruption) started to decrease in intensity slowly.
2.2. Methods
2.2.1. Permutation Entropy
Permutation entropy (PE) is a statistical tool to measure the dynamical complexity of a time series [10]. PE is maximum when the system is utterly chaotic and minimized when the order is stochastic up, or the characteristics changed to a more deterministic nature [11]. PE calculation encodes the time series into one symbol, which reflects the rank of the order of successive in sequences of length [25, 26]. PE is defined as follows.
The number of possible permutation patterns is . represent the relative frequencies of possible patterns of symbol sequences, termed permutations. is usually normalized to produce a value between . The common denominator for normalization is or () [27]. An example of the permutation entropy approach after normalization () written in Equation (2).
The common approach to handle equal values in a time series which might, e.g., be caused by low digitalization or oversampling is to lengthen the time interval between the values of considered sequence and the step size by means of and , respectively. The term is known from the phase space reconstruction employing time delay embedding where is the time delay, and is the embedding dimension.
2.2.2. Power Spectral Density.
The power spectral density (PSD) or is the Fourier transform of the autocorrelation function >, as shown in Equation (3) [14].
The < > symbol indicates the average over the time . Depending on whether is a displacement (), velocity (), or acceleration () record. is given in unit of m2/Hz, (m/s)2/Hz, or (m/s2)2/Hz for , , and . Power spectral density can also be expressed in units of dB, such as power spectral density for acceleration as in Equation (4).
In the PSD calculation, the seismic data in specific segment durations divided into many subsegments. For example, the seismic data with a one-hour duration and 40 Hz sampling rate can be divided into 13 subsegments. PSD is then calculated for each of these subsegments [28]. PSD can be displayed in time series for several frequency bands.
2.2.3. Data Processing
Processing begins with data corrections and conditioning. Data corrections include the detrend process to eliminate DC offset and the merger process to convert data that was initially in daily to a monthly format. Data conditioning consists of the downsampling process to change the number of data samples from 100 samples per second to 20 samples per second, resulting in a deterioration of frequency resolution above 5 Hz. For our study, we choose 10 minutes window length, embedding dimension , and time delay . These parameters are adapted from the study of the dynamical complexity of the 1996 Gjálp eruption [11]. The processing is then continued to the calculation of the temporal variation in permutation entropy and power spectral density. We also calculated the moving standard deviation (MOVSTD, ) of the PE values for further analysis. A time frame from 7 September to 11 November was considered suitable since it covered the period before and after the eruption with no data gaps. For processing the PSD (time and frequency domains), we used the ObsPy [29], the Python-based Seismology data analysis platform, and MATLAB 2018. The MATLAB is used for permutation entropy calculations with the same code as used in processing the data from Iceland [11].
3. Results
The temporal variations of PE and PSD for PUS, DEL, and PLA stations are shown in Figure 3. There was no significant reduction in signal randomness preceded the 2010 Merapi eruption as stated in the preliminary study [30]. The variations of PE are in the range of 0.9 to 1, which means that the character of the events in the preeruptive period and the crisis tends to be chaotic. PSD fluctuations in the low-frequency band (0.13-0.28 Hz and 0.36-0.7 Hz) almost coincide with the variation in the high-frequency band (0.71-1.43 Hz and 3.4-6.8 Hz). Two primary eruption onsets on 26 October 10 : 02 UTC and 4 November 17 : 01 UTC correspond to the two highest PSD peaks, both in the low and high-frequency bands. PSD fluctuations in the low-frequency band are more visible at DEL station. These fluctuations show that PUS station, which is closer to the summit than the rest of other stations, record more high-frequency events. Based on the PSD amplitudes, events with the frequency range 3.4-6.8 Hz appear to have the most dominant intensity.
Figure 4 shows a more detailed temporal variation of PE for PUS station. PE tended to increase towards the first eruption onset on 26 October. This increase was confirmed by the calculation result of moving standard deviation (MOVSTD, ) against PE values, which showed a decreasing trend for four days before the eruption onset. The decrease indicates the randomness of the signal in the time range was stable at the highest values. After the first eruption onset, PE values gradually returned to the original fluctuation pattern. The calculation of moving standard deviation also performs on the PE parameters of DEL and PLA stations, which is located further from the summit. Figure 5 shows the comparison of the standard deviation for three short-period stations (PUS, DEL, and PLA). The trend of decreasing standard deviation is also observed at DEL and PLA stations but not as clear as the decrease at PUS station.
We are using the same technique to calculate the temporal variations of PSD on the seismic data of HOT14 and HOT23 broadband stations that record the 1996 Gjálp eruption. Both stations were equipped with the broadband Seismometer Güralp CMG-3ESP/30s. We requested from IRIS (using the online BREQ_FAST request form) the 3-component data for the time period of 31 July 1996 to 11 October 1996. PSD fluctuations in the low-frequency band (0.13-0.28 Hz and 0.36-0.70 Hz) and high-frequency band (0.71-1.43 Hz and 3.4-6.80 Hz) can be clearly distinguished (Figure 6). Decreased PSD in the high-frequency band and increase in the low-frequency band started occurring eight days before the Mw ~5.6 earthquake or ten days before the onset of the subglacial eruption coincided with a significant decrease in the PE value. After the earthquake, the PSD value in the high-frequency band started to increase. The increase of PSD value in the high-frequency band could be associated with the crisis period of the Gjálp eruption. Low-frequency events (0.13-0.28 Hz) appear to be more dominant in the period of volcanic unrest. The results of moving standard deviation calculations against PE values for HOT14 and HOT23 stations are shown in Figure 7. The decrease in standard deviation before the eruption is seen in the results. Unlike the case with the Merapi eruption in 2010, the randomness of the signal about eight days before the Mw ~5.6 earthquake was stable at low values.
We use the data from short-period seismometers (Mark L-22/2 Hz and Mark L-4C/1 Hz) mainly due to data availability. The previous 1996 Gjálp eruption study [11] used data from Güralp CMG-3/30s broadband seismometer. The type of broadband seismometer used at Merapi (Güralp CMG-40T/30s) has a frequency response that is almost the same as the type of seismometer used in Iceland monitoring. The comparison of the frequency responses for four seismometers is shown in Figure 8. Both short-period seismometers have limitations in responding to low-frequency events compared to broadband seismometers. Short-period seismometers will attenuate signal components with low-frequency bands (<1 Hz) at certain amplitude limits. In comparison, Güralp CMG-3ESP and Güralp CMG-40T seismometers can still respond well to signal components with frequencies of up to 0.03 Hz. Figure 9 shows the comparison of background seismic noise in Iceland before the 1996 Gjálp eruption and Merapi before the 2010 eruption. In this comparison, we use data from two broadband stations (LBH and PAS) in June to calculate background seismic noise at Merapi. In general, the area in the flank of Merapi is relatively noisier compared to the Icelandic region.
Figure 10 shows the comparison of the analysis results using RSAM, temporal variation of PSD, and the daily number of events calculation of volcano-tectonics (VT) A+B, multiphase (MP), low-frequency (LF), and rockfalls (RF) events during the Merapi eruption in 2010. The RSAM result and the daily event observation were carried out by BPPTKG [2]. The PSD calculation looks more able to show fluctuations in volcanic activity mainly related to the amount of VT, MP, and RF events compared to RSAM calculation.
4. Discussion
The seismic activity of Merapi is characterized by a large variety of events that correspond to different locations and physical processes of the sources [7]. The seismicity at Merapi volcano during the 2010 crisis shows that all types of events previously identified at Merapi [23]. A large number and magnitude of events accompanied the eruption, including VT, MP, episodes of tremor, as well as LF and VLP events. The period from 20 September 2010 until the initial explosive eruption on 26 October 2010 was marked by a dramatic increase in seismicity rate [24]. The number of both volcano-tectonic (VT) events corresponding to shear fracturing in the edifice and multiphase events (MP, also called “hybrid” events) corresponding to magma movement increased exponentially. Besides the sharp increase of VT and MP events, the number and magnitude of rockfalls (RF) also intensified before the eruption. From 1 to 18 October, more than 200 very long-period (VLP) signals were recorded at summit stations, with some large VLP events recorded at all broadband stations. These seismic data indicate the transport of large volumes of magma and fluids. Volcano-Tectonics (VT) events, volcanic tremors, and surface processes, e.g., pyroclastic flows, lahars (volcanic debris flows), and rockfalls (RF) are associated with high-frequency content (>1 Hz) [13].
The temporal variation of PSD in several frequency bands can describe the intensity of events that occurred in the 2010 Merapi eruption. The dominance of events associated with high-frequency content (VT, MP, and RF) is estimated to be the cause of PE values tending to be high (>0.9). During the four days leading up to the first eruption on October 26, PE was stable at the highest values. This increase expected to be associated with an exponential increase leading up to the first eruption onset. The fluctuations in PSD values at the PAS station illustrate the reality of events associated with high-frequency content better because of its location close to the summit of Merapi. PSD calculations can also provide an overview of the intensity of events that occurred in the initial explosive phase (26 October-1 November) and the magmatic phase (1 November-7 November) [1].
The results of polarization analysis in the previous 1996 Gjálp eruption study determined that a very long period of tremor (<1 Hz) was not the cause of the significant reduction in signal randomness [31]. The reduction was most likely the result of decreased scattering caused by the strongly pressurized crust due to doming in the upper mantle. CF coincides with DF as the crustal heterogeneity was changed, and the higher frequency scattering was decreased. Therefore, PE variations in the ambient noise, for any given volcanic regime, are estimated likely localized to the areas undergoing strong pressurization, which alters the scattering properties of the crust. Monitoring and petrologic data show that the 2010 eruption was fed by the rapid ascent of magma from depths ranging from 5 to 30 km. Magma reached the surface with variable gas content resulting in alternating explosive and rapid effusive eruptions and released a total of ~0.44Tg of SO2 [1]. It is also possible to be another cause of the absence of a significant decrease in PE before the eruption.
Background seismic noise in an area will also affect the results of PE calculations. The natural sources and anthropogenic sources (also known as cultural noise) responsible for the generation of ambient seismic noise differ in frequency content. This difference made way for distinguishing natural sources (relatively low frequency) as ambient seismic noise (microseisms) and anthropogenic sources (relatively high frequency) as microtremors [32]. Ambient seismic noise and microtremors can overlap, resulting in the seismic noise signal becoming highly contaminated from cultural activities. The contribution of cultural noise (>1 Hz) is a factor that needs to be considered in processing ambient seismic noise at Merapi due to the density of the surrounding population. This condition is different from the conditions in Iceland. The contributions from cultural noise sources in the central part of Iceland are expected to be almost zero since the area (including Vatnajokull Glacier) is uninhabited [11]. PE captured changes in the stochasticity of the noise wavefield (at stations located up to 120 km away) 8.57 days before the Bárðarbunga earthquake and 10.76 days before the onset of the eruption [31].
Some external factors that also influence the results of PE calculations are the type of seismometer used and the completeness of the data. Short-period seismometers have limitations in responding to low-frequency events compared to broadband seismometers. Short-period seismometers will attenuate signal components with low-frequency bands (<1 Hz) at certain amplitude limits. The results of PE calculations that reflect the system stochasticity may not experience a significant decrease if events with low-frequency content dominate the volcanic unrest period. Data availability and minimum data gap also have an impact on calculations. A large number of data gaps often becomes a problem in programming.
The types of seismic phenomena associated with the eruption in Iceland are not as many as Merapi. Some of these events are Low- Frequency events, mixed frequency events, and volcanic tremor [33]. The tremor was interpreted as evidence of magma migration from a deep source to a shallow magma chamber underneath the Bárdarbunga volcano. The possible eruption mechanism in Iceland’s volcanic system is also different from the eruption mechanism at Merapi. The possible mechanism for the generation of earthquakes on a ring fault underlying the caldera is the inflation that occurs in a shallow magma chamber caused by magma injected from a deeper source. The shallow chamber experiences higher pressure leading to increased stresses on the rocks underneath the chamber, which could eventually initiate rupture along the ring fault [34]. Events with low-frequency content appear to be more dominant in the 1996 Gjálp eruption.
5. Conclusions
A large number and magnitude of events that accompanied the 2010 Merapi eruption are estimated to be the cause of PE values tended to be high. The highest PE values that were stable during the four days leading up to the first eruption on October 26 are estimated to be associated with an exponential increase of all types of events. The temporal variation of PSD can describe the intensity of events prior to the eruption and the dominance of the events related to high-frequency content. The rapid ascent of magma is also possible to be another cause of the absence of a significant reduction in signal randomness before the eruption because the strongly pressurized crust due to doming in the upper mantle has not yet dominated the system so that the crustal heterogeneity has not changed, and the higher frequency scattering has not decreased. The contribution of cultural noise is a factor that needs to be considered in processing ambient seismic noise at Merapi due to the density of the surrounding population. Some external factors that also influence the results of PE calculations are the type of seismometer used and the completeness of the data. Types of seismic phenomena and the possible mechanism associated with the previous eruption can be considered in the application of PE calculations as an analytical method.
Data Availability
The data is maintained by the Center of Volcanology and Geological Hazard, Bandung Indonesia.
Conflicts of Interest
The authors declare no conflict of interest.
Authors’ Contributions
A.R. wrote the original draft preparation; A.R., A.B.S., and W.S. wrote review and editing; W., A.B.S., and H.H. conceptualized; W. did the supervision.
Acknowledgments
We would like to thank the Ministry of Research and Technology-National Research and Innovation Agency of Indonesia for a doctoral scholarship granted to A.R. We thank Balai Penyelidikan dan Pengembangan Teknologi Kebencanaan Geologi (BPPTKG), Yogyakarta, for providing data needed in the research. We wish to thank Prof. Kirbani S.B. (Rahimahullah) for many helpful discussions and for suggesting several improvements. We are grateful to Moch. Nukman who also read the manuscript and helped us very much with his comments and suggestions to improve the paper. Special thanks to Ade Anggraini and Kuwat Triayana for fruitful discussions. We acknowledge the staff of the Physics Dept. Universitas Gadjah Mada for their support. Department of Physics through Hibah Dosen Fisika UGM.