© Author(s) 2013. CC Attribution 3.0 License. Nonlinear Processes in Geophysics

Abstract. We analyse the seismic noise recorded at the Colima Volcano (Mexico) in the period December 2005–May 2006 by four broadband three-component seismic stations. Specifically, we characterize the spectral content of the signal and follow its time evolution along all the data set. Moreover, we infer the properties of the attractor in the phase space by false nearest neighbours analysis and Grassberger–Procaccia algorithm, and adopt a time-domain decomposition method (independent component analysis) to find the basic constituents (independent components) of the system. Constraints on the seismic wavefield are inferred by the polarization analysis. We find two states of the background seismicity visible in different time-intervals that are Phase A and Phase B. Phase A has a spectrum with two peaks at 0.15 Hz and 0.3 Hz, with the latter dominating, an attractor of correlation dimension close to 3, three quasi-monochromatic independent components, and a relevant fraction of crater-pointing polarization solutions in the near-field. In Phase B, the spectrum is preserved but with the highest peak at 0.15 Hz, the attractor has a correlation dimension close to 2, two independent components are extracted, and the polarization solutions are dominated by Rayleigh waves incoming from the southwest direction. We depict two sources acting on the background seismicity that are the microseismic noise loading on the Pacific coastline and a low-energy volcanic tremor. A change in the amplitude of the microseismic noise can induce the switching from a state of the system to the other.


