Robust analysis for the characterization of the seismo-electromagnetic signals observed in Southern Italy

Seismo−electromagnetic signals (SES) are anomalous electromagnetic signals generated as a response to the propagation of a mechan− ical perturbation within the subsoil. Fluid presence plays a key role in determining SES generation and characteristics, therefore SES study could be useful for subsoil characterization. In a more general framework, it can give insight about the role of fluids in the earthquake generation and seismic wave propagation. A systematic study on SES and the related data analysis techniques is fundamental to define the characteristics of these signals which are superimposed to the natural electromagnetic field induced by the external variable magnetic field. To this aim, the Pollino seismic swarm was a great opportunity because continuous MagnetoTelluric (MT) data were recorded in a period in which numerous seismic events of various magnitudes occurred. During the observational period, SES have also been recorded in correspondence to earthquakes far from the MT stations over 800 km. In this paper, we present a procedure aimed to improve the SES detectability and gather as much information as possible on these sig− nals. The procedure is especially tuned for the analysis of MT time series and is based on the application of the Continuous Wavelet Trans− form (CWT) and frequency filters. As will be shown, the operational scheme allows minimizing the background variability of the MT signal facilitating SES detection and the characterization in terms of amplitude and duration.


INTRODUCTION
The occurrence of seismo-electromagnetic signals (SES) associated with the generation and propagation of mechanical perturbation in the subsoil is a matter of fact. Starting from the pioneering work by Frenkel (1944), in which the author postulated equations that estimated the amount of relative fluid motion induced by a seismic wave, several authors developed physical models to explain the growing number of experimental evidence [i.e. Fitterman, 1979;Gershenzon, 1992 and references therein, Pride, 1994;Pride and Haartsen, 1996;Honkura et al., 2009].
The most reliable hypothesis concerning the origin of the coupled seismic and electromagnetic wave propagation is the electrokinetic theory that assumes the existence of an electrochemical double layer along the solidgrain/fluid-electrolyte boundaries of porous media. Within the double layer the charge is redistributed, creating an excess of electrical charge in the fluid along the boundary.
One of the most complete and rigorous physical and mathematical formulations of the equations controlling such phenomena can be found in Pride [1994], in which Biot's poro-elastodynamic equations are coupled to Maxwell's equations of electrodynamics in fully saturated porous media.
Following this theory, seismic waves disturb the fluid excess charge creating an electric streaming current (seismoelectric effect) involving the existence of three kinds of response field.
The first and most common to observe is a coseismic EM anomaly produced by the seismic wave passage in a fluid-saturated porous media. EM and seismic signals are coupled and travel in the medium at the same velocity propagation. The generation of this signal is known as seismoelectric conversions.
The second is a pure EM signal created when the seismic wave crosses an interface with a contrast in electrical and/or mechanical properties and is often referred to as the electromagnetic Interface Response (IR). In this case the EM signals propagate outside the support of the seismic waves at a higher speed.
The third is a pure EM signal directly generated by the seismic source and is related to the relative motion between the fluid and the solid phase in the focal area.
According to the formulation developed by Pride [1994], electric and magnetic coseismic signals are mainly related to different seismic phases. In a homogeneous porous medium, within the seismic wave, an electric field is generated that travels with the P waves, which cause pressure gradients, and a magnetic field is generated by the S waves which cause grain accelerations and set up current sheets [Haartsen and Pride, 1997;Gao and Hu, 2010;Hu and Gao, 2011]. This aspect is not predicted by all the formulations of the SES generation. For example, Gershenzon [1992] provided an analytic expression for the electric and magnetic fields generated by a Rayleigh wave passing through the heterogeneous crust. He concluded that a magnetic anomaly may be generated only when seismic waves cross non-horizontally placed heterogeneities.
For instance, the recognition of pure EM signal generated in earthquake focal area could represents in perspective an important topic as it can be a tool for seis-mic early warning assessment and, in a more general framework, can give insight about the role of fluids in the earthquake generation. From the viewpoint of geophysical sounding, starting from the first attempt to define a SES-based exploration technique [Thompson, 1936], coseismic and IR can now be used, under certain conditions, for the subsoil characterization [e.g. Warden et al., 2012;Revil et al 2015].
In this work, we focus on earthquake-related SES in the attempt to improve the detectability, and thus the characterization, of the coseismic SES. These signals, potentially usable to gather information on the seismic source and the propagating medium, are usually fortuitously recorded during MagnetoTelluric (MT) surveys or monitoring [e.g. Balasco et al, 2014, and references therein] hence, a systematic observational database is still missing. Related to the way in which they are recorded, coseismic SES detection is not a trivial task due to the natural variability of the MT signals and to the unpredictable nature of their seismic source (earthquakes).
To this aim, the present paper focuses on the assessment of a reliable and, possibly, robust method for SES extraction of and characterization. As will be shown, the implemented procedures allow the recovering of SES characteristics associated to both local and regional earthquakes. This work will be used in the future to analyse, in a semi-automatic way, 8 years of MT continuous monitoring data acquired in southern Italy from 2006 to 2014.

