Maximal Lyapunov exponent variations of volcanic tremor recorded during explosive and effusive activity at Mt Semeru volcano , Indonesia

We analyze 25 episodes of volcanic tremor recorded from 22 November until 31 December 2009 at Mt Semeru volcano in order to investigate their spectral and dynamical properties. The overtone frequencies for most of the tremor events indicate a pattern of period-doubling, which is one possible route that can lead a system to chaotic behavior. Exponential divergence of the phase space orbits is a strong indicator of chaos and was quantified by estimating the maximal Lyapunov exponent (MLE) for all tremor events. MLEs were found to vary linearly with the number of frequency overtones present in the tremor signals. This implies that the tremor source at Semeru fluctuates between a quasi-periodic state with few overtone frequencies (2–3) and small MLEs (~0.013), and a chaotic one with more overtones (up to 8) and larger MLEs (up to 0.039). These results agree well with the tremor generation model suggested previously by Julian (1994), which describes wall oscillations of a crack excited by unsteady fluid flow. In this model, as fluid pressure increases, a period-doubling cascade leads to numerous new frequencies and a chaotic tremor signal. The temporal variation of MLEs exhibited significant fluctuations from 23 until 31 December when the eruptive activity shifted from explosive to effusive. Such a situation may reflect variable fluid pressure conditions inside the conduit, where at first magma is accumulated and subsequently is erupted, releasing the buildup of pressure. Our results give further evidence for the role of nonlinear deterministic processes in generating volcanic tremor and call for similar investigations to be conducted in other volcanoes.