Introduction
In recent years, the attention of the scientific community to the seismic noise has been increasing and many evidences are emerging about the existence of coherent signals within it.The analysis of the seismic noise is mostly based on the estimate of the cross-correlation function between stations of seismic arrays.This approach has been applied both for detecting low-energy travelling waves and for tomography purposes; in particular, it has led to correlate modifications of rheological properties of the medium to activation of volcanic phenomena or local and teleseismic earthquakes (e.g., Friedrich et al., 1998;Obara, 2002;Shapiro and Campillo, 2004;Almendros et al., 2007;Brenguier et al., 2008;Ibañez et al., 2008;Harmon et al., 2012;Minato et al., 2012).
Here we propose a different approach to analyse the seismic noise recorded at the Colima Volcano, Mexico, in order to constrain the active source processes and to face the open issue on the existence (or not) of volcanic tremor.We look for the basic constituents of the signals by a statistical technique working in time-domain and characterize the underlying dynamical system.Our approach is suitable for nonlinear systems, such as those displaying self-sustained oscillations, for which source and medium actively interact and cannot be distinguished.Active volcanic areas belong to this class of systems, as the seismicity is induced by the complex interaction between the flow of the magmatic fluids and the rock, developing feedback mechanisms (Julian, 1994;De Lauro et al., 2008).
Despite the complexity of the source mechanisms, volcanic signals can be usually described as the effect of Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.
Volcanic activity at the Colima Volcano is characterized by eruptive cycles consisting in the extrusion of lava flow followed by the growth of a lava dome and its subsequent destruction accompanied by large explosions.It shows a wide range of eruption styles, like dome growth, lava flow accompanied by frequent small block avalanches, phreatic explosions, explosive events and pyroclastic flows (Luhr, 2002;Petrosino et al., 2011).
We analyse a period of about five months starting about two months after the end of an eruptive cycle (September, 2005) and dominated by the occurrence of vulcanian-style explosions.Most vulcanian-style explosions are preceded by a high-energy transient seismic event with principal frequency peak at about 0.4 Hz and over-modes.This transient event is generated in the crater area and is induced by strong and rapid degassing flux interacting with the volcanic conduit and leading to the eruption (Zobin et al., 2006(Zobin et al., , 2010;;Palo et al., 2009;Petrosino et al., 2011;De Lauro et al., 2012a).Such pressurization (or long-period) events dominate the seismicity of active volcanic areas and share with the volcanic tremor very similar spectral content and wavefield properties, as the effect of a common source mechanism (e.g., Chouet, 1996;De Lauro et al., 2008).
In particular, De Lauro et al. (2012a) decomposed longperiod events at the Colima Volcano in few nonlinear fundamental signals associated with the vibration modes of two nearly orthogonal conduits displaying self-sustained oscillations.Here we neglect this high energy volcanic activity together with other transient events and focus on the background seismicity looking at its basic constituents and enlightening the associated dynamical system.
The paper is organized as follows.First, we will describe our data set and its spectral properties.In Sect. 3 we will reconstruct the trajectories in the phase space and constrain the asymptotic dynamics by false nearest neighbours method and by the Grassberger-Procaccia algorithm.In Sect. 4 we will adopt an entropy-based decomposition method (independent component analysis) to isolate the fundamental oscillations of the system, and in Sect. 5 we will perform the polarization analysis.Discussion and conclusions will follow.

Data set and spectral analysis
The data set is composed of continuous seismic signals recorded during a survey performed on the Colima Volcano from November 2005 to May 2006 by the Osservatorio Vesuviano (Istituto Nazionale di Geofisica e Vulcanologiasezione di Napoli, Italy), Observatorio Vulcanologico de Colima (Colima University, Mexico), and Instituto Andaluz de Geofisica (University of Granada, Spain).The signals were locally logged by four Lennartz Marslite stations, labelled as COCA, COME, COBA and COTE, and digitized at 62.5 samples/s.Each station was equipped with a threecomponent broadband seismometer (Lennartz LE-3D/20s).The actual operating time-intervals of each seismic station and their location with respect to the crater are shown in Table 1.
In Fig. 1 the topographic contours of the volcano and the position of the seismic stations are displayed.The raw signals have been filtered with a Butterworth four-pole highpass frequency filter with corner frequency equal to 0.05 Hz, in order to select the linear response band of the instruments.Examples of acquired waveforms are shown in Fig. 2.
As a preliminary characterization, we estimate the basic spectral content of the recordings by performing the Fourier analysis in a time-window of 60 s sliding along all the data set without overlap and convoluted with a hanning-type function.We exclude from the analysis the time-windows in which volcanic transient events, local earthquakes and teleseisms occur.Figure 3a contains the time-evolution of the power spectra averaged over non-overlapping time-intervals of 20 min.The basic spectral content is concentrated in the range 0.1-0.4Hz.In particular, two main spectral peaks can be revealed respectively around 0.15 Hz and 0.3 Hz, as shown in Fig. 3b by the mean normalized spectrum.The frequencies of these peaks may be slightly variable, and we can define two frequency intervals in which they range, that are respectively, 0.1-0.2Hz (frequency band 1, FB1) and 0.2-0.4Hz (frequency band 2, FB2).
The amplitudes of these peaks vary with time, as shown in Fig. 4a which displays the time-evolution of the standard deviation (STD) of the signal in FB1 and FB2.Both frequency bands show STD fluctuations spanning about one order of magnitude.STD of FB1 is prevalently higher than that of FB2 along all the data set except during a time-interval in the early stage of the recordings in which this behaviour is reversed as an effect of the reduction of STD of FB1.In order to well identify the edges of this time-interval, we compute the ratio between STD of FB1 and FB2 (R, hereafter) averaged over sliding windows of 4 days.R(t) is plotted in Fig. 4b.We label as Phase A the time-interval in which R is mainly lower than 1, whereas we define as Phase B the remaining time-interval, in which R is prevalently higher than 1.Phase A lasts approximately from 25/11/2005 to 15/12/2005, while Phase B begins on 15/12/2005 and lasts until the end of the data set.In Fig. 4 Phase A and Phase B are separated by a   vertical dotted line.In the behaviour of R(t), we note also an increasing linear trend superimposed to the fluctuations, which suggests that the STD of FB1 and FB2 evolve over very long time scales.

Phase space dynamics
In this section we analyse the dynamical system underlying the observed seismicity.For this aim, we recover the points in the phase space starting from a scalar series, i.e. our seismic recordings.We use the time delay method (Tak-ens, 1981) to reconstruct the asymptotic dynamics in phase space, i.e. the attractor.Being {s i } i=1,...,n our signal composed of a series of n scalar observations sampled at regular time-intervals, the reconstructed attractor consists of the vectors x i = (s i , s i+t ,..., s i+(m−1) * t ).In this formulation, the unknowns are m, the embedding dimension, and t, the time delay.We set t as the first minimum of the mutual information, which can be considered the nonlinear extension of the autocorrelation function (Fraser and Swinney, 1986).The parameter m is fixed by the false nearest neighbours (FNN) technique (Kennel et al., 1992).In this technique, the points of the trajectory in the phase space are projected onto spaces with  increasing dimensions, where they appear nearest as neighbours of some other points until the dimension of the spaces matches that of the embedding space.The procedure is complete when the fraction of these false nearest neighbours with respect to the total number of points is zero.The obtained value of m represents a lower limit for the embedding dimension and an upper bound for the attractor dimension d a .
To estimate d a we use the technique of Grassberger and Procaccia (1983), which is based on the correlation integral: where x i are the vectors above defined, and θ is the Heaviside function.The slope of C(l) in the scaling region (where it has a power law behaviour) gives the correlation dimension, which is an estimate of the dimension of the attractor (Bergé et al., 1986).We estimate the percentage of FNN and the correlation dimension of several samples of signal along all the data set.All the stations and the three ground motion components have been analysed.Moreover, three different lengths of the samples are adopted to verify the stability of the solutions.The results of FNN analysis are plotted in Fig. 5.Each column of symbols represent the fraction of FNN for the corresponding station and ground motion component relative to the embedding dimension equal to 1 (black), 2 (blue), and 3 (green).Different symbols are used for different lengths of the samples.The results are distinguished for Phase A and Phase B. We note that the differences between the different lengths of the samples are negligible, confirming that a stable process acts on the analysed time scales.The fraction of FNN is mostly zero for m = 4 in Phase A and for m = 3 in Phase B. We note also an increasing trend of the fraction of FNN estimated for the embedding dimension equal to 2 and 3 moving from COCA to COTE in Phase A, which suggests an increase of the complexity of the system going away from the crater area.
Starting from the FNN results, we apply the Grassberger and Procaccia procedure to refine the estimate of the dimension of the attractor (Albano et al., 1988).Figure 6 shows    the correlation dimension as a function of the embedding dimension, separating the results of Phase A and Phase B. In all the cases we find a plateau at low and non-integer dimensions (< 3) indicating a non-stochastic system and suggest-ing the presence of relevant elements of nonlinearity.Moreover, the two phases are associated with different dimensions of the attractor, which range 2.7-3 in Phase A and 2-2.2 in Phase B, confirming the results of FNN analysis.Thus, there

Independent component analysis
Previous analyses suggest the existence of two states of the system having different spectral properties and dynamical systems.In this section we apply a decomposition technique to find the basic constituents of the system, i.e. the independent component analysis (ICA).ICA is an entropy-based technique working in time-domain, which can find underlying factors or components from multivariate (multidimensional) statistical data on the basis of their statistical independence (Hyvärinen et al., 2001).Its goal is to find a linear representation of non-Gaussian data so that the components are statistically independent, or as independent as possible.ICA is based on the central limit theorem which states that the sum of two independent random variables has a distribution that is closer to Gaussian than any of the two original random variables.Hence, ICA finds one independent component maximizing its non-Gaussianity.
ICA is based on the following mathematical setting.We can suppose to have m different time-series x (mixtures) that we hypothesize to be the linear superposition of n mutually independent unknown signals s (independent components-ICs).How s are mixed to compose x is described by a constant unknown matrix A (mixing matrix).This mixing can be due to path, noise, instrumental transfer-functions, etc.The hypothesis is to have linear mixtures of some independent dynamical systems.We remark that in the ICA model, whereas the mixing is supposed to be a linear combination in time domain of the sources, nothing is assumed with respect to the independent dynamical systems, which can be both linear or not.Formally, the mixing model is written as x = As.In addition to the signal independence requirement, ICA assumes m ≥ n.Under these hypotheses, the goal of ICA is to obtain a matrix of separation W ≈ A −1 , in such a way that the vector y = Wx is an estimate y ∼ s of the original independent signals.
ICA cannot determine the variance (proportional to the energy) of the independent components because both s and A are unknown, implying that any scalar multiplier in one of the independent components could always be cancelled by dividing the corresponding column of A by the same scalar.For the same reason, ICA cannot establish an order among the independent components.Some heuristic approaches have been proposed in literature to achieve the separation of the independent components.The classical measure of non-Gaussianity is kurtosis, i.e. the fourth-order cumulant.For a random variable y it is defined as E(y 4 ) −3E(y 2 ) 2 , where E is the expectation value.However, kurtosis is not a robust estimator of non-Gaussianity as it can be very sensitive to outliers, whereas a good measure of non-Gaussianity is given by the contrast function negentropy J (Hyvärinen and Oja, 2000).Negentropy is defined as where y gauss is a gaussian variable with the same covariance matrix as y and H is the differential entropy H (y) = − f (y) logf (y) dy, being f (y) the density function of the random variable.The estimate of negentropy is difficult, and in practice, some approximations must be introduced.We use the fixed-point FastICA algorithm: a highly efficient computational method for performing ICA (Hyvärinen et al., 2001).Independence is here measured by the following approximation of the negentropy: where G is a suitable contrast non-quadratic function, w is an m-dimensional (weight) vector constituting one column of W , x represents our mixtures, E{w T x} = 1, ν is a standardized Gaussian random variable.Maximizing J G allows to find the weight vectors w 1 , . .., w n and hence the independent components.
ICA has been already successfully applied to a variety of experimental signals from biology to geophysics (for a review see, e.g., Fiori, 2003).Here, we apply ICA separately to samples extracted from Phase A and Phase B. Our mixture matrices are composed of 120 s of the simultaneous recordings of the background seismic signal of the three ground motion components of all the stations.In Fig. 7 we show the results of ICA for the two phases displaying the extracted time series and the relative amplitude spectra.In Phase A, ICA extracts three independent components (ICs) with a quasi-monochromatic frequency content (subplots sa.1-3 and a.1-3).In detail, an IC has a peak in the spectrum close to 0.15 Hz, whereas the other two ICs have a frequency peak around 0.3 Hz.In Phase B ICA extracts two ICs (subplots sb.1-2 and b.1-2).One of them has a quasi-monochromatic spectrum with a peak close to 0.15 Hz, sharing the same properties of one of the ICs extracted in Phase A. The other IC, instead, has a broadband spectrum with two main peaks at about 0.15 Hz and 0.3 Hz.We do not find higher-frequency ICs, confirming that the background signal is dominated by frequencies below 1 Hz.

Polarization analysis
In this section we analyse the polarization properties of the seismic noise by applying the Kanasewich algorithm (Kanasewich, 1981).This algorithm involves the diagonalization of the covariance matrix constructed with the three ground motion components.The best estimate for the polarization vector is given by the eigenvector corresponding to the highest eigenvalue.The algorithm returns three parameters: Rectilinearity (RL), azimuth (θ ), and dip (φ).RL is a measure of linearity of polarization.RL = 0 means that all the eigenvalues of the covariance matrix are equal and a preferential direction of polarization cannot be distinguished; RL = 1 means linearly polarized waves.As polarization solutions are reliable only for elliptical or linear waves, a threshold for RL can be fixed.θ is the clockwise angle between the north direction and the polarization vector.It indicates the direction of oscillation, but not the sense, therefore θ and θ +π correspond to the same polarization.ϕ is the angle between the z-axis and the polarization vector.the wave type and/or source positions will be recovered by the particle motion analysis.
We apply the Kanasewich algorithm to samples of one hour of background seismic signal extracted from the data set.Signals have been low-pass frequency filtered under 1 Hz to highlight the main spectral properties.We perform the polarization analysis within time windows of 10 s sliding along all the duration of the samples and overlapping of 50 % at each step.The lengths of the time windows are chosen to contain about two cycles of oscillations to guarantee the stability of the solutions.
In Fig. 8 we display the distribution of the polarization parameters separately for the two phases.RL solutions are on average lower during Phase A, suggesting less elongated and/or scattered oscillations, with a decrease of the mean RL approaching the crater.Moreover, in Phase A, the dip values are higher, indicating shallower oscillations.The azimuth directions do not provide a simple pattern.Indeed, the preferential oscillation direction is in the range 30 • -60 • at COBA and COTE (Figs. 8c and d, respectively); at COME it is about 150 • while at COCA it is about 100 • or 160 • depending on the phase (Figs.8b and a, respectively).We find a sharp change at COCA with an increase of the crater-pointing solutions from about 20 % to 50 %, moving from Phase B to Phase A. On the other hand, the differences in the azimuth distribution between the two phases are negligible for the other stations.
Particle motions of the two phases are plotted in Fig. 9.In Phase B the motion is prevalently elliptical with a clockwise direction in both EW-Z and NS-Z planes.In contrast, the particle motions of Phase A are more scattered, with a clockwise elliptical motion occasionally appearing.These plots also reveal that the decrease of RL inferred by the polarization analysis during Phase A is mainly the effect of a scattered wavefield and not of a less elongated particle motion.Note also the very shallow oscillations during Phase A, as inferred also by the polarization analysis.

Discussion and conclusions
We have studied the dynamical system and the wavefield properties of the background seismic signals recorded at the Colima Volcano in the period December 2005-May 2006 by four broadband three-component seismometers labelled as COCA, COME, COBA and COTE (sorted by crater-station distance).We have excluded from the analysis transient signals, such as those associated with volcanic explosions, local earthquakes, teleseisms, rockfalls, etc. Below we resume the principal results: -Fourier analysis indicates that the basic spectral content is contained in the range 0.1-0.4Hz and the main peaks are in the range 0.1-0.2Hz and 0.2-0.4Hz. -The amplitude of the main spectral peaks is timevariable.We define as Phase A the time-interval in which the frequency peak in the range 0.2-0.4Hz is dominant, and Phase B as the time-interval in which the peak in 0.1-0.2Hz is dominant.
-False nearest neighbours method and the Grassberger-Procaccia algorithm indicate different correlation dimensions of the attractor in the phase space for the two phases.This dimension is 2.7-3 in Phase A and 2-2.2 in Phase B.
-Independent component analysis -a time-domain decomposition technique -separates the seismic wavefield into three independent signals (independent components-ICs) in Phase A and two ICs in Phase B.
In Phase A, all the ICs have quasi-monochromatic spectral content with frequency peak at 0.15 Hz and 0.3 Hz.
In Phase B, one IC has quasi-monochromatic spectrum with peak at 0.15 Hz, while the other IC has a broad spectrum with two main peaks at 0.15 Hz and 0.3 Hz.
-Polarization analysis shows an increase of the craterpointing azimuths at the station COCA during Phase A. Moreover, at all the stations the dip is higher and the RL is lower (on average) in Phase A.
-Particle motion is elliptical with clockwise direction during Phase B, while during Phase A it is scattered with very shallow oscillations.
All the results suggest the existence of two states of the background seismicity, which appear during Phases A and B. These two dynamical states share a common basic constituent, that is the quasi-monochromatic IC with frequency peak at 0.15 Hz.It is straightforward to identify this signal with the microseismic noise, which shows a spectral peak at about 0.15 Hz (secondary microseisms; Bäth, 1974).The polarization direction at COBA and COTE, which are the stations farthest from the crater, identifies the SW-NE oscillation direction, and similar results appear for COCA in Phase B. In addition, the particle motion displays a predominance of oscillations with mainly anticlockwise direction, which means a retrograde motion for a source located along the SW direction, as occurs for the Rayleigh waves dominating the microseismic noise (e.g.Friederich et al., 1998).In this scheme, the source would be located along the SW direction.An imprinting of the sea waves-ground coupling could be visible in the oscillation directions nearly perpendicular to the coast line, suggesting that the source can be identified with the oceanic load on the Pacific Ocean coastline, which is about 50 km from our seismic network (Braun et al., 1996;Almendros et al., 2000;Stehly et al., 2006), whereas COME would be affected by local effects modifying the polarization direction (Hellweg, 2003).
On the other hand, the other ICs depend on the phase.Indeed, in Phase A two further ICs are extracted having quasimonochromatic spectral content close to 0.3 Hz, whereas in Phase B only one broadband further IC is extracted.This increase of the number of ICs in Phase A is in line with the results of the phase space reconstruction, which provides higher dimensions in this phase.
Many models of the microseismic noise generation involve the interaction between a periodic external forcing and medium (e.g.Correig and Urquizu, 2002;Ryabov et al., 2003;Correig et al., 2007).In this framework, our results suggest to associate the quasi-monochromatic IC with frequency peak at 0.15 Hz to the external forcing, which acts with different intensities in the two phases, as inferred by the time-evolution of R(t).In our case, the different strength of the external forcing would induce different responses of the medium.Specifically, for high-energy external force (Phase B) the medium and the source effects would be entangled and then a broadband IC extracted.Nevertheless, the correlation dimension during Phase B is very low (about 2) suggesting that the system is associated with a low-dimension dynamical system, as occurs for the systems organized in global coherent oscillations.Urquizu and Correig (1998) found a very similar fractal dimension (2.3) for the microseismic noise in the Pyrenees, supporting the hypothesis that this source is dominating in Phase B. Accordingly, in this case we find that the external forcing-induced polarization dominates.
For low-energy external forcing (Phase A), the medium/local effects would be visible as quasimonochromatic ICs with frequencies close to 0.3 Hz, with the effect of a more complex (higher dimension) system.Moreover, in this case the azimuth from the polarization analysis suggests the existence of an internal source, possibly located in the crater area, visible in the near field.The decreasing of RL in Phase A supports the existence of more than one source simultaneously acting, inducing more scattered oscillations especially at the stations closest to the crater, accordingly to the hypothesis of an internal source.Moreover, looking at the dip distributions, this additional source should induce very shallow oscillations.Also, the results of FNN analyses are in line with the existence of an internal source in Phase A as the fraction of the FNN increases moving away from the crater area, compatible with a lowering of the signal-to-noise ratio enhancing the complexity of the system.
We propose a volcanic source producing low-energy volcanic tremor in the frequency band 0.2-0.4Hz located in the crater area and affecting the near-field seismic wavefield.In this scheme, the two ICs with a frequency close to 0.3 Hz would represent the medium response to the external forcing and the volcanic tremor, which would become detectable for a decrease of the external source.The existence of volcanic tremor in the analysed frequency band is supported by the similar frequencies of the explosion-related transient seismic events, which display radial and shallow polarization in the     near field as well (Petrosino et al., 2011;De Lauro et al., 2012a).In addition, volcanic tremor with frequency content matching that of FB2 has been already observed in other areas (De Lauro et al., 2005;Haney, 2010) and a competition between volcanic and non-volcanic sources has been already detected at Stromboli, with the evidence of volcanic tremor during a phase of low-energy microseismic noise (De Lauro et al., 2005, 2006).
The volcanic source can be imputed to the interaction between the plumbing system and the low-intensity continuous degassing, which has been experimentally detected during the analysed period (Stevenson and Varley, 2008).Specifically, degassing permanently can put in vibration the main conduit of the plumbing system, exciting its fundamental vibration mode.During vulcanian-style explosions, the degassing flux is strongly enhanced due to the coherent rapid movement of aggregated magma-gas mixtures nucleated at a depth of 1-3 km, producing very energetic conduit oscillations with the appearance of higher vibration modes along the entire event and a broad spectrum just during the nucleation (Palo et al., 2009;Petrosino et al., 2011;De Lauro et al., 2012a).Our scheme follows the behaviour of many volcanoes worldwide for which explosion-related seismic events and background signals are both manifestations of a common source, i.e. the degassing process, but with different intensities.
Hence, our scheme involves two sources: an external source inducing the microseismic noise and a volcanic (internal) source.In this framework, it is straightforward to wonder if the two states of the background seismicity correlate with the volcanic activity and the atmospheric conditions, which are considered as the basis of the microseismic excitation.With this aim, we look at the simultaneous time-evolution of the density and the amplitude of the volcanic explosions, of R, and of the atmospheric pressure (Fig. 10).Both volcanic parameters do not show statistically relevant patterns of activity in the analysed time-interval, whereas the standard deviation of the atmospheric pressure during Phase A (0.7 hPa) is lower than in Phase B (1.7 hPa).Moreover, the highest values of R correspond to a lowering of the atmospheric pressure (April, 2006).Hence, an imprint of the atmospheric conditions appears in R, suggesting that the external forcing is induced by varying meteorological conditions.The volcanic source would appear for a lowering of the microseismic noise influence, which, in turn, is induced by changes of the meteorological conditions with effects on the oceanic load intensity.
Therefore, the intensity of the microseismic noise discriminates between state A and state B of the system, which can be interpreted as the superposition of two sub-systems: the volcanic source and the interaction between the external forcing and the medium.In this framework, the difference between the correlation dimensions in the two phases would be due to the contribution in Phase A of the volcanic source, whose dynamics thus should evolve on an attractor of dimension ∼ 1, whereas the other sub-system would display an attractor of dimension ∼ 2.
The volcanic seismicity is induced by instability in the magmatic flowing fluids inducing the coupling with the rock and a nonlinear behaviour (see, e.g., Chouet, 1996;Konstantinou and Schlindwein, 2003, and references therein).Such systems produce self-excited vibrations, whose phase space dynamics evolve on limit cycles that are closed trajectories asymptotically approached, developed in nonlinear systems and displaying a dimension ∼ 1, in line with our ob-servations (Julian, 1994;De Martino et al., 2002;Balmforth et al., 2005;De Lauro et al., 2008;Fujita, 2008).
On the other hand, the characteristics of the other subsystem depend on the phase.Indeed, in our framework, in Phase A the medium response is visible separately from the source, whereas in Phase B it cannot be distinguished.This change can be due to different attractors in the two phases induced by different behaviour of this sub-system depending on the amplitude of the forcing, as occurs for the systems able to develop chaotic attractors.In our case, the amplitude of the forcing is the system parameter tuning the transition to chaos.It is well known that, as the parameter is varied, before fully developing chaos, the systems commonly display two-frequency quasi-periodicity and then a periodic attractor.Both of these behaviours imply an evolution on a torus in the phase space, with the latter displaying closed trajectories, and the former evolving over all the surface (Bergé et al., 1986;Ott, 1993;Godano and Capuano, 1999).
We propose that in Phase A, i.e. for a low value of the system parameter, the sub-system displays a two-frequency quasi-periodicity, and the dynamics can be decomposed along the two periodicities of the torus, thus separating the source and the medium imprints and fitting the dimensionality of this sub-system.In Phase B, the system is moving towards chaos and the attractor changes to a closed trajectory, which can no longer be decomposed and the medium and source become "syncronized" and entangled, with the effect of extracting a two-frequency independent component (Bergé et al., 1986;Ott, 1993).In this phase, an imprint of the periodic forcing is however visible in the other nearly monochromatic independent component, making the total dimension ∼ 2 also in this case, whereas the volcanic effects are masked by the increase of the energy of the external source.We remark that a system with elements of chaoticity and with nonlinear interactions between the forcing effects and the medium can explain the evidence that microseismic noise can be not coherent in space (e.g., Correig and Urquizu, 2002).Moreover, this scheme implies that this sub-system can be described (in a first approximation) by a nonlinear second-order differential equation with a harmonic forcing in Phase A (e.g., Godano and Capuano, 1999;Correig and Urquizu, 2002), whereas in Phase B such synchronization phenomena must be associated with the interaction between two nonlinear oscillators (Ott, 1993).
Our results remark also the importance of the seismicity in the near field, i.e. within 2-3 km from the crater.The deployment of a dense seismic network in this region can enlighten the features of the degassing process and of the plumbing system of the Colima Volcano, limiting the effects of nonvolcanic sources.Indeed, strong constraints on the comprehension of the volcanic source mechanisms have been produced by the large distance between the crater and the stations installed in the seismic surveys so far performed at the Colima Volcano.
At the same time, for future investigations at the Colima Volcano, it will be needed to extend the period of monitoring; indeed, the linear trend emerging in the time-evolution of R indicates the existence of long-term variations, suggesting that the dynamics of the volcano and of the external sources act over very long time scales (years).

Figure 1 .
Figure 1.Topographic contours of Volcàn de Colima with the location of the four seismic stations (blue diamond) and of the crater area (black triangle).

Fig. 1 .
Fig. 1.Topographic contours of Volcàn de Colima with the location of the four seismic stations (blue diamonds) and of the crater area (black triangle).

Figure 3 .
Figure 3. Power spectra (computed by Fast Fourier Transform and averaged over 20-minutes time-intervals) of the North-South component of COCA as function of time.Horizontal colourbar indicates the correspondence between the colours and the power spectra.(b) Mean normalized power spectra obtained by averaging the spectra computed in the 60s time-windows over all the data-set and dividing by the maximum value.The principal results do not depend on the station nor on the ground motion component.

Fig. 3 .
Fig. 3. (a) Power spectra (computed by the fast Fourier transform and averaged over 20-min time-intervals) of the North-South component of COCA as function of time.Horizontal colour bar indicates the correspondence between the colours and the power spectra.(b) Mean normalized power spectra obtained by averaging the spectra computed in the 60 s time-windows over all the data set and dividing by the maximum value.The principal results do not depend on the station nor on the ground motion component.

Figure 4 .
Figure 4. a) Time-evolution of the standard deviations of the two FBs estimated within time-windows of 60 seconds and averaged over the three ground motion components.The curves are plotted in log-scale and transient seismic events are discarded from the analysis.b) Ratio between the standard deviation of FB1 and FB2 averaged over 4 days to get a smoothed curve.We define as Phase A the time interval 25/11/2005-15/12/2005 in which this ratio is mostly lower than one, whereas in the rest of the data-set (Phase B) there is a predominance of higher ratios.