OBSERVATIONAL DATABASE: EM MONITORING IN THE SOUTHERN APENNINE, ITALY (HIGH AGRI VALLEY AND POLLINO MOUNTAINS)
To assess the efficiency of the analysis method that we are going to propose, it is necessary to have a statistically relevant database on which to test it. For our purpose, we will use the electromagnetic time series continuously recorded by the MT stations installed in Southern Italy starting from 2007. Originally, the MT monitoring was not focused on SES studies, it was aimed to investigate the stability of the magnetotelluric transfer function [Chave and Jones, 2012] and eventually seismo-related variations in the estimated apparent resistivity. In 2007, a permanent MT station (Balasco et al. 2008) was installed at Tramutola (TRAM-MT: lat. 40.297° N, lon. 15.805° E, elevation 890 m), a very quiet site from the cultural electromagnetic noise point of view (Southern Italy, Agri Valley). By analysing about 4 yr of MT data characterized by a low seismic activity, long-term systematic variations of robust single station MT transfer function estimates were observed in two different sounding period ranges. First, a signifi-

ROMANO ET AL.
cant seasonal component of variability for short periods was noted; these short periods were up to 16 s and were linked to variations in wetting/drying of soil moisture in the shallower layers. Second, a connection between the monitored estimates and global geomagnetic activity, Ap index, was found, particularly in the [20-100 s] period range .
The first evidence of SES presence in the MT time series was noted on 28 th May 2012. Significant changes in amplitude and frequency content of the recorded EM time series were observed in the time series recorded on 28 th May 2012 at 01:06:27. Balasco et al. 2014 successfully linked these anomalies to the M w 4.3 earthquake occurred in the Pollino area at 55.5 km from the TRAM station. The temporal lag between the earthquake occurrence and SES detection at the MT stations allows the identification of this last as coseismic SES. This earthquake was part of a seismic sequence of thousands of small to moderate earthquakes occurred in the Pollino area starting from 2010. The increase of seismicity and the possibility to detect other SES, induced us to install a second MT monitoring station (Campotenese, CAMP-MT station, lat. 39.894° N, lon. 16.085° E, elevation 1487 m) in September 2012 exactly in the epicentral area of the seismic crisis. Both stations are equipped with an MT24LF receiver (EMI Schlumberger), two orthogonal induction coils (BF-4, EMI Schlumberger) that measure the time-varying magnetic field (frequency range 10 -4 -10 3 Hz) in NS and EW direction (H x-NS and H y-EW ), and two 50 m electrical dipoles to measure the electric field in the surface plane (E x-NS and E y-EW ). The TRAM station is also equipped with a vertical induction coil (H z-vert ) and a second dipole E y-EW . It is worth noting that x, y, and z directions correspond to northward, eastward, downward directions in the geomagnetic coordinate system, respectively. The declination angle is about 3° at our observation sites. The horizontal coils are buried in the trenches 0.5 m deep, whereas the vertical coil and the Pb_PbCl2 electrodes are placed in drilled holes 0.5 m in depth. The frequency of electric and magnetic data recording is set to 6.25 Hz with a 24-bit A/D system.
During the 2012-2014 period the large number of earthquakes occurred provided us with the opportunity to observe SES in the 2 MT stations: one is located at 50 km away from the epicentral area and the other in proximity to the swarm. Furthermore, a careful analysis of the EM time series revealed the presence of SES linked to earthquake occurrences in the whole Mediterranean area.

METHODOLOGY
As stated above, this paper deals with the study of earthquake-related SES associated to the passage of seismic waves. To this aim, the first and maybe the most crucial task is their detection within the EM time series. Contrary to what happens for the seismic time series, where even a low magnitude earthquake-related signal is clearly visible in the recorded time series, the same in not necessarily true for SES. EM time series are, in fact, characterized by a natural variability induced by source fluctuations (e.g. variation in of the solar activity) or by the effects of cultural noise, whose amplitude may totally mask SES presence. For SES characterization, it is, thus, necessary to define an operational scheme useful to detect SES presence in the EM time series and perform their characterization.
The one that we propose ( Figure 1) is based on the recovering of the SES characteristics by removing/minimizing the instrumental response and the contribution of all the non-SES-related signals from the EM time series. The proposed scheme is particularly suitable for analysing EM time series recorded by means of MT instruments and takes advantage of both the characteristic of the MT and SES. In what follows, a detailed description and justification of all the extraction steps is given.

