Multistationary Geomagnetic Vertical Intensity Polarization Anomalies for Predicting M ≥ 6 Earthquakes in Qinghai, China

Abstract

Single-stationed geomagnetic vertical intensity polarization (GVIP) anomalies have demonstrated good predictions of the occurrence of large earthquakes in Japan. Nonetheless, due to the lack of a previously densified geomagnetic network, how the multistationary GVIP anomaly (MGVIPA) corresponds to impending earthquakes remains poorly understood. Based on the newly constructed geomagnetic network from 2014 in Qinghai, China, which is composed of 23 electromagnetic stations, we suggested an MGVIPA method to analyze the correlation with large earthquakes since 2015. The results show that (1) the occurrence of MGVIPA is characterized by clusters in time that continue in a short period; (2) the spatial distribution of MGVIPA usually occurs with high values synchronously at several places over the same period; and (3) the Mw ≥ 6 earthquakes occurred in the regions indicated by MGVIPA within a period ranging from 3 months to 1 year from 2015 to 2021 in Qinghai, China.

1. Introduction

Previous studies on geomagnetic numerical simulations have shown that in magnetic field signals with frequency bands of 10⁻²~10² Hz, the amplitude of the vertical component of the primary source from the magnetic signal inside the crust observed at the surface is greater than or close to the amplitude of the horizontal component (e.g., [1]). The ratio of the amplitude of the vertical component to the horizontal component of the magnetic field from the Earth’s crust is greater than or close to one. However, based on electromagnetic induction theory, the ratio of the amplitude of the vertical component to the horizontal component of the magnetic signal from the outer space of the ionosphere is less than one. Hayakawa et al. [2] developed the method of geomagnetic vertical intensity polarization (GVIP) to extract the high amplitude ratio of the geomagnetic vertical component to the horizontal component and obtain the magnetic field signal from the subsurface. By analyzing the changes in polarization values before and after the Mw 8.0 earthquake in Guam on 8 August 1993, and using magnetic fluxgate magnetometer observations at the Guam Geomagnetic Station, which is located 65 km away from Guam, studies found that the polarization value in the 0.01–0.05 Hz (100–20 s) band gradually increased in the two months before the earthquake, reaching a maximum at the time of the earthquake, and gradually recovered after the earthquake (e.g., [2]). Numerous seismic events have since been identified anomalously through high polarization values in the days to three months prior to the earthquake (e.g., [3,4,5,6,7]). Hobara et al. [8] subjected the polarization results to a more refined frequency segmentation, and more significant anomalous signals of polarization values were extracted in the frequency bands of 0.02–0.022 Hz (50–45 s) and 0.05–0.1 Hz (20–10 s) before the 1997 Guam Mw 8.0 earthquake and the 2000 Izu swarm earthquakes. Based on fluctuation theory, a three-dimensional anomalous body model was constructed, and the seismic geomagnetic disturbance signal coinciding with the anomalous polarization value was obtained by numerical simulation (e.g., [8]). Feng et al. [9] found that geomagnetic vertical intensity polarization has an annual variation characteristic and eliminated it before extracting the high value of the polarization anomaly. The high value of the polarization anomaly at the Kashi station after eliminating the annual variation corresponds well to the earthquakes around the station that occur 2 months after it, with most of the earthquakes occurring within 1 month [9].