Fig. 4 .
Fig. 4. (a) Time-evolution of the standard deviations of the two FBs estimated within time-windows of 60 s and averaged over the three ground motion components.The curves are plotted in log-scale and transient seismic events are discarded from the analysis.(b) Ratio between the standard deviation of FB1 and FB2 averaged over 4 days to get a smoothed curve.We define as Phase A the time-interval 25/11/2005-15/12/2005 in which this ratio is mostly lower than one, whereas in the rest of the data set (Phase B) there is a predominance of higher ratios.

Figure5.
Figure5.Fractions of FNN vs. ground motion components.The symbols in each panel represent the fraction of FNN for the corresponding station and ground motion component relative to the embedding dimension equal to 1 (black), 2 (blue), 3 (green).Different symbols are used for different sample lengths: circles for 256s, squares for 320s and crosses for 400s.The upper panels (a.1, a.2 and a.3) are associated to the Phase A and the bottoms (b.1, b.2 and b.3) to the Phase B. The legend at the bottom shows the correspondence between the numbers indicated in the horizontal axes and the ground motion components.

Fig. 5 .
Fig. 5. Fractions of FNN vs. ground motion components.The symbols in each panel represent the fraction of FNN for the corresponding station and ground motion component relative to the embedding dimension equal to 1 (black), 2 (blue), 3 (green).Different symbols are used for different sample lengths: circles for 256 s, squares for 320 s and crosses for 400 s.The upper panels (a.1, a.2 and a.3) are associated to Phase A and the bottom (b.1, b.2 and b.3) to Phase B. The legend at the bottom shows the correspondence between the numbers indicated in the horizontal axes and the ground motion components.