SIGNAL PRECONDITIONING
Ad hoc instruments for SES recording are still not available; magnetic and electric sensors, as well as data loggers, are usually borrowed from other geophysical techniques (i.e. MT). Focusing on the MT equipment, both the magnetic induction coils and the electric sensors are similar to a band pass filter (Figure 2) in which the extension and the position of the flat area is a function of the investigation depth to achieve. To simplify, the higher the investigation depth is, the more the central frequency of the flat area is shifted toward lower frequencies. For a more detailed explanation, the reader is invited to read a general MT textbook [e.g. Chave and Jones, 2012].
The deconvolution of the instrumental response from the recorded time series is necessary, prior to any attempt of SES characterization, to avoid that the amplitude of any SES is artificially modified by falling in the not-flat area of the sensor response. This is particularly true for the magnetic field since the instrumental response has a relatively narrow (from 0.1 Hz to 100 Hz) flat area (Fig. 2, lower panel) but it may be not mandatory for the electric fields for which we have a flat response in a much wider frequency range (Figure 2, upper panel).
The deconvolution is easily performed in the frequency domain. Starting from the unprocessed time series, the transition to the frequency domain is performed by applying a standard FFT algorithm. The deconvolution is then made rescaling the FFT output by sensor response. The deconvolved signals are then converted in the time domain by applying an inverse FFT algorithm.
The above-mentioned procedures for MT data cannot be performed on the whole-time series. A segmentation of the input time series is required by the characteristics of the sensor response behaviour that is evaluated only in a limited frequency range.
For a posteriori SES detection and characterization, the segmentation of the time series does not represent a problem. Being linked to the earthquake generation and the passage of the seismic waves in the area where the MT station is installed, the only necessary condition is that the recording segment must include the time of the earthquake triggering and must be long enough to cover the seismic wave passage.

WAVELET ANALYSIS
The output signal of the preconditioning block is the superposition of the background EM fluctuations and the SES where and when these are present. For SES characterization (amplitude, duration, etc), background fluctuations have to be considered noise and removed as efficiently as possible. To this aim, any standard filtering procedure can be adopted as long as it preserves SES characteristics. For instance, if information on SES frequency content is known or acquired, a simple frequency-based filter can be used.
To infer hints on the spectral content of the earthquake-related SES, we use the Continuous Wavelet Transform (CWT) which allows analysing localized variation of power (as the SES are) within a time series. By decomposing a time series into time-frequency space, CWT identifies both the dominant modes of variability and how they vary in time. The basics of the CWT are described in a number of papers and books [Mallat, 1998;Torrence and Compo, 1998 and references therein]. Here only the CWT aspects relevant to the analyses are detailed. As mother wavelet employed, we utilize the Morlet wavelet to estimate the local and global frequency content of the series. This choice is classical in harmonic time-frequency analyses as the Morlet is sine-like shaped and limited in time [Morlet et al., 1982;Torrence and Compo, 1998].
The natural electromagnetic field can be described with a background red noise spectrum (increasing power with decreasing frequency); conversely SES are non-stationary events in a limited range of frequency, which coincides with the highest frequency range we are recording. When SES occur, thus a local maximum in the CWT is found. If this maximum is significantly above this background spectrum, then it can be assumed to be a true feature with a certain percent confidence [Torrence and Compo, 1998]. In the present paper, when CWT results are presented, the regions contoured in black individuate the significant ones at 95% confidence level.
CWT has been applied to all the segments of the EM time series in which the presence of earthquake-related SES was visually recognizable. CWT results for all the analysed segments allowed us to assess that the earthquake-related SEs are characterized by a frequency content similar to the one of the seismic waves responsible for their generation and with a lower frequency limit, for the P and S seismic phase related SE signals, of ~ 0.5 Hz.
This result confirms what shown in Balasco et al, [2014] and Balasco et al., [2015] and provides the information needed for a correct choice of the parameters of the frequency-based filter (i.e. corner frequency). Moreover, it clarifies why SEs are relatively easy to identify (even without deconvolving the instrumental response) in the MT time-series. SES frequency content falls, in fact, in the MT dead band [Simpson and Bahr, 2005], the frequency range in which the natural EM field has the minimum amplitude and in which MT instruments have the maximum sensitivity (Figure 2).