Introduction
A key issue that is highlighted during the analysis of time series is what kind of process is generating this observed signal.Linear stochastic processes with infinite degrees of freedom have been used in the past to explain irregular or broadband signals, and their Fourier spectrum was believed to be the only tool needed in order to characterize them (Hilborn, 1994).Nonlinear deterministic processes famously documented by Lorenz (1963) for the first time offer a plausible alternative to linear models and suggest that a much richer structure may exist in the data.This structure takes the form of orbits moving in an m-dimensional Euclidean space that represent the evolution of the states of the system under study.For systems that exhibit energy dissipation the orbits always move within a bounded region of the phase space forming an attractor.If sensitivity to initial conditions exists, orbits diverge exponentially from each other and the attractor becomes a fractal geometrical object referred to as strange or chaotic attractor.During the last two decades a whole new discipline of time series analysis has been established that provides methods to detect and exploit such deterministic structures in the data (for recent overviews see Kantz andSchreiber, 2004, andSmall, 2005).
Volcanic tremor is a persistent seismic signal that is observed near active volcanoes, lasting from several minutes to several days preceding and/or accompanying most volcanic eruptions (Konstantinou and Schlindwein, 2003).Julian (1994Julian ( , 2000) ) has suggested that any tremor generation process must necessarily be nonlinear, basing this claim on the fact that a time-invariant linear system cannot oscillate spontaneously, as its output cannot contain frequencies that were not present in its input.In this sense Julian argued that tremor is most likely an oscillatory response to a steady, nearzero frequency input such as stress, fluid pressure or heat.Observational evidence for the existence of such a nonlinear process can be found in the systematic relationship between the frequency and amplitude of tremor signals as documented in Karimsky and Sangay volcanoes (Lees et al., 2004;Lees and Ruiz, 2008).On the other hand, several studies have applied nonlinear time series analysis methods to tremor recorded at volcanoes worldwide in order to detect the presence of deterministic structure (Chouet and Shaw, 1991;Godano et al., 1996;Carniel, 1996;Carniel and Di Cecca, 1999;Konstantinou, 2002;Konstantinou and Lin, 2004;De Lauro et al., 2008;Cannata et al., 2010).Their results revealed the existence of a strange attractor that describes the tremor source process in phase space and is characterized by a relatively small number of degrees of freedom.
Mt Semeru is a 3676 m-high stratovolcano that lies on the eastern part of Java island in Indonesia (Fig. 1) and is one of the most active volcanic centers on the island.The composite andesitic cone of Semeru consists of two smaller edifices, namely the Mahameru (old) and the Seloko (younger) edifices (Thouret et al., 2007).Based on the Smithsonian Institution global volcanism database (http://www.volcano.si.edu),Semeru has produced numerous eruptions during the last 200 yr, with the majority of them having a Volcanic Explosivity Index (VEI) between 2 and 3. Eruptive activity is usually of vulcanian type generating several short-lived eruption columns within one day, as well as less frequent lava dome growth and collapse resulting in block-and-ash flows.The densely populated area around Mt Semeru is exposed, due to this activity, to a number of volcanic hazards such as pyroclastic flows and rainfall-triggered lahars.In this work, we analyze volcanic tremor signals recorded at Mt Semeru during late November-December 2009 in order to derive dynamical properties and their variations.First, we give an overview of the available waveform data and their spectral characteristics in terms of a time-frequency representation.Embedding parameters (delay time, embedding dimension) are then determined for each tremor signal and the phase space is reconstructed.The maximal Lyapunov exponent (thereafter called MLE) is a strong indicator of chaos and is estimated with the aim of investigating its possible variations and relationship with the physical mechanism that is causing tremor.

Data and spectral analysis
The Center for Volcanology and Geological Hazard Mitigation (CVGHM) of Indonesia has been monitoring the seismic activity at Mt Semeru using five stations (Fig. 1).These stations are equipped with short-period sensors having only vertical components recording at a sampling inter- val of 100 samples per second, while absolute timing is provided by GPS receivers.Apart from the instrumental recordings, CVGHM has also set up the Gunung Sawur volcano observatory about 10 km from the volcano summit for the purpose of obtaining visual observations of the ongoing activity.We analyze volcanic tremor recorded during a period when a plug of solidified viscous lava had developed at the summit vent and was followed by phreatic and/or steam explosion activity.We select tremor waveforms for analysis only if they conform to the following criteria: (a) there are no gaps due to transmission problems, and (b) the duration of each tremor episode is long enough (> 200 s), since the methods applied in this study are quite sensitive to the total number of samples in the time series.In total 25 volcanic tremor episodes spanning the period from 22 November until 31 December 2009 were found to fulfill these criteria.Additionally, in order to focus only on source characteristics and avoid as much as possible propagation effects through the highly heterogeneous volcanic edifice, we chose to analyze waveforms recorded only at the closest station to the active vent.Unfortunately, station PCK, which was installed next to the active Seloko crater, was hit by lightning and was disabled, therefore we performed our analysis using data from station KLP installed 1 km away from the vent.
Volcanic tremor usually starts with a high-amplitude explosion signal that gradually becomes quasi-periodic, decaying slowly with time (Fig. 2).We constructed spectrograms of each tremor episode using the fast Fourier transform in order to investigate the frequency content variation as a function of time by applying an overlapping window with a length of 10.24 s shifted by 0.5 s.The high-amplitude explosion signal at the beginning exhibits high frequency content that later switches to well-defined peaks consisting of a fundamental  1).First and last parts of the sequence are signals related to eruptive activity having a broader spectrum.The inset on top shows an enlarged section of the tremor waveform delineated by the red lines.
frequency and numerous overtones.Significant shifting of these peaks to higher or lower frequencies (often termed as "gliding") has not been observed in the tremor events under study, probably implying a rather stationary source process.For the purpose of quantifying the number of overtones in each tremor event we count frequency peaks having heights larger than 5 % of the maximum amplitude in the spectrum (Maryanto et al., 2008).Figure 3 shows the number of overtones defined in this way for each of the 25 tremor events selected, where it can be seen that this number varies between 2 and 8 while the fundamental frequency takes values between 0.5 and 2 Hz.These observations are similar to the ones reported earlier by Schlindwein et al. (1995) for the Semeru tremor where the authors noted 2-11 overtones and fundamental frequencies in the range of 0.5-1.8Hz.