Although some earthquake predictions have been obtained by the GVIP anomaly method, the reliability of this method for extracting information on seismic geomagnetic anomalies is still unknown. First, selecting data at midnight does not fully prevent the influence of the magnetospheric current system on the elevated polarization values, even though many studies have demonstrated that high polarization values are negatively correlated with the geomagnetic activity index (e.g., [9,10,11,12,13]). Second, due to the limitations of station observation data quality and station spacing, pre-seismic anomalies usually exist only at the single station with the closest epicenter, and the reliability of pre-seismic anomalies is often questioned when data interference cannot be excluded. However, based on the small-aperture array observation data established in northwest Yunnan, Li et al. [12] extracted the pre-seismic polarization anomalies at several stations using the polarization method. Their results showed that the polarization values of each station had abnormally high fluctuations half a month before the earthquakes occurred and returned to normal one week before the earthquake, and the amplitude of the anomalies gradually decreased with increasing distance from the epicenter, which was consistent with the attenuation characteristics of seismic electromagnetic disturbance signals (e.g., [14]). Feng et al. [9] obtained magnetic flux gate second sampling data from nearly 100 geomagnetic stations nationwide to study the relationship between Mw > 6 events and polarization value anomalies in western mainland China (west of 110° E) from 2015 to 2018 and found that multiple stations in a certain range near the epicenter before the earthquake would have anomalies at the same time, and the anomaly area was proportional to the magnitude. The larger the area is, the higher the magnitude. The seismic correspondence rate (seismic anomaly/all anomalies) reaches 65% in six months and reaches 80% within one year [15]; these works illustrated that pre-earthquake anomalies do not only exist at the single station closest to the epicenter. In addition, Fan et al. [16] also found groups of GVIP high anomalies at multiple stations in western China during the period of 2019–2020.

Based on the above work, we aim to investigate the multistation geomagnetic vertical intensity polarization high-value anomalies before the occurrences of several strong earthquakes, such as the Yangbi Mw 6.4 earthquake in Yunnan, the Mado Mw 7.4 earthquake in Qinghai, and the Luxian Mw 6.0 earthquake in Sichuan during 2020–2021. The correlations between the simultaneous GVIP anomaly and the occurrence of subsequent earthquakes in western China from 2015–2021 are discussed.

2. Data and Methods

The definition of geomagnetic vertical intensity polarization [2,9] is shown mathematically as:

$$Y_{zh}=\left|\frac{Z(\omega )}{H(\omega )}\right|$$

(1)

$$H(\omega )=\sqrt{H_x^2(\omega )+H_y^2(\omega )}$$

(2)

where Z(ω) is the spectral amplitude value of the geomagnetic vertical component, H(ω) is the spectral amplitude value of the full vector of the geomagnetic horizontal component, H_x(ω) is the north-south spectral value of the geomagnetic horizontal component, H_y(ω) is the east–west spectral value of the geomagnetic horizontal component, and ω is the circular frequency.

If the instrument frequency response needs to be considered, then

$$Z^{\prime }(\omega )=G_Z(\omega )\cdot Z(\omega )$$

(3)

$${H^{\prime }}_x(\omega )=G_{H_x}(\omega )\cdot H_x(\omega )$$

(4)

$${H^{\prime }}_y(\omega )=G_{H_y}(\omega )\cdot H_y(\omega )$$

(5)

where $Z^{\prime }(\omega )$ is the spectral amplitude value of the vertical component of the geomagnetic field observed by the instrument, ${H^{\prime }}_x(\omega )$ is the north-south spectrum of the horizontal component of the geomagnetic field observed by the instrument, ${H^{\prime }}_y(\omega )$ is the east–west spectrum of the horizontal component of the geomagnetic field observed by the instrument, and $G(\omega )$ is the instrument frequency response curve.

The fluxgate magnetometer G value is generally approximated as 1, but induction magnetometers need to be calibrated using the instrument. Induction magnetometers consist of a magnetic sensor and a recorder, so the product of the two is used in practice. In addition, attention should be given to the long-term stability of the G value.