SIGNAL FILTERING
The information coming from the application of the CWT is crucial for a correct definition of the characteristics of the frequency-based filter that will be used to minimize the EM background contribution in the timeseries and consequently improve the SES detectability. From all the analysed cases, the SES have a bounded content of frequencies that spam from 0.5 Hz (lower limit as retrieved from the CWT) to ~ 3 Hz. It is worth noting that the high frequency limit is imposed by the recording settings (the time series were acquired with a sampling frequency of 6.25 Hz) and doesn't necessary represent a real bound for the SES frequency content. It cannot be excluded, thus, that the SES frequency content is wider than the one here indicated and that it extends over the 3 Hz. Moreover, as shown by combined analysis of co-located seismic and EM signals, the electrical signals have higher frequency content above 1 Hz than seismics [Balasco et al., 2015] confirming the hypothesis of the proportionality of the electrical field with the acceleration of the seismic wave, as expected during coseismic effects [Garambois and Dietrich, 2001].
In the light of what stated above, we decided to apply a Finite Impulse Response (FIR) filter [Shenoi, 2005] with a corner frequency usually of0.5 Hz (but that can be changed on the basis of CWT results) for high pass filtering the EM time series. The filter was applied to the time series by using the Matlab© function "filtfilt" which performs zero-phase digital filtering by processing the input data, in both the forward and reverse directions, thus giving an output with zero phase distortion.

SES DETECTION AND CHARACTERIZATION
The filtered time series can then be examined for the identification of SES and the extraction of their characteristics. In case of a negative detection of SES, the original time series are reprocessed by applying the same scheme but reconsidering the CWT results. When the previously used filter parameters are recognized to be inappropriate for this specific case, a modification of the filter block is performed and the time series filtered. In the case of both a further negative detection and positive one, the time series is flagged as "anomalous" for further studies.
The output of the presented operational scheme is a "clean" EM trace dominated by SES. It can be hence used for a precise estimation of the SE arrival time and their amplitude. As will be shown by the examples presented for the SES generated by small amplitude earthquake, the operational scheme allows the extraction of the SES characteristic which are practically undetectable in the unprocessed time series.

RESULTS
To test the efficiency of the operational scheme, we report two cases study.
The first is related to SE signals generated by a regional earthquake and for which a certain ambiguity in the estimation of the arrival time and amplitude is present. The second is referred to a SE signal occurred as a consequence of a very small magnitude local event and which is difficult to detect/characterize in the unprocessed time series.

REGIONAL EARTHQUAKE
On 12 th October 2013, a M w 6.4 earthquake occurred near Crete, Greece, (Lat:35.52 N; Long: 23.28 E: depth 52 km, ISIDE database -Italian Seismological Instrumental and parametric Data-BasE). Several seconds after the earthquake triggering. both the TRAM and CAMP MT stations recorded anomalous EM signals. Figure 3 shows the deconvolved electric field (N-S component) time series recorded by the TRAM station. The origin time (O.T) is set on the earthquake occurrence time (13:11:53 UTC). The anomalous fields, visible starting from ~ 110 seconds after the O.T., have a variable amplitude difficult to estimate due to the fluctuations of the natural background signal. The SES delayed occurrence is comparable with the travel-time of the seismic waves to cover the ~850 km which separate the MT station from the hypocentral area.
Being the SES presence clearly visible, the application of the presented operational scheme is aimed to the background removal/minimization and a better definition of the SES characteristics.
Following the flow chart of Figure 1, the detrended time series were analysed by applying the CWT and on the basis of the results obtained, high-pass filtered. These operations are summarized in Figure 4 in which the upper panel shows the input time series (derived by those shown in figure 3 but standardized), the middle panel the results of the CWT and the lower panel the filtered time series.
As it is possible to appreciate, the CWT clearly delineates SES features both in the time and frequency domain. The SES is characterized by two main wave trains (respectively between ~110 s and ~160 s and between ~190 s and ~250 s after the O.T.) with most of the spectral power for frequencies above 0.5 Hz (2 seconds on the vertical scale of fig. 4b). The EM signal background fluctuation is generally characterized by a lower frequency content. Considering this, the high-pass filter has been applied. The resulting time series (Figure 4) does not show any long period variation anymore allowing a correct estimation of the SES amplitude. For this specific case, the presence of two distinct wave trains can be associate to the propagation of two different types of seismic waves: the body waves and the surface waves.
The effectiveness of the adopted processing scheme is better appreciable on the magnetic time series ( Figure  5). In this case, starting from time series in which the SES presence could only be guessed (Figure 5a) due to a higher variability of the background signal, we retrieved a "clean" time series (Figure 5c) that can be used for further analysis.

