Permutation entropy variation for 2007 effusive dome-forming eruption period of Kelud Volcano, Indonesia

The complexity of a system recorded in time series data can be measured statistically using permutation entropy (PE). The state of a system (e.g. regular, chaotic, random, etc.) that underlies the appearance of variations in time series can be determined with PE. Since volcanoes are considered as the complex dynamical system controlled by interactions of many processes. Permutation entropy can be applied to study the system mechanism of volcano. We utilized PE to study system mechanism of Kelud volcano in 2007 dome-forming eruption period, from 3 (KWH; KLD; UMBK) seismic stations with different distances from the crater lake. Then, we want to compare the results. The result of study shows that the PE pattern for each station is different. The unique PE pattern that can be used as an eruption precursor is only shown at KWH and KLD stations. This pattern began to appear 2.7 days before the eruption on 3 November 2007. Data from UMBK station doesn’t show unique PE pattern. The factors such as sensor distance from magmatic activity center, size, and type of eruption probably influenced the final PE result. Using PE as the addition to volcano monitoring can maximize efforts in mitigation activities.


Introduction
One of the Quaternary active and hazardous volcanoes in Java Island, Indonesia is Kelud volcano (1731 m above sea level) [1]. The volcano has a crater lake that contains water at its summit. The tectonic setting of the volcano is the subduction zone between the Indo-Australia plate beneath the Eurasian plate [2]. Position of Kelud volcano is surrounded by Kawi and Butak volcanoes situated at the east of Kelud volcano, and also Arjuno-Welirang volcano that located at the northeast of Kelud volcano (Fig. 1).
Based on Center for Volcanology and Geological Hazard Mitigation (CVHGM) records, Kelud volcano has history of explosive and effusive eruptions. Explosive eruption of Kelud volcano usually was shortlived period and initiated by phreatomagmatic eruption, then followed by explosive magmatic eruptions [3]. According to eruptive records after 1000, explosive eruptions of the volcano have occurred in 1586; 1901; 1919; 1951; 1966; 1990; 2014 with explosive events of VEI between 3 to 5 [3]. Meanwhile for effusive eruptions of the volcano formed a lava dome. Effusive dome-forming eruption events have occurred in 1376; 1920; 2007 [3]. The VEI scale of effusive eruption was not bigger compared to explosive eruption.
Volcanoes are complex dynamical system [4]. The system is controlled by interaction of various processes, Seismic waves like ambient noise seismic which propagates through the system of volcano can be considered as complex time series [5]. The randomness amount contained in complex time series can be measured statistically by permutation entropy (PE). Glynn and Konstantinou (2016) introduced PE as one of methods to approach volcanic eruption forecasting. They applied PE to forecast the 1996 Gjálp fissure eruption at Iceland. Position of Gjálp fissure is between Grimsvotn volcano and Bardarbunga volcano with Mid Atlantic Ridge as its tectonic setting. The result of their study represents similar temporal pattern PE variation from 3 different seismic stations, namely 7 days before the eruption, PE reached minimum value [6]. This pattern can be considered as the precursor signal of the 1996 Gjálp eruption.
In this study, we applied PE on the 2007 domeforming eruption event of Kelud volcano. We used data from 3 seismic stations with different distance from the crater lake. The aims of the study are to know the PE patterns of Kelud volcano especially for the 2007 effusive eruption event from seismic stations with different distance and to compare them. Furthermore, we want to know the effectiveness of PE for forecast effusive eruption at Kelud volcano.