Estimation of embedding parameters
The delay embedding theorem first suggested by Takens (1981) and later mathematically proven by Sauer et al. (1991) is applied for the phase reconstruction from the tremor time series.This theorem states that a series of scalar measurements s(t) can be used in order to define the orbits describing the evolution of the states of the system in an m-dimensional Euclidean space.The orbits will then consist of points x with co-ordinates as follows: where τ is referred to as the delay time and for a digitized time series (such as a seismogram) a multiple of the sampling interval is used, while m is the embedding dimension.
There are two methods for estimating the delay time required by the embedding theorem from an observed time series.The first is to calculate the autocorrelation function of the data and select as τ the time of its first zero crossing.The second method makes use of a nonlinear autocorrelation function called average mutual information (AMI) that is defined as (Fraser and Swinney, 1986) where pi is the probability that the signal s(t) will take a value inside the ith bin of a histogram and pij is the probability that s(t) is in bin i and s(t + τ ) is in bin j .The first minimum in the AMI graph is considered as the most suitable choice for delay time, since this is the time when s(t +τ ) adds maximum information to the knowledge we have from s(t).
For the Semeru volcanic tremor events we calculated both the autocorrelation function and AMI using time lags in the range of 0-50 (0-0.5 s), and an example can be seen in Fig. 4.Both methods give consistently the same values of the delay time (between 7 and 11) for all tremor events, which is in accordance with similar studies performed previously (e.g., Konstantinou and Lin, 2004).In order to cross-check this result we also constructed two-dimensional phase portraits for different delay times with the aim of examining the resulting orbit geometry.It is expected that if the delay time is suitable for embedding, the orbits should not be aligned along the diagonal of the phase portrait.We find that the delay time values determined earlier result in phase portraits with orbits that are spread around the diagonal (Fig. 5).
The sufficient embedding dimension is estimated using the false nearest neighbors method developed by Kennel et al. (1992).The basic principle that the method is founded on is that the orbits of a chaotic attractor should not intersect or overlap with each other.Such an intersection or overlap may result when the attractor is embedded in a dimension lower than the sufficient one needed by the delay embedding theorem.In this way, two points are considered as true neighbors only if the dynamic evolution of the orbits brought them  close, and they are termed as false neighbors if they are close because of the projection of the attractor to a lower dimension.By comparing the Euclidean distance of two points in two consecutive embedding dimensions d and d +1 it is possible to decipher between these two cases.If the ratio of these distances is greater than a predefined value s, then the two points are false neighbors.An additional criterion for characterizing two points as false neighbors is that their distance in dimension d + 1 should be larger than σ/s where σ is the standard deviation of the data around the mean.Since σ can be thought of as a representative measure of the size of the attractor, the latter criterion reflects the fact that if two points are false neighbors they will be stretched to the extremities of the attractor when they are separated at dimension d + 1.This procedure is repeated for all pairs of points at higher dimensions until the percentage of false neighbors becomes zero, and then the attractor is unfolded.We use values for s in the range of 9-15 and find that there are no significant differences in the resulting embedding dimensions, in agreement with previous studies (e.g., Konstantinou, 2002).An example of the results obtained can be seen in Fig. 6, while a summary of the embedding parameters is given in Table 1.The embedding dimension for the tremor events considered varies from a minimum value of 5 to a maximum value of 7.