Figure 6 .
Figure 6.Correlation dimension versus embedding dimension for two representative samples from the Phase A (09/12/2005 13:00) and the Phase B (14/02/2006 13:00).All the available stations and all the ground motion components have been analyzed.The correlation dimensions have been estimated as the slopes of the correlation integral as function of the scaling parameter l (see eq. 1).The correlation dimension converges to different limit values for the two Phases, that are in the range (2.7-3) for the Phase A and (2-2.2) for the Phase B.

Fig. 6 .
Fig. 6.Correlation dimension versus embedding dimension for two representative samples from Phase A (09/12/2005 13:00) and Phase B (14/02/2006 13:00).All the available stations and all the ground motion components have been analysed.The correlation dimensions have been estimated as the slopes of the correlation integral as function of the scaling parameter l (see Eq. 1).The correlation dimension converges to different limit values the two phases, which are in the range of 2.7-3 for Phase A and 2-2.2 for Phase B.
Figure 7 ICA results in time and frequency domain.The time-series extracted in the Phase A (sa.1, sa.2 and sa.3) show one quasi-monochromatic peak close to 0.15 Hz (a.1), and two close to 0.3 Hz (a.2 and a.3).On the other hand, the waveforms of the two ICs extracted in the Phase B (sb.1 and sb.2) display one quasi-monochromatic peak close to 0.15 Hz (b.1), and a broad spectrum with two main peaks at about 0.15 Hz and 0.3 Hz (b.2).

