Fractal dimension analysis of the magnetic time series associated with the volcanic activity of Popocat́epetl

Fractal analysis of the total magnetic field (TMF) time series from 1997 to 2003 at Popocat épetl Volcano is performed and compared with the TMF-series of the Teoloyucan Magnetic Observatory, 100 km away. Using Higuchi’s fractal dimension method ( D). TheD changes over time for both series were computed. It was observed, when the time windows used to compute D increase in length, both series show nearly the same behavior. Some criteria of comparison were employed to discriminate the local effects inherent to volcano-magnetism. The simultaneous maximum in D (1.8) of the TMF series at Popocat épetl Volcano and the recovered volcanic activity indicates a scaling relation of the TMF at Popocat́ epetl Volcano and demonstrates a link between the magnetic field and volcanic activity.


Introduction
In this work, we study fluctuations of the total magnetic field (TMF) near Popocatépetl volcano (Tlamacas site) that could be related to triggering its activity in terms of Higuchi's fractal dimension.To assess the use of magnetic time series as an indicator of precursory volcanic activity, the fundamental issue to address is if these parameters are able to reveal dynamical characteristics of volcanic activity.Some previous methods have been used to find changes in TMF linked to volcanic activity at Popocatépetl (Martin-del-Pozo et al., 2002, 2003;Cifuentes-Nava, 2009).These methods are based on the Rikitake (1968) approach, which consists of a simple direct comparison between data recorded at different locations.Recently, the use of monofractal and multifractal methods in investigating the temporal fluctuations of geoelectrical and geomagnetic signals has been shown to be a good candidate for being used to better understand and characterize the complexity of these fields (Telesca et al., 2004;de Santis et al., 2004;Hongre et al., 1999;Uritsky et al., 2006;Telesca and Hattori, 2007;Varotsos et al., 2002Varotsos et al., , 2003)).
The study of volcanic triggering has received special attention in recent years through both direct field observations and historical descriptions of eruptions and/or earthquake activity.The direct field observations comprise gravity, seismic, and magnetic data, geochemical measurements of the ejecta, including direct flight observations and changes in crater or dome morphology.Some reports of clustered eruptions and earthquakes may indicate that interactions exist in some regions.Therefore, the interactions between volcanic edifices and their tectonic surroundings should also be studied.According to Eggert and Walter (2009), volcanoes are systems that interact with their environments at different scales.Assuming this point of view, it is necessary to explore different approaches of data analysis.
The geophysical phenomenon underlying magnetic field variations is complex and the physical laws that govern the process are not completely known.The use of monofractal and multifractal methods in investigating temporal fluctuations of geoelectrical signals (Telesca et al., 2003) can also be employed to investigate fluctuations in magnetic field and to better understand such complexity.
The stress field of the lithosphere is an expression of a hierarchy of volumes, or blocks, that move in relation to each other and are divided into smaller blocks, such as shields or mountain chains.This stress field determines the collisions between these blocks, and the accumulated energy is radiated as seismic waves waves (Telesca et al., 2003).A system consisting of porous media could develop due to the decomposition of elements on smaller scales or due to the largerscale structure.Both of these processes could cause scaling between geoelectrical or geomagnetical signals and seismic emissions from microcurrents along cracks (Telesca et al., 2003).
The evolution of the earth's crust toward a self-organized critical state takes place not only seismically, as a result of the block-spring fractal structure in the fault zone, but also electromagnetically, as a result of the fractal conductor-dielectric structure.
Volcanoes are also manifestations of the stress field of the lithosphere, and they are considered complex systems.Magnetic time series recorded at volcanoes can be useful for extracting quantitative information about their complexity (Uyeda et al., 2002(Uyeda et al., , 2009)).Monofractal analyzes could show the presence of scaling behavior, which would imply the existence of correlated time structures.In this context, we have selected Higuchi's fractal analysis to characterize our data set and to find time structures in the magnetic time series.Higuchi (1988) showed that this method provides a precise estimation of the fractal dimension, even for a small number of data.Higuchi developed this method as an alternative to spectral analysis because there is a relation between D and β (the slope in plots of the spectral analysis).
The aim of this work is fractal analysis of the TMF data series at Popocatépetl Volcano (Tlamacas magnetic station, TLA), the second highest point in Mexico (5450 m above sea level), and the Teoloyucan Magnetic Observatory (TEO), 100 km away.The goal of this analysis is discriminating the local effects inherent to volcano magnetism and inferring if some time structures of the TMF are related to volcanic activity at fumaroles or domes, seismic activity, ash emissions, etc.Our results in terms of the Higuchi's fractal dimension show an important characterization of the behavior of these series and the scaling between them.

