Introduction

In the last few decades on the ground-based observations, more data regarding the existence of the electromagnetic precursors related to the earthquakes have been accumulated (Biagi et al. 2011; Fraser-Smith et al. 1990; Fenoglio et al. 1995; Han et al. 2015; Hattori et al. 2013; Hayakawa and Fujinawa 1994; Huang 2011; Johnston 1997; Kopytenko et al. 1994; Nagao et al. 2002; Stanica and Stanica 2010, 2011, 2012; Stanica et al. 2015; Uyeda et al. 2011; Uyeda 2013, 2015; Varotsos 2005). The satellite observations over the seismic areas indicate changes of the ionosphere state before the earthquakes (Błęcki et al. 2010, 2011; Hayakawa et al. 2011; Parrot 1995; Parrot et al. 2015). Even if all the above works have presented a posteriori observations related to the different pre-seismic signatures, their very optimistic results, together with the theoretical models of fracture induced electromagnetic phenomena (Rabinovitch et al. 2007; Hunt et al. 2007), laboratory experiments on the rock samples (Freund 2000; Freund et al. 2006), ionosphere state observation by the vertical magnetic component (Ernst et al. 2010), the correlation of resistivity variations and geodynamic processes (Bataleva et al. 2013) and very long time electromagnetic data acquisition, processing and analysis realized for a unique and very active seismic zone in Romania, offer a great advantage in emphasizing useful information regarding the possible interrelations between the Vrancea’s earthquakes generation and geomagnetic signature. Thus, in this work, we focused on the geomagnetic methodology able to identify in the ultra-low frequency (ULF) range a pre-seismic anomalous signal related to the Mw5.7 earthquake triggered on September 24, 2016, and to bring strong arguments in favor of the need to asses such information, daily, on the institute website (http://www.geodin.ro), as it is shown in Fig. 1, with the aim to promote an international scientific co-operation.

Fig. 1
figure 1

Daily mean distribution of the geomagnetic parameters Bzn and Bzn* presented on the Institute of Geodynamics of the Romanian Academy website (photo image captured on September 28, 2016). White vertical bar is Bzn; red vertical bar is standard deviation (SD); red full circle is earthquake; the ratio 5.7/92 km is earthquake magnitude/hypocentre depth; blue vertical bar is Bzn*abs (abs means absolute value); red dashed line is threshold for anomaly; red vertical arrow indicates the pre-seismic signal occurrence; the explanation of the parameter Bzn and Bzn* will be done below

Full size image

Methodology, data processing and analyzing

An earthquake of Mw5.7 struck the Vrancea zone, Romania, on September 24, 2016 at 02:11:20 local time (LT). The epicenter was located at the geographic coordinates 45.710N and 26.620E, with a focal depth at 92 km and, about 150 km NE of Bucharest, as it was determined by the Euro-Mediterranean Seismic Centre (http://www.emsc-csem.org). Although this earthquake was moderate, it was felt in Bucharest and generated panic to a certain extent through its inhabitants, being intensively mediated on the TV “Breaking News”, as seismologists systematically claimed a possible imminence of a large earthquake, based on the recurrence rate of about 40 years, on the analogy of the previous two earthquakes incidence of Mw7.7 on November 10, 1940 and Mw7.5 on March 4, 1977, respectively.

Basic concept of the geomagnetic precursor

It was demonstrated by Word et al. (1970) that for a 2D geoelectric structure, the vertical component (Bz) of the geomagnetic field is produced mainly by the horizontal magnetic component perpendicular to strike (\({\text{B}} \bot\)) and, as a consequence, the normalized function (Bzn) having the form:

$${\text{Bzn(}}f ) { = }\frac{{{\text{Bz(}}f )}}{{{\text{B}} \bot (f )}}$$
(1)

should be time invariant and it becomes unstable due to the geodynamic processes related to the intermediate-depth seismicity, being associated with the resistivity changes along the good conducting pass in the lithosphere (Stanica and Stanica 2011), according to the relation (2):

$$\left| {{\text{Bzn(}}f )} \right|{ = }\sqrt {\frac{{\rho_{\parallel } (f )}}{\rho z (f )}} ,$$
(2)

where \(\rho_{\parallel }\) is resistivity parallel [Ωm] to the geoelectric strike, ρz is vertical resistivity [Ωm] and f is frequency [Hz].

Consequently, the existence of a 2-D structure gives rise to normalized function Bzn that has the magnitude proportional to the intensity of the geomagnetic component perpendicular to the geoelectric strike, which is in turn determined by the electric field concentrations released in pre-seismic conditions and propagated along the strike, due to the torsion process of the Vrancea’s seismogenic volume (Stanica et al. 2004; Stanica and Stanica 2012).

To fulfill the conditions imposed by relation (1), the Geodynamic Observatory Provita de Sus (GOPS) is placed on the Carpathian electrical conductivity anomaly (CECA), at about 100 km westwards of the seismic active Vrancea zone (Fig. 2). This conductivity anomaly is generated by a 2D geoelectric structure (Pinna et al. 1993; Stanica et al. 1999) and it was identified in the ULF range (0.001–0.0083 Hz), on the basis of the magnetotelluric dimensionality parameters skewness and strike (Stanica and Stanica 2010).

Fig. 2
figure 2

Map with crustal (black dots) and intermediate (read dots) earthquakes in Vrancea zone. Carpathian Electrical Conductivity Anomaly (yellow dashed line); epicenter of M5.7 earthquake (yellow star); Geodynamic Observatory Provita de Sus (blue star)

Full size image

According to relation (3), the ULF range (0.001–0.0083 Hz) may be associated with the intermediate-depth earthquakes interval (60–180 km) where possible geomagnetic signals are generated. This supposition is based on the electromagnetic skin depth relation, which in this study is equated with penetration depth of the electromagnetic field into the Earth and can be approximated as:

$$p(f) \approx 500\frac{1}{{\sqrt {\sigma f} }},$$
(3)

where p is the penetration depth (m), σ is conductivity (S/m) and f is frequency (Hz).

If in relation (3) it is assumed that the Earth’s lithosphere has an average conductivity of about 10−2 S/m and f = 0.001 Hz (minimum value in ULF range), then the geomagnetic signal generated in Vrancea’s seismogenic volume, at about 160 km depth, can reach to the GOPS.

To identify the distance for pre-seismic signal detection, depending on the earthquake magnitude, we used the Morgunov and Malzev (2007) relation:

$$R^*({\text{km}}) \, = \, 10^{0.5M - 0.27} ,$$
(4)

where R* is epicentral distance and M is earthquake magnitude.

In particular, for the Mw5.7 earthquake that occurred on September 24, the epicentral distance R* ≈ 380 km. As the distance between GOPS and the earthquake epicenter is about 100 km (Fig. 2), the condition for pre-seismic signal detection imposed by relation (4) is fulfilled.

Satellite observations

The registrations performed by the Swarm satellites over Vrancea zone, before the M5.7 earthquake, are used as an additional information about variations in the ionosphere. Swarm is a constellation of three satellites to measure the Earth’s magnetic field and identify the sources of its variations originated from core, ionosphere, magnetosphere, mantle, crust and as well as the oceans. This mission consists of the three identical Swarm satellites (A, B, and C), which were launched on 22 November 2013 into a near-polar orbit. Swarm A and C form the lower pair of satellites flying side by side (1.4° separation in longitude) at an altitude of about 470 km (inclination angle is equal to 87.30°), whereas Swarm B is cruising at higher orbit of about 520 km (inclination angle is equal to 87.75°). They are equipped with a set of six identical instruments—absolute scalar magnetometer, vector field magnetometer, star tracker, electric field instrument, GPS receiver, and accelerometer. Further, we use data from vector field magnetometer to study magnetic field variations and Langmuir probe (being a part of electric field instrument) to provide the representation of the main filed.

The magnetometers installed on-board of Swarm satellites measure main magnetic field with sampling rate 50 Hz (3 components and absolute value). On this way, we are looking for the ionosphere effects associated with seismic activity, that are mainly seen in the variations of the electromagnetic fields in very broad frequency range from fraction of Hz up to several Hz (Olsen et al. 2013).

Geomagnetic and satellite data processing and analysis

Geomagnetic data

The ground-based monitoring system used for this study is installed at GOPS (Fig. 2) and consists of: (1) data logger (MAG-03 DAM) and a three-axes fluxgate magnetic sensor (MAG-03MS), both used for the geomagnetic time series (B⊥ and Bz) collection with a sampling rate of 60 s; (2) Computer with programs dedicated to the data acquisition, storage and daily transfer, via internet cable, to Bucharest at the Institute of Geodynamics of the Romanian Academy (IG-RA).

The automatic system installed at IG-RA is composed by: (1) Computer (server) used to receive and storage daily the geomagnetic data records; (2) work-station with specific programs for data processing, analysis and display on the institute webpage.

In this study, the daily mean distributions of the Bzn and Bzn* with their standard deviation (SD) are analyzed on the span of time 01–30 September 2016, by the following procedures:

  • FFT Band-pass filtering (BPF) analysis in the ULF range was applied to Bzn time series, for two successive time windows of 1024 samples, with about 30% overlapping on 1440 data acquired each day (Table 1; Fig. 3);

    Table 1 B⊥, Bz, Bzn and Bzn (FFT-BPF) time series obtained for 20 min on September 22, 2016, Bzn mean and Bzn (FFT-BPF) mean are daily averaged values
    Full size table
    Fig. 3
    figure 3

    FFT Band-pass filtering (red line) applied on Bzn time series (blue line) for a time windows of 1024 samples recorded on September 22, 2016

    Full size image

  • Statistical analysis based on the standardized random variable has been used to discriminate the pre-seismic anomalous intervals as follows:

$${\text{Bzn}}^* = \frac{{X - \overline{X} }}{{\overline{Y} }},$$
(5)

where \(X\) is daily mean value of Bzn, \(\overline{X}\) is mean value of Bzn obtained for 30 consecutive days before \(X\), \(\overline{Y}\) is mean value of SD obtained for 30 consecutive days before \(X\), \({\text{Bzn*}}\) is threshold for anomaly using SD.

On this way, it is expected to eliminate seasonal variation of the normalized function Bzn and to discriminate with high accuracy the pre-seismic geomagnetic signature related to the seismic event.

A very important element of the initial stage of interpretation is to ascertain that the geomagnetic component Bz is not of external origin, i.e., it is not an over-ground effect of ionosphere currents. We verified it on the basis of records from two magnetic observatories belonging to the INTERMAGNET network: Surlari (Romania) and Tihany (Hungary).

We divided the Bz component variations into the internal (induction) and external parts (Ernst and Jankowski 2005; Ernst et al. 2010). It is worth noting that the external part of Bz variations is often much bigger than the induction (internal) one.

Such a situation is shown in Fig. 4, from which it is also evident how important it is to take into account and analyze the external part Bze in the study of pre-seismic processes.

Fig. 4
figure 4

Twelve-hour time series from Magnetic Observatory Surlari. Bz is total vertical component; Bzi is internal (induced) part of vertical component; Bze is external part of vertical component

Full size image

In Fig. 5, we present 5-day simultaneous time series of the horizontal geomagnetic components B x and B y and the external part of Bz (denoted Bze) recorded at the observatories Surlari and Tihany, where it is clearly emphasized that changes of Bze in the both observatories are negligible (< 3 nT), which makes it sure that Bz used in relation (1) is not of external origin, and the anomalous behavior of the Bzn may be interpreted as the earthquake precursor.

Fig. 5
figure 5

Five-day simultaneous time series from Magnetic Observatories Surlari and Tihany. B x S and B y S are horizontal components for at Surlari; B x T and B y T are horizontal components from Tihany; Bze S is external part of vertical component from Surlari; Bze T is external part of vertical component from Tihany

Full size image

Satellite data

The registrations performed by the Swarm satellites limits our studies only to very low frequencies called ULF (ultra-low frequencies up to single Hz) and ELF (extra low frequencies). The magnetic vector is measured with a frequency of 50 Hz, what gives theoretically the possibility to study variations up to 25 Hz, but the filter inside the instrument limits this range to 15 Hz.

Next problem of our studies is related to the extremely small values of these variations (typically fraction of nT) in comparison with the main geomagnetic field measured in tens thousands of nT. Their observation requires sensitive magnetometers and carefully selected methods of signal analysis. The use of the vector field magnetometer on board Swarm, intended for observation of the main geomagnetic field, requires application of special filtering techniques.

Steps distinguished in the data processing chain are as follows:

  1. 1.

    δBi residuals retrieval from the measured signal for three B components;

  2. 2.

    FFT transformation applied to the δBi residuals;

  3. 3.

    Generation of time–frequency spectrograms for δBi along the orbit in the frequency range up to 25 Hz;

  4. 4.

    Integration of lightning’s database with derived Swarm δBi spectra.

In the δBi residuals retrieval procedure, a second-order polynomial approximation is applied to set of 1024 samples of the VFM waveform measurements. Also, non-linear least squares method and the Levenberg–Marquardt algorithm are used to provide the representation of the main filed. Thus, δBi is obtained as a difference between measured signal and computed best fit. In the next step FFT is applied, giving power spectrum for each component and intensity of the B field residuals, in the frequency range up to 25 Hz. The FFT window is 1024 samples wide, advanced by 256 samples for overlapping.

Results and discussion

The continuous monitoring of the Bzn and Bzn* time series realized in the last 10 years at GOPS has shown small variability related to their normal trend observed in non-geodynamic conditions and has emphasized consistent anomalies some days before the onset of the seismic event.

To have a comprehensive view on the applied methodology, in this work the daily mean distributions of the Bzn and Bzn* obtained in September 2016, in correlation with Mw5.7 earthquake, are presented in Figs. 6 and 7.

Fig. 6
figure 6

Daily mean distribution of Bzn and SD on September 2016. Blue line is Bzn; red vertical bar is SD; red star is earthquake; the ratio 5.7/92 km is earthquake magnitude/hypocentre depth, in km

Full size image
Fig. 7
figure 7

Daily mean distribution of the Bzn* and earthquakes magnitude obtained on September 2016. Blue bar is Bzn*abs; red full circle is earthquake; the ratio 5.7/92 km is earthquake magnitude/hypocentre depth in km; red dashed line is threshold for anomaly using SD

Full size image

In Fig. 6, there is a significant anomalous domain of minimum, extended on the interval September 19–September 26, with values ranging from 1.933 to 1.942, easily identifiable on the Bzn distribution. This abnormal change of Bzn may be associated with the increasing magnitude of B⊥ (see relation 1) which is basically controlled by the electric charges released and propagated along the CECA, prior the onset of the seismic event, as the result of high stress reached in seismogenic volume (Stanica and Stanica 2012).

The Bzn* time series, displayed in absolute value (Bzn*abs) on the same time interval as Bzn, is shown in Fig. 7. On the Bzn* anomaly, extended on the interval September 21–25, a pre-seismic geomagnetic signature corresponding to a magnitude greater than 5·SD (red dashed line) is emphasized starting with 3 days before the onset of the Mw5.7 earthquake on September 24.

The analysis of the ULF/ELF emissions and electron density in the ionosphere in the vicinity of the earthquake epicenter has been performed for 4 days, on the time interval September 20–23 and 1 day before the Vrancea earthquake. A very weak effect of increasing of the intensity and variations of the electron density has been registered close to the epicenter. For these days, the disturbances registered by Swarm are presented below in details.

Figure 8 shows the part of orbit of the Swarm B satellite during observations on September 20. The color of the dots corresponds to the intensity of the magnetic field variations in the determined point. The enhancement of the intensity is seen in the closest vicinity part to the epicenter.

Fig. 8
figure 8

The part of Swarm B orbit related to the flight over Vrancea zone on September 20. The color of the dots on the orbit line corresponds to the intensity of the magnetic field variations

Full size image

In Fig. 9 are shown the spectra of the magnetic field variations and changes of the electron concentration in the vicinity (around 600 km) of the Vrancea earthquake epicenter 4 days before this event. The variations of the electron density are seen in time interval between 08:17 and 08:19UT. The satellite flown by, in this time, a distance of about 900 km. The variations of the magnetic field were present during this entire interval. The effect shown in this case is very weak in comparison to the effect reported from DEMETER registration (Parrot et al. 2006; Błęcki et al. 2010), but clearly seen. The spectra have a maximum in the lowest part of the frequency range and correspond to the range below the oxygen ions O+ gyro-frequency which is of the order of 34 Hz. It can be associated with Alfven waves.

Fig. 9
figure 9

The spectra of the magnetic field variations (upper panel) and electron concentration registered by Swarm B satellite in the vicinity of the Vrancea zone on September 20

Full size image

Figure 10 presents spectra of the magnetic field variations and electron concentration taken on September 23 (1 day before the discussed earthquake) from half orbit crossing the Vrancea zone. One can see strong effects associated with high latitude regions (auroral oval, ionospheric through) in comparison with the weak variations in the vicinity (distance about 500 km) of the epicenter. The closest vicinity was at 16:04UT. The disturbances are seen around 16:04 to 16:05UT.

Fig. 10
figure 10

The same as in Fig. 9, but for the September 23 and taken for entire half orbit crossing area in vicinity of the Vrancea zone. White line represents the position of the satellite

Full size image

The weakness of the variations seen by Swarm satellites in comparison to reported after DEMETER mission can be related to the smaller effects in the magnetic field than in electric measured by DEMETER, the earthquake was also weaker and may be the local conditions not amplified the electromagnetic effect.

Conclusions

In this paper, we have investigated the ULF geomagnetic data recorded on September 2016 to find a possible pre-seismic signal associated with Mw5.7 earthquake occurred on September 24. On this interval of analysis, a very clear and unusual anomaly of minimum, extended from September 19 to September 26, has been detected on the Bzn distribution. The new Bzn* time series, obtained after applying a statistical analysis based on Eq. 5, indicates that the variability of Bzn is not at random one, this being a significant and reliable anomalous pre-seismic effect associated with Mw5.7 earthquake. The both results suggest that the pre-seismic anomalous behavior of the Bzn and Bzn*, emphasized by Figs. 6 and 7, was triggered with 3 days before the onset of the seismic event. Complimentarily, the analysis of external part of the vertical geomagnetic field, recorded in the observatories Surlari and Tihany, demonstrated that vertical component (Bz) recorded at GOPS is of internal origin and, consequently, the pre-seismic anomalous signatures observed in the Bzn time series may be interpreted as the earthquake precursor. The measurements originated from Swarm satellites also indicate some disturbances of the magnetic field and electron concentration over Vrancea zone with 3 days before the Mw5.7 earthquake, but the effect was rather small because the earthquake was weaker.