Fig. 7 .
Fig. 7. ICA results in time and frequency domain.The time series extracted in Phase A (sa.1, sa.2, and sa.3) show one quasi-monochromatic peak close to 0.15 Hz (a.1), and two close to 0.3 Hz (a.2 and a.3).On the other hand, the waveforms of the two ICs extracted in Phase B (sb.1 and sb.2) display one quasi-monochromatic peak close to 0.15 Hz (b.1), and a broad spectrum with two main peaks at about 0.15 Hz and 0.3 Hz (b.2).

Figure 8 .
Figure 8. Examples of polarization analysis for (a) COCA, (b) COME, (c) COBA, (d) COTE in Phase A (redFig. 8. Examples of polarization analysis for (a) COCA, (b) COME, (c) COBA, (d) COTE, in Phase A (red histograms) and in Phase B (blue histograms).For each station, the panel on the left contains the distribution of RL fraction; the central panel displays the azimuth distribution; the right graphic shows the dip distribution.
histograms) and in Phase B (blue histograms).For each station, the panel on the left contains the distribution of RL fraction; the central panel displays the azimuth distribution; the right graphic shows the dip distribution.

Figure 9 .
Figure 9. Example of particle motions at COCA, COME, COBA and COTE for a sample of 16 seconds of the Phase A (09/12/2005 13:00, 30/11/2005 03:00 for COTE) and the Phase B (14/02/2006 13:00, 22/01/2006 03:00 for COTE).The motion is projected into EW-Z an NS-Z planes.Red circles indicate the initial point.Units are counts.Clockwise motion dominates the Phase B. In the Phase A, the motion is scattered and a stable oscillation pattern does not appear.Clockwise motion is occasionally detectable, especially at the stations farthest from the crater; it is visible in the plots as trajectories escaping from the core of the motion.