Sequence of 2007 eruption event
The information on sequence of 2007 eruption event is summarized from CVHGM records in Hidayati, et.al. According to CVHGM records, intense degassing at the lake floor in July 2007 was observed as the first sign of volcanic activity increased after 17 years in the rest period [3]. After that in August, the sign was followed by an increase of the crater lake water temperature and CO2 concentration. The colour change of crater lake water also recorded from typical green to yellowish in August [1,3,7]. Meanwhile from deformation monitoring showed that continuous uplift occurred at the crater rim from July until early September [7].
The increase of seismic activity was marked by first volcano-tectonic (VT) swarm that occurred on 9 September. CVHGM raised the level of Kelud activity from level I (normal) to level II (watch) [3]. Two days later on 11 September, VT earthquakes reached 13 events in 5 hours with hypocenters distribution at 2 -4 km below the crater [1]. After that day, VT earthquake event occurred on average one each day. A remarkable VT earthquake occurred on 26 to 29 September. In this period, second VT swarm that contained larger amplitude occurred and amount of VT reached 61 events with no significant change of hypocenters distribution [1]. These events made CVHGM raised volcanic activity to level III (standby) [3]. The occurrence of VT earthquakes decreased after 29 September, with VT hypocenters distribution shifted to 0.5 -1 km below the crater [1]. Low frequency (LF) earthquakes began to occur as VT earthquakes decreased. The occurrence of LF earthquake continued until 13 October. The first volcano tectonic type B (VB) swarm that contained 306 events occurred on 16 October for 7 hours [1]. This caused CVHGM raised volcanic activity again to level IV (warning) [1,3]. The VB earthquake continued until 21 October with decreasing in number. The event was followed by seismicity quiescence that occurred from 22 to 23 October [1].
According to the record on 24 October, the occurrence of VT earthquakes was followed by occurrence of VB earthquakes [1]. Then on 31 October, VT swarm with larger amplitude reached 31 events occurred, and followed by 89 VB earthquake events. On the following day recorded 55 VT earthquake events were followed by 1437 VB events. Furthermore, on 1 November, there was also tremor earthquakes with 274 events (marked as line a in Fig. 3, Fig. 4, Fig 5) [1]. After that on 2 November, seismicity activity was signed by the occurrence of VT that reached 52 events and VB earthquakes that attained 690 events [1]. That event was followed by 8 hours of quiescent period [7]. The quiescence was ended by continuous tremor that began at 11:07 local time (marked as line b in Fig. 3, Fig. 4, Fig 5) [1]. The amplitude of continuous tremor was getting larger and at 14:43 continuous tremor was replaced by spasmodic tremor for 70 minutes (marked as line c in Fig. 3, Fig. 4, Fig 5) [3,7]. The occurrences of tremor were interspersed by emission earthquakes. On 19:47 local time, VT earthquake was occurred again with the number and amplitude increasingly until 3 November 06:31 local time [3]. On the same day at 06:00, the temperature measurement of water crater lake showed the temperature has increased quite significantly [3].
On 3 November at 12:32 from seismogram showed continuous LF tremor began as rapid series of shallow LF and the event was followed by spasmodic LF tremor at 13:11 [3]. At 15:55 (marked as line d in Fig. 3, Fig. 4, Fig 5), the amplitude was getting larger and became saturated [1,3,7]. At the same time, there was sharp increase of temperature and sudden deflation of radial tilt. This moment was estimated as the time where extrusion lava occurred [3]. Unfortunately, the weather was heavy rain so it was too dangerous for field observation. A small fuming black dot in crater lake was observed from CCTV at 4 November [1, 3, 7].

Methodology
The activity of Kelud volcano is monitored by Kelud Volcano Observatory and operated by CVHGM. Since April 2007, the seismic activities are monitored by 5 seismic stations and deformation of the volcano is recorded by 2 tiltmeter stations. Each seismic station is equipped by L4-C vertical short-period seismometer [1]. The location of seismic stations is showed at Fig. 2. The format data for input in PE calculation is MSEED. Because the original format data is WIN, so it must be converted to MSEED. After that, we cut the time series with window length of 10 minutes to generate several time series. These data were used as input in PE calculation.

Calculation of permutation entropy
Permutation entropy (PE) was introduced by Bandt and Pompe on 2002 as statistical measurement method of complexity amount contained in the time series data [8]. This method can also be used to determine the state of system that underlies the appearance of variations in time series [9]. There are 2 important parameters in PE calculation, namely the embedding dimension ( ) and the delay time ( ). Embedding dimension describes the amount of information contained in each vector. Meanwhile, represents delay time on successive of vectors.
PE transforms time series data into graph of permutation entropy value. A time series is arranged become several vectors based on value input in and parameters. Then the values contained in vector are mapped and substituted with unique permutation pattern π (0,1,.., -1). Variation amount of the pattern are determined by factorial ( !). After that, the probability distribution of each pattern is calculated as the input value in PE Equation (1) [9].
Normalizing PE( ) by ( !) we obtain normalized quantity with interval value from 0 to 1. The value of PE normalization can describe how far the data from maximum complexity (PE = 1) or how close the data from minimum complexity (PE = 0). The zero value of PE indicates time series data are decreasing or increasing monotonously [6]. This condition shows the presences of regular system mechanism.
The values of PE depend on the combination of chosen parameter and [9]. For practical purposes, Bandt and Pompe recommend = 3,…,7 [8]. Riedl, et.al. in 2013 summarized several studies of PE on various fields showing that most of them used = 3-8, from tens to hundreds for time series from measurements of dynamical processes being observed at a fixed rate per unit time, = 1-5 for time series from measurements of processes that possess an inherent cycle represented by triggering events [10]. Previous studies such as [6] used = 5 and = 5 for monitor the 1996 Gjálp eruption; the same values for m and L were adapted by [11] to study the 2010 Merapi eruption; [5] applied = 5,6 and = 1,2,3,4,5 for monitor the January 2011 Shinmoedake eruptions sequence; [12] utilized = 3,4,5 and = 4,5,6 for study the 2014 Kelud eruption. Thus, we used variation of = 3,5 and = 2,3,4,5 to perform permutation entropy in this study.