Higuchi's fractal dimension method
New methods based on nonlinear dynamics have become important tools for obtaining relevant information from time series.Higuchi (1988Higuchi ( , 1990) ) proposed a method to calculate the fractal dimension of self-affine curves in terms of the slope of the straight line that best fits the curve length versus the time interval (the lag) on a log-log graph.This method consists of considering a finite set of data taken at a regular interval ν(1), ν(2), ν(3), • • • , ν(N).From this series, we construct a new series ν m k , defined as with m = 1, 2, 3, ..., k, where denotes Gauss' notation and k and m are integers that indicate the initial time and the time interval, respectively.For a time interval k, there are k new sets of time series.For instance, in the case of k = 3 and N = 100, the three series that are obtained with this process are The length of the curve ν m k is defined as where the term (N − 1) (N − k) k k is a normalization factor.Then, the length of the curve for the interval k is L (k) , the average over k sets L m (k).Finally, if then there is a scaling relation with a scaling exponent of D, which is Higuchi's fractal dimension.Higuchi showed that this method provides a precise estimate of the fractal dimension, even for a small number of data.This method was developed as an alternative to spectral analysis because although there is a relation between D and β (the slope in plots of the spectral analysis).For auto-affine curves, this fractal dimension is related to the exponent β by means of β = 5-2-D.If D is in the range 1 < D < 2, then 1 < β < 3 ( Van Ness, 1968).Higuchi's method allows the clear definition of two or more regions in which the graph of log L m (k) versus log k is divided by crossovers, i.e. the points that divide different scaling regions with different values of the fractal dimension D (Higuchi, 1988(Higuchi, , 1990)).When D is 1.5, the dynamics of the systems are Brownian, while D less than 2 corresponds to pink noise and D = 2 is white noise (Gálvez-Coyt et al., 2011).

Popocatépetl Volcano
Popocatépetl (19.02  2008).It is part of the Sierra Nevadas, a north-south-trending mountain range (Fig. 1), and its geologic past clearly indicates that it is capable of producing major catastrophic eruptions: three Plinian events have occurred within the past 5000 yr, well within the period of human settlement in central Mexico (Siebe et al., 1996;Siebe and Macías, 2004).
Popocatépetl's history of fumarolic activity, ash emissions and Plinian events has been recorded by stratigraphy and by oral and written documents of both pre-Columbian and Spanish colonial times.Throughout its history, Popocatépetl has been repeatedly active.According to Robin (1984) and Martin Del Pozzo et al. (1997), the most recent large Plinian eruptions occurred approximately 3000 BC, 200 BC and 800 AC, but small eruptions have been reported in each century, the last of which occurred in the 1920s.However, Siebe and Macias (2004) reported 450 yr between each of the last three Plinian eruptions.
Fumarolic and seismic activity increased on 21 December 1994.Popocatépetl began erupting and it has since had fluctuating amounts of activity, principally ash emissions and seismicity.Subsequent to this activity, volcano monitoring was implemented by the Instituto de Geofísica, UNAM, which consisted of seismic stations, four global positioning system stations, two meteorological stations, surface temperature measurements, periodical sampling of water, ashes and gas (SO 2 ) for chemical analysis, gravity field measurements, and measurements of the total magnetic field with a Geometrics G856 protonic precession magnetometer.
In the present study, we focused on the TMF recorded from 1997 to 2003 at Tlamacas station.