Fig. 9 .
Fig. 9. Example of particle motions at COCA, COME, COBA and COTE for a sample of 16 seconds of Phase A (09/12/2005 13:00, 30/11/2005 03:00 for COTE) and Phase B (14/02/2006 13:00, 22/01/2006 03:00 for COTE).The motion is projected into EW-Z and NS-Z planes.Red circles indicate the initial point.Units are counts.Clockwise motion dominates Phase B.In Phase A, the motion is scattered and a stable oscillation pattern does not appear.Clockwise motion is occasionally detectable, especially at the stations farthest from the crater; it is visible in the plots as trajectories escaping from the core of the motion.

Figure 10 .
Figure 10.Time-evolution of the density (upper panel), and amplitude (second panel from the top) of volcanic explosion, R (third panel from the top) and atmospheric pressure (bottom panel).Dotted vertical line separates the Phase A and the Phase B. We do not observe statistically relevant modifications of the volcanic parameters along the data-set, whereas correlations between R and atmospheric pressure appear.See the text for details.

Fig. 10 .
Fig. 10.Time-evolution of the density (upper panel), and amplitude (second panel from the top) of volcanic explosions, R (third panel from the top) and atmospheric pressure (bottom panel).Dotted vertical line separates Phase A and Phase B. We do not observe statistically relevant modifications of the volcanic parameters along the data set, whereas correlations between R and atmospheric pressure appear.See the text for details.

Table 1 .
(a-c) Position of the stations with respect to the crater in polar coordinates.The azimuth is measured from north direction clockwise.(d) Operating time-interval of each station.