Estimation and characteristics of MLEs
An important indicator of chaotic behavior is the exponential divergence of nearby orbits in the phase space, which is expressed through the Lyapunov exponents.Chaotic processes therefore have at least one finite and positive Lyapunov exponent: for linear processes this exponent is zero, while it becomes infinite for stochastic processes.In general, for an m-dimensional phase space the rate of expansion or contraction of orbits is described for each dimension by one Lyapunov exponent, resulting in m different Lyapunov exponents of which some are zero or negative.It should be noted that for practical applications the main interest is focused on the largest of these exponents, as it can be readily estimated from time series data.Towards this end, we apply the method suggested by Rosenstein et al. (1993) for calculating the MLE for our tremor events.After reconstructing the phase space using suitable values for τ and m, a point x n0 is chosen, all neighbor points x n closer than a distance ε are found and their average distance from that point is calculated.This is repeated for N number of points along the orbit so as to calculate an average quantity S known as the stretching factor where u xno is the number of neighbors found around point x n0 .A plot of the stretching factor S versus the number of points N (or the equivalent time t = N δt) will yield a curve with a linear increase at the beginning, followed by a flat region.The first part of the curve represents the exponential increase of S as more points from the orbit are included, while the flat region signifies the saturation effect of exponential divergence due to the finite size of the attractor.A least-squares fit for the slope of the linear part of the curve would then give an estimate of the MLE.We applied the method described above to the 25 tremor events using the embedding parameters determined previously.The minimum value of ε was taken as the data interval divided by 1000, while in order to avoid temporal correlations among consecutive points in the phase space, we use a "Theiler window" (Theiler, 1990) of 500 as recommended by Kantz and Schreiber (2004).The curves for all events exhibited clear linear increase regions with some fluctuations superimposed on the linear part (Fig. 7).These fluctuations can be explained based by the fact that the stretching factor is just an average of the local stretching or shrinking rates in the attractor, therefore these different rates may not be always smoothed by the averaging process of the algorithm.For each curve of the stretching factor we first define the linear part by including only the portion of the curve that maximizes the correlation coefficient r; then we obtain the slope of the line after a least-squares fit using a maximum likelihood estimator (Press et al., 1992).The correlation coefficient was found in most cases higher than 0.85, except from two events that also exhibit the smallest MLEs.For the purpose of cross-checking whether some curves may exhibit a power law rather than exponential behavior, we plotted each curve in a doubly logarithmic scale.Again, we find clear linear increase for all tremor events and we conclude that the exponential divergence is a real feature of the observed time series.Table 2 lists the estimated MLE for each tremor event along with the corresponding error of the least-squares fit.
The relationship of the MLE and the number of overtones present in each of the tremor events can be seen in Fig. 8.In general, there seems to be a linear relationship between the two quantities, implying that the addition of more overtone frequencies results in an increase in the MLE.However, more observations are clearly needed in order to validate this relationship.A more careful look at the spectral analysis results (cf.Fig. 3) reveals that for every three consecutive frequencies the one in the middle lies midway from the other two (except from events #4 and #7, which have only two frequencies).This pattern resembles the well-known perioddoubling process that is one of the routes a system may take in order to reach chaotic behavior (Hilborn, 1994).It should be noted that Julian (2000) had recognized a similar pattern in the spectral analysis results of Schlindwein et al. (1995) for the Semeru tremor.We also constructed the temporal variation of MLEs within the period spanned by  the selected tremor events (Fig. 9).From 22 November (Julian day 326) until 23 December (Julian day 357) the MLEs exhibit fluctuating values ranging between 0.013 and 0.030, and the highest MLE value (∼ 0.039) is registered in the afternoon of 23 December.For the next four days a significant drop is observed, only to be followed by a sharp increase again on 31 December (Julian day 365).Figure 9 also shows the temporal distribution of the number of eruptions per day during the same period when tremor events were occurring.
It should be noted that visual observations start on 4 December, hence there are no observations before this date.Eruptive activity seems to be most intense at the time when the MLE values exhibit fluctuations, while it almost ceases after 23 December when MLE values have a sharp reduction.
All of these eruptions produced visible white smoke with an estimated column height of 100-200 m.Unfortunately, low visibility due to the fog that had covered the Semeru edifice did not allow the obtaining of further visual observations of the eruptive activity during early January 2010.

Discussion
Volcanic tremor signals similar to the ones recorded at Mt Semeru have been observed at several other volcanoes worldwide such as Sangay in Ecuador (Johnsson and Lees, 2000), Arenal in Costa Rica (Hagerty et al., 2000), Karimsky in Kamtchatka (Lees et al., 2004) or Sakurajima in Japan (Maryanto et al., 2008).A common characteristic of all these volcanoes is the development in the top part of their conduit of a plug, which is a solidified portion of magma that was not erupted.The plug effectively seals the conduit, prohibiting the escape of gas that is progressively accumulating below it.Once the gas pressure inside the conduit exceeds a critical threshold, an explosion occurs that may destroy the plug completely, propelling material out of the active crater.Alternatively, the gas explosion may only damage the plug by creating fractures that can play the role of pathways for the gas flow.Volcanic tremor signals are in this respect the oscillatory response of these fractures as the gas flows through them.After the gas has escaped and the pressure inside the conduit has been lowered, the fractures may close until gas pressure in the conduit increases again and the cycle begins anew.Many authors have pointed out the similarity of the fractured plug with the reed of a wind instrument like the clarinet (Lesage et al., 2006), or with the valve that intermittently releases steam in a pressure cooker (Johnson and Lees, 2000).The visual observations provided by the CVGHM observatory, which lies a few kilometers away from the active vent, seem to corroborate a tremor source mechanism such as the one described above.The formation of a plug closing the active conduit at Seloko crater was observed during November 2009.After that several explosions occurred every day, followed by visible white smoke and the recording of volcanic tremor without any audible sound being heard at the CVGHM observatory site.The white smoke implies that gas rather than ash escaped into the atmosphere, therefore the magma flow rate in the conduit did not change significantly.Until the end of December 2009 there were no observations of pyroclastic flows, lava incandesence or rockfalls originating near the crater.The number of explosions and white smoke columns decreased considerably after 23 December to only one or two in one day.Finally, in the morning hours of 31 December a lava tongue was observed on the flanks of the volcano at a distance of 200 m from the Seloko crater rim.
The existence of finite and positive MLEs in volcanic tremor signifies that a nonlinear deterministic process may play a dominant role in the generation of these signals.The variations of MLEs also demonstrate that the system essentially fluctuates between two states, an almost quasi-periodic (with few overtones, small MLE values) and a chaotic one (more overtones added, larger MLE values).Such variations agree rather well with the predictions of the tremor generation model suggested by Julian (1994)  planar crack excited by unsteady flow.The crack may connect two fluid filled reservoirs, or alternatively one fluid reservoir with the Earth's surface.Numerical simulations of this model show that as the pressure inside the reservoir gradually increases, the character of the wall oscillations changes through a period-doubling cascade from periodic to quasiperiodic and finally to chaotic by the progressive addition of frequencies to the tremor signal.At Semeru it is possible that a fracture in the plug could have functioned as the fluidtransporting crack, while the conduit below the plug would represent the reservoir.In this context, we interpret the fluctuations observed prior to 23 December as fluid pressure variations indicating cycles of fluid accumulation beneath the plug and its partial escape into the atmosphere.Accumulation of fluid is probably faster than its removal, therefore there is a peak in MLE values during 23 December followed by a significant decrease on 26 December.This sudden change can be interpreted as fluid pressure decrease when the activity shifted to effusive mode; lava slowly filled the Seloko crater and finally became visible in the early morning of 31 December.In the meantime, the lava covering the crater was partially solidified, prohibiting any further escape of fluid.This leads once more to a buildup of fluid pressure inside the conduit and a new increase in the value of MLEs, setting the ground for a new period of explosive activity.

Conclusions
The main conclusions of this work can be summarized as follows: 1. Volcanic tremor signals recorded at Mt Semeru during November-December 2009 consist of 2-8 overtone frequencies and a fundamental frequency between 0.5 and 2 Hz without exhibiting gliding effects.The overtone frequencies for some of the tremor events indicate a pattern of period-doubling, which is one possible route to chaotic behavior.
2. The tremor signals can be embedded in a Euclidean phase space with a dimension that varies between 5 and 7, and exhibit a variable degree of exponential divergence of nearby orbits as quantified by the MLEs.The MLEs increase as more frequencies are included in the tremor signal and show a significant fluctuation when the eruptive activity shifts from explosive to effusive.
3. The tremor source model suggested by Julian (1994) can be used to explain these observations, since it predicts that increasing fluid pressure within a conduit will result in tremor signals with more overtone frequencies, eventually leading to chaotic behavior through a period-doubling cascade.The physical observations at Semeru indicate the existence of a plug where fluid can escape from the conduit through a fracture.In this context, as fluid pressure increases below the plug, the character of the tremor signal changes by the inclusion of more frequencies and the shifting from quasi-periodic to chaotic behavior.This process may be reversed when there is a reduction in fluid pressure when fluids escape into the atmosphere.
Finally, it should be mentioned that the results presented here point to the necessity for conducting a similar analysis to other volcanic tremor data sets that exhibit similar transitions from quasi-periodic to broadband behavior and vice versa.This would help to further constrain the relationship between the number of overtones and MLEs, as well as the temporal variation of MLEs with respect to the eruptive activity.

Fig. 1 .
Fig. 1.Topographic map of Mt Semeru volcano showing the locations of seismic stations.The inset at the bottom of the image shows the location of Semeru on the island of Java.

Fig. 2 .
Fig. 2. Vertical velocity waveform and corresponding spectrogram for event #20 that was recorded on 23 December 2009 (see Table1).First and last parts of the sequence are signals related to eruptive activity having a broader spectrum.The inset on top shows an enlarged section of the tremor waveform delineated by the red lines.

Fig. 3 .
Fig. 3. Diagram showing the number of different frequency overtones for each volcanic tremor event under study.

Fig. 4 .
Fig. 4. Example of delay time selection for event #25.(a) shows the variation of average mutual information (AMI) as a function of time lag.The arrow highlights the first minimum of AMI.(b) shows the variation of the autocorrelation function as the time lag varies.The arrow shows the first zero crossing.

Fig. 5 .
Fig. 5. Two-dimensional phase portraits of the volcanic tremor attractor for event #25 reconstructed using different delay times.

Fig. 6 .
Fig. 6.Diagram showing the distribution of the percentage of false nearest neighbors as a function of the embedding dimension for event #25.
Fig. 7. (a) Example of the estimation of MLE using the method of Rosenstein et al. (1993) for event #5.(b) The same plot as above; only the x axis is logarithmic (see text for more details).The portion of the curve used for the least-squares line fit starts from the beginning until the saturation point shown by the blue arrow.The red line represents the resulting least-squares line fit.

Table 2 .
Summary of the results of the estimation of maximal Lyapunov exponent.The correlation coefficient r refers to the linear part of the stretching factor curve.

Fig. 8 .
Fig. 8. Diagram showing the relationship between the number of frequency overtones and the values of MLE.The red line indicates the result of the linear regression that was applied to the data.The correlation coefficient R, the slope and intercept along with their error estimates are shown in the upper left hand corner of the plot.
Fig. 9. (a) Diagram showing the temporal variation of MLEs and its relationship with the number of overtones for the tremor events under study.The vertical bars indicate the corresponding error of each MLE estimate.Number of frequency overtones for each tremor event is also shown using a color scale.(b) Temporal variations of number of eruptions per day for the same time period as registered by the Gunung Sawur volcano observatory.Visual observations of eruptions start from 4 December (Julian day 338) until 31 December (see text for more details).

Table 1 .
Estimated delay time (τ ) and embedding dimension (m) for each of the 25 volcanic tremor episodes under study.