Total Magnetic Field data
The TMF data were recorded from 1997 to 2003 at two sites: Tlamacas station and Teoloyucan Magnetic Observatory of Mexico, located at 19.745 • N, 99.193 • W, approximately 100 km away from the volcano.The Popocatépetl magnetic station (Tlamacas) is 4 km NNW of the volcano's crater (19.057 • N, 98.637 • W, 4029 m a.s.l.) and 300 m from the Tlamacas hiking shelter; this is the first permanent magnetic station on a Mexican volcano.At both sites, the magnetometric data were recorded with a Geometrics G856 protonic precession magnetometer with a sampling period of one minute.The resolution of the instruments is 0.1 nT, and the precision is determined according to established standards at magnetic observatories (Wienert, 1970;Jankowski and Sucksdor, 1996).For TLA, the precision is 0.4 nT, following accepted practices of inter-instrumental comparison proposed by Martin-del Pozo et al. (2003).
The TMF measurements at TLA (Fig. 2) were transmitted by radio modem with a FreeWave 900-MHz spread-spectrum to the Geophysics Institute of the National Autonomous University of Mexico (UNAM) in Mexico City during some periods of the monitoring.During other periods, the data were stored on floppy disks at a local computer due to difficultly in accessing the abandoned hiking shelter at Tlamacas to maintain the equipment, which resulted in a lack of some important information.This lack of data was treated with moving average (MA), as shown in Fig. 4.
Following the geodynamo theory, the main contribution to the TMF is produced in the earth's interior, with smaller contributions produced by currents in the ionosphere and magnetosphere and geomagnetic storms that contribute to all other electromagnetic interactions.Therefore, some actions are necessary to isolate the anomalous TMF produced by volcano-magnetic activity.First, we assume that the measurements of TMF include the aforementioned contributions and that the differences between series of measurements from TLA and TEO can be attributed to the volcanomagnetic activity of Popocatépetl with some other minor local contributions.We also assume that the anomalous magnetic field produced by magmatic activity may be due to three phenomena: (1) thermal fluctuations in the vicinity of a magma chamber affect surrounding rocks, shifting the Curie point of the host rocks and decreasing the magnetic field when temperature is increased (e.g., Dzurisin et al., 1990;Zlotnicki and Bof, 1998).This effect depends on the thermal inertia and the volume of newly intruded magma.Other transients may be caused by ( 2) piezomagnetism (e.g., Johnston and Stacey, 1969) and ( 3) electrokinetic (e.g.Zlotnicki and Le Mouël, 1988) effects on rocks.
Following Rikitake (1968) and Martin-del-Pozo et al. (2002), the weighted differences method can be applied to cancel the effects of non-volcanic external sources.This method assumes that every measurement F can be expressed as where F c is the earth's magnetic field, F e represents the external sources in the magnetosphere and ionosphere, and F n represents the local subsurface anomalies, such as those originating from magmatic processes.By subtracting the corresponding component of the TMF, the effect of F c will be canceled over short time spans.Some previous works   (Martin-Del Pozzo et al., 2002Pozzo et al., , 2008) ) have focused on the analysis of TMF fluctuations, including normalized differences (Rikitake, 1968), spectral analysis and correlation with synthetic Those studies show some evidence of correlation between these fluctuations and volcanic events, such as ash emission and dome creation.
While the geomagnetic field changes on time scales from milliseconds to millions of years, it is important to focus on short-term instabilities of the magnetic field due to volcanic activity.We will statistically characterize the magnetic series at TEO and TLA independently and will then compute differences based on the Rikitake method and using Higuchi's fractal analysis over different time windows to find possible scaling relations.

