Chaos, Solitons and Fractals 21 (2004) 491–500 www.elsevier.com/locate/chaos
Long-range correlation analysis of earthquake-related geochemical variations recorded in Central Italy Vincenzo Lapenna a, Giovanni Martinelli b, Luciano Telesca
a,*
a
b
Istituto di Metodologie per l’Analisi Ambientale, Consiglio Nazionale delle Ricerche, C.da S.Loja 5, 85050 Tito Scalo (PZ), Italy ARPA, Regione Emilia Romagna, Sezione Reggio Emilia, Centro Eccellenza Acque, Via Amendola 2, 42100 Reggio Emilia, Italy Accepted 9 December 2003 Communicated by I. Antoniou
Abstract The long-range correlation properties in the hourly time variability of geochemical signals measured in a 70 m depth well at Triponzo (Umbria region), are investigated by the detrended fluctuation analysis (DFA). DFA is a data processing method that allows for the detection of scaling behaviors in observational time series even in the presence of non-stationarities. The procedure adopted has allowed for the unambiguous identification of possible correlations among the recorded signals and local earthquakes. Ó 2004 Elsevier Ltd. All rights reserved.
1. Introduction Geochemical signals can often bear witness to deep geological processes such as seismic, volcanic or geothermal activities. In particular, in seismically active areas, the study of chemical variations can represent an important tool in understanding seismological phenomena [21,22]. Co-seismic geochemical variations have been detected in some gaseous vents and natural springs during the last seismic crisis occurring in the Umbria region (Central Apennines, Italy), which started on September 26th, 1997 with several moderate earthquakes (up to ML 5:8) [12,14]. Temporal variations were also recorded in the physico-chemical parameters of three thermal springs within two and half hydrologic years including the 1997–1998 Colfiorito (Umbria-Marche) seismic sequence [7]. Manual weekly sampling utilized by Favara et al. [7] and by Italiano et al. [14] was not capable of unambiguously identifying the shortterm geochemical variations occurring in the thermal waters of the monitored area. Thus, an automatic probe was let down into a 70 m deep well located in a thermal spring area, previously selected by Martinelli and Albarello [17].
2. Data The sampling site, called Triponzo, is located in the Central Apennines (Fig. 1) that are affected by frequent seismicity associated with an extensional strain field. In this area, earthquakes generally occur in the sedimentary cover at a
*
Corresponding author. Tel.: +39-0971-427206; fax: +39-0971-427271. E-mail address:
[email protected] (L. Telesca).
0960-0779/$ - see front matter Ó 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2003.12.008
492
V. Lapenna et al. / Chaos, Solitons and Fractals 21 (2004) 491–500
Fig. 1. Location of the Triponzo station in central Italy.
shallow depth [5,11]. At least 22 (ML > 5:0) local earthquakes had occurred in the studied area between 1279 and 1984 [4]. The focal mechanisms of the mainshocks occurred in 1979, 1984 and 1997, also confirmed by further stress indicators [18], highlighting the active extension processes in the crust characterized by NE–SW and E–W directions [1,8]. Travertine rocks occur in the vicinity of CO2 gas emissions of the Umbria Apennine, evidencing deep-fluid ascent toward surface [3]. Crossing the data of seismic hazard evaluation, the thermal water data-base and historical information on local deep fluid behavior in concomitance with past seismic events [17] have allowed us to previously select the thermal spring located at Triponzo for geochemical monitoring. The Triponzo spring area is located in a limestone formation and is characterized by high permeability and important local karst phenomena. The monitored well is located beside the River Nera at the bottom of a canyon generated by a deep regional fault where the river itself flows [16]. It is worth noting that in the 13th century the Triponzo springs were located a few kilometers away from their present site, and that after the 1328 earthquake, they shifted to their present setting, thereby providing evidence of a close link between local seismicity and thermal water uprising due to faulting. At the Triponzo well site the main chemico-physical water parameter data were recorded over a time-interval of more than two years, from May 14, 1998 to July 7, 2000. As a result of the presence of gaps generated by power supply failures, we considered three time periods for continuous hourly measurements: (1) from May 14, 1998 to September 9, 1998, (2) from July 29, 1999 to October 14, 1999, and (3) from February 24, 2000 to July 21, 2000. The well is 70 m deep and the instrumental equipment (Greenspan multiparametric probe) was let down to a 40 m depth. Electrical conductivity, temperature and pH data were collected in all three periods (CONDi , TEMPi , pHi , with i ¼ 1, 2, 3), while level data (LEV3 ) were also recorded in the last period. In these periods several mild earthquakes, characterized by magnitudes ranging from 3.0 to 4.5, occurred in the investigated area. According to the theoretical findings of Dobrovolsky [6], the Triponzo well site might be sensitive to local seismic events. Recorded signals seem to be influenced by non-stationarities, because data clearly show the presence of trends that could mask possible correlations with earthquakes (Fig. 2).
3. Methods In order to analyze temporal data fluctuations, we performed a robust analysis based on detrended fluctuation analysis (DFA). Detrended fluctuation analysis (DFA) was proposed by Peng et al. [19], and it avoids the spurious detection of correlations that are artifacts of non-stationarity which often affects experimental data. Such trends have to be well distinguished from the intrinsic fluctuations of the system in order to find the correct scaling behavior of the fluctuations. Very often the reasons for underlying trends in collected data are unknown and the scales of the underlying
V. Lapenna et al. / Chaos, Solitons and Fractals 21 (2004) 491–500
493
trends are also unknown. DFA is a well-established method for determining data scaling behavior in the presence of possible trends without knowing their origin and shape [13].
Fig. 2. Hourly geochemical signals recorded at the Triponzo station: (a) from May 14, 1998 to September 15, 1998; (b) from July 29, 1999 to October 14, 1999; (c) from February 24, 2000 to July 21, 2000. The triangles indicate the earthquakes ðM P 3:0Þ occurred during the period of observation and satisfying Dobrovolsky’s rule.
494
V. Lapenna et al. / Chaos, Solitons and Fractals 21 (2004) 491–500
The methodology operates on the time series xðiÞ, where i ¼ 1; 2; . . . ; N and N is the length of the series. With xave we indicate the mean value xave ¼
N 1 X xðkÞ: N k¼1
ð1Þ
The signal is first integrated yðkÞ ¼
k X ½xðiÞ xave :
ð2Þ
i¼1
Next, the integrated time series is divided into boxes of equal length, n. In each box a least-squares line is fitted with data, representing the trend in that box. The y coordinate of the straight line segments is denoted by yn ðkÞ. Next, we detrend the integrated time series yðkÞ by subtracting the local trend yn ðkÞ in each box. The root-mean-square fluctuation of this integrated and detrended time series is calculated by vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N u1 X F ðnÞ ¼ t ½yðkÞ yn ðkÞ2 : ð3Þ N k¼1 By repeating this calculation over all the box sizes, we obtain a relationship between F ðnÞ, that represents the average fluctuation as a function of box size, and the box size n. If F ðnÞ behaves as a power-law function of n, data present scaling: F ðnÞ / na :
ð4Þ
Under these conditions fluctuations can be described by the scaling exponent a, representing the slope of the line fitting log F ðnÞ to log n. For a white noise process, a ¼ 0:5. If there are only short-range correlations, the initial slope may be different from 0.5 but will approach 0.5 for large window sizes. 0:5 < a < 1:0 indicates the presence of persistent long-range correlations, meaning that a large (compared to the average) value is more likely to be followed by large value and vice-versa. 0 < a < 0:5 indicates the presence of anti-persistent long-range correlations, meaning that a large (compared to the average) value is more likely to be followed by small value, and vice-versa. a ¼ 1 indicates flickernoise dynamics, typical of systems in a self-organized critical state. a ¼ 1:5 characterizes processes with Brownian-like dynamics.
4. Results Fig. 3a shows the DFA results concerning the first period; two scaling regions are clearly visible, with a crossover around the diurnal timescale (24 h). The scaling exponents are estimated from the slope fitting the curves by a leastsquares method. We have examined the second region, corresponding to long timescales. The numerical value of a indicate that the signals are not realizations of purely random stochastic processes, but they are characterized by longrange correlation features. Analogous considerations regard the DFA results for the other two periods; only the crossover timescale changes from 24 h in the second period (Fig. 3b) to almost 40 h in the third period (Fig. 3c). As the area under investigation is uninhabited, the crossover 24-h-timescale could be put into relation with meteorologicallyinduced diurnal cycle and the 40-h-timescale could be tentatively attributed to as yet uninvestigated possible earth tidal effects. The numerical values of the scaling exponents a show a slight variation among the different periods; in particular, acond varies from 1.50 to 1.85, atemp ranges from 1.64 to 1.96, apH lies between 1.38 to 1.60, and alev is almost 1.37. In order to probe the presence of locally-correlated sequences possibly connected with the earthquakes occurring during the periods of the measurements, we have constructed a so-called observation box (a mathematical probe) of 103 h in size placed at the beginning of the data, and we calculated a for the data in the box. The box size was in some sense chosen arbitrarily, but sufficiently to evaluate the slope of the local F ðnÞ n. However, the box should be larger than 2 h in order to avoid finite size effects [2,15,23]. The box size has been selected considering that: (i) it has to contain a sufficient number of signal values to perform reliable estimate of the DFA scaling exponents; (ii) it has to be larger than the crossover timescales observed in the DFA results plotted in Fig. 3. Furthermore, we obtained similar results using different box sizes, varying from 750 to 1500 h. Of course, the size must not be larger than the whole observation period. For each period the involved timescales range from the crossover timescale (24 or 40 h) to 1/4 of the box-size. The lower timescale in the record was chosen to take into account the breakpoint of the F ðnÞ n (Fig. 3a–c). Then, we
V. Lapenna et al. / Chaos, Solitons and Fractals 21 (2004) 491–500
495
COND 1 TEMP 1 PH 1
24h
3
α =1. 50+ − 0.04
2
log (F(n))
1 α =1. 96+ 0. 01
0
-1 α =1. 60+ 0.02
-2
-3
-4
0.5
1.0
(a)
1.5
2.0
2.5
3.0
log (n) (hour) 24 h
2.0
COND 2
1.5
TEMP 2
α =1.7 1+ 0.04
PH 2
1.0
log (F(n))
0.5 α =1. 64+ 0.03
0.0 -0.5
α =1. 54+ 0.03
-1.0 -1.5 -2.0 -2.5 -3.0 0.5
1.0
(b)
1.5
2.0
2.5
3.0
log (n) (hour) COND3 TEMP3 PH3 40 h
2
LEV3 α =1. 85+ 0.02
1
α =1.3 7+ 0.03
log (F(n))
α =1. 38+ 0.01
0
-1 α =1.7 2+ 0.03
-2
-3
0.5
(c)
1.0
1.5
2.0
2.5
3.0
log (n) (hour)
Fig. 3. DFA results regarding the geochemical time series recorded at the Triponzo station during the periods: (a) May 14, 1998 to September 15, 1998; (b) July 29, 1999 to October 14, 1999; (c) February 24, 2000 to July 21, 2000. It is clearly visible the power-law behavior with scaling index a varying from 1.37 to 1.96.
496
V. Lapenna et al. / Chaos, Solitons and Fractals 21 (2004) 491–500
Fig. 4. Time variation of the DFA scaling exponents a. The vertical bars indicate the earthquakes with M P 3:4, selected by Dobrovolsky’s rule, that can be responsible for decreasing variations (indicated by arrows in the plots) of the parameters a.
V. Lapenna et al. / Chaos, Solitons and Fractals 21 (2004) 491–500
497
moved this box by 1 h rightwards along the signal sequences and again calculated a. By iterating this procedure for the data sequences, we obtained the measure of the degree of the local long-range correlations a.
Fig. 4 (continued)
498
V. Lapenna et al. / Chaos, Solitons and Fractals 21 (2004) 491–500
Fig. 4 (continued)
V. Lapenna et al. / Chaos, Solitons and Fractals 21 (2004) 491–500
499
Fig. 4 (continued)
Table 1 Space-time-magnitude parameters of the earthquakes plotted in Fig. 4 Latitude (deg)
Longitude (deg)
Depth (km)
Time
Magnitude
From May 14, 1998 to September 9, 1998 42.803 43.171 42.445
13.010 12.815 13.058
5.000 7.758 5.000
Jul 1 09h 10m 18s Aug 11 05h 22m 59s Aug 15 05h 18m 10s
3.4 3.7 4.5
From July 29, 1999 to October 14, 1999 43.018 42.685
12.840 13.108
5.000 5.000
Oct 2 19h 27m 23s Oct 10 15h 35m 52s
3.4 4.0
From February 24, 2000 to July 21, 2000 43.151
12.763
5.000
Jun 11 22h 55m 31s
4.0
Fig. 4 shows the results of this analysis. The plots also show the earthquakes with magnitude M P 3:4 satisfying Dobrovolsky’s rule (Table 1) Although the values of the local a seem to depend on the period and on the parameter, the evolutionary trends are very similar to one another. The first striking feature in the plots is the variability of a, meaning clear temporal variability in the strength of the long-range correlations of the analyzed signals. In almost all the cases we can observe that the local a is characterized by a decreasing behavior corresponding to the occurrence of an earthquake with magnitude M P 3:4. In the first period (Fig. 4a–c), significant decreasing trends characterize all the three parameters, in correspondence with the occurrence of the first two earthquakes; only the last stronger seismic event seems to be not correlated with any negative variation of the signals. In the second period (Fig. 4d–f), it is clearly evident the presence of the 24-h oscillation, that modulates the avariability for all the measured signals. Although acond2 seems much more influenced by the diurnal oscillation, we can observe a decreasing behaviour of this parameter possibly related with the occurrence of the two earthquakes. More evident sharp negative variation characterizes the evolution of atemp2 and aph2 in correspondence with the two seismic events. In the third period (Fig. 4g–j), only one moderate earthquake ðM ¼ 4:0Þ occurred during the period of observation; acond3 , atemp3 and alev3 present a clear minimum in relation with the event; aph3 shows more than one minimum located close to the time of occurrence of the earthquake, superimposed to a negative evolutionary trend.
500
V. Lapenna et al. / Chaos, Solitons and Fractals 21 (2004) 491–500
5. Conclusions Data processing of recorded signals has allowed us to address the following issues: (i) the scaling behavior of the geochemical signals investigated by means of DFA, a method insensitive to the presence of non-stationarities and used to reveal long-range correlations in observational data; (ii) the possible correlation between the temporal evolution of the local degree of correlation, given by the local DFA scaling exponent a, and the earthquakes occurring in the investigated area, not otherwise unambiguously revealed by the observation of the pure variations in the recorded signals. In most cases, the DFA scaling exponent a shows a decreasing behavior in correspondence with the occurrence of moderate earthquakes, suggesting that the temporal signal fluctuations tend to become less smooth during the earthquake preparation process. Similar behavior has been observed in the time evolution of the spectral power-law index of ULF emissions prior to strong earthquakes [9,10] and in the temporal variations of the spectral exponent of geoelectrical signals in correspondence with mild earthquakes [20].
References [1] Amato A, Azzara R, Chiarabba C, Cimini GB, Cocco M, Di Bona M, et al. The 1997 Umbria-Marche, Italy, earthquake sequene: a first look at the main shocks and aftershocks. Geophys Res Lett 1998;25:2861–4. [2] Ausloos M, Ivanova M. Introducing false EUR and false EUR exchange rates. Physica A 2000;286:353–66. [3] Barbier E, Fanelli M. Main fractures of Italy from ERTS satellite images and correlation with thermal springs, volcanoes and earthquakes. In: Aoki H, Iizuka S, editors. Volcanoes and Tectonosphere. Tokai University Press; 1976. p. 321–31. [4] Boschi E, Guidoboni E, Ferrari G, Valensise G. I terremoti dell’ Appennino Umbro-Marchigiano. Area sud orientale dal 98 a.C. al 1984. I.N.G.-S.G.A., Bologna, 1998. p. 267. [5] Deschamps A, Iannaccone G, Scarpa R. The Umbrian earthquake (Italy) of 19 September 1979. Ann Geophys 1984;2:19–36. [6] Dobrovolsky IP, Zubkov SI, Miachkin VI. Estimation of the size of earthquake preparation zone. Pure Appl Geophys 1979;117:1025–44. [7] Favara R, Italiano F, Martinelli G. Earthquake-induced chemical changes in the thermal waters of the umbria region durino the 1997–1998 seismic swarms. Terra Nova 2001;13:227–33. [8] Frepoli A, Amato A. Contemporaneous extension and compression in the Northern Apennines from earthquake fault-plane solution. Geophys J Int 1997;129:368–88. [9] Hayakawa M, Ito T, Smirnova T. Fractal analysis of ULF geomagnetic data associated with the Guam earthquake on August 8, 1993. Geophys Res Lett 1999;26:2797–800. [10] Hayakawa M, Ito T, Hattori K, Yumoto K. ULF electromagnetic precursors for an earthquake at Biak, Indonesia on February 17, 1996. Geophys Res Lett 2000;27:1531–4. [11] Haessler H, Gaulon R, Rivera L, Console R, Forogneux C, Gasparini L, et al. The Perugia (Italy) Earthquakes of 29 April 1984: a microearthquake survey. Bull Seismol Soc Am 1988;78(6):1948–64. [12] Heinicke J, Italiano F, Lapenna V, Martinelli G, Nuccio PM. Coseismic geochemical variations in some gas emissions of Umbria region (Central Italy). Phys Chem Earth 2000;25:289–93. [13] Kantelhardt JW, Konscienly-Bunde E, Rego HHA, Havlin S, Bunde A. Detecting long-range correlations with detrended fluctuation analysis. Physica A 2001;295:441–54. [14] Italiano F, Martinelli G, Nuccio PM. Anomalies of mantle-derived helium durino the 1997–98 seismic swarms of Umbria-Marche, Italy. Geophys Res Lett 2001;28:839–42. [15] Ivanova K, Ausloos M. Application of the detrended fluctuation analysis (DFA) method for describing cloud breaking. Physica A 1999;274:349–54. [16] Lotti B. Descrizione Geologica Dell’Umbria. In: Memorie descrittive della Carta Geologica d’Italia, vol. XXI. Roma, 1926. [17] Martinelli G, Albarello D. Main constraints for siting monitorino networks devoted to the study of earthquake related hydrogeochemical phenomena in Italy. Annali di Geofisica 1997;XL(6):1505–25. [18] Montone P, Amato A, Frepoli A, Mariucci MT, Cesaro M. Crustal stress regime in Italy. Annali di Geofisica 1997;XL(3):741–57. [19] Peng C-K, Havlin S, Stanley HE, Goldberger AL. Quantification of scaling exponents and crossover phenomena in nostationary heartbeat time series. CHAOS 1995;5:82–7. [20] Telesca L, Cuomo V, Lapenna V, Macchiato M. A new approach to investigate the correlation between geoelectrical time fluctuations and earthquakes in a seismic area of southern Italy. Geophys Res Lett 2001;28:4375–8. [21] Thomas D. Geochemical precursors to seismic activity. Pageoph 1988;126:241–65. [22] Toutain JP, Baubron JC. Gas geochemistry and seismotectonics: a review. Tectonophysics 1999;304:1–27. [23] Vanderwalle N, Ausloos M. Coherent and random sequences in financial fluctuations. Physica A 1997;246:454–9.