The calculation process extracts the high ratio of polarization values to obtain the magnetic field signal from the crust of the primary source. The specific calculation process includes the following steps: first, obtaining the daily polarization value series of stations, fitting the daily polarization value series, and establishing the threshold curve that reflects the annual variation of polarization value; second, removing the daily polarization value below the threshold to obtain the daily polarization value series of high ratio, and using the fitting method to further eliminate the residual annual variation in the daily polarization value series of high ratio; third, the polarization values with higher than 2 times the mean-variance are considered anomalies. To reduce the influence of factors such as the geoelectrical structure, instruments, and the environment on the polarization value, we normalize the polarization value by 2 times the mean squared difference of the one-year window length of each station before drawing the contour of the spatial distribution of the polarization value and then subtract 1 after normalization. The high-value anomalies above 0 are called normalized zero-polarized values. Take the day with the most anomalous stations as the anomaly day and plot the normalized zero polarization values of all stations on the anomaly day. If the conditions of the grouping anomalies are met, the anomaly spatial contour map is drawn after superimposing the multiple anomalies in a group.

The distribution of fluxgate second sampling observation stations in western China in 2019 is given in Figure 1. Before 2014, the distribution of geomagnetic stations in western China was relatively sparse, and there were large blind monitoring areas in Xinjiang, Tibet, and Qinghai. Since 2015, geomagnetic stations have been rapidly constructed, and the current scale was formed. Therefore, the second sampling observation from the fluxgate magnetometer output of the geomagnetic observation stations west of 110° E in mainland China from 2015 to 2021 is adopted in this study.

3. Results

3.1. Multistationary Geomagnetic Vertical Intensity Polarization Anomalies

From January 2015 to December 2021, there were 18 geomagnetic vertical intensity polarization high-value anomalies in western mainland China [13], including nine anomalies over two or more anomalous regions. The spatial distribution of the nine anomalies is given in Figure 2 in the order of their occurrence.

The yellow and red areas are the anomalous regions, and their mathematical significance is that the polarization value exceeds two times the mean squared difference. Moreover, the area above the contour 0.2 is the predicted seismogenic region according to the statistical results of the example. The first four times all had two anomalous regions, then 20161202 and 20190924 had three anomalous regions, 20201015 had four anomalous regions, while 20170318 and 20191112 each formed five anomalous regions. Among them, the most significant anomaly is the anomaly of 18 January 2015, in Tibet, Qinghai, Xinjiang, and Gansu (Figure 2a), where a total of 16 high-value anomaly stations appeared, forming a very large anomaly area named “1-1” with an area of 288 × 104 km². The anomaly of 15 September 2016 formed two larger anomaly areas with a total area of more than 288 × 104 km². The anomaly of 15 October 2020 has the largest number of anomalous stations in the dense area of stations in the north-south seismic zone among the nine grouped anomalies, with 11 high-value anomalous stations and an anomalous area of 58 × 104 km².

3.2. Characteristics of Anomaly Spatial Groupings

Geomagnetic vertical intensity polarization anomalies have the characteristic of occurring in groups in space [17], i.e., geomagnetic vertical intensity polarization usually appears high in multiple locations synchronously in time, and a single anomaly can be spatially divided into several anomaly areas, each with an area larger than 3 × 10⁴ km² (e.g., [15]). In this study, such anomalies with several anomaly areas occurring spatially and synchronously in time are referred to as anomaly grouping or anomaly groups. After analyzing the spatial distribution map of the nine groups of anomalies in

Table 1. Geomagnetic vertical intensity polarization anomalies and spatial distributions of earthquakes in the following year.

Number	Date	Origin Time	Magnitude	Epicenter	Abnormal Lead Time (Day)	Abnormal Area Number	Abnormal Area (104·km2)	Earthquake Time Difference
1	20150118	20150415	5.8	Alxa Zuoqi, Inner Mongolia	87	1-1	288	-
		20150425	8.2	Nepal	97	1-1	288	10
		20150425	5.9	Tingri, Tibet	97	1-1	288	10
		20150703	6.5	Pishan, Xinjiang	166	1-1	288	79
		20151207	7.4	Tajikistan	323	1-1	288	157
		20160121	6.4	Menyuan, Qinghai	368	1-1	288	45