Data analysis
Figure 2 shows the TMF recorded at one-minute intervals at the Tlamacas and Teoloyucan sites, where the data series are well correlated.Gaps in data are shown as lines connecting neighboring data.Some periodic features can be observed in Fig. 2b, which shows data recorded from 15-21 December 1999.There, the scale of observation clearly shows different behavior, scales of the order of minutes show local effects at TLA, scales of the order of some hours show similar behavior in both series.
Two approaches were followed to analyze the data, the first one was focused in examined the global series and the second one was devoted to evaluate the evolution of the series in time by using different window sizes.
We try to find the global fractal dimension for TLA and TEO, and we look for the fractal dimension of the differences between them (D-diff).Using the Higuchi method and k = 500, we obtained the plot shown in Fig. 3.The linear fits for D-TLA and D-Diff have similar slopes, implying the same fractal dimension (D-TLA = 1.66 k ± 0.0022 and D-diff = 1.66 k ± 0.0018).The slope for the D-TEO line (1.46k ± 0.0021) suggests a fractal dimension < 1.5 for the TEO series.We assume the differences in D-TLA and D-TEO are associated with magnetic local effects at TLA and D-diff preserves these local effects.Autocorrelations of TEO-TMF and TLA-TMF series for different window sizes were performed.We observed that a long-range correlation is present for all window sizes.Periodic features corresponding to diurnal and seasonal signatures are present in the autocorrelation, but as these features are not the focus of this study, we do not explore them further.
To independently characterize both series in terms of their fractal properties, first a moving average (MA) was used in windows to account for some gaps in data and weakly smooth them.Steps of 1, 10, 20, 50, and 60 min glided the sliding window.Second, we compute the fractal dimension using Higuchi's method (1988Higuchi's method ( , 1990) ) and analyze the data with sliding windows with lengths of 10 to 720 samples, should be noted that there are two types of windows: one for MA and the other one for Higuchi method, as is showed in Fig. 4. The fractalities of the D-TEO and D-TLA series are largely independent of window length, however when windows are shorter than six hours, changes in the TMF reflect the local effects produced by external sources, including ionospheric effects.
As a result of the effect of smoothing in MA, a log-log plot (Fig. 5) was made for a MA of non-overlapping windows of an hour (approximately 60 samples).For points 0 0.2 0.4 0.6 0. 8 1 1.2 1.4 1.6 1.8 2  1 corresponding to short windows of 1 to 10 data points, the fitted slopes of the linear regressions to these data provide the fractal dimension, which is D = 1.4736 for the TEO-TMF series and D = 1.4747 for the TLA-TMF series.In this window, the global fractal behavior is similar for both series.These D-values correspond to Brownian motion (D = 1.5), similar to the results in Fig. 3 (D = 1.46 and 1.47).However, the value for D-TLA is smaller than in Fig. 3, which could be due to local magnetic effects, particularly volcanism.The behavior for k > 10 shows the existence of periodic components, for which k = 10 1.38 corresponds to a time of 24 h.This feature is present in both series and corresponds to diurnal variations.The fact that no differences exist between the observed TMF-values at both sites in terms of their global fractality is unsurprising because the most important component in the recorded TMF corresponds to the earth's field.Nevertheless, the data recorded at TLA are expected to reveal some evidence of volcanic activity in the TMF.Del Pozo et al. (2002) reported that some features of the TMF weighted difference (Rikitake, 1968;Zeng et al., 1998) between TLA and TEO should be linked to volcanic activity, such as dome creation, ash emission and volcanic tremor.
The evolution of D-values for TEO and TLA for overlapping sliding windows glided by steps of 10 to 60 samples was computed for both series to characterize their fractality.Figure 6 shows the behavior of D for the recovered data with three different windows.Many trials were completed to choose which scale clearly shows changes in D-values without a loss of generality; we present three of these trials for each site.The pattern of fluctuations throughout the evolution of D-values in time is the same for any window size less than 360 samples (∼ 6 h) and for any step, showing scale invariance.However, windows longer than 360 samples produce the same pattern of fractal evolution in both series.Independent of window length, D-values are always greater than 1 for both series.The D-TLA values fluctuate approximately 1.5 with a maximum of 1.8, while the D-TEO values are always approximately 1.3, corresponding to Brownian motion.This implies an organization of the system and the existence of a scaling relation for specific periods.Figure 6c shows that for both series, the fractal dimension is approximately D = 1.5 for windows that are 360 samples in length (about six hours) or shorter.Following the philosophy of the Rikitake method, we compute the numerical differences (without the weight factor) between TLA and TEO for some windows.As shown by Eq. ( 5), the effect of F c will be canceled, but the addition of F e and F n for both series will remain in the computed differences.Even when we know that the external sources of the magnetic field are not necessarily the same for both sites, we presume that the contribution of this magnetic source is small and that these computed differences correspond to the local subsurface anomalies, which are the volcano-magnetic effects.We compute the MA for selected windows between 10 and 720 samples in length (from 10 min to 12 h) and steps from 1 to 60 samples, and then compute D-values for windows from 10 to 360 samples and several combinations of them.The fractal dimension is then computed using Higuchi's method for the resulting differences using the same windows and steps.Figure 7 shows the computed D-values, which range between 1.2 and 1.8 and are extraordinarily similar to those observed for D-TLA; the maximum D-values are occur at the same time in both series.The behavior of the D-Differences is independent of the size of the window for lengths shorter than 720 samples.We observe the squares of the differences in fractal dimension for both data series.Because D-values fluctuate between 1 and 2 for both series, these differences D should be between 0 and 1.They are computed as Even as Higuchi's method reveals the scale invariance of the data, we expect the influence of volcanic activity on the magnetic field to be reflected in these D.
Figure 8 shows the computed D for different MA windows, steps and windows for D-values of both data series.The pattern of the D evolution follows the same that of the evolution of D for TLA for windows shorter than 360 samples (about six hours).We conclude that the organization shown in these differences is due to the organization of data at TLA.This organization could be due to a volcanic mechanism, while the D-TEO pattern is quite regular.However, Fig. 8d shows that as the window length increases, D approaches zero except for during gaps in the data (green line).For windows greater than 360 samples, both series have the same D-value.The value of D increases during data gaps, especially for windows longer than 360 samples.