LOCAL EARTHQUAKE
On 6 th June 2014, a M w 4.0 earthquake occurred in the Pollino area (Lat:39.90 N; Long: 16.09 E: depth 8 km, 13:43:34 UTC, ISIDE). This earthquake was part of a seismic sequence in the area and was followed by a M L 2.6 event (Lat:39.90 N; Long: 16.12 E: depth 8 km, 13:41:38 UTC). Due to the fact that the CAMP MT was located only few kilometres far away from the epicentral area, the larger earthquake generated visible EM anomalies ( Figure 6) easily detectable also without any particular processing procedure. Conversely, the smaller event did not generate clear SES (in terms of occurring time and amplitude).
In what follows, in order to prove the effectiveness of the processing scheme in enhancing the SES detectability, we analysed a segment of EM time series which contains the SES anomalies generated by the M L 2.6 event.
By focusing on EM time series segment shorter than the one shown in fig 6, it is possible to observe how SES represents a relatively small anomaly whose amplitude is difficult to estimate due to the presence of a decreasing trend of the back-ground signal (Figure 7a). The CWT (Figure 7b) clearly highlights the SES presence and gives insights about its duration and frequency content. The SES lasts for a relatively short time (less than ten seconds) and as for the case of the regional earthquake, most of its frequency content can be found for frequencies higher than 0.5 Hz. Based on this, the time series were filtered ( Figure 8) and the SE features extracted. It should be noted the absence of the SES in the Hy magnetic component (recorded along the E-W direction, second panel in Figure 8). This result is also confirmed by the CWT.

EXTENSION OF THE APPLICATION OF THE PRO-CESSING SCHEME TO A LARGER DATASET
Verified the effectiveness of the processing procedure, we used it to analyse part of the MT time series recorded during the Pollino seismic swarm. From the ISIDE database, we randomly selected 36 seismic events occurred in the Pollino area from November 2011 to June 2014, with a magnitude ranging from 4 to 2, and we looked for SES in the MT recording. According to the empirical relationship between SES amplitude, the seismic event magnitude and the distance from the observatory site , we analyzed only CAMP MT data since the lowest mag-nitude events should be not observable at the TRAM MT (~ 50 km far away from the source area).
The results of this analysis are shown in Table 1. The first 4 columns report the information on the seismic events (occurrence time, magnitude, depth and distance from the observatory site); the remaining columns report the semi-amplitude (peak to peak) of the SES, when observed, in all the recorded channels. For each channel, it is also reported the difference (expressed in percentage) between the amplitude of the SES extrapolated respectively from the unfiltered and filtered time series. In these columns, the label "detected" indicates the detection in the filtered ROMANO ET AL.  time series of SES not visible in the non-processed data. As can be observed, relevant differences (larger than 10%) in the semi-amplitude of the SES are reported especially for the magnetic recordings for which the minimization of the background signal (performed in the proposed processing scheme) is particularly effective.
The results of Table 1 also confirm the directionality of the SES; the "n. d." label (which stands for "not detected") is mainly retrieved for the EM signals recorded along the E-W direction.
The observed directionality can be explained or with dependence of the SES recorder waveform to the observation orientation [Gao and Hu, 2010;Hu and Gao, 2011] or by the characteristic of the electrical structure of the subsoil beneath the observing station (CAMP MT) [Gershenzon, 1992].

CONCLUSIONS
The seismo-electromagnetic signals are generated as a response to the propagation of a mechanical pertur-bation within the subsoil. Performing a recording of the natural electromagnetic fields (as it is usually done in the MT method), SES occurrence can be recognized as a transient anomaly with very specific features (in terms of amplitude, frequency content and duration). SES amplitude is generally connected with the one of the generation machanisms (i.e. the magnitude of the earthquake) and with the distance between the recording station and the source. Excluding the cases when the source is either very energetic or close to the recording station, SES amplitude is usually small and of the same magnitude order as the natural background signal. Nevertheless, thanks to the fortuitous coincidence that places the SES characteristic frequency range in the so-called "magnetotelluric Dead-Band", their detection and characterization is generally possible.
The interest in the SES characterization is linked to their generation mechanisms. By studying the SES it should be possible to get new insights for the understanding of the role of fluids in the earthquake generation and seismic waves propagation and the characterizing of the subsoil in terms of fluid presence.