2	20150213	20150415	5.8	Alxa Zuoqi, Inner Mongolia	61	2-2	5	-
		20150703	6.5	Pishan, Xinjiang	140	2-1	88	79
		20151207	7.4	Tajikistan	166	2-1	88	157
3	20150718	20151207	7.4	Tajikistan	142	3-1	12	-
		20160121	6.4	Menyuan, Qinghai	178	3-2	7	36
4	20160915	20161017	6.3	Zadoi, Qinghai	32	4-2	123	-
		20161208	6.2	Hutubi, Xinjiang	84	4-1	80	52
		20161220	5.8	Qiemo, Xinjiang	96	4-1	80	12
		20170808	7.0	Jiuzhaigou, Sichuan	327	4-2	123	231
		20170809	6.6	Jinghe, Xinjiang	328	4-1	80	1
		20170916	5.8	Kuqa, Xinjiang	366	4-1	80	38
5	20161202	20161208	6.2	Hutubi, Xinjiang	7	5-1	18	-
		20170808	7.0	Jiuzhaigou, Sichuan	249	5-2	43	242
		20170809	6.6	Jinghe, Xinjiang	250	5-1	18	1
		20170916	5.8	Kuqa, Xinjiang	288	5-1	18	38
		20171118	6.9	Milin, Tibet	351	5-2	43	61
6	20170318	20170808	7.0	Jiuzhaigou, Sichuan	143	6-4	3	-
		20170809	6.6	Jinghe, Xinjiang	144	6-1	106	1
		20170916	5.8	Kuqa, Xinjiang	182	6-1	106	38
		20171118	6.9	Milin, Tibet	243	6-2	20	61
7	20190924	20200119	6.4	Payzawat, Xinjiang	117	7-1	63	-
		20200320	5.9	Tingri, Tibet	178	7-1	63	61
		20200626	6.4	Yutian, Xinjiang	276	7-1	63	98
		20200723	6.6	Nima, Tibet	303	7-1	63	27
8	20191112	20200626	6.4	Yutian, Xinjiang	227	8-1	8	-
		20200723	6.6	Nima, Tibet	254	8-1	8	27
9	20201015	20210521	6.4	Yangbi, Yunnan	218	10-3	5	-
		20210522	7.4	Maduo, Qinghai	219	10-2	58	1
		20210813	5.8	Maduo, Qinghai	299	10-2	58	80
		20210916	6.0	Luzhou, Sichuan	336	10-2	58	117

, all of these anomalies have the characteristic of anomaly grouping. Figure 2 shows that the earthquakes basically occurred in or near an area higher than 0.2. Each of the nine anomalies in Figure 2 occurred on the same day, and each anomaly produced more than two high-value anomaly areas. For example, two large and one small high-value anomaly zones were formed on 18 January 2015, 13 February 2015, 18 July 2015, 15 September 2016, and 2 December 2016, and multiple anomaly zones were formed on 18 March 2017, 24 September 2019, and 12 November 2019.

The MGVIPA formation events with multiple strong earthquakes occurring in and around the high-value anomaly area within 1 year are listed in Table 1, with a total of 10 events; in addition, the high-value anomaly on 15 October 2020, subsequently corresponds to the 22 May 2021, Mado 7.4 magnitude earthquake, and the corresponding period of this anomaly has not yet ended. According to the study of previous earthquake cases, there is still a possibility of strong earthquakes in this high-value area. In fact, on 16 September 2021, the Luzhou and Sichuan Mw 6.0 earthquakes occurred at the edge of this high-value region.

3.3. Earthquake Grouping Characteristics

To facilitate the discussion of the grouping characteristics of earthquakes occurring after anomalies, strong earthquakes with a magnitude Mw ≥ 5.8 that occurred within the high-value anomaly area and near its edges within 1 year after a group of anomalies are called a group of earthquakes. Earthquakes of magnitude 5.8 or more that occurred within 1 year after the anomaly appeared are marked on the anomaly distribution map (Figure 2).