Discussion
The significant parameters in the determination of fractal dimension are the size of the window for the MA, the step, and the number of samples used to find the temporal evolution of D. The numerous trials in this study have helped us certify that our findings are consistent and invariant of scale below a six-hour threshold (360 samples).A clear temporal similarity is observed in volcanic activity and the maxima in computed fractal dimension of D-TLA, D-Diffs and D. Despite these similarities, gaps in data for TLA-TMF and/or TEO-TMF influence the computed values of D-Diffs and D, as is obvious in Fig. 8.The computed differences, D, show the same organization as TLA for windows shorter than six hours.This organization may be due to volcanic mechanisms, while the evolution of D-TEO is quite regular, reflecting the global magnetic effects for all analyzed windows.
The pattern of D-TLA can be analyzed to find a robust relationship between fractal organization and volcanic activity.1).
The recovered volcanic activity and the maxima in fractal dimension of TMF-TLA are concurrent in time, and this fractal dimension corresponds to D ≈ 1.8.The change of fractal dimension in time denotes the existence of correlated time structures, characterized by scaling laws distributions.Such properties are features of phenomena closely to systems in a critical point (Varotsos et al., 2009).
To find evidence of a relationship between the pattern of D, D-TLA, D-diffs and volcanic activity, we compare the corresponding data series from December 1997 to April 2003 (Fig. 9).The maximum D-values for all studied series are concurrent with recorded volcanic activity (Martin Del Pozzo et al., 1997, 2002, 2003, 2008and Martin Del Pozzo, 2012) consisting of fumaroles of varying height (1000 to 4000 m), ash emissions, incandescent scoria, ballistics, volcano-tectonic seismicity, and dome growth.These signals were cataloged with an arbitrary code as described in Table 1.
From these comparisons, we observe clear coincidences between maxima in computed fractal dimension for all series and volcanic activity.However, we cannot conclude that these coincidences show evidence of precursors to volcano-magnetism, although changes in fractal dimension occur fairly simultaneously with other signals of volcanic activity.
Table 1.Codes used to catalog the volcanic activity reported for Popocatépetl Volcano.

