TY - JOUR
AU1 - Natsume,, Yuki
AU2 - Ichihara,, Mie
AU3 - Takeo,, Minoru
AB - SUMMARY We perform non-linear time-series analysis on a harmonic tremor seismogram recorded at 830 m away from the centre of the crater during the 2011 eruption at Shinmoedake, Japan. We found features suggesting the existence of period doubling bifurcation in the harmonic tremor signal, implying that the harmonic tremor might be generated by a non-linear process. In order to quantify the non-linearity in the harmonic tremor signal, we measure the correlation dimension D and the maximal Lyapunov exponent λ. For one short but stable segment of the harmonic tremor seismogram, we obtained D = 1.12 and λ = 0.03 s−1. This result implies that the stable oscillation of the harmonic tremor is predominantly a limit cycle with small amounts of chaos present. We then use surrogate data analysis to check that our measurements of D and λ do not include any false positive detection of non-linearity. Limit cycles imply that the harmonic tremor is generated by self-sustained oscillations. We show that the autonomous Julian tremor model is able to exhibit period doubling bifurcation. We also show that the non-autonomous Julian tremor model with a step like increase and decrease in input pressure is able to exhibit oscillations of varying amplitude while keeping a constant frequency spectrum. Both phenomena were observed at Shinmoedake. We also demonstrate that the non-autonomous Julian tremor model with a transient input pressure is able to exhibit long period-like events in addition to harmonic tremor-like events, implying that the same non-linear mechanism could be responsible for the generation of both type of events. Japan, Time-series analysis, Volcano seismology 1 INTRODUCTION Harmonic tremor is a continuous seismic signal that often accompanies volcanic eruptions. As its name implies unlike other forms of volcanic tremor, harmonic tremor exhibits a set of narrow, evenly spaced frequency peaks with a fundamental frequency of f1 and additional harmonic peaks at nf1 where n is some positive integer. Harmonic tremor has been detected at various active volcanoes worldwide such as Karymsky volcano, Russia (Lees et al.2004), Mt Semeru volcano, Indonesia (Schlindwein et al.1995), Lascar volcano, Chile (Hellweg 1999, 2000), Sakurajima volcano, Japan (Maryanto et al.2008) and Shinmoedake volcano, Japan (Ichihara et al.2013). In some cases the harmonic tremor exhibits only a seismic wavefield (seismic harmonic tremor) while in other cases both a seismic and an accompanying acoustic wavefield (seismic acoustic harmonic tremor) have been observed. The presence of an accompanying harmonic acoustic wavefield implies that the observed harmonics in the tremor signal must arise from source effects rather than path or site effects as the acoustic wavefield is presumably transmitted directly from the source to the receiver through the atmosphere. Despite numerous observations, the source mechanism of harmonic tremor remains poorly understood although various models have been proposed. Early models include the resonance of a fluid filled buried cavity such as a spherical magma chamber (Kubotera 1974), a magmatic pipe (Chouet 1985) or a fluid filled crack (Chouet 1988). Julian (1994, 2000) proposed a lumped parameter volcanic tremor model based on non-linear fluid flow through an elastic channel buried within the volcanic edifice. Balmforth et al. (2005) and Sakuraba & Yamauchi (2014) further elaborated on Julian’s ideas by conducting further investigation on the types of instability which might arise in fluid flowing through an elastic conduit, and how those instabilities might be related to volcanic tremor generation. Ripepe & Gordeev (1999) proposed that shallow volcanic tremor at Stromboli could be caused by the the coalescence and bursting of bubbles in magma. In a study on the tremor at Karymsky and Sangay volcanoes (Johnson & Lees 2000), the authors proposed a pressure cooker type mechanism where harmonic tremor is excited when gas periodically escapes through a porous plug situated in the upper part of the conduit. The harmonic tremor at Arenal volcano was proposed to be generated by intermittent gas flow through a fracture located at the top of the conduit, and the harmonic peaks are stabilized by the resonance of the conduit in a similar fashion to how a clarinet produces musical notes (Lesage et al.2006). Hellweg (2000) proposed three possible physical models for the harmonic tremor observed at Lascar volcano, Chile: eddy shedding in a flow of fluid near the volcano’s surface, slug flow within the volcanic conduit and gradual degassing similar to a slightly opened soda bottle. Recently, a magma wagging model wherein a stiff, vertical magma column oscillates within a surrounding annulus of foamy magma was proposed to explain the origin of volcanic tremor in silicic systems (Jellinek & Bercovici 2011; Bercovici et al.2013; Liao et al.2018). Evidence for period doubling bifurcation in the harmonic tremor seismograms has been found at Mt Semeru, Indonesia and Arenal volcano, Costa Rica (Julian 2000; Konstantinou et al.2013). Rust et al. (2008) designed a experiment reminiscent of the Julian tremor model in which gas was forced through a channel cut through a block of gelatine, and showed that under certain conditions period doubling bifurcation in the acoustic oscillations can occur. Assuming autonomous systems, period doubling bifurcation can only occur in a non-linear system of at least three dimensions, implying that harmonic tremor must be generated by a non-linear process. Previous non-linear time-series studies on seismograms of the the volcanic tremor and gas piston events at Kilauea volcano, Hawaii showed that both event types have non-integer correlation dimensions lying in the range of 3.2–4.1 with a mean value of about 3.7 which is highly indicative of non-linear chaos. The similarity of the correlation dimensions of both event types suggest that the tremor has a non-linear source mechanism similar to that of the gas piston events (Chouet & Shaw 1997). Similar studies by Godano et al. (1996) showed that volcanic tremor observed at Mt Etna has a correlation dimension of 1.6, and a study by Godano & Capuano (1999) showed that low frequency events observed at Stromboli and Vulcano have correlation dimensions of 2.68 and 2.83–2.92, respectively. Additionally, studies on the volcanic tremor at Sangay volcano, Ecuador showed that 12 observed tremor events exhibited maximal Lyapunov exponents over a range of 0.013–0.043 and correlation dimensions of relatively low values of 1.2–3.5 from which the authors suggested that volcanic tremors might arise from a system described by a second order non-autonomous non-linear differential equation driven by a sinousoidal forcing term (Konstantinou & Lin 2004). This is in contrast to the autonomous non-linear fluid flow system proposed by Julian (1994). Also, a study regarding of the volcanic tremors recorded at Mt Semeru, Indonesia showed that they exhibited maximal Lyapunov exponents over the range of 0.013–0.039. Additionally based on their results the authors suggested that a linear relationship might exist between the value of the maximal Lyapunov exponent and the number of overtones present in the tremor seismogram (Konstantinou et al.2013). Our aim in this study is to quantify the non-linearity present in the harmonic tremor signal observed during the 2011 Shinmoedake eruption in order to gain further insight into the physical processes driving the harmonic tremor oscillations. In Section 4 of this paper we present evidence for the presence of period doubling bifurcation in the seismograms of the harmonic tremor. In Section 5, we make detailed non-linear time-series analyses of the harmonic tremor seismograms which include attractor reconstruction techniques to measure the correlation dimension and Lyapunov exponent of the harmonic tremor so as to quantify the type of non-linearity of the reconstructed attractor of the harmonic tremor. We then employ a technique called surrogate data analysis (Theiler et al.1992) to rule out the possibility of any false positive detections of non-linearity. Finally, we qualitatively consider the Julian volcanic tremor model (Julian 1994) as a plausible source mechanism for the harmonic tremor at Shinmoedake. 2 PREVIOUS STUDIES OF THE HARMONIC TREMOR DURING THE 2011 SHINMOEDAKE ERUPTION Shinmoedake is an active volcano in the Kirishima volcano group in Kyushu, Japan (Fig. 1). The 2011 eruption at Shinmoedake started with a sub-Plinian stage characterized by three sub-Plinian eruptions on 2011 January 26–27. This was followed by a magma effusive stage from 2011 January 28–31 where a magma pancake formed at the top of the volcano whilst large amounts of degassing was observed together with Vulcanian events. During 2011 February 1–10, Vulcanian eruptions occurred frequently. In short the eruptive behaviour of Shinmoedake from 2011 January 26 to February 10 can be broken down into three distinct stages: a sub-Plinian stage from January 26 to 27, a magma effusive stage from January 28 to 31, and a Vulcanian stage from February 1 to 10 (Nakada et al.2013; Takeo et al.2013). Figure 1. View largeDownload slide Topographic map of the Kirishima volcano group and the emplaced seismic stations. Shinmoedake’s summit crater is marked with a red triangle. In our study we perform non-linear time-series analysis on the harmonic tremor waveform recorded at SMN which is the station closest to the crater of Shinmoedake (830 m away from the centre of the crater). The inset shows the location of the Kirishima volcano group in Kyushu, Japan. Figure 1. View largeDownload slide Topographic map of the Kirishima volcano group and the emplaced seismic stations. Shinmoedake’s summit crater is marked with a red triangle. In our study we perform non-linear time-series analysis on the harmonic tremor waveform recorded at SMN which is the station closest to the crater of Shinmoedake (830 m away from the centre of the crater). The inset shows the location of the Kirishima volcano group in Kyushu, Japan. The Earthquake Research Institute of the University of Tokyo maintains a permanent seismic network at the Kirishima volcano group (Fig. 1). Nanometrics Trillium 120 broad-band seismometers with a sampling rate of 100 Hz were set up at SMN and KRS2 while Guralp CMG-3T broad-band seismometers with a sampling rate of 100 Hz were set up at SMW and TKW. A Hakusan SI102 infrasonic microphone was also set up at SMN together with the seismometer. The Japan Meteorological Agency installed a borehole short period seismometer at KITK at a depth of 100 m below the surface, located about 3  km southeast from the crater of Shinmoedake. In addition to these stations, the Japan Meteorological Agency and the National Research Institute for Earth Science and Disaster Resilience maintain other seismic stations whose locations are marked in Fig. 1 as unnamed points. Seismic harmonic tremor at Shinmoedake was first observed during the magma effusive stage on January 30 at around 11:30, and was then observed to occur intermittently until January 31. During this period, the seismic harmonic tremor had clear but unstable modes in the frequency spectra and the fundamental frequency was noted to range from 0.85 to 1.34 Hz (Table 1). Seismic acoustic harmonic tremor was first observed on January 31 at 21:49 after magma effusion stopped, and by February 6 a total of six seismic acoustic harmonic tremor events were recorded simultaneously on both the seismometers and microphone. All the seismic acoustic harmonic tremor events displayed more distinct and stable modes as compared to the preceding seismic harmonic tremors. It has been proposed that seismic acoustic harmonic tremor is likely to be generated by gas flow rather than magma flow, and since both the seismic harmonic tremor and seismic acoustic harmonic tremor share common harmonic modes and occur within the same time frame, it is highly plausible that both types of events share the same common gas driven source mechanism (Ichihara et al.2013). Additionally, earlier studies by Takeo et al. (2013) also demonstrated synchronization of harmonic tremor events with the inflation-deflation cycles of the volcanic edifice. The start of a deflation cycle coincided with the start of a harmonic tremor event, and the amplitude of the tremor gradually decreased at the start of the subsequent inflation cycle. This synchronization was interpreted by the authors to mean that the harmonic tremor events were caused by the self-sustained oscillation of fluid flow through the conduit with changing pressure or volume. Table 1. List of harmonic tremor events at Shinmoedake. Event Date Time (total duration (s)) f1 (Hz) Type Period type 1 2011 January 30 11:36:04–11:45:22 (558) 0.85 SHT – 2 2011 January 30 12:17:08–12:19:18 (130) 0.88 SHT Period 2 3 2011 January 30 13:28:46–13:35:24 (398) 0.95 SHT Period 2 4 2011 January 30 14:00:45–14:14:18 (853) 0.98 SHT Period 2 5 2011 January 30 14:45:18–14:58:12 (774) 0.98 SHT Period 2 6 2011 January 30 15:24:19–15:28:09 (230) 1.03 SHT Period 2 7 2011 January 30 17:04:06–17:22:36 (1110) 1.03 SHT Period 2 8 2011 January 30 18:27:12–18:28:24 (72) 1.05 SHT Period 2 9 2011 January 30 19:46:48–19:54:54 (486) 0.88 SHT Period 2 10 2011 January 31 00:03:40–00:20:04 (984) 0.85 SHT Period 2 11 2011 January 31 01:00:24–01:20:48 (1224) 0.85 SHT Period 2 12 2011 January 31 01:51:24–02:17:06 (1542) 1.34 SHT Period 1 13 2011 January 31 02:56:12–03:12:30 (978) 1.32 SHT Period 1 14 2011 January 31 03:54:57–04:06:33 (696) 1.29 SHT Period 1 15 2011 January 31 04:46:42–04:56:09 (567) 1.07 SHT Period 2 16 2011 January 31 06:34:48–06:41:18 (390) 1.02 SHT Period 2 17 2011 January 31 21:49:05–21:53:10 (245) 1.54 SAHT Period 1 18 2011 February 01 08:07:55–08:10:01 (126) 0.63 SAHT Period 2 19 2011 February 02 05:23:30–05:25:00 (90) 1.78 SHT Period 1 20 2011 February 02 06:16:48–06:27:24 (636) 1.03 SHT Period 2 21 2011 February 02 07:10:50–07:33:00 (1330) 0.88 SHT Period 2 22 2011 February 02 20:43:42–21:21:42 (2280) 0.93 SAHT Period 2 23 2011 February 03 06:08:30–06:10:48 (138) 1.12 SAHT Period 2 24 2011 February 03 13:24:05–13:58:10 (2045) 1.71 SAHT Period 1 25 2011 February 06 11:59:48–12:01:53 (125) 0.85 SAHT Period 2 Event Date Time (total duration (s)) f1 (Hz) Type Period type 1 2011 January 30 11:36:04–11:45:22 (558) 0.85 SHT – 2 2011 January 30 12:17:08–12:19:18 (130) 0.88 SHT Period 2 3 2011 January 30 13:28:46–13:35:24 (398) 0.95 SHT Period 2 4 2011 January 30 14:00:45–14:14:18 (853) 0.98 SHT Period 2 5 2011 January 30 14:45:18–14:58:12 (774) 0.98 SHT Period 2 6 2011 January 30 15:24:19–15:28:09 (230) 1.03 SHT Period 2 7 2011 January 30 17:04:06–17:22:36 (1110) 1.03 SHT Period 2 8 2011 January 30 18:27:12–18:28:24 (72) 1.05 SHT Period 2 9 2011 January 30 19:46:48–19:54:54 (486) 0.88 SHT Period 2 10 2011 January 31 00:03:40–00:20:04 (984) 0.85 SHT Period 2 11 2011 January 31 01:00:24–01:20:48 (1224) 0.85 SHT Period 2 12 2011 January 31 01:51:24–02:17:06 (1542) 1.34 SHT Period 1 13 2011 January 31 02:56:12–03:12:30 (978) 1.32 SHT Period 1 14 2011 January 31 03:54:57–04:06:33 (696) 1.29 SHT Period 1 15 2011 January 31 04:46:42–04:56:09 (567) 1.07 SHT Period 2 16 2011 January 31 06:34:48–06:41:18 (390) 1.02 SHT Period 2 17 2011 January 31 21:49:05–21:53:10 (245) 1.54 SAHT Period 1 18 2011 February 01 08:07:55–08:10:01 (126) 0.63 SAHT Period 2 19 2011 February 02 05:23:30–05:25:00 (90) 1.78 SHT Period 1 20 2011 February 02 06:16:48–06:27:24 (636) 1.03 SHT Period 2 21 2011 February 02 07:10:50–07:33:00 (1330) 0.88 SHT Period 2 22 2011 February 02 20:43:42–21:21:42 (2280) 0.93 SAHT Period 2 23 2011 February 03 06:08:30–06:10:48 (138) 1.12 SAHT Period 2 24 2011 February 03 13:24:05–13:58:10 (2045) 1.71 SAHT Period 1 25 2011 February 06 11:59:48–12:01:53 (125) 0.85 SAHT Period 2 View Large Table 1. List of harmonic tremor events at Shinmoedake. Event Date Time (total duration (s)) f1 (Hz) Type Period type 1 2011 January 30 11:36:04–11:45:22 (558) 0.85 SHT – 2 2011 January 30 12:17:08–12:19:18 (130) 0.88 SHT Period 2 3 2011 January 30 13:28:46–13:35:24 (398) 0.95 SHT Period 2 4 2011 January 30 14:00:45–14:14:18 (853) 0.98 SHT Period 2 5 2011 January 30 14:45:18–14:58:12 (774) 0.98 SHT Period 2 6 2011 January 30 15:24:19–15:28:09 (230) 1.03 SHT Period 2 7 2011 January 30 17:04:06–17:22:36 (1110) 1.03 SHT Period 2 8 2011 January 30 18:27:12–18:28:24 (72) 1.05 SHT Period 2 9 2011 January 30 19:46:48–19:54:54 (486) 0.88 SHT Period 2 10 2011 January 31 00:03:40–00:20:04 (984) 0.85 SHT Period 2 11 2011 January 31 01:00:24–01:20:48 (1224) 0.85 SHT Period 2 12 2011 January 31 01:51:24–02:17:06 (1542) 1.34 SHT Period 1 13 2011 January 31 02:56:12–03:12:30 (978) 1.32 SHT Period 1 14 2011 January 31 03:54:57–04:06:33 (696) 1.29 SHT Period 1 15 2011 January 31 04:46:42–04:56:09 (567) 1.07 SHT Period 2 16 2011 January 31 06:34:48–06:41:18 (390) 1.02 SHT Period 2 17 2011 January 31 21:49:05–21:53:10 (245) 1.54 SAHT Period 1 18 2011 February 01 08:07:55–08:10:01 (126) 0.63 SAHT Period 2 19 2011 February 02 05:23:30–05:25:00 (90) 1.78 SHT Period 1 20 2011 February 02 06:16:48–06:27:24 (636) 1.03 SHT Period 2 21 2011 February 02 07:10:50–07:33:00 (1330) 0.88 SHT Period 2 22 2011 February 02 20:43:42–21:21:42 (2280) 0.93 SAHT Period 2 23 2011 February 03 06:08:30–06:10:48 (138) 1.12 SAHT Period 2 24 2011 February 03 13:24:05–13:58:10 (2045) 1.71 SAHT Period 1 25 2011 February 06 11:59:48–12:01:53 (125) 0.85 SAHT Period 2 Event Date Time (total duration (s)) f1 (Hz) Type Period type 1 2011 January 30 11:36:04–11:45:22 (558) 0.85 SHT – 2 2011 January 30 12:17:08–12:19:18 (130) 0.88 SHT Period 2 3 2011 January 30 13:28:46–13:35:24 (398) 0.95 SHT Period 2 4 2011 January 30 14:00:45–14:14:18 (853) 0.98 SHT Period 2 5 2011 January 30 14:45:18–14:58:12 (774) 0.98 SHT Period 2 6 2011 January 30 15:24:19–15:28:09 (230) 1.03 SHT Period 2 7 2011 January 30 17:04:06–17:22:36 (1110) 1.03 SHT Period 2 8 2011 January 30 18:27:12–18:28:24 (72) 1.05 SHT Period 2 9 2011 January 30 19:46:48–19:54:54 (486) 0.88 SHT Period 2 10 2011 January 31 00:03:40–00:20:04 (984) 0.85 SHT Period 2 11 2011 January 31 01:00:24–01:20:48 (1224) 0.85 SHT Period 2 12 2011 January 31 01:51:24–02:17:06 (1542) 1.34 SHT Period 1 13 2011 January 31 02:56:12–03:12:30 (978) 1.32 SHT Period 1 14 2011 January 31 03:54:57–04:06:33 (696) 1.29 SHT Period 1 15 2011 January 31 04:46:42–04:56:09 (567) 1.07 SHT Period 2 16 2011 January 31 06:34:48–06:41:18 (390) 1.02 SHT Period 2 17 2011 January 31 21:49:05–21:53:10 (245) 1.54 SAHT Period 1 18 2011 February 01 08:07:55–08:10:01 (126) 0.63 SAHT Period 2 19 2011 February 02 05:23:30–05:25:00 (90) 1.78 SHT Period 1 20 2011 February 02 06:16:48–06:27:24 (636) 1.03 SHT Period 2 21 2011 February 02 07:10:50–07:33:00 (1330) 0.88 SHT Period 2 22 2011 February 02 20:43:42–21:21:42 (2280) 0.93 SAHT Period 2 23 2011 February 03 06:08:30–06:10:48 (138) 1.12 SAHT Period 2 24 2011 February 03 13:24:05–13:58:10 (2045) 1.71 SAHT Period 1 25 2011 February 06 11:59:48–12:01:53 (125) 0.85 SAHT Period 2 View Large A previous experimental study by Lyons et al. (2013) showed that the switching between seismic harmonic tremor and seismic acoustic harmonic tremor might arise from the increase in viscosity of the magma pancake situated on top of the crater of Shinmoedake. The results from this experimental study was compared with observation evidence and it was proposed that during the course of the magma effusive stage the magma pancake became more viscous, allowing a permanent conduit to eventually form through the accumulated magma. The formation of this permanent conduit allowed gases to escape directly from the harmonic tremor source into the atmosphere, therefore allowing the tremor to exhibit both seismic and acoustic wavefields (Ichihara et al.2013). However these studies did not attempt to explain the possible mechanism behind tremor generation. A study using MUSIC analysis was also done on the spatial change of the source location of one harmonic tremor event on 2011 February 02, 20:57 (Matsumoto et al.2013). The authors concluded that the harmonic tremor source is likely to be located in a region just beneath the crater, although for several short intervals the direction of the seismic wave propagation was observed to suddenly change to the geodetic pressure source determined by Nakao et al. (2013) and Ueda et al. (2013). This result was interpreted to indicate the possible presence of at least two major tremor sources within the volcanic edifice, situated at opposite ends of the conduit linking the magma source to a shallow depth beneath the crater. When tremor occurred at one end of the conduit, some sort of instability might arise within the conduit due to pressure or flow gradients leading to the excitation of tremor at the other end (Matsumoto et al.2013). 3 RELATIONSHIP BETWEEN LONG PERIOD EVENTS AND HARMONIC TREMOR AT SHINMOEDAKE Swarms of long period events were also detected during the same time frame as the harmonic tremor events during the 2011 Shinmoedake eruption. Long period events resemble small tectonic earthquakes in duration but have a much longer period range of about 0.2–2 s compared to tectonic events of similar magnitude, while harmonic tremor is characterized by a sustained harmonic signal, generally lasting from minutes to days (Chouet 1996). The long period events at Shinmoedake share the same frequency band as that of the harmonic tremor despite having a broader peak due to a shorter signal in the time domain (Fig. 2), hinting that both types of events share a common source mechanism. In such a case the harmonic tremor is the system’s response to a prolonged excitation while the long period event is the corresponding response to an impulsive excitation (Chouet 1996). In Section 9 of this paper we show that the Julian volcanic tremor model (Julian 1994) is able to exhibit long period event-like waveforms to support this common source assumption. Figure 2. View largeDownload slide Harmonic tremor waveform (a) and power spectral density (b) observed on 2011 January 31 at 21:49:05 and long period event waveform (c) and power spectral density (d) observed at SMN on 2011 February 01 at 00:32:35. Both events share the same dominant peak frequency at 1.5 Hz indicating a possible common source mechanism. Figure 2. View largeDownload slide Harmonic tremor waveform (a) and power spectral density (b) observed on 2011 January 31 at 21:49:05 and long period event waveform (c) and power spectral density (d) observed at SMN on 2011 February 01 at 00:32:35. Both events share the same dominant peak frequency at 1.5 Hz indicating a possible common source mechanism. We emphasize at this point that we have currently not considered other seismological information and data such as the source location and depth, as well as the moment tensor solutions of both event types. In addition, we note that in a recent study it was demonstrated that slow-rupture failure of unconsolidated material in the volcanic edifice is able to produce simple, short duration waveforms similar to the long period events observed at Shinmoedake (Bean et al.2014). Therefore, we stress that further analyses such as moment tensor inversion must be done on both event types before a definitive conclusion can be reached on whether the long period events and harmonic tremor truly share the same source mechanism. In this study, we are essentially proposing that the advantage of the Julian volcanic tremor model (Julian 1994) is its ability to generate both long period-like events and harmonic tremor-like events with the same peak frequency under different conditions. 4 PERIOD DOUBLING BIFURCATION OF HARMONIC TREMOR 4.1 Characterization of harmonic tremor at Shinmoedake The seismic acoustic harmonic tremor waveform recorded on 2011 February 02, 20:52:20 across 5 seismic stations (SMN, SMW, KRS2, TKW and KITK) and 1 microphone (same location as SMN) exhibits the same set of frequency peaks (Fig. 3), confirming that the harmonic peaks cannot arise from path or site effects, and must arise from the source. We also note that in Fig. 3 the dominant mode of oscillation of the seismic wavefield recorded at all 5 seismic stations is the second mode at 1.86 Hz, while the dominant mode of oscillation of the acoustic wavefield is the fundamental mode at 0.93 Hz. Figure 3. View largeDownload slide Vertical component waveforms (left-hand panels) and corresponding power spectral densities (right-hand panels) of a harmonic tremor event recorded on 2011 February 02, 20:52:20 across five seismic stations (SMN, SMW, KRS2, TKW and KITK) and 1 microphone (SMN). All waveforms share the same harmonic peaks implying that they arise from source effects rather than path or site effects. Figure 3. View largeDownload slide Vertical component waveforms (left-hand panels) and corresponding power spectral densities (right-hand panels) of a harmonic tremor event recorded on 2011 February 02, 20:52:20 across five seismic stations (SMN, SMW, KRS2, TKW and KITK) and 1 microphone (SMN). All waveforms share the same harmonic peaks implying that they arise from source effects rather than path or site effects. We observed a total of 25 harmonic tremor events from 2011 January 30 to February 6 (Fig. 4, Table 1). Only continuous events which had a duration of at least 1 min, and exhibited at least two harmonic peaks were classified to be harmonic tremor events in this study. We also note that not all harmonic tremor events occurred as one long continuous event. In such cases if the time interval between individual harmonic tremor events was less than 30 min long, and if the structure of the frequency spectrum does not change across the time interval, we considered all the individual events to be subevents of one long event. In such cases we took the difference of the end time of the last individual event and the start time of the first individual event as the duration of the event. On the other hand if the time interval exceeded 30 min, even if the frequency spectrum was constant we classified the events as separate events. In order to identify the fundamental mode of a particular event, we checked the spectrograms for regions of constant and stable harmonic modes, and recorded the fundamental frequency of this stable region. We admit that in many events spectral gliding was present, and such a classification method captures only the frequency content of the event at a particular time. Figure 4. View largeDownload slide Harmonic modes for the (a) east–west, (b) north–south and (c) up–down seismogram (SMN) components of 21 harmonic tremor events observed during the 2011 Shinmoedake eruption. Black dots indicate period 2 events while red dots indicate period 1 events (see text for definitions). Some amount of spectral gliding of the modes exist. In some events some of the modes were buried within noise and precise identification of those peaks was difficult, so we have left them out in the figure. Figure 4. View largeDownload slide Harmonic modes for the (a) east–west, (b) north–south and (c) up–down seismogram (SMN) components of 21 harmonic tremor events observed during the 2011 Shinmoedake eruption. Black dots indicate period 2 events while red dots indicate period 1 events (see text for definitions). Some amount of spectral gliding of the modes exist. In some events some of the modes were buried within noise and precise identification of those peaks was difficult, so we have left them out in the figure. Analysis of the frequency spectrum shows that the harmonic tremor exhibits a relatively stable set of frequency modes throughout the entire sequence despite the presence of spectral gliding, implying that the harmonic tremor mechanism remained relatively stable in time. We also note the presence of mode hopping of the dominant oscillation mode. Referring to Fig. 5, it can be seen that for the up–down seismic component for the time duration of 4–6 min the fundamental mode is the dominant mode, while for the time duration of 9–10 min the second mode is dominant. In contrast no mode hopping was observed in the east–west and north–south seismic components where the second mode is dominant or in the acoustic wavefield where the fundamental mode is dominant. Mode hopping implies that the harmonic tremor source does not oscillate consistently at a particular dominant mode, indicating that a transitional process is occurring at the source region despite generating a continuous signal with harmonic peaks. We note that mode hopping is a phenomenon frequently encountered in lasers, where the dominant mode of laser light oscillation hops to other nearby modes due to the optical gain spectrum of the laser shifting to other wavelengths (for an example, see fig. 5 in Saulys et al.2010). Figure 5. View largeDownload slide Harmonic tremor recorded at SMN (3 seismic components + infrasound) on 2011 February 02, 20:43:00. Despite showing relatively stable frequency peaks, the amplitude of the seismic signal changes significantly over time, and mode hopping occurs in the up–down seismic component. The dominant oscillation mode can be different between the different wavefields. Figure 5. View largeDownload slide Harmonic tremor recorded at SMN (3 seismic components + infrasound) on 2011 February 02, 20:43:00. Despite showing relatively stable frequency peaks, the amplitude of the seismic signal changes significantly over time, and mode hopping occurs in the up–down seismic component. The dominant oscillation mode can be different between the different wavefields. 4.2 Phase portrait analysis and period doubling bifurcation A common way to analyse a non-linear time-series is to study its phase portrait. Phase portraits show the relationship between two different state variables of a non-linear system. For example, a common textbook example of a phase portrait is the analysis of the relationship between the velocity and displacement of an oscillating pendulum by plotting its velocity against its displacement (Strogatz 1994a). In the same way, we can study the relationship between the velocity and displacement of the harmonic tremor signal through a velocity–displacement phase portrait by plotting the velocity seismogram against the displacement seismogram. Analysis of the phase portrait can reveal inherently non-linear features such as period doubling bifurcation. Period doubling bifurcation occurs when a variation in a parameter of the non-linear system causes the original periodic trajectory to take twice as long to repeat itself. When this happens, the original periodic trajectory in phase space changes from a single loop (period 1 state) to a double nested loop (period 2 state). According to the Poincare–Bendixson theorem, any trajectory that is confined to a closed, bounded region of a two dimensional plane that does not contain any fixed points must eventually form a closed orbit. The consequence of this theorem is that for non-linear autonomous systems, behaviour such as period doubling bifurcation and chaotic motion can only occur in systems of at least three dimensions (Strogatz 1994b; Jordan & Smith 2007). Trajectories cannot cross each other on a two dimensional plane, however the presence of a third (or higher) dimension allows the trajectories to go over/under each other without crossing to form the characteristic double nested loop when period doubling bifurcation occurs. Therefore the presence of period doubling bifurcation means that the underlying non-linear system of the harmonic tremor must be of at least 3 dimensions. Also, when period doubling bifurcation occurs, additional harmonic peaks appear exactly in between the original harmonic peaks of the time-series frequency spectrum. The appearance of these new harmonic peaks corresponds to the doubling of the period. Analysis of the harmonic tremor’s frequency spectrum reveals the presence of period doubling bifurcation in the harmonic tremor oscillation as follows. On January 31 (Fig. 6, 110131), the frequency bands present are 1.8, 3.6 and 5.4 Hz (period 1 state). Three days later on February 2 (Fig. 6, 110202), frequency bands of 0.9, 2.7 and 4.5 Hz appear exactly in between the original three frequency bands of 1.8, 3.6 and 5.4 Hz (period 2 state). One day later on February 3 (Fig. 6, 110203), the three new frequency bands of 0.9, 2.7 and 4.5 Hz are absent, exhibiting only the original frequency bands (period 1 state). Three days later on February 6 (Fig. 6, 110206), the three new bands appear again (period 2 state). This corresponds to a period doubling bifurcation occurring during January 31–February 2, a reverse period doubling bifurcation during February 2–February 3, and another period doubling bifurcation during February 3–February 6. Figure 6. View largeDownload slide Spectrograms of the harmonic tremor at Shinmoedake in the order of 2011 January 31 (110131), February 2 (110202), February 3 (110203) and February 6 (110206), corresponding to events 17, 22, 24 and 25 in Fig. 4. Note how the frequency bands of 0.9, 2.7 and 4.5 Hz are absent on January 31, present on February 2, absent on February 3 and present again on February 6. Figure 6. View largeDownload slide Spectrograms of the harmonic tremor at Shinmoedake in the order of 2011 January 31 (110131), February 2 (110202), February 3 (110203) and February 6 (110206), corresponding to events 17, 22, 24 and 25 in Fig. 4. Note how the frequency bands of 0.9, 2.7 and 4.5 Hz are absent on January 31, present on February 2, absent on February 3 and present again on February 6. Plotting the velocity–displacement phase portraits of the north–south motion of the harmonic tremor recorded at SMN (Fig. 1) shows the following changes: on January 31 (Fig. 7, 110131) the phase portrait is a single loop (period 1 state), three days later on February 2 (Fig. 7, 110202) the phase portrait is a double nested loop (period 2 state). On February 3 (Fig. 7, 110203), the phase portrait returns to a single loop (period 1 state) and on February 6 (Fig. 7, 110206), the phase portrait is again a double nested loop (period 2 state). The appearance of the nested second loop is equivalent to the appearance of the new frequency bands in between the original frequency bands in the spectra in Fig. 6 and is indicative of period doubling bifurcation (period 1 → period 2), while the disappearance of the nested second loop is equivalent to the disappearance of the new frequency bands and is indicative of reverse period doubling bifurcation (period 2 → period 1). Figure 7. View largeDownload slide Phase portraits of the harmonic tremor at Shinmoedake corresponding to Fig. 6. Note how the phase portrait is a single loop on January 31, a double nested loop on February 2, a single loop on February 3 and a double nested loop again on February 6. This corresponds to the appearance and disappearance of the alternate, evenly spaced frequency bands in Fig. 6. Such behaviour is indicative of period doubling bifurcation. Figure 7. View largeDownload slide Phase portraits of the harmonic tremor at Shinmoedake corresponding to Fig. 6. Note how the phase portrait is a single loop on January 31, a double nested loop on February 2, a single loop on February 3 and a double nested loop again on February 6. This corresponds to the appearance and disappearance of the alternate, evenly spaced frequency bands in Fig. 6. Such behaviour is indicative of period doubling bifurcation. In additon, we discovered a harmonic tremor event on 2011 February 03, 06:08:38 that appears to transition from the period 2 state to the period 1 state (reverse period doubling bifurcation) within the same time-series (Fig. 8). This particular transition occurred towards the end of the event, just as the oscillatory motion was decreasing in amplitude (Fig. 8a). The velocity–displacement phase portrait during the period 2 region shows a double nested loop (Fig. 8b) while the phase portrait during the period 1 region shows a single loop (Fig. 8c) as with the case in Fig. 7. From a non-linear dynamic system standpoint, this hints that the activation/deactivation of a harmonic tremor event is a direct result of a changing parameter in the harmonic tremor system. Figure 8. View largeDownload slide (a) Waveform of the north–south motion at SMN showing transition from period 2 (from 0 to 9 s) to period 1 (from 9 to 18 s) oscillation within the same harmonic tremor event detected on 2011 February 03, 06:08:38. Phase portraits of velocity versus displacement (b) in the period 2 region and (c) in the period 1 region. Figure 8. View largeDownload slide (a) Waveform of the north–south motion at SMN showing transition from period 2 (from 0 to 9 s) to period 1 (from 9 to 18 s) oscillation within the same harmonic tremor event detected on 2011 February 03, 06:08:38. Phase portraits of velocity versus displacement (b) in the period 2 region and (c) in the period 1 region. In order to track how many times period doubling bifurcation of the harmonic tremor occurred during the 2011 Shinmoedake eruption, we then attempted to classify each of the observed harmonic tremor events (Fig. 4, Table 1) as either 'period 1' or 'period 2' events based on the structure of the velocity–displacement phase portraits. If the amount of noise was so large that it was difficult to resolve the structure of the phase portraits accurately, we did not attempt to clarify the presence of period doubling bifurcation. Period doubling bifurcation implies that the harmonic tremor must be generated by an autonomous non-linear system of at least three dimensions. This allows us to rule out all linear resonant models or linear superposition models for the source mechanism of the harmonic tremor and gives support to the hypothesis that harmonic tremor is an excitation arising from a non-linear process. In the following sections we characterize the non-linearity of the harmonic tremor in order to elucidate its non-linear dynamics. 5 RECONSTRUCTING THE ATTRACTOR OF THE HARMONIC TREMOR The non-linear time-series analysis methods used in this study is applicable to signals generated by autonomous systems. We extract a 50-s-long signal with relatively stable amplitude, no spectral gliding, and no mode hopping and assume that this signal is generated by an autonomous sytem (Fig. 9). This signal was extracted from the longest and most stable event which was observed on 2011 February 2, 20:51:50. We have chosen not to use the other portions of the signal as certain characteristics such as the changes in the signal amplitude, spectral gliding or mode hopping appear to arise from a non-autonomous system (Fig. 5). The signal used was recorded at the station nearest to the crater (830 m from the centre of the crater, SMN NS component, 2011 February 2, 20:51:50) in order to minimize any non-linear transformations due to wave propagation. The signal was detrended to remove any linear trends, and was filtered with a 200-tap FIR filter for the frequency range of 0.5–4 Hz to remove high frequency components which are more susceptible to scattering. Data from non-linear systems should not be filtered with IIR filters which might corrupt the original non-linear dynamics of the data. FIR filters do not exhibit this issue (Broomhead et al.1992). Figure 9. View largeDownload slide Harmonic tremor recorded at SMN (north–south component) on 2011 February 02, 20:51:50. The entire length of the time-series is 5000 data points long sampled at a rate of 100 Hz for a total time of 50 seconds. This seismogram is used in the non-linear time-series analysis. Figure 9. View largeDownload slide Harmonic tremor recorded at SMN (north–south component) on 2011 February 02, 20:51:50. The entire length of the time-series is 5000 data points long sampled at a rate of 100 Hz for a total time of 50 seconds. This seismogram is used in the non-linear time-series analysis. We use Taken’s embedding theorem to create from a single discretized time-series yt, where t = 0, 1, ..., M − 1, the following N dimensional vector in time delay coordinate space using some time delay T \begin{eqnarray*} x_t(N) = [y_t,y_{t+T},y_{t+2T}, ..., y_{t+(N-1)T}], \end{eqnarray*} (1) which is a reconstructed attractor that is topologically equivalent to the original attractor of the harmonic tremor (Takens 1981). We need to determine the appropriate value of the time delay T and the number of embedding dimensions N to be used in the reconstruction. Both values of T and N used to reconstruct the attractor must be carefully chosen; if N is too small the reconstructed attractor is not fully unfolded and the calculated topological properties of the attractor are most likely to be wrong, on the other hand if N is too large the calculated topological properties of the attractor become spurious and imprecise. During reconstruction we normalized the size of the entire reconstructed attractor by its maximum norm. We estimated the appropriate value of the time delay T using two different methods. In the first method T was estimated from the point nearest to the first zero crossing of the autocorrelation function A(τ) of the time-series yt where τ is the relative discretized time delay in order to minimize the correlation between the different dimensions of the reconstructed attractor (Theiler 1990; Albano et al.1991). For the harmonic tremor the point at τ = 13 lies closest to the first zero crossing (Fig. 10). In the second method T was estimated from the point of the first minimum of the mutual information function defined as (Fraser & Swinney 1986): \begin{eqnarray*} MI(\tau )=\sum _{ij}p_{ij}(\tau ) \log p_{ij}(\tau ) - 2 \sum _k p_k(\tau ) \log p_k(\tau ), \end{eqnarray*} (2) where pi is the probability of some data point from the time-series yt falling within the ith bin of a histogram, while pij is the probability of yt falling within the ith bin and yt + τ falling within the jth bin. For the harmonic tremor, the first minimum occurs at the point at τ = 12 (Fig. 11). As the values of τ from both methods are very close to each other, we used a delay time value of T = 13 in the reconstruction of the harmonic tremor attractor. Figure 10. View largeDownload slide Autocorrelation function A(τ) for the signal in Fig. 9. The point closest to the first zero crossing is located at τ = 13. Figure 10. View largeDownload slide Autocorrelation function A(τ) for the signal in Fig. 9. The point closest to the first zero crossing is located at τ = 13. Figure 11. View largeDownload slide Graph of the mutual information for the signal in Fig. 9. The first minimum occurs at τ = 12. This value is very close to the value of 13 obtained through the autocorrelation function in Fig. 10. Figure 11. View largeDownload slide Graph of the mutual information for the signal in Fig. 9. The first minimum occurs at τ = 12. This value is very close to the value of 13 obtained through the autocorrelation function in Fig. 10. We then used a modified nearest neighbours method (Cao 1997) to determine the value of the minimum embedding dimension N to be used in the reconstruction (the details of Cao’s method is given explicitly in the Appendix). In order to overcome the problems involved with choosing an arbitrary threshold used in traditional nearest neighbours methods, Cao calculated the relative change of the mean of the ratio of the distance between 2 nearest neighbours as the embedding dimension is increased from N to N + 1. The distance between two true nearest neighbours do not change significantly when N is increased while the distance between two false nearest neighbours will increase significantly as the attractor is unfolded. The value of N where the relative change reaches some saturation value is the minimum embedding dimension needed to fully unravel the reconstructed attractor where no crossing of trajectories occur. Additionally, Cao’s method is also able to distinguish deterministic signals from stochastic noise. Our calculations of the two Cao’s coefficients (Cao 1997) E1(N) and E2(N) for the harmonic tremor signal show that for a time delay of T = 13, an embedding dimension of N = 7 is needed to fully reconstruct the attractor (Fig. 12). The resulting reconstructed attractor (Fig. 13) shows that the harmonic tremor resembles a limit cycle, however, the trajectories do not repeat exactly but diverge slowly over time leading to some amount of ‘thickness’ in the trajectories. Figure 12. View largeDownload slide Cao coefficients E1(N) and E2(N) calculated for the signal in Fig. 9. Both E1(N) and E2(N) saturate at embedding dimension N = 7. Figure 12. View largeDownload slide Cao coefficients E1(N) and E2(N) calculated for the signal in Fig. 9. Both E1(N) and E2(N) saturate at embedding dimension N = 7. Figure 13. View largeDownload slide Three dimensional plots of 3 out of the 7 dimensions of the reconstructed attractor of the signal in Fig. 9, showing that the time-series of the harmonic tremor is not a true limit cycle as the orbits do not repeat exactly, but diverge gradually over time leading to some amount of thickness in the trajectories. Figure 13. View largeDownload slide Three dimensional plots of 3 out of the 7 dimensions of the reconstructed attractor of the signal in Fig. 9, showing that the time-series of the harmonic tremor is not a true limit cycle as the orbits do not repeat exactly, but diverge gradually over time leading to some amount of thickness in the trajectories. 6 CORRELATION DIMENSION OF THE HARMONIC TREMOR The fractal dimension of an attractor is a representation of the minimum number of degrees of freedom of the underlying non-linear system. Amongst the various possible methods available we utilize the algorithm proposed by Grassberger & Procaccia (1983) to measure the correlation dimension of the harmonic tremor given by the equation \begin{eqnarray*} C(r,N_p) &=& \frac{2}{(N_p+1-W)(N_p-W)} \nonumber \\ && \times \, \sum _{n=W}^{N_p-1} \sum _{i=0}^{N_p-1-n} \Theta (r - || \vec{x}_i - \vec{x}_{i+n} ||), \end{eqnarray*} (3) where r is the Euclidean distance in N dimensional phase space, Np is the number of data points in the reconstructed attractor, W is the Theiler window (W = 500 is used in this study) and Θ(r) is the step function. This algorithm calculates the average distances between all pairs of points as long as they are at least W data points apart. Using this Theiler window W removes the effects of autocorrelation which might give rise to spurious estimates of the correlation dimension (Theiler 1990). The correlation dimension D of the harmonic tremor is then estimated using |$\frac{d \log C}{d \log r}$| of the constant gradient region of the graph of log C versus log r (Grassberger & Procaccia 1983). We calculated the values of |$\frac{d \log C}{d \log r}$| for a range of embedding dimensions N ∈ [2, 14] to observe how the value of the correlation dimension D changes with N. Fig. 14 shows that the value of |$\frac{d \log C}{d \log r}$| increases rapidly from N = 2 and saturates to a constant value from N = 7 with mean value of 1.1205 ± 0.0004, where the error was obtained from the standard deviation of the 8 individual slopes. This result is consistent with the calculations above using Cao’s method which determined the minimum embedding dimension to be N = 7. The value of |$\frac{d \log C}{d \log r}$| for N = 7 was measured to be D = 1.12. The important result is that the measured value of D = 1.12 is close to the value of 1 which is the correlation dimension of a limit cycle, implying that the dynamics of the harmonic tremor is predominantly a limit cycle with some small amounts of chaotic dynamics present. Figure 14. View largeDownload slide Left-hand panel: plots of log C versus log r for embedding dimensions of 2 to 14 for the harmonic tremor signal in Fig. 9. Dotted lines show the region used in the fit; for embedding dimensions N = 2, 3 the fitting region lies within log r ∈ [ − 1, −0.4] while for the other values of N the fitting regions lie within log r ∈ [ − 1, 0]. Right-hand panel: the values of the correlation dimension for N = 7 is D = 1.12. Figure 14. View largeDownload slide Left-hand panel: plots of log C versus log r for embedding dimensions of 2 to 14 for the harmonic tremor signal in Fig. 9. Dotted lines show the region used in the fit; for embedding dimensions N = 2, 3 the fitting region lies within log r ∈ [ − 1, −0.4] while for the other values of N the fitting regions lie within log r ∈ [ − 1, 0]. Right-hand panel: the values of the correlation dimension for N = 7 is D = 1.12. 7 LYAPUNOV EXPONENT OF THE HARMONIC TREMOR Trajectories on chaotic attractors exhibit extremely sensitive dependence on initial conditions, those that were originally very close together diverge quickly from each other and the futures of both trajectories become unpredictable. For low dimensional deterministic chaotic systems the Lyapunov exponent λ represents the rate of divergence δ(t) = δ(0)eλt, λ > 0 of two trajectories that were initially very close to each other. Strictly speaking an N dimensional system will have N independent Lyapunov exponents, many of which are either zero or negative. However the maximal Lyapunov exponent effectively determines the dynamics of the trajectories and we will measure only this exponent of the harmonic tremor. We utilize an algorithm proposed by Kantz (1994) to measure the maximal Lyapunov exponent of a reconstructed attractor as follows. First a reference point in phase space |$\mathbf {x}_t$| is chosen on the reconstructed attractor and the space around this point is searched within an ε-neighbourhood of Ut for all neighbours |$\mathbf {x}_i$|⁠. The average distance from the reference point |$\mathbf {x}_t$| to all the neighbours |$\mathbf {x}_i$| is then measured for some relative shifted time τ. This measurement is then repeated for all points on the trajectory Np. As with the correlation dimension we used a Theiler window of 500 points in the calculations in order to minimize the effects of autocorrelation (Theiler 1990). Mathematically this algorithm is written as \begin{eqnarray*} S(\tau )=\frac{1}{N_p} \sum _{t=1}^{N_p} \log \left( \frac{1}{V} \sum _{i \in U_t} |\mathbf {x}_{t+\tau } - \mathbf {x}_{i+\tau }| \right). \end{eqnarray*} (4) where V is the total number of neighbours |$\mathbf {x}_i$| in the ε-neighbourhood of Ut and S(τ) is termed the stretching factor. The plot of S(τ) against τ should increase with constant gradient before saturating at some value. The region of constant gradient arises from the exponential divergence of the distance between two nearest neighbours while the saturation arises from the finite size of the attractor. The value of this constant gradient |$\frac{d S(\tau )}{d \tau }$| gives an estimate of the maximal Lyapunov exponent λ. Using the results from the previous section we defined the ε-neighbourhood Ut to be points that are located within a distance of ε = e−1 (with regards to the size of the normalized reconstructed attractor) from the reference point |$\mathbf {x}_t$| in N dimensional phase space. The value of |$\frac{d S(\tau )}{d \tau }$| of the region of constant gradient for the harmonic tremor gives the estimate for the maximal Lyapunov exponent to be λ = 0.03 s−1 (Fig. 15). This value of λ is thirty times smaller than the frequency of the fundamental mode (0.9 Hz) of the harmonic tremor and is close to the value of zero, which is the value of the maximal Lyapunov exponent of a limit cycle. This result is consistent with the results for the correlation dimension that the dynamics of the harmonic tremor is predominantly a limit cycle with some small amounts of chaos present. Figure 15. View largeDownload slide Estimation of the maximal Lyapunov exponent of the harmonic tremor in Fig. 9. The plot of the stretching factor S(τ) against relative time τ exhibits a region of constant gradient λ = 0.03 s−1. The vertical dotted lines demarcate the region τ ∈ [300, 1200] used in the fit. Figure 15. View largeDownload slide Estimation of the maximal Lyapunov exponent of the harmonic tremor in Fig. 9. The plot of the stretching factor S(τ) against relative time τ exhibits a region of constant gradient λ = 0.03 s−1. The vertical dotted lines demarcate the region τ ∈ [300, 1200] used in the fit. 8 SURROGATE DATA ANALYSIS OF THE HARMONIC TREMOR If the data used to calculate the correlation dimension and maximal Lyapunov exponent is too short or too noisy, the calculations may result in false positive detection of chaos and non-linearity. For example in previous studies on the climate, when different methodologies were applied to the data, different conclusions were made regarding the positive presence of non-linearity (Nicolis & Nicolis 1984; Grassberger 1986; Grassberger 1987; Nicolis & Nicolis 1987). This ambiguity is due to the fact that fluctuations and oscillations in supposedly chaotic data sets are due to combinations of several factors such as chaos, nonchaotic non-linear determinism as well as noise and linear correlations in both the dynamics as well as in the instruments used (Theiler et al.1992). As the seismogram used to reconstruct the attractor of the harmonic tremor is only 50 s long, we used the method of surrogate data analysis proposed by Theiler et al. (1992) to confirm the presence of non-linearity in the data. As the harmonic tremor time-series displays harmonic spectral peaks, we used the null hypothesis that the observed data come from linearly autocorrelated Gaussian noise, and that all the structure in the data is given by its autocorrelation function or by its Fourier spectrum (Theiler et al.1992). This null hypothesis assumes that the data can be described by an autoregressive-moving-average model, and that the perturbations in the trajectories of the original data arise from Gaussian noise. Using this null hypothesis we created an ensemble of 39 surrogate data time-series by randomizing the phases of the Fourier transform of the original data (Kantz & Schreiber 1997; Henry et al.2012). The phases are randomized by multiplying each complex amplitude of the Fourier transform by eiϕ where ϕ is some random value in the range [0, 2π). The phase randomized Fourier transform is then symmetrized about the Nyquist frequency to ensure that the inverse Fourier transform is real. The surrogate data is the inverse Fourier transform of this phase randomized data, and has the same mean, variance, and Fourier power spectra as the original data. The phase randomization process introduces random fluctuations into the original time-series, making the reconstructed attractors of the surrogate data much more space filling than the original (compare Figs 13 and 16). If both the original and the surrogate data time-series exhibit the same topological values, the null hypothesis that the original data is a result of some well-defined linear process cannot be ruled out and we cannot conclude positively that the data is non-linear. On the other hand if the surrogate data have topological values which are different from that of the original data, we can rule out the null hypothesis and conclude that the original data is non-linear in nature. Figure 16. View largeDownload slide Three dimensional plots of 3 out of the 7 dimensions of the reconstructed attractor for one of the surrogate data sets obtained through Fourier phase randomization. Compared to the reconstructed attractor of the original time-series (Fig. 13), the reconstructed attractor for the surrogate data appears to be more space filling due to the effect of Gaussian noise. Figure 16. View largeDownload slide Three dimensional plots of 3 out of the 7 dimensions of the reconstructed attractor for one of the surrogate data sets obtained through Fourier phase randomization. Compared to the reconstructed attractor of the original time-series (Fig. 13), the reconstructed attractor for the surrogate data appears to be more space filling due to the effect of Gaussian noise. The correlation dimension of each reconstructed attractor in the ensemble of surrogate data was measured and compared to the original. Fig. 17 shows that only the original data exhibits the region of constant gradient characteristic of a deterministic non-linear system in the plots of log C versus log r. Therefore, we can reject the null hypothesis that the data arises from autocorrelated Gaussian noise with a confidence level of |$\frac{39}{39+1}=97.5\hbox{ per cent}$|⁠. We note however, that the rejection of this null hypothesis does not automatically mean that the data is proven to be non-linear in nature (Kantz & Schreiber 1997). Therefore we show in the next section that a non-linear system can generate waveforms similar to those observed at Shinmoedake in order to support the hypothesis that it is possible for long period events and harmonic tremor to share the same non-linear source mechanism. Figure 17. View largeDownload slide Plots of log C versus log r (top figure) and the corresponding gradients (bottom figure) for both the original data (black curve) and the 39 surrogate data (cyan curves). Only the original data exhibits the constant gradient region characteristic of non-linear systems, while all 39 surrogate data do not. Figure 17. View largeDownload slide Plots of log C versus log r (top figure) and the corresponding gradients (bottom figure) for both the original data (black curve) and the 39 surrogate data (cyan curves). Only the original data exhibits the constant gradient region characteristic of non-linear systems, while all 39 surrogate data do not. 9 PLAUSIBILITY OF THE JULIAN TREMOR MODEL AS THE HARMONIC TREMOR SOURCE MECHANISM AND DISCUSSIONS As our current measurements above imply that the harmonic tremor events at Shinmoedake arise from some sort of non-linear process, and since earlier studies have hinted that the source mechanism might involve fluid flow (Ichihara et al.2013; Lyons et al.2013), we qualitatively consider the Julian tremor model (Julian 1994) as a plausible mechanism for the harmonic tremor. The Julian tremor model is a three dimensional lumped parameter model invoking non-linear incompressible fluid flow through an elastic channel given by the following set of coupled equations \begin{eqnarray*} \frac{dv}{dt} &=& \frac{p_1-p_2}{\rho L} - \frac{12 \eta v}{\rho h^2} \nonumber \\ \frac{dh}{dt} &=& \dot{h} \nonumber \\ \frac{d \dot{h}}{dt} &=& \left[ L \left( \frac{p_1+p_2}{2} - \frac{\rho v^2}{2} \right) - k(h-h_0) \right. \nonumber \\ &&\left. - A \dot{h} - \frac{L^3}{12h} \left( \frac{12 \eta }{h^2} - \frac{\rho \dot{h}}{2 h} \right) \dot{h} \right] \left( M + \frac{\rho L^3}{12 h} \right)^{-1} , \end{eqnarray*} (5) where |$\dot{h}$| represents the time derivative of the channel wall displacement h, v is the fluid flow velocity, p1 and p2 are the input and output fluid pressures, ρ is the fluid density, L is the channel length, η is the fluid viscosity, h0 is the channel width equilibrium value, k, A and M are, respectively the stiffness, damping constant and mass of the channel wall. The parameter values used here are listed in Table 2. Table 2. Parameter values used in the Julian Tremor Model. Parameter name Parameter symbol Parameter value Input pressure p1 As stated in text Output pressure p2 105 Pa Fluid density ρ 2500 kg m−3 Channel length L 10 m Fluid viscosity η 50 Pa s Equilibrium channel width h0 1 m Rock stiffness k 0.6 GPa Damping constant A 107 kg m−1 s−1 Rock mass M 3 × 105 kg Parameter name Parameter symbol Parameter value Input pressure p1 As stated in text Output pressure p2 105 Pa Fluid density ρ 2500 kg m−3 Channel length L 10 m Fluid viscosity η 50 Pa s Equilibrium channel width h0 1 m Rock stiffness k 0.6 GPa Damping constant A 107 kg m−1 s−1 Rock mass M 3 × 105 kg View Large Table 2. Parameter values used in the Julian Tremor Model. Parameter name Parameter symbol Parameter value Input pressure p1 As stated in text Output pressure p2 105 Pa Fluid density ρ 2500 kg m−3 Channel length L 10 m Fluid viscosity η 50 Pa s Equilibrium channel width h0 1 m Rock stiffness k 0.6 GPa Damping constant A 107 kg m−1 s−1 Rock mass M 3 × 105 kg Parameter name Parameter symbol Parameter value Input pressure p1 As stated in text Output pressure p2 105 Pa Fluid density ρ 2500 kg m−3 Channel length L 10 m Fluid viscosity η 50 Pa s Equilibrium channel width h0 1 m Rock stiffness k 0.6 GPa Damping constant A 107 kg m−1 s−1 Rock mass M 3 × 105 kg View Large The original autonomous Julian tremor model is able to exhibit one of the main features of the harmonic tremor observed at Shinmoedake: period doubling bifurcation. Figs 18(a) and (b) show respectively the time-series for |$\dot{h}$| before and after period doubling bifurcation. The corresponding phase portraits of h versus |$\dot{h}$| show the transition from a single to a double nested loop (Fig. 19) as with the harmonic tremor seismograms presented earlier (Fig. 7). We note that period doubling bifurcation was also discovered in an experiment where air was forced through a vertical slit cut in a gelatine block (Rust et al.2008). Figure 18. View largeDownload slide Time-series of the channel wall velocity for the Julian tremor model for (a) p1 = 10 MPa, (b) p1 = 14 MPa and (c) |$p_1(t)=30 \cdot 10^6 \cdot e^{-(t-30)^2/0.1}+p_2$|⁠. Between (a) and (b) period doubling bifurcation has occured due to the increase in p1. In (c) the time-series now depicts a long period event like waveform caused by a transient input pressure. The normalized Power Spectral Densities of (a) and (c) are plotted in (d) as solid and dotted lines, respectively, showing that the waveforms of (a) and (c) share a common frequency peak at 4.25 Hz. Figure 18. View largeDownload slide Time-series of the channel wall velocity for the Julian tremor model for (a) p1 = 10 MPa, (b) p1 = 14 MPa and (c) |$p_1(t)=30 \cdot 10^6 \cdot e^{-(t-30)^2/0.1}+p_2$|⁠. Between (a) and (b) period doubling bifurcation has occured due to the increase in p1. In (c) the time-series now depicts a long period event like waveform caused by a transient input pressure. The normalized Power Spectral Densities of (a) and (c) are plotted in (d) as solid and dotted lines, respectively, showing that the waveforms of (a) and (c) share a common frequency peak at 4.25 Hz. Figure 19. View largeDownload slide Phase portraits of channel wall displacement versus channel wall velocity for the Julian tremor model showing period doubling bifurcation (see Fig. 18 for time-series waveforms). (a) When p1 = 10 MPa the system is in the period 1 state exhibiting a single loop, and (b) when p1 = 14 MPa the system is in the period 2 state exhibiting a double nested loop. Figure 19. View largeDownload slide Phase portraits of channel wall displacement versus channel wall velocity for the Julian tremor model showing period doubling bifurcation (see Fig. 18 for time-series waveforms). (a) When p1 = 10 MPa the system is in the period 1 state exhibiting a single loop, and (b) when p1 = 14 MPa the system is in the period 2 state exhibiting a double nested loop. Next we change the original autonomous Julian tremor model into a non-autonomous model by changing p1 from a constant to a time dependent parameter. Firstly, the original autonomous Julian model is unable to exhibit changes in the signal amplitude while keeping the frequency spectra constant. We modify the Julian model by changing the input pressure p1 to a step like increase and decrease, and the resulting channel wall velocity time-series exhibited a corresponding increase and decrease in signal amplitude while keeping a constant frequency spectrum (Fig. 20). This simulated result qualitatively matches the observations made at Shinmoedake where the seismogram amplitude fluctuates despite a constant set of frequency peaks (Fig. 5). Figure 20. View largeDownload slide (a) Time-series for the channel wall velocity of the Julian model, given a step like increase and decrease in input pressure p1 shown in (b). (c) Corresponding spectrogram of the time-series in (a) showing a constant set of frequency peaks despite changing waveform amplitude and parameter p1. Figure 20. View largeDownload slide (a) Time-series for the channel wall velocity of the Julian model, given a step like increase and decrease in input pressure p1 shown in (b). (c) Corresponding spectrogram of the time-series in (a) showing a constant set of frequency peaks despite changing waveform amplitude and parameter p1. Second, we change p1 to the form: |$p_1(t)=30 \cdot 10^6 \cdot e^{-(t-30)^2/0.1}+p_2$|⁠. This is able to produce a long period-like event (Fig. 18c), implying that the same physical non-linear fluid flow source mechanism is able to exhibit both harmonic tremor-like events and long period-like events just by changing the behaviour of p1. Fig. 18(d) shows that both the harmonic tremor-like waveform and the long period-like waveforms from the Julian tremor model share the same frequency peak, similar to the observations made at Shinmoedake (Fig. 2). We also make p1 a linearly decreasing function in time of the form: p1(t) = −6 · 104t + 14 · 106, and the resulting time-series from the Julian model undergoes period doubling bifurcation within the same event (Fig. 21). This is very similar to the observations made at Shinmoedake (Fig. 8) and hints that a change in a parameter of the underlying non-linear system is responsible for the activation/deactivation of the harmonic tremor events. However, we emphasize that these results do not prove that long period events and harmonic tremor at Shinmoedake share the same non-linear source mechanism, and stress again that more analyses of the observed waveforms must be completed before a definitive conclusion can be reached. While the original autonomous Julian tremor model has been studied extensively, to the best of our knowledge the non-autonomous Julian model has not and remains an open field for further research. Figure 21. View largeDownload slide (a) Time-series for the channel wall velocity of the Julian model, given a linearly decreasing input pressure of the form: p1(t) = −6 · 104t + 14 · 106. At the start of the time-series the system is in the period 2 state (15–25 s), and at the end the system has transitioned to the period 1 state (30–35 s). This can be seen in the phase portraits of velocity versus displacement in the (b) period 2 region and in the (c) period 1 region. This behaviour is very similar to the observations of the harmonic tremor made at Shinmoedake as depicted in Fig. 8. Figure 21. View largeDownload slide (a) Time-series for the channel wall velocity of the Julian model, given a linearly decreasing input pressure of the form: p1(t) = −6 · 104t + 14 · 106. At the start of the time-series the system is in the period 2 state (15–25 s), and at the end the system has transitioned to the period 1 state (30–35 s). This can be seen in the phase portraits of velocity versus displacement in the (b) period 2 region and in the (c) period 1 region. This behaviour is very similar to the observations of the harmonic tremor made at Shinmoedake as depicted in Fig. 8. We also note that at Shinmoedake the harmonic tremor was observed to only exist in the period 1 and 2 limit cycle states, whereas according to the Julian tremor model, higher order states (period 4, 8 ... chaos) are possible (see figs 12 and 13 in (Julian 1994)). This apparent lack of higher order states at Shinmoedake could be a result of a few possible reasons. For example certain parameter values in the Shinmoedake harmonic tremor driving system could have simply not been large/small enough to force the limit cycle oscillations into higher order periods. Also in the original Julian tremor model the damping response was defined to be the linear term |$A\dot{h}$| (see eq. 11 in Julian 1994). We note that in studies of vocal fold oscillations of the human body, a non-linear damping term was introduced ad hoc in the vocal fold response in order to account for various non-linear effects such as pressure losses in the glottis, vocal fold collision, formation of a downstream air jet, tissue viscoelasticity, etc. (Titze 1988; Laje et al.2001; Lucero 2005). Using these studies as a guide, we found that if the damping response in the original Julian tremor model |$A\dot{h}$| is changed to the non-linear form |$A(h-h_0)^2\dot{h}$|⁠, the system’s dynamics becomes more restricted and the system progresses from period 1 state to period 2 state and back to period 1 state as p1 is increased (Fig. 22). This shows that the elastic response of the channel wall greatly influences the dynamics of the system and has to be investigated further. Figure 22. View largeDownload slide Bifurcation diagram for the channel wall velocity of the Julian tremor model with non-linear damping term |$A(h-h_0)^2 \dot{h}$| instead of the term |$A\dot{h}$| (eq. 4). Unlike the original model which exhibits higher order period limit cycles and even chaotic motion (see fig. 16 in Julian 1994), this modified system is only able to exhibit up to period 2 limit cycle oscillations for the same range of parameters used. Figure 22. View largeDownload slide Bifurcation diagram for the channel wall velocity of the Julian tremor model with non-linear damping term |$A(h-h_0)^2 \dot{h}$| instead of the term |$A\dot{h}$| (eq. 4). Unlike the original model which exhibits higher order period limit cycles and even chaotic motion (see fig. 16 in Julian 1994), this modified system is only able to exhibit up to period 2 limit cycle oscillations for the same range of parameters used. Based on the results and comparisons above it is our opinion that a non-linear fluid flow model is a more likely explanation of both harmonic tremor and long period events at Shinmoedake compared to those invoking linear resonance or repetitive triggering mechanisms. It is not necessary to invoke resonant models to explain the presence of the harmonic spectral peaks (Julian 2000), resonance in cavities seldom produce continuously more than three or four overtones (Hellweg 1999, 2000), and resonant models require an arbitrary initial condition while the Julian tremor model produces self-sustained oscillations (Julian 1994, 2000). In addition using the same physical model to explain both harmonic tremor and long period events allows us to reduce the entire complexity of the volcanic plumbing system without necessitating different sources for different event types. While we have been able to reproduce some features of the harmonic tremor as well as long period events using the Julian tremor model, we have not been able to reproduce other features such as mode hopping. Also, the Julian tremor model makes certain assumptions in order to simplify the physical problem. For example in the model incompressible fluid flow is assumed, while earlier studies have hinted that the harmonic tremor at Shinmoedake might be driven by gas flow (Ichihara et al.2013; Lyons et al.2013). This means that we have to consider the effects of compressibility, which could drastically affect the dynamics of the system. The model also assumes that the channel walls are perfectly rigid, so the model does not take into account the dynamics of crack waves along the interface of the channel wall and the fluid which might result in certain contributions to the observed waveform and frequency spectra (Chouet 1988). A more refined non-linear fluid flow model is therefore required to fully elucidate the physical processes behind harmonic tremor generation. Additionally, analogue experiments such as those designed by Rust et al. (2008) will allow us to explore phenomena which might be difficult to simulate accurately, and are more accessible than observations. Also, further inspiration for the development of source mechanism models could possibly come from other fields involving non-linear phenomena such as laser physics (Haken 1975), glaciology (Müller et al.2005), the physics of musical instruments (Lesage et al.2006) or voice generation (Titze 1988; Laje et al.2001; Lucero 2005). 10 CONCLUDING REMARKS Multiple episodes of harmonic tremor as well as long period events were observed during the 2011 eruption at Shinmoedake, Kyushu, Japan. Both the harmonic tremor and long period events share a common set of frequency peaks, indicating a common source mechanism for both types of events. We discovered several occurrences of period doubling bifurcation in the harmonic tremor seismogram, indicating that the harmonic tremor at Shinmoedake must arise from some sort of non-linear process. In order to quantify the non-linearity of the harmonic tremor seismogram we then selected a 50-s-long segment of the harmonic tremor seismogram and measured the correlation dimension D and maximal Lyapunov exponent λ of this segment to be D = 1.12 and λ = 0.03 s−1 respectively. Both values indicate clearly that the harmonic tremor is predominantly a limit cycle with some small amounts of chaos present. Due to the short 50-s-long data set used in the non-linear time-series analyses, in order to ensure that we did not have false positive detections of non-linearity we performed surrogate data analysis on the data. As limit cycle oscillations model systems that exhibit self-sustained oscillations, we considered the plausibility of the Julian tremor model as a source mechanism for the harmonic tremor. We showed that the original autonomous model is able to exhibit period doubling bifurcation. When the parameter p1 in the model was changed to a step like increase and decrease, we were able to generate oscillations which had varying amplitudes while maintaining a constant frequency spectrum. Additionally, when p1 was changed to a transient pulse we were able to generate long period-like events. We also showed that non-linear response of the channel walls can drastically affect the dynamics of the system. We propose that future studies should investigate the effects of compressibility, crack waves and non-linear response of the channel walls on the fluid flow. Non-linear time-series analyses and forward simulations of a lumped parameter model will not elucidate fully the source mechanism of the observed harmonic tremor, and more detailed studies will be required. For example more detailed waveform analysis is required to study the seismological source mechanism and to constrain the source geometry. In addition further advancements could be made by combining the study of analogue experiments and non-linear systems in other fields. ACKNOWLEDGEMENTS We thank K. I. Konstantinou for discussions regarding his previous work on non-linear time-series analysis. We thank the SGH foundation for generously providing YN with a scholarship to pursue a Masters degree at the University of Tokyo, during which this study was conducted. We are grateful to the Japan Meteorological Agency for providing data from KITK. Data from the other stations listed in the text were obtained from the Earthquake Research Institute, the University of Tokyo, for which we thank T. Ohminato, J. Oikawa and A. Watanabe. The calculations of the correlation dimensions were run on the EIC system at the Earthquake Research Institute, the University of Tokyo. This study was supported by MEXT (KAKENHI Grant Nos. 16K13872 and 15K13561) and MEXT under its Observation and Research Program for Prediction of Earthquake and Volcanic Eruptions. YN performed the non-linear analyses of the harmonic tremor seismogram, proposed the various modifications to the Julian tremor model and drafted the manuscript. YN and MT discussed about the implementation of the various non-linear time-series analysis algorithms. All three authors YN, MT and MI discussed about the contents of the manuscript, making it more readable to nonspecialists in non-linear dynamics. REFERENCES Albano A.M. , Passamante A. , Farrell M.E. , 1991 . Using higher-order correlations to define an embedding window , Phys. D , 54 , 85 – 97 .. https://doi.org/10.1016/0167-2789(91)90110-U Google Scholar Crossref Search ADS Balmforth N.J. , Craster R.V. , Rust A.C. , 2005 . Instability in flow through elastic conduits and volcanic tremor , J. Fluid Mech. , 527 , 353 – 377 .. https://doi.org/10.1017/S0022112004002800 Google Scholar Crossref Search ADS Bean C.J. , De Barros L. , Lokmer I. , Métaxian J.P. , O’ Brien G. , Murphy S. , 2014 . Long-period seismicity in the shallow volcanic edifice formed from slow-rupture earthquakes , Nat. Geosci. , 7 , 71 – 75 .. https://doi.org/10.1038/ngeo2027 Google Scholar Crossref Search ADS Bercovici D. , Jellinek A.M. , Michaut C. , Roman D.C. , Morse R. , 2013 . Volcanic tremors and magma wagging: gas flux interactions and forcing mechanism , Geophys. J. Int. , 195 ( 2 ), 1001 – 1022 .. https://doi.org/10.1093/gji/ggt277 Google Scholar Crossref Search ADS Broomhead D.S. , Huke J.P. , Muldoon M.R. , 1992 . Linear filters and non-linear systems , J. R. Soc., B (Methodological) , 54 ( 2 ), 373 – 382 .. https://doi.org/10.1111/j.2517-6161.1992.tb01887.x Google Scholar Crossref Search ADS Cao L. , 1997 . Practical method for determining the minimum embedding dimension of a scalar time series , Phys. D , 110 , 43 – 50 .. https://doi.org/10.1016/S0167-2789(97)00118-8 Google Scholar Crossref Search ADS Chouet B. , Shaw H.R. , 1991 . Fractal properties of tremor and gas piston events observed at Kilauea volcano, Hawaii , J. geophys. Res. , 96 ( B6 ), 10177 – 10189 .. https://doi.org/10.1029/91JB00772 Google Scholar Crossref Search ADS Chouet B. , 1985 . Excitation of a buried magmatic pipe: a seismic source model for volcanic tremor , J. geophys. Res. , 90 , 1881 – 1893 .. https://doi.org/10.1029/JB090iB02p01881 Google Scholar Crossref Search ADS 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 ( B5 ), 4375 – 4400 .. https://doi.org/10.1029/JB093iB05p04375 Google Scholar Crossref Search ADS Chouet B.A. , 1996 . Long-period volcano seismicity: its source and use in eruption forecasting , Nature , 380 , 309 – 316 .. https://doi.org/10.1038/380309a0 Google Scholar Crossref Search ADS Fraser A.M. , Swinney H.L. , 1986 . Independent coordinates for strange attractors from mutual information , Phys. Rev. A , 33 ( 2 ), 1134 – 1140 .. https://doi.org/10.1103/PhysRevA.33.1134 Google Scholar Crossref Search ADS Godano C. , Cardaci C. , Privitera E. , 1996 . Intermittent behaviour of volcanic tremor at Mt. Etna , Pure appl. Geophys. , 147 , 729 – 744 .. https://doi.org/10.1007/BF01089699 Google Scholar Crossref Search ADS Godano C. , Capuano P. , 1999 . Source characterization of low frequency events at Stromboli and Vulcano islands , J. Seismol. , 3 , 393 – 408 .. https://doi.org/10.1023/A:1009842411153 Google Scholar Crossref Search ADS Grassberger P. , Procaccia I. , 1983 . Characterization of strange attractors , Phys. Rev. Lett. , 50 , 346 – 349 .. https://doi.org/10.1103/PhysRevLett.50.346 Google Scholar Crossref Search ADS Grassberger P. , 1986 . Do climatic attractors exist? . Nature , 323 , 609 – 612 .. https://doi.org/10.1038/323609a0 Google Scholar Crossref Search ADS Grassberger P. , 1987 . Evidence for climatic attractors , Nature , 326 , 524 . https://doi.org/10.1038/326524a0 Google Scholar Crossref Search ADS Haken H. , 1975 . Analogy between higher instabilities in fluids and lasers , Phys. Lett. A , 53 ( 1 ), 77 – 78 .. https://doi.org/10.1016/0375-9601(75)90353-9 Google Scholar Crossref Search ADS Hellweg M. , 1999 . Listening carefully: unique observations of harmonic tremor at Lascar volcano, Chile , Ann. Geof. , 42 , 451 – 464 . Hellweg M. , 2000 . Physical models for the source of Lascar’s harmonic tremor , J. Volc. Geotherm. Res. , 101 , 183 – 198 .. https://doi.org/10.1016/S0377-0273(00)00163-3 Google Scholar Crossref Search ADS Henry B. , Lovell N. , Camacho F. , 2012 . Nonlinear dynamics time series analysis , in Nonlinear Biomedical Signal Processing, Volume 2: Dynamic Analysis and Modeling , pp. 1 – 39 ., Wiley-IEEE Press . Ichihara M. , Lyons J.J. , Yokoo A. , 2013 . Switching from seismic to seismo-acoustic harmonic tremor at a transition of eruptive activity during the Shinmoe-dake 2011 eruption , Earth Planets Space , 65 , 633 – 643 .. https://doi.org/10.5047/eps.2013.05.003 Google Scholar Crossref Search ADS Jellinek A.M. , Bercovici D. , 2011 . Seismic tremors and magma wagging during explosive volcanism , Nature , 470 ( 7335 ), 522 – 525 .. https://doi.org/10.1038/nature09828 Google Scholar Crossref Search ADS PubMed Johnson J.B. , Lees J.M. , 2000 . Plugs and chugs - seismic and acoustic observations of degassing explosions at Karymsky, Russia and Sangay, Ecuador , J. Volc. Geotherm. Res. , 101 , 67 – 82 .. https://doi.org/10.1016/S0377-0273(00)00164-5 Google Scholar Crossref Search ADS Jordan D.W. , Smith P. , 2007 . Nonlinear Ordinary Differential Equations: An Introduction for Scientists and Engineers , 4th edn , pp. 383 – 385 ., p. 453 , Oxford University Press . Julian B.R. , 1994 . Volcanic tremor: Nonlinear excitation by fluid flow , J. geophys. Res. , 99 ( B6 ), 11859 – 11877 .. https://doi.org/10.1029/93JB03129 Google Scholar Crossref Search ADS Julian B.R. , 2000 . Period doubling and other nonlinear phenomena in volcanic earthquakes and tremor , J. Volc. Geotherm. Res. , 101 , 19 – 26 .. https://doi.org/10.1016/S0377-0273(00)00165-7 Google Scholar Crossref Search ADS Kantz H. , Schreiber T. , 1997 . Nonlinear Time Series Analysis , pp. 91 – 111 ., Cambridge University Press . Kantz H. , 1994 . A robust method to estimate the maximal Lyapunov exponent of a time series , Phys. Lett. A , 185 , 77 – 87 .. https://doi.org/10.1016/0375-9601(94)90991-1 Google Scholar Crossref Search ADS Konstantinou K.I. , Lin C.H. , 2004 . Nonlinear time series analysis of volcanic tremor events recorded at Sangay Volcano, Ecuador , Pure appl. Geophys. , 161 , 145 – 163 .. https://doi.org/10.1007/s00024-003-2432-y Google Scholar Crossref Search ADS Konstantinou K.I. , Perwita C.A. , Maryanto S. , Surono , Budianto A. , Hendrasto M. , 2013 . Maximal Lyapunov exponent variations of volcanic tremor recorded during explosive and effusive activity at Mt. Semeru volcano, Indonesia , Nonlin. Process. Geophys. , 20 , 1137 – 1145 .. https://doi.org/10.5194/npg-20-1137-2013 Google Scholar Crossref Search ADS Kubotera A. , 1974 . Volcanic tremors at Aso volcano , in Physical Volcanology , eds Civetta L. , Gasparini P. , Luongo G. , Rapolla A. , pp. 29 – 47 ., Elsevier Sci Laje R. , Gardner T. , Mindlin G.B. , 2001 . Continuous model for vocal fold oscillations to study the effect of feedback , Phys. Rev. E , 64 , 056201 . https://doi.org/10.1103/PhysRevE.64.056201 Google Scholar Crossref Search ADS Lees J.M. , Gordeev E.I. , Ripepe M. , 2004 . Explosions and periodic tremor at Karymsky volcano, Kamchatka, Russia , Geophys. J. Int. , 158 , 1151 – 1167 .. https://doi.org/10.1111/j.1365-246X.2004.02239.x Google Scholar Crossref Search ADS Lesage P. , Mora M.M. , Alvarado G.E. , Pacheco J. , Métaxian J.P. , 2006 . Complex behavior and source model of the tremor at Arenal volcano, Costa Rica , J. Volc. Geotherm. Res. , 157 , 49 – 59 .. https://doi.org/10.1016/j.jvolgeores.2006.03.047 Google Scholar Crossref Search ADS Liao Y. , Bercovici D. , Jellinek M. , 2018 . Magma wagging and whirling in volcanic conduits , J. Volc. Geotherm. Res. , 351 , 57 – 74 .. https://doi.org/10.1016/j.jvolgeores.2017.12.016 Google Scholar Crossref Search ADS Lucero J.C. , 2005 . Bifurcations and limit cycles in a model for a vocal fold oscillator , Commun. Math. Sci. , 3 ( 4 ), 517 – 529 .. https://doi.org/10.4310/CMS.2005.v3.n4.a3 Google Scholar Crossref Search ADS Lyons J.J. , Ichihara M. , Kurokawa A. , Lees J.M. , 2013 . Switching between seismic and seismo-acoustic harmonic tremor simulated in the laboratory: insights into the role of open degassing channels and magma viscosity , J. geophys. Res.: Solid Earth , 118 , 277 – 289 .. https://doi.org/10.1002/jgrb.50067 Google Scholar Crossref Search ADS Maryanto S. , Iguchi M. , Tameguri T. , 2008 . Constraints on the source mechanism of harmonic tremors based on seismological, ground deformation and visual observations at Sakurajima volcano, Japan , J. Volc. Geotherm. Res. , 170 , 198 – 217 .. https://doi.org/10.1016/j.jvolgeores.2007.10.004 Google Scholar Crossref Search ADS Matsumoto S. , Shimizu H. , Matsushima T. , Uehira K. , Yamashita Y. , Nakamoto M. , Miyazaki M. , Chikura H. , 2013 . Short-term spatial change in a volcanic tremor source during the 2011 Kirishima eruption , Earth Planets Space , 65 , 323 – 329 .. https://doi.org/10.5047/eps.2012.09.002 Google Scholar Crossref Search ADS Müller C. , Schlindwein V. , Eckstaller A. , Miller H. , 2005 . Singing icebergs , Science , 310 , 1299 . Google Scholar Crossref Search ADS PubMed Nakada S. , Nagai M. , Kaneko T. , Suzuki Y. , Maeno F. , 2013 . The outline of the 2011 eruption at Shinmoe-dake (Kirishima), Japan , Earth Planets Space , 65 , 475 – 488 .. https://doi.org/10.5047/eps.2013.03.016 Google Scholar Crossref Search ADS Nakao S. et al. , 2013 . Volume change of the magma reservoir relating to the 2011 Kirishima Shinmoe-dake eruption--charging, discharging and recharging process inferred from GPS measurements , Earth Planets Space , 65 , 505 – 515 . Google Scholar Crossref Search ADS Nicolis C. , Nicolis G. , 1984 . Is there a climatic attractor? . Nature , 311 , 529 – 532 .. https://doi.org/10.1038/311529a0 Google Scholar Crossref Search ADS Nicolis C. , Nicolis G. , 1987 . Evidence for climatic attractors , Nature , 326 , 523 . https://doi.org/10.1038/326523b0 Google Scholar Crossref Search ADS Ripepe M. , Gordeev E. , 1999 . Gas bubble dynamics model for shallow volcanic tremor at Stromboli , J. geophys. Res. , 104 ( B5 ), 10639 – 10654 .. https://doi.org/10.1029/98JB02734 Google Scholar Crossref Search ADS Rust A.C. , Balmforth N.J. , Mandre S. , 2008 . The feasibility of generating low-frequency volcano seismicity by flow through a deformable channel , in Fluid Motions in Volcanic Conduits: A Source of Seismic and Acoustic Signals , Vol. 307 , pp. 45 – 56 ., eds Lane S.J. , Gilbert J.S. , Geological Society , Special Publications , London . Sakuraba A. , Yamauchi H. , 2014 . Linear stability of plane Poiseuille flow in an infinite elastic medium and volcanic tremors , Earth, Planets Space , 66 , 19 . https://doi.org/10.1186/1880-5981-66-19 Google Scholar Crossref Search ADS Saulys B. , Matukas J. , Palenskis V. , Pralgauskaite S. , Vysniauskas J. , 2010 . Analysis of mode-hopping effect in Fabry-Perot Laser Diodes , Paper presented at 18th International Conference on Microwaves , Radar and Wireless Communications , Vilnius, Lithuania . Schlindwein V. , Wassermann J. , Scherbaum F. , 1995 . Spectral analysis of harmonic tremor signals at Mt. Semeru volcano, Indonesia , Geophys. Res. Lett. , 22 ( 13 ), 1685 – 1688 .. https://doi.org/10.1029/95GL01433 Google Scholar Crossref Search ADS Strogatz S.H. , 1994a . in Nonlinear Dynamics And Chaos: With Applications To Physics, Biology, Chemistry And Engineering , pp. 168 – 171 ., Perseus Books . Strogatz S.H. , 1994b . in Nonlinear Dynamics And Chaos: With Applications To Physics, Biology, Chemistry And Engineering , pp. 203 – 210 ., pp. 377–378 , Perseus Books . Takens F. , 1981 . Detecting Strange Attractors in Turbulence , Lecture Notes Math. , 898 , 366 – 381 . Google Scholar Crossref Search ADS Takeo M. , Maehara Y. , Ichihara M. , Ohminato T. , Kamata R. , Oikawa J. , 2013 . Ground deformation cycles in a magma-effusive stage, and sub-Plinian and Vulcanian eruptions at Kirishima volcanoes, Japan , J. geophys. Res.: Solid Earth , 118 , 4758 – 4773 .. https://doi.org/10.1002/jgrb.50278 Google Scholar Crossref Search ADS Theiler J. , Eubank S. , Longtin A. , Galdrikian B. , Farmer J.D. , 1992 . Testing for nonlinearity in time series: the method of surrogate data , Phys. D , 58 , 77 – 94 .. https://doi.org/10.1016/0167-2789(92)90102-S Google Scholar Crossref Search ADS Theiler J. , 1990 . Estimating the fractal dimension of chaotic time series , Lincoln Lab. J. , 3 ( 1 ), 63 – 86 . Titze I.R. , 1988 . The physics of small-amplitude oscillation of the vocal folds , J. acoust. Soc. Am. , 83 , 1536 – 1552 .. https://doi.org/10.1121/1.395910 Google Scholar Crossref Search ADS PubMed Ueda H. , Kozono T. , Fujita E. , Kohno Y. , Nagai M. , Miyagi Y. , Tanada T. , 2013 . Crustal deformation associated with the 2011 Shinmoe-dake eruption as observed by tiltmeters and GPS , Earth Planets Space , 65 , 517 – 525 .. https://doi.org/10.5047/eps.2013.03.001 Google Scholar Crossref Search ADS APPENDIX: CAO’s NEAREST NEIGHBOURS METHOD According to Cao’s nearest neighbours method (Cao 1997) , the ratio of the ‘nearest neighbour distance’ between dimensions N and N + 1 is defined as \begin{eqnarray*} a(t,N)= \frac{||x_t(N+1) - x_{n(t,N)}(N+1) ||}{|| x_t(N) - x_{n(t,N)}(N) ||} , \end{eqnarray*} (A1) where || · || is the Euclidean distance given by the maximum norm \begin{eqnarray*} ||x_k(m)-x_l(m)|| = \max _{0 \le j \le m-1}|y_{k+j \tau } - y_{l+j \tau }| , \end{eqnarray*} (A2) where xt(N + 1) is the tth reconstructed vector with dimension N + 1, n(t, N) is an integer such that xn(t, N)(N) is the nearest neighbour of xt(N) in the N-dimensional reconstructed phase space. By this construction, for some embedding dimension N, two points which are close together in N-dimensional phase space will still remain close together in (N + 1)-dimensional phase space. Points that behave in such a manner are termed true nearest neighbours, while all other points that become far away when the dimension is increased are termed false nearest neighbours. For minimum embedding dimension N, no more false neighbours exist and the reconstructed attractor is fully unfolded in phase space. The presence of false nearest neighbours is determined by comparing the calculated value of a(t, N) with some threshold value which introduces certain complications such as the correct choice for the threshold. This threshold value depends on the dimension N as well as on time t along the trajectory of the reconstructed attractor. Additionally, threshold values may differ from different time-series arising from the same non-linear system. Therefore choosing an appropriate threshold value is arbitrary and difficult. In order to overcome this problem, Cao introduced the method of using the mean of all the values of a(t, N) \begin{eqnarray*} E(N) = \frac{1}{M-NT} \sum _{i=1}^{M-NT} a(t,N), \end{eqnarray*} (A3) where E(N) is dependent only on the time delay T and the embedding dimension N. The determination of true from false nearest neighbours is done by comparing the change of E(N) as the dimension is increased from N to N + 1 \begin{eqnarray*} E1(N)=\frac{E(N+1)}{E(N)}. \end{eqnarray*} (A4) Cao found that the value of E1(N) starts to increase from N = 0 and then saturates after some value N0, where N0 + 1 is the minimum embedding dimension to be used for the reconstruction of the attractor. In order to distinguish deterministic signals from stochastic noise, Cao defined another quantity \begin{eqnarray*} E^*(N)= \frac{1}{M-NT} \sum _{i=1}^{M-NT}|y(t+NT) - y(n(t,N)+NT)| , \end{eqnarray*} (A5) where as defined above n(t, N) is an integer where the nearest neighbour of xt(N) is xn(t, N)(N). Again we are interested in how this quantity changes from N to N + 1 dimensions so we evaluate \begin{eqnarray*} E2(N)=\frac{E^*(N+1)}{E^*(N)}. \end{eqnarray*} (A6) Ideally a stochastic signal will not result in E1(N) reaching a saturation value as N is increased. However, in reality if the gradient is very small it might be difficult to determine whether E1(N) is still increasing or if it has already reached its saturation value. Additionally, for a stochastic time-series with a limited number of points, the value of E1(N) might stop changing after a certain value of N. E2(N) is on the other hand able to distinguish between a stochastic and deterministic time-series. For stochastic data E2(N) saturates at value 1 for all values of N. In constrast if the data contains a deterministic signal, E2(N) will depend on N and there exists a range of N such that E2(N) is smaller than 1. © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - A non-linear time-series analysis of the harmonic tremor observed at Shinmoedake volcano, Japan
JF - Geophysical Journal International
DO - 10.1093/gji/ggy522
DA - 2019-03-01
UR - https://www.deepdyve.com/lp/oxford-university-press/a-non-linear-time-series-analysis-of-the-harmonic-tremor-observed-at-cKVrHP4Ah2
SP - 1768
VL - 216
IS - 3
DP - DeepDyve
ER -