From our analysis, 100% of the grouped earthquakes occurred within 1 year after the appearance of the grouped anomaly, and the epicenter location is still located inside or near the boundary of the high-value zone. However, there are two cases of earthquake groupings: first, all the earthquakes occurred within the same anomaly area, and no earthquakes occurred in other anomaly areas (e.g., [18]), such as Anomalies 7 and 8 in Table 1; second, earthquakes occurred in multiple anomaly areas of grouping anomalies, except for Anomalies 7 and 8, and the remaining seven groups of anomalies belong to this case, accounting for 70% of grouping anomalies.

In addition, as for the subsequent year of anomaly occurrences, multiple earthquakes of magnitude 5.8 or greater occurred in the subsequent year of five anomalies of 18 January 2015, 15 September 2016, 2 December 2016, 24 September 2019, and 12 November 2019. Specifically, on 15 September 2016 (Anomaly 4) and 24 September 2019 (Anomaly 7) had four earthquakes of magnitude 5.8 or greater in the anomaly area (e.g., [19,20]).

Among the 28 anomalous zones discussed in this section, 10 had more than two strong earthquakes in the follow-up, six had only one strong earthquake in the follow-up, and 12 had no strong earthquake in the follow-up.

3.4. Temporal Distribution of Earthquake Occurrences

The time difference between anomalies and subsequent corresponding earthquakes, as well as the seismic interval between sets of earthquakes corresponding to the same set of anomalies, are also presented in Table 1.

The subsequent first strong earthquakes after the appearance of the geomagnetic vertical intensity polarization anomaly ranged from 7–227 days. Of these, 1/9 occurred within 30 days, 4/9 occurred within 90 days, and 7/9 occurred within 180 days, with the longest being the Yutian 6.4 magnitude earthquake that occurred 227 days after the anomaly on 12 November 2019. In addition, the double earthquakes or series earthquakes have the longest corresponding time of 368 days.

The time difference between the two earthquakes in the group was 1–242 days, with twenty-two of the next earthquakes occurring within 100 days and five occurring more than 100 days after (e.g., [21]), including two events that occurred on 7 December 2015 (Tajikistan earthquake), two events on 8 August 2017 (Jiuzhaigou earthquake), and one on 16 September 2021 (Luzhou earthquake).

4. Discussion

Geomagnetic polarization anomalies are generally considered to be caused by electromagnetic radiation due to fault creep or rock rupture in the pre-earthquake source areas [22], which implies that multiple areas may experience fault creep or rock rupture simultaneously prior to an earthquake. According to the results of Ma et al. [23], the pre-earthquake regional stresses will appear to be concentrated at multiple points, and there is likely to be a similar connection between multiple strong earthquakes across the region or a strong earthquake sequence (e.g., [23]). The strong earthquakes in western mainland China discussed above are mainly distributed in the Qinghai-Tibet Plateau (i.e., Tibetan Plateau) and surrounding areas, and the occurrence of multiple earthquakes within a short period of time in the aforementioned areas may have the same tectonic genesis. The collision between the Indian and Eurasian plates led to the rapid uplift of the Tibetan Plateau, and the strong tectonic shifts have continued until the present time, resulting in strong seismic activity within and around the Tibetan Plateau [24]. The subsequent seismogenesis of multiple polarized high-value zones occurring simultaneously depends on the development of each gestation terrestrial body. Based on the solid-body earthquake gestation model [25], the areas can gradually accumulate a large amount of elastic potential energy during the slow relative plate motion. The high-speed block zones with solid-body properties and the ring-shaped distribution of a certain range of seismic activity enhancement zones on their periphery are caused by the higher rupture intensity of the source zone than the surrounding medium (e.g., [26,27]), and this difference also leads to the difference in electromagnetic radiation intensities.