Conclusions
A clear relationship and a long-range correlation between the TEO-TMF and TLA-TMF series are observed when classical statistics are applied to both series.
The computed D-TLA and D-diffs show clear changes in organization, with values ranging between 1.4 and 1.8, while D-TEO values are always between 1.2 and 1.4.
The aforementioned volcanic activity signals and maxima in fractal dimension of TMF-TLA are clearly concurrent in time, and they correspond to D ≈ 1.8, indicating a complex system and the existence of correlated time structures, a scaling relation, which denotes some type of dynamic organization.
Changes of the TMF for windows shorter than six hours reflect the local effects produced by external sources, including ionospheric effects.For windows greater than six hours, the local effects are minimized and the global effects of the geomagnetic field dominate.The volcano-magnetism effects should be observed in short time intervals from seconds up to six hours.
Due to the coincidence of maxima in D-TLA, D-Diffs and D, we conclude that it is only necessary to analyze the evolution of the fractal dimension corresponding to the TMF signal monitored near the volcano of interest while a second station, 100 km away, just provides global information.This is important for volcanic monitoring, as volcanologists are always in search of new ways to examine precursors to volcanic activity.Changes in the fractal dimension of the TMF measured in volcanoes could contribute to the study of precursors.
The observed maxima in fractal dimension occur fairly simultaneously with volcanic activity, and they could be associated with volcano-magnetism.Thus, the investigation of fractal changes over short time lags should be a new focus in volcanic research.

Fig. 2 .
Fig. 2. Upper: total magnetic field recorded at one-minute intervals at Tlamacas station and Teoloyucan Magnetic Observatory from 1997 to 2003.Data gaps are filled with local averages.Lower: data from 15-21 December 1999 showing periodic features.

Fig. 3 .Fig. 4 .
Fig.3.The length L m (k) for the magnetic field data on a doubly logarithmic scale for the fractal analysis of the total TMF series at TLA and TEO and the differences between them.

Fig. 5 .
Fig. 5.The length L m (k) for the magnetic field data on a doubly logarithmic scale, (a) analysis for the TEO-TMF series; (b) analysis for the TLA-TMF series.

Fig. 6 .
Fig. 6.Three examples of windows used to compute the moving average (MA) and D over time, with different steps used to scan the TMF series at TLA and TEO, as indicated in each subplot (a, b, and c).

Fig. 7 .
Fig. 7. Three examples of the D-diffs for TLA and TEO, (a) the computed TMF differences for TLA-TMF and TEO-TMF with windows of 20 samples and a step of 1; (b, c, and d) D-diffs computed with different windows of MA and scan steps, as indicated in each subplot.

Fig. 8 .
Fig. 8.The square of the differences in fractal dimension, D, obtained for different windows, (a) an example of the computed Dvalues for TLA-TMF and TEO-TMF for MA-windows of 30 samples and a step of 1; (b, c, and d) D computed for different windows of MA and steps, as indicated in each subplot.

Fig. 9 .
Fig. 9. Comparison of the computed values for (a) the computed D-values for TLA-TMF (red) and TEO-TMF (blue); (b) D-TLA; (c) D-diffs: differences between TLA-TMF and TEO-TMF; (d) D, the square of the D difference; and (e) the recorded volcanic activity from December 1997 to April 2003 (Table1).