Results and discussion
Based on the PE results of KWH, KLD, and UMBK stations show varying result. The PE graph of KLD station shows clear pattern at combination of parameter = 5 and = 2 (Fig. 3). PE value from 270 Julian date or 27 September to 304 Julian date end presents the state of volcano system on the maximum complexity during that period. It can be seen from the high PE value, which is close to 1. The differences PE pattern trend began to appear on 305 Julian date or 1 November, namely 2.7 days before eruption. It indicates that there has been a change in the system of volcano. PE value began to decrease, this trend corresponds with recorded tremor events (marked as a in Fig. 3). The properties of tremor which have monotonous wave pattern and lasting for long-duration indicated that the mechanism in volcano system began well-ordered. Then, the quiescent period that lasted for 8 hours caused PE value increased again (marked as pink highlight). This condition stated system on maximum complexity for a moment.
A sudden drop in the PE value (marked as b in Fig.  3) associated with the continuous tremor event that ended the quiescent period. Vertical line c in PE graph associated with amplitude of tremor became larger and the VT event occurred again. The trend of PE graph was increasing slowly but the maximum value was not high as before. Then, a sudden drop of PE value (marked as d in Fig. 3) coincided with the time that was interpreted as the period of first lava extrusion. At that time, the results of deformation and temperature monitoring also showed the precursor of eruption.
The PE graph of KWH station shows unclear pattern at all combination of parameter and . So, we tried to filter time series before it was used again in PE calculation. We filtered the data using LPF in 3 stages. Firstly, we tried to apply LPF 3 Hz. Then, the filtered data were processed using PE but the PE graph still showed unclear pattern. So, we tried to filter data using LPF 2 Hz and processed the filtered data again using PE.  The PE graph of filtered data LPF 2 Hz shows clear pattern at combination = 5 and = 5 (Fig. 4). We also tried to filter data using LPF 1 Hz, but the final PE result doesn't show any trend.
The unique pattern of PE began appear at 305 Julian date (2.7 days before eruption) same as the result from KLD station (Fig 3 ; Fig. 4). But the both trend are not similar like the result of study [6]. The spectrogram of KWH showed high intensity at the dominant frequency <10 Hz from 270 -304 Julian dates. However, KLD did not show same plot spectrogram at the same period. This may be due to small-scale activity in the Kelud center. Thus, signal from the activity had low energy. As a result, the signal can only be recorded by the nearest seismometer to crater lake (KWH). The volcano tectonic setting is interpreted to have influence on PE result. The tectonic setting of Gjálp fissure is mid ocean ridge and location of the fissure at neovolcanic zones. Neovolcanic zones are the active magmatic zones within divergent plate boundaries [13]. The magmatic system with mid ocean ridge tectonic setting is estimated to have larger scale than subduction tectonic setting. Based on [14], mid ocean ridges are the largest and most voluminous magmatic environment on Earth. So, it is possible that the installed seismometers at around neovolcanic zone will give the similar response.
The PE graphs from UMBK station that located 4.5 km from the crater presents unclear pattern before eruption (Fig. 5). This unclear pattern appear at all combination parameter, without and includes a filtering stage in processing. It is probably due to the location of sensor which is far from the center of magmatic activity (4.5 km). Thus, PE graph cannot describe the pattern that useful as one of eruption precursor. The type and size of eruption are interpreted will give influence because the amount of released seismic energy will different. Effusive eruption will release less amount of energy in comparison with explosive eruption. Based on the result of this study, PE graph will show the pattern that useful as one of sign Kelud eruption precursor, especially dome-forming eruption, if the location of the sensor (short-period seismometer) is not too far, namely about 1 km from magmatic activity center. If the installed instrument is broadband seismometer probably data from sensor at UMBK can show PE pattern as dome-forming eruption precursor. Because the recorded frequency range is getting wider. From the result of this study, we know that using PE as an addition to existing monitoring method can maximize efforts to mitigate volcanic eruption.

Conclusion
PE patterns of dome-forming Kelud eruption period from 3 station are different. The unique PE pattern that appeared before eruption was only shown at Kawah (KWH) and Kelud (KLD) stations. According to this study, PE is effective to forecast the dome-forming eruption of Kelud with a note the sensor (short-period seismometer) distance from magmatic activity center is about 1 km. Using PE as addition volcano monitoring method can help mitigation activities.