From 1 January 2018 to 21 May 2021, only one Mw 4.9 earthquake occurred on 2 October 2019, in Tongren, Guizhou within the region of high geomagnetic vertical intensity polarization on 15 October 2020, and no Mw > 4 magnitude earthquake occurred in the rest of the region; from 3 October 2019 to 21 May 2021, no > Mw 4.0 earthquake occurred in the region, while the periphery of the high-value region, especially the periphery of the southwest boundary, had a higher concentration of seismic activity. In fact, 10 earthquakes of Mw > 4.0 occurred in the region during the 3 years from 1 January 2015 to 31 December 2017, including the Jiuzhaigou Mw 7.0 earthquake and two Mw > 5.0 earthquakes, which indicates that a more significant seismic calm occurred in the region from 3.5 to 1.5 years before the earthquake, which is in line with the characteristics of the solid-body that the peripheral areas rupture earlier than the solid-body interior (e.g., [28,29,30]). Xue et al. [31] summarized the characteristics of seismic activity before strong earthquakes occurred and concluded that the long axis of the gestation void area appeared from 370 to 780 km and 1–7 years before the earthquake, which agrees with the anomalous area counted in this study of 3~288 × 10⁴ km² and the circle diameter of 196~1970 km.

The results of geomagnetic sounding indicated that the main source and aftershock body of the Mado earthquake is in the junction area between the high conductivity of the middle and lower crust and the high-resistance body in the upper part [32]; similar results were obtained from studies in the eastern section of the East Kunlun and Longmenshan Fault Zones. These studies imply that strong earthquakes mostly occur near the significant electrical boundary zone at the block boundary, and this region is also the area where remarkable differences in electromagnetic radiation characteristics occur. The difference in the medium physical properties between and outside the seismogenic zones can last for a certain period. The reason why electromagnetic anomalies can be extracted only at specific times may be due to the limitation of the data processing level. The electromagnetic data using the geomagnetic vertical intensity polarization method still contain annual variations, monthly variations, and semimonthly variations, making the true high values of electromagnetic radiation submerged by noise or appear only under specific conditions or excited by external factors.

In addition to the electromagnetic signals [33,34], various forms of energy can be used as precursors to regular earthquakes for short-term forecasting, such as the acoustic (e.g., [35,36,37,38]) and neutron emission (e.g., [39,40,41,42]). The geological structure and acoustic emission (AE) can imply earthquake forecasting, for example, a large abrupt increase in the AE signals was observed at the Assisi earthquake [43]. The neutron flux variations can correspond to seismic activity in the Pamir region [44], and the relationship between neutron flux variations and the processes of the Earth’s crust can be used for the short-term prediction and monitoring of earthquakes (e.g., [41]). Furthermore, the geothermal variation near the surface was suggested to greatly affect the geomagnetic anomalies and indicates a potential for achieving specific understanding in earthquake prediction (e.g., [45,46]). Considering the spatial gradient of thermal sources beneath the active tectonic settings, geothermal variation including Joule heat is a candidate to explain the spatiotemporal distribution of geomagnetic changes and the time delay of the earthquake precursors (e.g., [47]). These approaches, from acoustic and neutron emission to the shallow geothermal investigation, are also important for understanding the observed earthquake precursors in Qinghai China, and these concerns and the potential mechanism (e.g., [48,49]) are left open for future multidisciplinary studies.

5. Conclusions

Through the above study, we draw the following conclusions:

(1)

The spatial distribution of seismic electromagnetic radiation anomalies is characterized by grouping, and the geomagnetic vertical intensity polarization usually has high values synchronized in time at multiple locations. Multiple regions of high-value anomalies in different anomalous regions are also synchronized in time.

(2)

Earthquakes occur in groups after the occurrence of anomalies within a period ranging from 3 months to 1 year. Multiple earthquakes occur in one or multiple anomalous zones within a group of anomalies.

(3)

The time interval of grouped earthquakes after anomalies is concentrated within 100~200 days.